DOI: 10.5061/dryad.mgqnk998g
title: Do fungi look like macroparasites? Quantifying the patterns and mechanisms of aggregation for host-fungal parasite relationships

description: Most hosts contain few parasites, whereas few hosts contain many. This
 pattern, known as aggregation, is well-documented in macroparasites where
 parasite intensity distribution among hosts affects host-parasite
 dynamics. Infection intensity also drives fungal disease dynamics, but we
 lack a basic understanding of host-fungal aggregation patterns, how they
 compare to macroparasites, and if they reflect biological processes. To
 begin addressing these gaps, we characterized aggregation of the fungal
 pathogen Batrachochytrium dendrobatidis (Bd) in amphibian hosts. Utilizing
 the slope of Taylor’s Power Law, we found Bd intensity distributions were
 more aggregated than many macroparasites, conforming closely to lognormal
 distributions. We observed that Bd aggregation patterns are strongly
 correlated with known biological processes operating in amphibian
 populations, such as epizoological phase (i.e., invasion, post-invasion,
 and enzootic), and intensity-dependent disease mortality. Using
 intensity-dependent mathematical models, we found evidence of evolution of
 host resistance based on aggregation shifts in systems persisting with Bd
 following disease-induced declines. Our results show that Bd aggregation
 is highly conserved across disparate systems and contains signatures of
 potential biological processes of amphibian-Bd systems. Our work can
 inform future modeling approaches and be extended to other fungal
 pathogens to elucidate host-fungal interactions and unite host-fungal
 dynamics under a common theoretical framework.
description: This dataset contains records of aggregated populations of
 amphibians and their fungal pathogen <em>Batrachochytrium
 dendrobatidis </em>(Bd) loads or intensities. Bd infection intensity
 obtained from amphibian skin swabs collected in the field. Bd loads were
 obtained through DNA extraction and qPCR, which detects the number of
 genomic equivalents or ITS1 copy number of Bd on amphibian skin. These
 methods were standardized across four dataset utilized for this study and
 published elsewhere. These included data taken from sites in São Paulo,
 Brazil [1; and manuscript in review], the East Bay region [2] and Sierra
 Nevada Mountains [3;4] in California, and across four states in the
 Eastern US [5]. Samples within each dataset were grouped based on host
 species, life stage (larva (i.e., tadpole in anuran species), subadult, or
 adult), research site, season (Brazil: Wet or Dry; East Bay and Sierra:
 Summer; Eastern US: winter, spring, summer, or fall), and year. 
 Martins RA, Greenspan SE, Medina D, Buttimer S,
 Marshall VM, Neely WJ, et al. Signatures of functional bacteriome
 structure in a tropical direct-developing amphibian species. Animal
 Microbiome [Internet]. 2022;4(1):40. Available from: https://doi.org/10.1186/s42523-022-00188-7 Wilber MQ, Johnson PTJ, Briggs CJ, Seabloom E. Disease hotspots or hot species? Infection dynamics in multi‐host metacommunities controlled by species identity, not source location. Ecology letters. 2020;23(8):1201–11. Vredenburg VT, Knapp RA, Tunstall TS, Briggs CJ. Dynamics of an emerging disease drive large-scale amphibian population extinctions. Proceedings of the National Academy of Sciences. 2010 May;107(21):9689–94. Briggs CJ, Knapp RA, Vredenburg VT. Enzootic and epizootic dynamics of the chytrid fungal pathogen of amphibians. Proceedings of the National Academy of Sciences of the United States of America. 2010 May;107(21):9695–700. Wilber MQ, Ohmer MEB, Altman KA, Brannelly LA, LeSage EH, LaBumbard BC, et al. Once a reservoir, always a reservoir? Seasonality affects the pathogen maintenance potential of amphibian hosts. Ecology. 2022;103:e3759.
