Leveraging machine learning to uncover multi-pathogen infection dynamics across co-distributed frog families

Sampling

We obtained tissue samples through field collections and museum loans (Table S1). Individuals were largely sampled during summer breeding months from 2009–2023 under appropriate state and local permits (see Field Study Permissions below) and were archived at the Museum of Southwestern Biology (MSB), University of New Mexico (UNM). Post-metamorphic frogs were captured by hand, often along public roadways or near lentic waterways, while tadpoles were collected via dipnet. Individuals were either toe clipped and released at the site of capture or euthanized shortly after capture following protocols approved by the UNM Institutional Animal Care and Use Committee (Protocols 20-201006-MC & 23-201375-MC) and Florida State University (FSU) Animal Care and Use Committee (Protocols 1017 & 1313). Specifically, we applied 20% benzocaine to the ventral side of the frog, as described in the Guidelines for Use of Live Amphibians and Reptiles in Field and Laboratory Research (Herpetological Animal Care and Use Committee (HACC), 2004). Once the frog was completely unresponsive, we removed the heart as a secondary means of euthanasia. Dissection tools were flame sterilized between individuals, but not between tissue types collected. Tissues (toe, tail, muscle, and/or liver) were preserved either in 95% ethanol or tissue buffer (20% DMSO, 0.25 M EDTA, salt-saturated), or were immediately flash frozen in liquid nitrogen and kept frozen until DNA extraction.

In total, we obtained 1,281 samples from 32 states across the central and eastern U.S. (Fig. 1). Samples were selected to represent three major anuran families (Bufonidae n = 320, Hylidae n = 456, Ranidae n = 505) with four taxa in each family (11 species and one species complex; Fig. 1). Specifically, we targeted common and abundant species with overlapping distributions to facilitate comparisons across families while controlling for geographic and environmental variables.

Molecular methods

Genomic DNA was extracted from tissue samples that consisted of either external only (toe, tail n = 259), internal only (liver, muscle n = 363), or a combination of both tissue types (n = 586; Fig. S1). We note that the use of internal tissues is not recommended for detecting Bd (World Organisation for Animal Health (WOAH), 2021a); however, during dissection, Bd zoospores are occasionally transferred from the skin to internal tissues. We therefore expected occasional, but not reliable detection of Bd from internal tissues (demonstrated in a comparison by Torres López et al., 2024), which is why extractions of these samples were included in select downstream analyses.

We used the E.Z.N.A. Tissue DNA Extraction Kit (Omega Bio-tek, Inc, Norcross, GA, USA) following the manufacturer’s recommendations. Whole genomic DNA concentrations were determined using the Broad Range Qubit Assay (Invitrogen, Eugene, OR, USA) and standardized to DNA concentrations of 9–15 ng/µL. Pathogen presence and quantity were determined via quantitative PCR (qPCR) following established protocols (Bd—Boyle et al., 2004; Rv—Allender, Bunick & Mitchell, 2013; Pr—Karwacki et al., 2018). Reactions were prepared in a UV-sterilized, clean air workstation. For each pathogen-specific qPCR, we used 25 μL reaction volumes, each containing 5 μL of extracted DNA template, 2 µL each of 10 μM forward and reverse primers, 5 µL of 1 μM probe (Eurofins Genomics, Louisville, KY, USA), and 8 µL SsoAdvanced Universal Probes Supermix (Bio-Rad, Hercules, CA, USA). Serial dilutions of pathogen-specific gBlocks (idtDNA, Coralville, IA, USA) were mixed with 0.1 ng/μL yeast tRNA carrier (Eurofins Genomics) to prevent small nucleic fragments from binding to the low-bind tubes and were used in duplicate as standards on each plate. Two positive controls (5 μL at 10 ng/μL concentration of known positive samples: Bd—MSB:Herp:104601; Rv—MSB:Herp:104600; Pr—MSB:Herp:104643) and two negative controls (molecular grade water) were included on each plate. Plates were run on a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA), with cycling conditions of 95 °C for 5 min, followed by 40 cycles of 95 °C for 15 s and either 60 °C for 1 min (Bd/Rv) or 59 °C for 1 min (Pr).

Initial screening included two independent runs on each sample for each pathogen type; samples that tested negative in both runs were deemed uninfected, requiring no further screening. Partial and consistent positives underwent re-testing on two more independent plates to rule out false positives and to obtain final infection intensity values, resulting in at least four independent qPCR runs. Infection intensity values were calculated using the amplification curves (Cq) converted into starting quantities of pathogen DNA (SQ), with averages derived from the final two qPCR runs to account for differences in reagent sensitivity between runs. Individuals testing positive half of the time over the four plates of screening were categorized as “low” infections and included in all prevalence evaluations, but their infection intensity values were unreliable and were not considered in downstream analyses.

Calculating pathogen prevalence and infection loads

Pathogen prevalence was determined by summing the number of individuals found positive for each pathogen within a family, species, age, and sex, divided by the total number of samples within that category (Tables S2–S4). Because of the relatively low proportion of infected individuals within groups, we calculated 95% binomial confidence intervals using the logit method available in the binom package in R (R Core Team, 2023; Sundar, 2006). All individuals were taxonomically assigned and identified to the species or species complex level (i.e., Hyla chrysoscelis-versicolor). Samples lacking age class or sex identification at the time of collection were categorized as an unknown and were excluded from prevalence or infection intensity comparisons for those factors.

Reducing dimensionality of environmental variables

We obtained 11 temperature and nine precipitation bioclimatic variables at 30 s (~1 km²) resolution from the WorldClim 2.1 database, which includes climate data from 1970–2000 (Fick & Hijmans, 2017). Variables were centered, scaled, and reduced via principal components analyses (PCAs) using the prcomp() function in R, with separate analyses conducted for temperature and precipitation. The first two axes for temperature (hereafter, TPC1-2) explained a combined 88.7% of the variation and the first two axes for precipitation (PPC1-2) explained 95.8% of the variation (Fig. S2). Variable loadings indicated TPC1 primarily represented year-round temperatures, with higher values relating to warmer, more stable overall temperatures, while TPC2 represented mean diurnal temperature ranges, with higher values relating to more extreme daily temperature fluctuations. Similarly, PPC1 described the overall precipitation amount, with higher values indicating wetter conditions, whereas PPC2 represented precipitation variability, with higher values representing more variability in rainfall (Tables S5 and S6).

Identifying factors related to infection status and intensity via random forests

We employed the machine learning technique, RF, to identify important predictor variables related to infection status and infection intensity. We created classification models to predict pathogen status (infected or uninfected) for datasets with sufficient infection counts, including models encompassing all frogs (Bd_all: n = 1,281, Rv_all: n = 1,187) and family-specific models for Bd only (Bd_Bufonidae: n = 320, Bd_Hylidae: n = 456, Bd_Ranidae: n = 505). We also created a Bd infection status model excluding internal only and unknown tissues to investigate the impacts of tissue type on our results (Bd_external: n = 845). Additionally, we created an RF regression model for Bd infections (n = 216) to examine variables important to infection intensity (average SQ values). Models were created and assessed using the randomForests (Liaw & Wiener, 2002) and caret (Kuhn, 2008) R packages.

To address the large imbalance between infected and uninfected cases across the data, classification models were first balanced by either downsampling or upsampling. For downsampling, majority cases (uninfected) were randomly subsampled with replacement to equal the number of minority cases (infected). For upsampling, minority cases (infected) were randomly duplicated with replacement to equal the number of majority cases (uninfected). RF models were then optimized to predict both majority and minority cases by using the “strata” and “sampsize” functions to pull equal counts of infected and uninfected cases when creating forests. We assessed model overfitting via cross-validation and selected the model with the highest performance without overfitting, ultimately using downsampling for all models.

To further mitigate effects of small sample sizes of infected cases, we employed regularization via limited maximum depth by using the optimal mtry with the tuneRF() function of the randomForests package. We also used feature selection methods to reduce over-fitting of the data and determined variables for inclusion in our final classification models by comparing average out-of-bag (OOB) error rates. Cross-validations were performed using the predict() function from the stats package (R Core Team, 2023) to create a confusion matrix to assess predictive accuracy. To assess regression models, we used mean squared error and R² metrics. Final models were constructed with 2,000 trees and 100 permutations, runs were averaged over 100 iterations, and mean error rates were recorded. Average relative importance for each independent variable in the final model was measured using the importance() function and visualized using the ggplot2 package (Wickham, 2016).

Lastly, we evaluated final model performance, assessed predictive accuracy, and documented misclassification by splitting the data into training (70%) and testing (30%) datasets with equal proportions of positive infections and passed testing data through final models 100 times. We applied cross-validation to the final models to assess performance using the same methods described above.