description: # Do fungi look like macroparasites? Quantifying the patterns and
 mechanisms of aggregation for host-fungal parasite relationships
 [https://doi.org/10.5061/dryad.mgqnk998g](https://doi.org/10.5061/dryad.mgqnk998g) ## Description of the data and file structure * Repository of code and data for the manuscript "Do fungi look like macroparasites? Quantifying the patterns and mechanisms of aggregation for host-fungal parasite relationships". ### Files and variables #### File: abund\_and\_prev.csv **Description:** Contains records of Sierra sites over multiple years and includes information on amphibian abundance, Bd prevalence, mean fungal load, and Poulin's D. Fields with a value of NA are defined as "**Not Available**" and indicates unknown or unrecorded data. ##### Variables * site_id: (numeric) unique identifier for specific sampling sites * year: (numeric) year in which samples were collected * capture_life_stage: (character) life stage of grouped amphibians sampled (three levels: larva, subadult, and adult) * bd_present: (numeric) indicates if amphibian group had Bd-positive individuals (1) or not (0) * num_infected: (numeric) the number of Bd-positive individuals in the sampled group * num_samps: (numeric) the number of individuals in the sampled group * adult_abundance: (numeric) the total adult abundance at the site at the time of sampling * prev: (numeric) the calculated prevalence based on number of Bd-positive individuals compared to total number of individuals sampled * mean_bd: (numeric) the average Bd load within the amphibian group (-1 indicates no Bd-positive individuals) * poulin_D: (numeric) Poulin's Discrepancy Index (see Poulin R. The disparity between observed and uniform distributions: A new look at parasite aggregation. International Journal for Parasitology. 1993;23(7):937–44.) #### File: analysis\_aggregation\_dataset.csv **Description:** Contains groups aggregated by dataset, site, species, life stage, and season. This dataset includes aggregated values of mean, variance, Poulin's D, coefficient of variation metrics, and AIC values on fitted distributions. Fields with a value of NA are defined as "**Not Applicable**" and indicate that the value could not be calculated based on the specified methods for analysis (see each variable description for details). There are exceptions for the last three variables (prev, phase, site_id) in which NA is defined as "**Not Available**". ##### Variables * site: (character) unique identifier for specific sampling sites * location: (character) study country (two levels: brazil and usa)  * region: (character) study region ( : santa_virginia, louisiana, tennessee, vermont, pennsylvania, east_bay, sierra_nevada) * species: (character) amphibian host species * life_stage: (character) life stage of grouped amphibians sampled (four levels: larva, subadult, adult, and unknown) * temporal_group: (numeric/factor) biologically relevant temporal group for each dataset [Sierra & East Bay: grouped by year; Eastern US: grouped by austral season; Brazil: grouped by Wet/Dry season] * season: (character) astronomical season for Sierra, East Bay, and Eastern US datasets and wet/dry for Brazil * year: (numeric) year in which samples were collected * N: (numeric) the number of individuals in the group * num_infected: (numeric) the number of Bd-positive individuals in the group * num_infected_discrete: (numeric) the number of Bd-positive individuals in the group after rounding fungal loads to a whole number (i.e. very low numbers were considered zero and therefore not infected) * mean_natural_with_zeros: (numeric) the arithmetic mean Bd load for the group, including zeros (Bd-negative individuals) * var_natural_with_zeros: (numeric) the variance Bd load for the group, including zeros (Bd-negative individuals) * mean_natural_without_zeros: (numeric) the arithmetic mean Bd load for the group, excluding zeros (Bd-negative individuals); NAs occur when the group has no Bd-positive individuals and therefore, Bd-positive mean cannot be calculated  * var_natural_without_zeros: (numeric) the variance Bd load for the group, excluding zeros (Bd-negative individuals); NAs occur when the group has no or only one Bd-positive individual and therefore, Bd-positive variance cannot be calculated  * mean_log10_without_zeros: (numeric) the arithmetic mean of log10 Bd load for the group, excluding zeros (Bd-negative individuals); NAs occur when the group has no Bd-positive individuals and therefore, Bd-positive log mean cannot be calculated  * var_log10_without_zeros: (numeric) the variance of log10 Bd load for the group, excluding zeros (Bd-negative individuals); NAs occur when the group has no or only one Bd-positive individual and therefore, Bd-positive log variance cannot be calculated  * mean_discrete_ceiling_with_zeros:  (numeric) the arithmetic mean of discretized Bd load (rounded up) for the group, including zeros (Bd-negative individuals) * var_discrete_ceiling_with_zeros:  (numeric) the variance of discretized Bd load (rounded up) for the group, including zeros (Bd-negative individuals) * mean_discrete_floor_with_zeros:  (numeric) the arithmetic mean of discretized Bd load (rounded down) for the group, including zeros (Bd-negative individuals) * var_discrete_floor_with_zeros:  (numeric) the variance of discretized Bd load (rounded down) for the group, including zeros (Bd-negative individuals) * nbd_k: (numeric) dispersion parameter for fitted  negative binomial distribution to discretized load data; NAs occur when discrete distributions were not fitted to the group because there were fewer than three Bd-positive individuals * nbd_mu: (numeric) mean parameter for fitted  negative binomial distribution to discretized load data; NAs occur when discrete distributions were not fitted to the group because there were fewer than three Bd-positive individuals * nbd_aic: (numeric) AIC value of modeling distribution as negative binomial; NAs occur when discrete distributions were not fitted to the group because there were fewer than three Bd-positive individuals * pois_lam:  (numeric) rate parameter for fitted  Poisson distribution to discretized load data; NAs occur when discrete distributions were not fitted to the group because there were fewer than three Bd-positive individuals * pois_aic: (numeric) AIC value of modeling distribution as Poisson; NAs occur when discrete distributions were not fitted to the group because there were fewer than three Bd-positive individuals * gamma_aic: (numeric) AIC value of modeling distribution as Gamma; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * gamma_shape:  (numeric) shape parameter for fitted  Gamma distribution to continuous load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * gamma_rate:  (numeric) rate parameter for fitted  Gamma distribution to continuous load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * lnorm_aic:  (numeric) AIC value of modeling distribution as Lognormal; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * lnorm_meanlog:  (numeric) location parameter for fitted  Lognormal distribution to continuous load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * lnorm_sdlog: (numeric) scale parameter for fitted  Lognormal distribution to continuous load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * weibull_aic: (numeric) AIC value of modeling distribution as Weibull; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * weibull_shape: (numeric) shape parameter for fitted  Weibull distribution to continuous load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * weibull_scale: (numeric) scale parameter for fitted  Weibull distribution to continuous load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * exp_aic: (numeric) AIC value of modeling distribution as Exponential; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * exp_rate: (numeric) rate parameter for fitted  Exponential distribution to continuous load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * norm_gof_pvalue_sw: (numeric) p-value from test for normality (Shapiro-Wilk) of log10-transformed load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * norm_gof_statistic_sw: (numeric) test statistic from test for normality (Shapiro-Wilk) of log10-transformed load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * norm_gof_pvalue_ks:  (numeric) p-value from test for normality (Kolmogorov-Smirnov) of log10-transformed load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * norm_gof_statistic_ks:  (numeric) test statistic from test for normality (Kolmogorov-Smirnov) of log10-transformed load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * norm_gof_pvalue_kslc: (numeric) p-value from test for normality (Lilliefors-corrected Kolmogorov-Smirnov) of log10-transformed load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * norm_gof_statistic_kslc:  (numeric) test statistic from test for normality (Lilliefors-corrected Kolmogorov-Smirnov) of log10-transformed load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * norm_gof_pvalue_ad: (numeric) p-value from test for normality (Anderson-Darling) of log10-transformed load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * norm_gof_statistic_ad:  (numeric) test statistic from test for normality (Anderson-Darling) of log10-transformed load data; NAs occur when continuous distributions were not fitted to the group because there were fewer than ten Bd-positive individuals * poulin_D:  (numeric) Poulin's Discrepancy Index calculated on Bd loads (see Poulin R. The disparity between observed and uniform distributions: A new look at parasite aggregation. International Journal for Parasitology. 1993;23(7):937–44.); NAs occur when groups had fewer than two Bd-positive individuals * poulin_D_log:  (numeric) Poulin's Discrepancy Index calculated on log10 loads (see Poulin R. The disparity between observed and uniform distributions: A new look at parasite aggregation. International Journal for Parasitology. 1993;23(7):937–44.); NAs occur when groups had fewer than two Bd-positive individuals * cv_log: (numeric) Coefficient of Variation (CV) calculated on log10-tranformed loads; NAs occur when groups had fewer than two Bd-positive individuals * log_cv: (numeric) Log Coefficient of Variation (CV) calculated on Bd loads; NAs occur when groups had fewer than two Bd-positive individuals * cv: (numeric) Coefficient of Variation (CV) calculated on Bd loads; NAs occur when groups had fewer than two Bd-positive individuals * dataset: (character) name of the dataset (four levels: brazil, east_bay, sierra_nevada, and serdp (i.e., eastern us) * prev:  (numeric) the calculated prevalence based on number of Bd-positive individuals compared to total number of individuals sampled for Sierra data (linked from abund_and_prev.csv); NAs occur for all datasets except for sierra_nevada, as the data was only available and used for sierra_nevada  * phase: epizoological phase based on putative knowledge and Bd prevalence [enzootic: Northern Sierra populations; invasion: Southern Sierra populations with less than 0.5 prevalence; epizootic: Southern Sierra populations with greater than 0.5 prevalence]; NAs occur for all datasets except for sierra_nevada, as the data was only available and used for sierra_nevada  * site_id: (numeric) unique identifier for specific sampling sites for Sierra data (linked from abund_and_prev.csv); NAs occur for all datasets except for sierra_nevada, as the data was only available and used for sierra_nevada  ## Code/software * Data was compiled and analyzed in R 4.3.2 * `code`: Contains the following R scripts to reproduce analyses and figures * `01_Aggregation_Groups.R`: Contains the code to aggregate raw data into groups and perform group-level analyses (calculating fungal intensity means and variances, fitting distributions, and calculating individual-group aggregation metrics (Poulin's D and variations of the Coefficient of Variation). -- Output: `data/formatted/summary_aggregation_data.csv` -- Note: Raw data is provided through previously-published studies: DOI 10.1073/pnas.0914111107; DOI 10.1073/pnas.0912886107; DOI 10.1186/s42523-022-00188-7; DOI 10.1111/ele.13518; DOI 10.1002/ecy.3759 * `02_Assign_Epi_Phase.R`: Contains code to classify Sierra sites as "Invasion", "Post-Invasion", or "Enzootic" phase -- Output: `data/formatted/abund_and_prev.csv` `data/formatted/summary_sierra_with_epi_phase.csv` `data/formatted/analysis_aggregation_data.csv` -- Note: Raw data is provided through previously-published studies: DOI 10.1073/pnas.0914111107; DOI 10.1073/pnas.0912886107 * `03_Fig.1_TPL_Plot.R`: Contains the code used to create Figure 1 and calculate the slope of Taylor's Power Law (TPL) for each dataset. * `04_TPL_Life_Stage.R`: Contains generalized linear mixed models, modeling the relationship between the log mean fungal load on the log variance of fungal load, incorporating region and life stage as possible predictor variables. * `05_Fig.2_Distribution_Fit_Boxplot.R`: Contains code to create Figure 2 and compare which distributions fit the data best. * `06_Fig.3_Trajectory_Plot.R`: Contains code to produce Figure 3. * `07_Mean_Agg_Quadratic_Model.R`: Contains code that fits linear and quadratic models to the Sierra dataset and compares the performance. Also compares if this relationship changes based on life stage or epizoological phase. * `08_IPM_Fig.S5-S8.R`: Contains the code for the Integral Projection Model (IPM), as well as different parameters used. Contains code to produce Supplemental Materials Figures S5-S8. -- Output: `data/model/eco_df.csv` -- `data/model/evo_df.csv` * `09_Fig.4_Empirical_Model_Agg_Plot.R`: Contains the code to produce Figure 4. * `10_Supp_Table.S1.R`: Contains code to extract values for Supplemental Materials Table S1. * `11_Supp_Fig.S2_Spec_TPL.R`: Contains generalized linear mixed models, modeling the relationship between the log mean fungal load on the log variance of fungal load, incorporating region and species as possible predictor variables. Includes code to create Figure 3. * `12_Supp_Fig.S3_GOF.R`: Contains code to create Supplemental Materials Figure S3. Also contains code to extract information on goodness of fit (GOF) of a normal distribution on log-transformed data and find proportion of log-transformed data that conforms to a normal distribution / proportion of natural data that conforms to a log-normal distribution. * `13_Supp_Fig.S4_PA_GreenFrogs.R`: Contains code to create Supplemental Materials Figure S4. ## Access information Data was derived from the following sources: * Vredenburg VT, Knapp RA, Tunstall TS, Briggs CJ. [Dynamics of an emerging disease drive large-scale amphibian population extinctions.](https://doi.org/10.1073/pnas.0914111107) Proceedings of the National Academy of Sciences. 2010 May;107(21):9689–94. * Briggs CJ, Knapp RA, Vredenburg VT. [Enzootic and epizootic dynamics of the chytrid fungal pathogen of amphibians](https://doi.org/10.1073/pnas.0912886107). Proceedings of the National Academy of Sciences of the United States of America. 2010 May;107(21):9695–700. * Martins RA, Greenspan SE, Medina D, Buttimer S, Marshall VM, Neely WJ, et al. Signatures of functional bacteriome structure in a tropical direct-developing amphibian species. Animal Microbiome [Internet]. 2022;4(1):40. Available from: [https://doi.org/10.1186/s42523-022-00188-7](https://doi.org/10.1186/s42523-022-00188-7) * Wilber MQ, Johnson PTJ, Briggs CJ, Seabloom E. Disease hotspots or hot species? Infection dynamics in multi‐host metacommunities controlled by species identity, not source location. Ecology letters. 2020;23(8):1201–11. * Wilber MQ, Ohmer MEB, Altman KA, Brannelly LA, LeSage EH, LaBumbard BC, et al. Once a reservoir, always a reservoir? Seasonality affects the pathogen maintenance potential of amphibian hosts. Ecology. 2022;103:e3759.