One-way analyses of infection status and intensity

We used parametric and non-parametric tests via the stats R package to examine relationships between the potential factors associated with infection prevalence and intensity. We assessed model assumptions using the plotNormalhistogram() function in the rcompanion package and the var.test() function in the stats package. To assess differences in infection prevalence among categorical variables such as family, species, age class, sex, and tissue type, we used Pearson’s chi-squared tests for each pathogen type. We compared aspects of geography (latitude and elevation) and environmental variables (TPC1-2, PPC1-2) between infected and uninfected frogs using Student’s two-tailed t-tests. We performed logistic regressions to examine odds of infection across continuous variables such as latitude, longitude, elevation, TPC1-2, and PPC1-2 using the stats package and assessed model assumptions using the rcompanion package.

To examine differences in infection intensity between categories (family, species, age class, sex), we performed one-way ANOVAs assuming equal variances or Welch’s ANOVAs for tests with unequal variances. When comparing mean infection intensity between two groups (e.g., Rv species infection intensity, age class, and sex), we performed two-tailed Student’s t-tests. When significant differences were detected in ANOVA results, post-hoc analyses were carried out as follows. For parametric data, multiple two-tailed t-tests with Bonferroni correction or Tukey’s Honest Significant Difference tests were employed, while for non-parametric data, Games-Howell post-hoc tests were conducted using the rstatix package (Kassambara, 2023). Additionally, given the possibility of coinfections acting synergistically to increase negative infection outcomes (reviewed in Herczeg et al., 2021), we also evaluated differences in infection intensity between coinfected individuals and individuals with single infections of each pathogen using two-tailed t-tests. Finally, we performed linear regressions to assess infection intensity across continuous variables such as latitude, longitude, elevation, TPC1-2, and PPC1-2 using the stats package and assessed model assumptions using the erikmisc package (Erhardt, 2024). All R code and datasets are available from Figshare (https://doi.org/10.6084/m9.figshare.26849554), specimens are searchable via Arctos (arctos.database.museum; Cicero et al., 2024), and pathogen data are available on the Amphibian Disease Portal (Expedition GUID: https://n2t.net/ark:/21547/FkC2).

Pathogen prevalence summary

Approximately 20% (n = 256) of the 1,281 individual frogs screened were infected with at least one pathogen, 231 with single infections and 25 with coinfections. We found substantial variation in prevalence among pathogens (Figs. 2 and S3), with Bd exhibiting the highest prevalence (16.9%; 95% CI [14.9–19%]), followed by Rv (4.38%; 95% CI [3.35–5.7%]) and Pr (1.06%; 95% CI [0.618–1.82%]). Significant variation in Bd and Rv infections was evident both between and within family groups (Fig. 2A; Table 1), with Pearson residual scores indicating higher-than-expected counts for both pathogens in the family Ranidae (Fig. S4). For Bd, we found higher counts of infection in the species R. catesbeiana, R. clamitans, and R. sphenocephala (Fig. S5). Within Bufonidae and Hylidae, Bd infections were primarily associated with a single species in each family, namely A. americanus and P. crucifer (Fig. 2A; Table S2). For Rv, we also observed higher-than-expected counts in Ranidae, specifically R. catesbeiana and R. clamitans (Table 1, Fig. S5). Pr prevalence was consistently low with no discernable pattern of infection across taxonomic groups.

Rv infection prevalence was notably higher in juveniles compared to adults (Tables 1 and S7). We found no significant difference in pathogen prevalence between sexes for any pathogen (Table S8). For Bd, as expected, there were significant differences based on tissue type, with fewer infections from internal only tissues compared to external and combination tissue types (Table 1; Fig S6; Table S9). When we re-ran Bd analyses excluding internal only tissues, family and species remained significant factors for Bd infection status (Table S10).

Bd coinfections represented more than half of the individuals infected with Pr (7 of 13, or 53%) and about a third of those infected with Rv (18 of 52, or 35%). Whereas only 12% of samples infected with Bd represented a coinfection (25 out of 216).

Infection intensity summary

Overall, infection intensities were low, with the highest values observed in Bd infections (Table S11). Bd and Rv intensity varied significantly among families, with the highest Bd intensity found in Ranidae, followed by Hylidae and Bufonidae (Table 1, Fig. 2B). For Bd, significant differences in infection intensity across families were primarily driven by Ranidae and Bufonidae, while for Rv, differences were driven by Hylidae and Ranidae (Figs. 2B and 2D; Table S11). No clear pattern was found in Pr infection intensities across families due to small sample sizes. Interspecific variation revealed differences in infection intensity among species, with significantly higher Bd intensity in R. clamitans compared to A. americanus, R. catesbeiana, and R. pipiens (Fig. 2C; Table S12). Rv intensity was somewhat higher in R. catesbeiana than R. clamitans, but the difference was not significant (Fig. 2E; Table S12).

We found slightly higher Bd infection intensity in juveniles compared to adults, but there was no significant variation in infection intensity for any pathogen screened related to sex (Fig. S7; Tables S13 and S14). When we re-ran Bd analyses excluding internal only tissues, family and species remained significant for Bd infection intensity, while age class became significant (Table S10). We found no significant differences between Bd infection intensity for individuals with single infections compared to individuals coinfected with Rv or Pr, nor were there differences in Rv or Pr infection intensity between single infections and coinfections with Bd (Table 1).

RF classification and regression models

The Bd_all RF classification model identified species-level taxonomic rank as the most important factor associated with Bd infection status (Fig. 3A). Overall annual temperature (TPC1), latitude, and daily temperature fluctuations (TPC2) were also important predictors (Fig. 3A). The most important variable differed depending on which infection status class was being considered, with species as the top predictor of uninfected Bd cases (Fig. 3B), while overall temperature (TPC1) and daily fluctuations (TPC2) showed greater accuracy in predicting infected Bd cases (Fig. 3C). The model exhibited a consistent error rate across classification types at 23.6%, achieving an overall accuracy of 74%. Notably, the model performed similarly at classifying both infected and uninfected cases with sensitivity and specificity scores of 0.70 and 0.79, respectively. The top predictors and model performance for Bd infection status were consistent when we excluded internal only tissues (Fig. S8).

Though there were fewer infection points to sample from, the Rv_all classification model also identified species as the most important predictor of Rv, followed by family, latitude, and TPC1 (Fig. 3D). Species was the most important predictor of uninfected status, while latitude and overall temperature (TPC1) were the top contributing factors to infected Rv cases (Figs. 3E and 3F). The model had an average error rate of 30.2%, with slightly lower error rates for predicting infected status (28.0%) compared to uninfected status (32.5%). Model validation showed moderate accuracy at 71%, with 0.75 sensitivity and 0.71 specificity.

Family-specific RF classification models highlighted differences in top predictors depending on the dataset (Fig. S9). Specifically, for the Bd_Bufonidae and Bd_Hylidae models, species was the most important factor (Figs. S9A–S9F). These results differed from the Bd_Ranidae model, where we observed that species had low importance for predicting infection status and that daily temperature fluctuation (TPC2) was instead the most important factor, followed by latitude and year-round temperature (TPC1; Figs. S9G–S9I). The Bd_Ranidae model had an average error rate of 27.7%, with slightly lower error rates at predicting infected status (27.4%) compared to uninfected (28.0%). Model validation showed moderate accuracy at 70%, with 0.62 sensitivity and 0.74 specificity. For Bd_Bufonidae and Bd_Hylidae, cross-validation revealed higher model accuracy (87% and 89%, respectively), but both models showed signs of model class imbalance, with low sensitivity values (0.5, and 0.55) and high specificity (0.9, 0.91; Figs. S9A and S9D).

The Bd infection intensity RF regression model indicated latitude, overall temperature (TPC1), and precipitation (PPC1) were the top predictors. The model only explained 12.6% of variation, however, with validation showing large mean absolute error (MAE = 1,173) and low R² (0.14; Fig. S10).

Geographic and environmental predictors of infection

We observed significant differences in latitude between infected and uninfected individuals across all pathogens (Table 1). For Bd and Rv, infected individuals tended to occur at higher mean latitudes than uninfected (Figs. 4A and 4C), while Pr infections were only found at lower latitudes (Fig. 4B). Within species, the difference in mean latitude between Bd infected and uninfected samples was only statistically significant for R. clamitans, with infected samples tending to occur at higher latitudes (Fig. S11). A similar trend was observed in A. americanus and R. catesbeiana, while P. crucifer showed a trend in the opposite direction, with infected individuals tending to occur at lower latitudes.

Elevation and overall precipitation (PPC1) did not exhibit significant differences between infected and uninfected individuals for any pathogen, but differences in the other environmental variables (TPC1, TPC2, PPC2) were significant for one or more of the three pathogens (Table 1; Fig. 4). For Bd and Rv, infected individuals were found at lower mean TPC1 than uninfected, suggesting an association with cooler year-round temperatures, while Pr infections were only found at higher TPC1 values (Figs. 4D–4F). For TPC2, Bd infections were associated with lower daily temperature fluctuations, Pr infections with higher daily temperature fluctuations, while the differences between infected and uninfected individuals were not significant for Rv (Figs. 4G–4I). PPC2 only differed for Bd, with slightly higher mean values for infected individuals compared to uninfected (Figs. 4J–4L), indicating Bd infections were associated with higher variability in rainfall. When we re-ran Bd analyses excluding internal only tissues, only differences in TPC2 remained significant (Table S10).

No significant relationships were identified between infection intensity and geographic and environmental factors for any pathogen, but there was a trend of higher infection intensities associated with lower TPC2, which only increased in significance once we removed internal-only samples.

Prevalence and geographic distribution vary among pathogens

Batrachochytrium dendrobatidis (Bd)—Our survey is consistent with previous reports (Li et al., 2021), distribution models (Rödder, Schulte & Toledo, 2013; Olson et al., 2021), and public database counts (https://amphibiandisease.org; Koo et al., 2021) that indicate the widespread occurrence of Bd across much of the central and eastern United States. We found at least one frog infected with Bd in 29 of the 32 states screened, reinforcing the conclusion that Bd is likely endemic throughout much of this region (Longcore et al., 2007). We also note that northern states such as Maine, Vermont, and Indiana exhibited the highest proportion of Bd infections, with nearly half of all frogs screened testing positive. This high relative prevalence at higher latitudes could be explained by lower overall summer temperatures, which support fungal growth while limiting exposure to above-optimal temperatures that slow disease progression (Lindauer, Maler & Voyles, 2020). Similar patterns have been observed across the U.S. (Petersen et al., 2016; Sonn, Utz & Richards-Zawacki, 2019; Belasen et al., 2024) and Australia (Kriger, Pereoglou & Hero, 2007). However, it is important to acknowledge that our assessment likely underestimates the prevalence of Bd within southeastern regions since many of these samples were derived from internal tissues (e.g., liver, muscle) which are not reliable when screening for Bd (Fig. S1; World Organisation for Animal Health (WOAH), 2021a).

When we excluded samples derived from internal and unknown tissue types, TPC2 (corresponding to mean diurnal temperature range) remained significantly correlated with Bd prevalence. Experimental studies have shown reduced growth and reproduction in Bd when exposed to heat pulses compared to constant temperature treatments (Greenspan et al., 2017; Lindauer, Maler & Voyles, 2020). Likewise, our results suggest that areas with more stable daily temperatures—such as humid, temperate regions near bodies of water—may serve as important refugia for Bd, while more fluctuations in daily temperatures result in decreased pathogen prevalence. Although extensive research has examined climate in relation to Bd infection patterns (e.g., reviewed in Sasso, McCallum & Grogan, 2021), the effect of daily temperature fluctuations in relation to pathogen prevalence and distribution is underreported.

Amphibian Perkinsea (Pr)—Our results support three key insights regarding the distribution of Pr. First, our study suggests the current known distribution of Pr is likely underestimated, as we document novel infections in Missouri and Oklahoma (Isidoro-Ayza et al., 2017; Garner & Ruden, 2019; Midwest Association of Fish and Wildlife Agencies (MAFWA), 2020). Second, the westward range expansion to these states, along with reported mortality events as far north as Wisconsin (Isidoro-Ayza, 2019) and Alaska (Isidoro-Ayza et al., 2017), indicate that Pr likely possesses broad climatic tolerances. Although researchers have yet to successfully culture Pr in the lab to explore these limits, previous surveys have identified seasonal Pr outbreaks throughout the southeastern U.S. (Atkinson & Savage, 2023b). Lastly, the low infection prevalence we document here supports the hypothesis that Pr infections are currently localized and rare in adult frogs during summer months, though not impossible (Jones et al., 2012; Karwacki et al., 2018). Nonetheless, with rising concern for species in recent decline, potentially due to the introduction of Pr, such as Rana capito (Crawford et al., 2022; Devitt et al., 2023) and Rana sevosa (Atkinson, 2016), we recommend targeted screening across seasons of abundant co-distributed adult frogs that do not currently show significant population declines, such as R. catesbeiana, R. sphenocephala, and R. clamitans, to assess how these pathogens are potentially spread and maintained in the wild.

Ranavirus (Rv)—Rv has been documented widely across North America, as shown in the Global Ranavirus Reporting System (brunnerlab.shinyapps.io/GRRS_Interactive/) and subsequent reviews (e.g., Duffus et al., 2015; Brunner et al., 2021). It is important to note that the GRRS database includes reports of infections in both wild and captive frog populations. Therefore, even though Rv distributions are recorded widely, some states only have Rv outbreaks documented in captivity (Duffus et al., 2015). Our study adds such a novel case of Rv infection in wild frogs of Mississippi. The majority of Rv infections, however, were detected in northern states such as Maine and New Jersey, where 28% and 13% of individuals tested were found positive, respectively. Similar to Bd research, and in line with previous observational studies (e.g., Youker-Smith et al., 2018), we observed significantly higher Rv prevalence in areas characterized by cooler year-round temperatures. We note, however, that the higher Rv prevalence we found in high-latitude areas may be influenced by the species sampled at northern sites (primarily Rana catesbeiana and R. clamitans) and the age class of those sampled (many juveniles), both of which are factors that could conceivably correspond with increased infection.

The relatively low overall prevalence of Pr and Rv in our study was likely influenced by a few factors. First, while both liver and toe clips can reliably detect Rv (Torres López et al., 2024) and are recommended diagnostic tissues (World Organisation for Animal Health (WOAH), 2021b), the use of toe clips for Pr detection has not yet been validated. Additionally, the absence of larval stages in our sampling likely had a strong impact on prevalence of both Pr and Rv. While infections (Gray, Miller & Hoverman, 2009; Karwacki et al., 2018) and mortality (Teacher, Cunningham & Garner, 2010) have been documented in adult frogs, mass die-off events for both pathogens primarily occur before and during metamorphosis (Rv: Green, Converse & Schrader, 2002; Pr: Isidoro-Ayza, 2019). This trend is at least partly explained by the underdeveloped immune response documented in hosts prior to and during larval development (Gray, Miller & Hoverman, 2009; Miller, Gray & Storfer, 2011; Herczeg et al., 2021), as well as the habitat overlap between these water-borne pathogens and larval stages (Miller, Gray & Storfer, 2011; Itoïz et al., 2022). We did, however, record higher than expected Rv infections in post-metamorphic juvenile frogs compared to adults, perhaps explained by the continued lag in host immune function between metamorphosis and sexual maturity (Rollins-Smith, 1998; Gantress et al., 2003). Lastly, the presence of both Rv and Pr infections in adult frogs in our study may either represent coincidental infections prior to clearance or be related to other mechanisms, such as coinfections, discussed below.

Coinfections account for many cases of rare pathogen infections

Even though we report no significant differences in infection intensity between single and coinfected individuals, we still document a surprising proportion of coinfections between Bd-Pr and Bd-Rv. This observation, along with findings from other recent multi-pathogen screening efforts (e.g., Landsberg et al., 2013; Karwacki, Martin & Savage, 2021; Atkinson & Savage, 2023b), support the idea that coinfections may facilitate heightened susceptibility in frogs (Herczeg et al., 2021), but more work is needed to understand how these coinfections are affecting seemingly robust species such as Rana catesbeiana and R. clamitans. Though we did not find coinfections between Pr and Rv, which primarily impact larval stages, infection with Bd may increase susceptibility of adult frogs to Pr and Rv. One concerning aspect of this finding is that given their increased immune function, adults can maintain chronic, low-level infections, thus providing a potential mechanism for pathogen spread and persistence in populations. Further investigation in a controlled setting is warranted to assess how coinfections at different life stages may impact the infection dynamics and outcomes we observe in natural environments and communities.

Infection prevalence differs across host taxonomy

Our study adds to a growing body of research that aims to understand differences in infection prevalence among hosts. Specifically, we found significantly higher Bd and Rv prevalence in four species in the family Ranidae compared to other families (Fig. 2; Fig. S2), while Pr counts were too low to draw significant conclusions. Higher pathogen prevalence in the family Ranidae has been reported from other surveys of wild communities (e.g., Olson et al., 2013; Rothermel et al., 2016; Karwacki et al., 2018), and increased susceptibility has been directly reported in controlled experiments (e.g., Hoverman, Gray & Miller, 2010; Gahl, Longcore & Houlahan, 2012). One potential explanation for increased pathogen prevalence in this family could be proximity and time spent in permanent and semi-permanent breeding ponds that act as pathogen refugia (Hoverman et al., 2011; Greenberg, Palen & Mooers, 2017).

Two species, Rana catesbeiana and R. clamitans, accounted for 38% of all infected frogs and 64% of all coinfections in our dataset. Both species have been classified as pathogen reservoirs which contribute to the maintenance of Bd and Rv across their native ranges over time (Rothermel et al., 2016; Yap et al., 2018; Brunner et al., 2019; Hoverman, Chislock & Gannon, 2019; Hossack et al., 2023). In addition, the global invasion of pathogen-positive R. catesbeiana has been linked to the spread of both Bd and Rv with primarily negative effects for native species (reviewed in Atkinson & Savage, 2023a). Other species-level differences observed in our study, such as increased Bd prevalence in Anaxyrus americanus and Pseudacris crucifer compared to other species in their respective families, also warrants further investigation to identify potential drivers of increased susceptibility outside of primary habitat type. Research into mechanisms of resistance and tolerance (e.g., Eskew et al., 2015), host-pathogen coevolutionary history (e.g., Carvalho et al., 2024) and immunogenetic adaptation (e.g., Trujillo et al., 2021) focused on additional host species could provide further important insights.

Infection intensity varies across host taxonomy and age class

In amphibians, infection intensity serves as a proxy for infection outcomes, with higher pathogen loads often corresponding with increased spread and host mortality (Vredenburg et al., 2010). Our survey of widespread species that are not currently documented as experiencing pathogen-induced decline, still showed significant associations between infection intensity, host species, and age class. Specifically, we found the highest Bd infection intensities in Rana clamitans (Fig. 2B) and significantly higher Bd intensity in juveniles compared to adults (Table S1). The majority of Bd-positive juvenile frogs in our dataset were R. clamitans, however, meaning that our data lack the power to distinguish infection intensity variation between age classes and species-level correlates. Importantly, multiple studies have shown that Bd infection intensities increase during winter months (Savage, Sredl & Zamudio, 2011), which our study also fails to capture. Therefore, more research is needed to assess the differences in Bd intensity across life stages within R. clamitans and across seasons to extrapolate the mechanisms responsible for this pattern.

Broadly, age class has been associated with pathogen-induced host mortality, but these results have been mixed. For example, Bd-induced mortality has been documented during metamorphosis (Humphries et al., 2024) and in juvenile stages (Abu Bakar et al., 2016), but also the inverse has been recorded, where juveniles harbor lower infection loads compared to adult frogs (Bradley et al., 2019). This inconsistency demonstrates the contextual nature of infection dynamics, where multiple aspects can influence infection outcomes. Factors such as length of exposure (Lips, 2016; Bielby et al., 2021), host skin peptides and microbiomes (Woodhams et al., 2007; Bernardo-Cravo et al., 2020), infection history (Greenberg, Palen & Mooers, 2017), and phylogenetic constraint (Hoverman et al., 2011; Longo, Lips & Zamudio, 2023) have all been shown to affect pathogen-specific infection intensities and/or host mortality. Continued work documenting species-specific mechanisms of pathogen tolerance and resistance is needed to better understand the drivers of infection intensity differences and host survival.

Application and limitations of RF models in multi-pathogen dynamics

Our application of balanced RF classification models builds on efforts that have used machine learning to address factors influencing both species-level (Murray et al., 2013) and site-specific (Atkinson & Savage, 2023b; Roth, Griffis-Kyle & Barnes, 2024) amphibian infection dynamics. Machine learning techniques have also been applied to examine projected distributions of Bd in the context of future climate change (Xie, Olson & Blaustein, 2016). We demonstrated the utility of RF models to quickly and efficiently provide valuable insights into the complex patterns encountered in observational disease ecology studies (Cutler et al., 2007). Though there are limitations such as over-fitting due to small sample sizes or amplifying noise in the data rather than true relationships, models built from balanced training datasets with ample sample sizes can achieve moderately low out-of-bag error rates and have high predictive accuracy, specificity, and sensitivity. The application of RF models to other wildlife disease systems such as avian malaria-causing parasites (e.g., Aželytė et al., 2023), Pseudogymnoascus destructans, which causes white-nose syndrome in bats, or Ophidiomyces ophiodiicola, which causes snake fungal disease, are other promising avenues for discovery.