Pathogen load predicts host functional disruption: A meta-analysis of an amphibian fungal panzootic

Supplementary Information

GitHub Repository
The Rmarkdown file can be downloaded from the Code drop down menu (top right).

This supplementary file contains the R workflow for processing and analysing the raw data, and creating figures for the manuscript ID: FE-2022-00608 titled “Pathogen load predicts host functional disruption: A meta-analysis of an amphibian fungal panzootic”.

---

Supplementary methods

Chytridiomycosis disease description

Both chytrid fungus species (Batrachochytrium dendrobatidis, Bd, and Batrachochytrium salamandrivorans, Bs) colonies the host skin via free-swimming motile zoospores and complete their life cycle within the epidermal layer of the skin (Berger et al., 2005; Martel et al., 2013). The chytrid fungus affects the host’s skin function through both the physical breakdown and chemical disruption of the epidermis (Brutyn et al., 2012; Voyles et al., 2009; Wu et al., 2019). Because amphibians rely heavily on their skin for maintaining physiological homeostasis such as electrolyte and fluid balance, and respiration (Hillyard et al., 2008), the resulting cutaneous disruption directly affects the host’s behaviour and physiology (Campbell et al., 2012; Kindermann et al., 2017; Peterson et al., 2013; Van Rooij et al., 2015; Voyles et al., 2007). The cause of death is proposed to be cardiac arrest due to low circulating ion levels (Voyles et al., 2009). The presence of the chytrid fungus also elicits a range of immunological responses (Grogan et al., 2018a; Rollins-Smith and Le Sage, 2021; Rollins-Smith et al., 2011) which in turn alters the host’s energy metabolism, behaviour, activity and movement patterns, and body condition (Grogan et al., 2018b; Richards-Zawacki, 2010; Wu et al., 2018).

Abstract inclusion and data exclusion criteria

Title and abstract screening of 1,038 records after duplicates were removed was conducted in Rayyan (Ouzzani et al., 2016), where the abstracts were screened for whether the study performed experimental infections, whether anurans were tested, and whether a functional trait or survival was measured. In the full-text screening, data were excluded for the following criteria:

• Non-English literature, books without links to primary source, and non-peer-reviewed articles.
• Studies that extracted Bd load from dead anurans in the wild, as death in the wild may not be directly caused by Bd directly.
• Studies that did not quantify Bd load, and only presented Bd prevalence or categorical detection of Bd (positive or negative detection).
• Studies that capped Bd load based on maximum standard runs, e.g. Brannelly et al. (2015).
• Studies that exposed amphibians to Bs, or co-exposed Bd with Bs.
• Exposed amphibians to Bd during the egg or larval life stage.
• Treatments with additional experimental stressors e.g., agrochemicals, pharmaceuticals, UV exposure, presence of predators.
• Data from isolated or cultured cells (in-vitro results), as they do not provide functionally relevant responses.
• Non-quantifiable traits such as histological sections, electron-microscope images, multivariate principal components with arbitrary units.
• Molecular responses not directly related to functional traits such as transcriptomics, proteomics, and mRNA expression.

Conversion statistics

The correlation coefficient (r) was calculated from inferential statistics with the following equations (1) from Lipsey and Wilson (2001), Nakagawa and Cuthill (2007), and Noble et al. (2017):

\[\begin{equation} r = \sqrt{[t^2 / (t^2 + d.f.)]} = \sqrt{[F / (F + d.f.)]} = \sqrt{(\chi^2 / n)}, \tag{1} \end{equation}\]

where \(t\) = t-statistics, \(F\) = F-statistics, \(d.f.\) = model denominator degrees of freedom, \(\chi^2\) = chi-squared statistics, and \(n\) = sample size. Effect sizes from inferential statistics were only retained when directions could be determined (e.g., \(F\) and \(\chi^2\) alone do not contain directional information). Studies that presented standardised mean differences, as Cohen’s d, were converted to r with the following equation (2) from Borenstein et al. (2009):

\[\begin{equation} r = \frac{d}{\sqrt{d^2 + a}}, \tag{2} \end{equation}\]

where \(a\) is a correction factor (3) for cases where \(n_1 \neq n_2\),

\[\begin{equation} a = \frac{(n_1 + n_2)^2}{n_1n_2}, \tag{3} \end{equation}\]

else, \(a\) = 4.

Keyword co-occurrence network analysis

During the systemetic review, a keyword co-occurrence network analysis was performed using the näive Boolean search from Web of Science and citations in Table S1 from Wu (2019) using the litsearchr R package (Grames et al., 2019).

# Import search results downloaded from Github
naive_results <- litsearchr::import_results(file = "/Users/nicholaswu/Library/CloudStorage/OneDrive-WesternSydneyUniversity/Chytrid meta-analysis/Initial search/Wu_table_ref.txt")
naive_results_2 <- litsearchr::import_results(file = "/Users/nicholaswu/Library/CloudStorage/OneDrive-WesternSydneyUniversity/Chytrid meta-analysis/Initial search/naive_search_WoS_211028.txt")

# Extract terms from titles
terms_title <- litsearchr::extract_terms(text = naive_results[, "title"], method = "fakerake",
    min_freq = 4, max_n = 2)

terms_title2 <- litsearchr::extract_terms(text = naive_results_2[, "title"], method = "fakerake",
    min_freq = 2, max_n = 2)

# Combine search terms and create document-feature matrix (DFM)
terms <- unique(c(terms_title, terms_title2))  # Combine both term files
docs <- unique(c(naive_results[, "title"], naive_results_2[, "title"]))  # extract titles
dfm <- litsearchr::create_dfm(elements = docs, features = terms)  # create dfm
network <- litsearchr::create_network(dfm, min_studies = 5)

# Pruning
strengths <- igraph::strength(network)
term_strengths <- data.frame(term = names(strengths), strength = strengths, row.names = NULL) %>%
    dplyr::mutate(rank = rank(strength, ties.method = "min")) %>%
    dplyr::arrange(strength)

Fig. S1. Output from the keyword co-occurrence network analysis. (a) Visualisation of the network where each node represents a term, and the line thickness represents the weight of linking nodes, i.e., the terms appearing in more articles together. (b) The strength of each term in ascending order from left to right.

PRISMA flow-diagram

Fig. S2. PRISMA flow diagram for the systematic data-collection process. n = number of papers remaining after each stage of selection. k = number of effect size after processing individual data, and n is number of observations/records after processing data. Dashed boxes indicate Boolean search string for Web of Science and Scopus.

---

Dataset

Calculate effect size, \(Z_r\), and sampling variance, \(v(Z_r)\)

The raw data for analysis is available on GitHub. This section is the workflow to import and clean the raw data for analysis. Data were separated by raw data, correlation data, and survival data. The correlation coefficient and sample size was calculated from the raw data and merged with the correlation data (es_all_dat).

Effect size as the Fisher z-transformation of the correlation coefficient r, \(Z_r\) (4), and the sampling variance, \(v(Z_r)\) (5), was calculated following Hedges and Olkin (2014):

\[\begin{equation} Z_r = \frac{1}{2}ln \left( \frac{1+r}{1-r}\right), \tag{4} \end{equation}\]

\[\begin{equation} v(Z_r) = \frac{1}{n - 3}, \tag{5} \end{equation}\]

where r is the correlation coefficient, and n is the sampling size. The standard error, \(se_i\) (6), was calculated as:

\[\begin{equation} se_{i} = \frac{1}{\sqrt{n - 3}}. \tag{6} \end{equation}\]

The inverse of \(se_i\) or precision (7) was calculated as:

\[\begin{equation} p_{i} = \frac{1}{se_i}. \tag{7} \end{equation}\]

The inverse of sampling variance or weight (8) for calculating heterogeneity:

\[\begin{equation} w_{i} = \frac{1}{v(Z_r)}. \tag{8} \end{equation}\]

All conversions for the standard error (\(se_i\)), sampling variance, precision (the inverse of \(se_i\)) and weight (the inverse of variance) followed Nakagawa et al. (2022).

# Load and clean raw data
trait_dat <- read.csv("https://raw.githubusercontent.com/nicholaswunz/chytrid-meta-analysis/main/data/trait_raw_data.csv")  %>% 
  dplyr::filter(trait != "Survival") %>% 
  dplyr::mutate(species_OTL = as.factor(species_OTL),
                trait_value = as.numeric(trait_value),
                trait       = as.factor(trait),
                PCR_meth    = as.factor(PCR_meth),
                response    = as.factor(response),
                lnDose      = log(dose_ZE_ml + 1),
                lnBd        = log(Bd + 1),
                temp_K      = trt_temp + 273.15, # convert temperature Celcius to Kelvin
                inv_temp    = 1 / 8.62e-5  * (1 / mean(temp_K, na.rm = TRUE) - 1 / temp_K) # standardise to mean temp, Boltzmann constant as eV/K
  ) %>%
  dplyr::select(-one_of(c("notes", "title", "ref")))

# Clean effect size data
es_dat <- read.csv("https://raw.githubusercontent.com/nicholaswunz/chytrid-meta-analysis/main/data/trait_corr_data.csv")  %>% 
  dplyr::mutate(study_ID    = ifelse(notes == "repeated", study_ID, study_ID + max(trait_dat$study_ID)),
                species_OTL = as.factor(species_OTL),
                lnDose      = log(dose_ZE_ml + 1),
                lnBd        = log(Bd + 1),
                resistance  = as.factor(resistance),
                trait       = as.factor(trait),
                response    = as.factor(response),
                Zr          = (1 / 2) * log((1 + corr_coeff) / (1 - corr_coeff)), # Fisher-transformed correlation coefficient
                Zr_v        = 1 / (sample_size - 3), # sampling variance (v)
                Zr_sei      = 1 / sqrt(sample_size - 3), # standard error (SE)
                Zr_inv      = 1 / Zr_sei, # precision (inverse of SE) 
                Zr_w        = 1 / Zr_v) %>% # weight (inverse of Zr_v))
  dplyr::select(-one_of(c("notes", "title", "ref")))

# Clean survival data
surv_dat <- read.csv("https://raw.githubusercontent.com/nicholaswunz/chytrid-meta-analysis/main/data/trait_raw_data.csv")  %>% 
  dplyr::filter(trait == "Survival") %>% 
  dplyr::mutate(species_OTL = as.factor(species_OTL),
                resistance  = as.factor(resistance),
                life_stage  = as.factor(life_stage),
                lnDose      = log(dose_ZE_ml + 1),
                lnBd        = log(Bd + 1),
                lnMass      = log(mass_g),
                lnTime      = log(time_post_days),
                temp_K      = trt_temp + 273.15, # convert temperature Celcius to Kelvin
                inv_temp    = 1 / 8.62e-5  * (1 / mean(temp_K, na.rm = TRUE) - 1 / temp_K), # standardise to mean temp, Boltzmann constant as eV/K
                alive       = ifelse(trait_value == "Alive", 1, 0)) %>%
  tidyr::unite("species_lifestage", c("species", "life_stage"), sep = "_", remove = FALSE) %>%
  dplyr::select(-one_of(c("notes", "title", "ref"))) 

# For raw data (trait_dat), calculate correlation coefficient, then calculate effect size
effect_size_dat <- as.data.frame(trait_dat %>%
                                    dplyr::group_by(study_ID, bio_hier, trait, species, response) %>%
                                    dplyr::summarise(corr_coeff     = cor(Bd, trait_value, method = "pearson"),
                                                     sample_size    = length(trait_value),
                                                     author_first   = unique(author_first),
                                                     year           = median(year),
                                                     extracted      = unique(extracted),
                                                     order          = unique(order),
                                                     family         = unique(family),
                                                     species        = unique(species),
                                                     species_OTL    = unique(species_OTL),
                                                     origin         = unique(origin),
                                                     setting        = unique(setting),
                                                     life_stage     = unique(life_stage),
                                                     strain         = unique(strain),
                                                     mean_mass_g    = mean(mass_g),
                                                     trt_temp       = mean(trt_temp),
                                                     dose_ZE_ml     = mean(dose_ZE_ml[dose_ZE_ml != 0]),
                                                     Bd             = mean(Bd[Bd != 0]),
                                                     time_post_days = median(time_post_days[Bd != 0]),
                                                     lnBd           = mean(lnBd[lnBd != 0]))) %>%
  dplyr::mutate(Zr     = (1 / 2) * log((1 + corr_coeff) / (1 - corr_coeff)), # Fisher-transformed correlation coefficient
                Zr_v   = 1 / (sample_size - 3), # sampling variance (v)
                Zr_sei = 1 / sqrt(sample_size - 3), # standard error (SE)
                Zr_inv = 1 / Zr_sei, # precision (inverse of SE) 0
                Zr_w   = 1 / Zr_v, # weight (inverse of Zr_v)
                Zr     = ifelse(response %in% c("Evaporative water loss", "AQP activity", "Skin permeability", 
                                                'Cutaneous ion loss', 'Muscle water content', "Caspase activity"), -abs(Zr), Zr)) 

# Combine effect size dataset
es_all_dat <- bind_rows(effect_size_dat, es_dat) %>%
  tibble::rowid_to_column("es_ID") %>% # add effect size ID
  dplyr::mutate(year_centre = year - mean(year)) # mean-centring year of publication)

Data summary

The es_all_dat used for the meta-analysis comprises of 307 individual effect sizes from 69 studies across 52 species (detailed in Supplementary Table S1).

Table S1 - Summary

Table S1 - Responses. Trait categories associated with Bd infection and the specific responses within these traits. Number of effect sizes, studies and species are shown. Responses with asterisks were corrected for direction. AQP = aquaporins, ENaC = epithelial sodium channel, NKA = sodium potassium pump, CT_max = critical thermal maxima, CT_min = critical thermal minima.

Categorised trait	Specific reponse	Effect size (k)	Studies (n)	Species (n)
Behaviour	Activity	2	2	2
Behaviour	Feeding rate	8	1	6
Behaviour	Sloughing	7	3	5
Behaviour	Breathing rate	2	1	1
Behaviour	Call duration	1	1	1
Behaviour	Call rate	3	3	3
Behaviour	Calling effort	1	1	1
Behaviour	Dominant frequency	1	1	1
Behaviour	Intercall interval	1	1	1
Behaviour	Note duration	1	1	1
Behaviour	Pulse rate	1	1	1
Behaviour	Water conserving posture	1	1	1
Body condition	Body condition (residual or index)	9	8	8
Body condition	Body size (body mass or snout-vent-length)	14	6	12
Body condition	Mass change	24	19	18
Body condition	Body mass index	4	2	4
Body condition	Fat mass	2	1	1
Cardiovascular	Contraction force	2	1	2
Cardiovascular	Haematocrit	1	1	1
Cardiovascular	Heart rate	2	1	2
Cardiovascular	Relative ventricular mass	2	1	2
Cardiovascular	Blood CO2	1	1	1
Cardiovascular	Haemoglobin	3	2	2
Cardiovascular	pH	1	1	1
Cardiovascular	RBC	3	2	2
Energy metabolism	RMR	8	4	3
Hormone level	Corticosterone	7	5	7
Hormone level	Testosterone	1	1	1
Immune	Bd-antibodies	2	1	2
Immune	Brevinin-1Sa	1	1	1
Immune	Lymphocyte/neutrophil ratio	1	1	1
Immune	Lymphocytes	15	6	6
Immune	Plasma protein	3	2	3
Immune	Serotonin	1	1	1
Immune	Spleen size	1	1	1
Immune	Total peptides	1	1	1
Immune	Bacterial killing ability	8	1	2
Immune	Basophils	4	3	3
Immune	Bone marrow	2	1	1
Immune	Eosinophils	4	3	4
Immune	Hematopoietic density	3	1	2
Immune	Hematopoietic tissue	6	1	3
Immune	Hepatosomatic index	1	1	1
Immune	Liver ratio	2	1	2
Immune	Monocytes	4	3	4
Immune	Neutrophil/lymphocyte ratio	8	1	2
Immune	Neutrophils	14	5	6
Immune	Peptide profiles	4	1	1
Immune	Spleen cells	3	1	3
Immune	Spleen ratio	2	1	2
Immune	Splenic cell viability	2	1	2
Immune	Splenosomatic index	1	1	1
Immune	Swelling response	3	1	1
Immune	White blood cells	6	3	6
Osmoregulation	Active sodium transport	2	2	1
Osmoregulation	AQP activity	1	1	1
Osmoregulation	Cutaneous ion loss*	1	1	1
Osmoregulation	ENaC abundance	1	1	1
Osmoregulation	Evaporative water loss*	9	3	9
Osmoregulation	Muscle water content	1	1	1
Osmoregulation	NKA abundance	1	1	1
Osmoregulation	Plasma chloride	4	4	2
Osmoregulation	Plasma osmolarity	3	3	2
Osmoregulation	Plasma potassium	6	5	4
Osmoregulation	Plasma sodium	7	6	4
Osmoregulation	Water uptake	6	4	5
Osmoregulation	Plasma calcium	3	2	2
Osmoregulation	Plasma glucose	4	3	2
Reproduction	Clutch size	1	1	1
Reproduction	Colour score	1	1	1
Reproduction	Egg size	1	1	1
Reproduction	Ovary area	1	1	1
Reproduction	Sperm count	1	1	1
Reproduction	Testis area	1	1	1
Reproduction	Tubule area	1	1	1
Reproduction	Egg numbers	1	1	1
Reproduction	Ovary mass	2	2	2
Reproduction	Oviduct size	1	1	1
Reproduction	Sexual display	1	1	1
Reproduction	Spermatic cyst diameter	1	1	1
Reproduction	Spermatocytes	2	1	2
Reproduction	Spermatogenesis-stage cycts	2	1	2
Reproduction	Spermatogonia	2	1	2
Reproduction	Testes size	4	3	4
Reproduction	Testes tubules	2	1	2
Skin integrity	Caspase activity*	1	1	1
Skin integrity	Skin permeability*	1	1	1
Skin integrity	Skin resistance	5	3	4
Skin integrity	Skin sheds*	2	1	1
Thermoregulation	Body temp	4	4	4
Thermoregulation	CTmax	1	1	1
Thermoregulation	CTmin	2	1	2
Locomotor capacity	Jumping acceleration	1	1	1
Locomotor capacity	Jumping velocity	1	1	1

Life stage

Table S1 - Life stage. Continuation of data summary with additional information on the number of effect size, studies, and species in the database by life stage (adult or Juvenile).

Life stage	Effect size (k)	Studies (n)	Species (n)
Adult	206	48	38
Juvenile	101	25	22

Species

Table S1 - Species. Continuation of data summary with additional information on the number of effect size, and studies in the database by species.

Species	Effect size (k)	Studies (n)
Litoria caerulea	58	13
Litoria verreauxii alpina	24	6
Rana cascadae	19	2
Pseudacris regilla	17	1
Lithobates pipiens	15	3
Litoria infrafrenata	15	1
Litoria aurea	13	3
Pseudophryne corroboree	12	2
Litoria rheocola	10	1
Anaxyrus boreas	8	3
Rana muscosa	8	2
Lissotriton helveticus	6	2
Litoria raniformis	6	1
Litoria serrata	6	2
Rhinella marina	6	3
Litoria chloris	5	1
Pseudacris triseriata	5	3
Atelopus zeteki	4	3
Xenopus laevis	4	2
Anaxyrus americanus	3	3
Bufo boreas	3	2
Dendropsophus�minutus	3	1
Hyla versicolor	3	2
Ischnocnema parva	3	1
Litoria wilcoxii	3	1
Physalaemus albonotatus	3	1
Pseudacris feriarum	3	1
Rana catesbeiana	3	2
Anaxyrus terrestris	2	1
Brachycephalus pitanga	2	1
Hyla japonica	2	1
Lechriodus fletcheri	2	2
Limnodynastes peronii	2	2
Limnodynastes tasmaniensis	2	2
Lithobates sylvaticus	2	1
Lithobates yavapaiensis	2	1
Osteopilus septentrionalis	2	1
Platyplectrum ornatum	2	2
Rana pipiens	2	1
Rana sphenocephala	2	1
Acris crepitans	1	1
Bufo bufo	1	1
Crinia signifera	1	1
Hyla chrysoscelis	1	1
Hyla cinerea	1	1
Leiopelma hochstetteri	1	1
Leiopelma pakeka	1	1
Lithobates palustris	1	1
Litoria citropa	1	1
Litoria ewingii	1	1
Litoria lesueurii	1	1
Pseudophryne pengilleyi	1	1
Rana boylii	1	1
Rana clamitans	1	1
Rana sylvatica	1	1

Traits

Table S1 - Traits Continuation of data summary with additional information on the number of effect size, studies, and species in the database traits.

Trait	Effect size (k)	Studies (k)	Species (n)
Immune	102	15	16
Body condition	53	34	34
Osmoregulation	49	12	13
Behaviour	29	10	16
Reproduction	25	6	6
Cardiovascular	15	4	5
Skin integrity	9	6	6
Energy metabolism	8	4	3
Hormone level	8	5	7
Thermoregulation	7	6	7
Locomotor capacity	2	1	1

---

Prepare phylogeny for analysis

This section provides the workflow to extract the phylogenetic tree from the Open Tree of Life (OTL; see phylogenetic reconstruction in main text), match the OTL names with the data set, and create a correlation matrix for subsequent analysis. The phylogenetic tree produced in R for the Supplementary Figure S3 was modified in Adobe Illustrator for clarity.

species_comb <- rbind(es_all_dat %>%
    dplyr::select(order, family, species, species_OTL), surv_dat %>%
    dplyr::select(order, family, species, species_OTL))  # 'surv_data_clean' created below
species_all <- sort(unique(as.character(species_comb$species_OTL)))  # generate list of species (as character format)
taxa <- rotl::tnrs_match_names(names = species_all)  # match taxonomic names to the OTL

# check if species list match OT identifier taxa[taxa$approximate_match ==
# TRUE,] # none so far

# retrieving phylogenetic relationships among taxa in the form of a trimmed
# sub-tree
tree <- rotl::tol_induced_subtree(ott_ids = rotl::ott_id(taxa), label_format = "name")

# Compute branch lengths
tree <- ape::compute.brlen(tree, method = "Grafen", power = 1)
tree <- ape::multi2di(tree, random = TRUE)  # use a randomization approach to deal with polytomies

# Check tree is ultrametric ape::is.ultrametric(tree) # TRUE

# Create correlation matrix for analysis
phylo_cor <- ape::vcv(tree, cor = T)

Fig. S3. Phylogenetic reconstruction of 60 amphibians used in the study from the Open Tree of Life. The phylogeny was converted to a correlation matrix for the meta-analysis.

---

Meta-analysis

Phylogenetically controlled, multi-level meta-analysis and meta-regression

Run analysis to first look at the null effect of Bd infection on all responses (\(Z_r\)) without moderators (null_model), then with the following moderators days post-exposure, life stage (adults or juveniles), origin (wild caught or lab raised), natural logarithm of inoculation dose, and treatment temperature (overall_model). A follow up analysis was performed to estimate of coefficients for every factor level to obtain posterior distributions of effect size between biological and functional traits (trait_model).

meta_prior <- c(set_prior("cauchy(0, 1)", class = "sd"))  # Cauchy on tau (random effect variance), normal on fixed effect

# Null model - Without moderators
set.seed(10)
null_model <- brms::brm(Zr | se(Zr_sei) ~ 1 + (1 | es_ID) + (1 | study_ID) + (1 |
    species_OTL) + (1 | gr(species, cov = phylo)), data = es_all_dat, family = gaussian,
    data2 = list(phylo = phylo_cor), prior = meta_prior, iter = 10000, warmup = 5000,
    cores = 4, chains = 4, control = list(adapt_delta = 0.999, max_treedepth = 18))

# Meta-regression model - With moderators
overall_model <- brms::brm(Zr | se(Zr_sei) ~ time_post_days + life_stage + origin +
    lnDose + trt_temp + (1 | es_ID) + (1 | study_ID) + (1 | species_OTL) + (1 | gr(species,
    cov = phylo)), data = es_all_dat, family = gaussian, data2 = list(phylo = phylo_cor),
    prior = meta_prior, iter = 10000, warmup = 5000, cores = 4, chains = 4, control = list(adapt_delta = 0.999,
        max_treedepth = 18))

# Trait-specific meta-regression model
trait_model <- brms::brm(Zr | se(Zr_sei) ~ -1 + trait + life_stage + year_centre +
    (1 | es_ID) + (1 | study_ID) + (1 | species_OTL) + (1 | gr(species, cov = phylo)),
    data = es_all_dat %>%
        dplyr::group_by(trait) %>%
        dplyr::filter(n() >= 5), family = gaussian, data2 = list(phylo = phylo_cor),
    prior = meta_prior, iter = 10000, warmup = 5000, cores = 4, chains = 4, control = list(adapt_delta = 0.999,
        max_treedepth = 18))

Publication bias test

Publication bias was tested by using the same formula as the null_model with the inclusion of the mean centring of year (test for time-lag bias) and the square root of the sampling variance (test for samm-sample bias) as moderators (bias_model).

bias_model <- brms::brm(Zr | se(Zr_sei) ~ 1 + year_centre + sqrt(Zr_v) + (1 | es_ID) +
    (1 | study_ID) + (1 | species_OTL) + (1 | gr(species, cov = phylo)), data = es_all_dat,
    family = gaussian, data2 = list(phylo = phylo_cor), prior = meta_prior, iter = 10000,
    warmup = 5000, cores = 4, chains = 4, control = list(adapt_delta = 0.999, max_treedepth = 18))

Fig. S4. Scatterplots of the observed data (y) vs the average simulated data (y_rep) from the posterior predictive distribution for the (a) null model, (b) the overall model, (c) the trait model, and (c) the bias model. Dashed line represents a slope of 1.

Model output

Table S2 - Null model

Table S2. Mean parameter estimates, estimate error, and 95% Bayesian credible intervals for the null model, which includes the intercept (\(\beta_0\)). Group-level effects include the standard deviations (\(\sigma\)) for individual-level observations (\(\sigma_{residual}^2\)), study-level observations (\(\sigma_{study}^2\)), species identity (\(\sigma_{species}^2\)), and phylogenetic relatedness (\(\sigma_{phylogeny}^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	-0.13	0.08	-0.28	0.05
Group-level effects
\(\sigma_{residual}^2\)	0.32	0.02	0.27	0.36
\(\sigma_{phylogeny}^2\)	0.09	0.08	0.00	0.32
\(\sigma_{species}^2\)	0.16	0.05	0.06	0.26
\(\sigma_{study}^2\)	0.19	0.05	0.09	0.29

Table S3 - Heterogeneity and variance

Table S3. Heterogeneity and \(R^2\) of the null and overall model. k = number of estimates, and \(I^2\) = heterogeneity. \(R_{marginal}^2\) represents the variance explained by moderators, while \(R_{conditional}^2\) represents the variance explained by both moderators and group-level effects. Estimates shown correspond to modes and 95% highest posterior density intervals.

Model	\(I_{residual}^2\)	\(I_{study}^2\)	\(I_{species}^2\)	\(I_{phylogeny}^2\)	\(I_{total}^2\)	\(R_{marginal}^2\)	\(R_{conditional}^2\)
Null	41.6	25.2	22	1.8	97.3	0	72.2
	[29.3–52.7]	[13.3–35.7]	[8.3–31.8]	[0–29.5]	[96.4–97.8]	[0–0]	[66.3–77.8]
Overall	31.2	28.8	22.5	9.1	98.1	12.81	78.3
	[19.6–41.7]	[11.8–45.5]	[1.3–38.9]	[0–38.5]	[97.3–98.6]	[2.96–24.89]	[72.3–83.8]

Table S4 - Bias model

Table S4. Mean parameter estimates, estimate error, and 95% Bayesian credible intervals for the bias model, which includes the intercept (\(\beta_0\)), publication year and sampling standard error as covariates to test and correct for publication bias. Group-level effects include the standard deviations (\(\sigma\)) for individual-level observations (\(\sigma_{residual}^2\)), study-level observations (\(\sigma_{study}^2\)), species identity (\(\sigma_{species}^2\)), and phylogenetic relatedness (\(\sigma_{phylogeny}^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	-0.06	0.12	-0.29	0.17
Publication year	0.01	0.01	-0.01	0.03
Sampling standard error	-0.37	0.42	-1.19	0.46
Group-level effects
\(\sigma_{residual}^2\)	0.31	0.02	0.27	0.36
\(\sigma_{phylogeny}^2\)	0.10	0.09	0.00	0.33
\(\sigma_{species}^2\)	0.15	0.05	0.04	0.25
\(\sigma_{study}^2\)	0.19	0.05	0.09	0.30

Fig S5 - Funnel plot

Fig. S5. Funnel plot on the effect size against the study standard error. The white area bordered by dashed lines represents the region of 95% pseudo confidence intervals where 95% of studies are expected to fall in the absence of bias and heterogeneity.

Table S5 - Overall model

Table S5. Mean parameter estimates, estimate error, and 95% Bayesian credible intervals for the overall model, which includes the intercept (\(\beta_0\)), days post-exposure, life stage (juvenile or adult), origin (wild caught or lab-reared), the natural logarithm of inoculation dose (lnDose), and treatment temperature as fixed effects. Group-level effects include the standard deviations (\(\sigma\)) for individual-level observations (\(\sigma_{residual}^2\)), study-level observations (\(\sigma_{study}^2\)), species identity (\(\sigma_{species}^2\)), and phylogenetic relatedness (\(\sigma_{phylogeny}^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	-1.19	0.80	-2.78	0.38
Days post-exposure	0.00	0.00	0.00	0.00
Life stage - Juveniles	0.08	0.17	-0.24	0.42
Origin - Wild	-0.04	0.26	-0.55	0.46
lnDose	0.01	0.06	-0.10	0.12
Treatment temperature	0.05	0.02	0.01	0.09
Group-level effects
\(\sigma_{residual}^2\)	0.33	0.03	0.28	0.40
\(\sigma_{phylogeny}^2\)	0.20	0.16	0.01	0.61
\(\sigma_{species}^2\)	0.24	0.11	0.03	0.48
\(\sigma_{study}^2\)	0.32	0.11	0.12	0.56

Table S6 - Trait model

Table S6. Mean parameter estimates, estimate error, and 95% Bayesian credible intervals for the trait model, which includes the biological and functional traits, life stage (juvenile or adult), and publication year as fixed effects. Group-level effects include the standard deviations (\(\sigma\)) for individual-level observations (\(\sigma_{residual}^2\)), study-level observations (\(\sigma_{study}^2\)), species identity (\(\sigma_{species}^2\)), and phylogenetic relatedness (\(\sigma_{phylogeny}^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
Behaviour	-0.18	0.13	-0.42	0.07
Body condition	-0.21	0.11	-0.41	0.03
Cardiovascular	0.04	0.17	-0.28	0.37
Energy metabolism	0.58	0.19	0.21	0.97
Hormon elevel	0.05	0.19	-0.32	0.43
Immune	-0.15	0.12	-0.37	0.09
Osmoregulation	-0.21	0.12	-0.43	0.04
Reproduction	-0.14	0.15	-0.44	0.17
Skin integrity	-0.06	0.17	-0.40	0.29
Thermoregulation	-0.22	0.19	-0.59	0.16
Life stage - Juveniles	0.03	0.09	-0.15	0.22
Publication year	0.01	0.01	-0.01	0.04
Group-level effects
\(\sigma_{residual}^2\)	0.30	0.02	0.25	0.34
\(\sigma_{phylogeny}^2\)	0.11	0.10	0.00	0.35
\(\sigma_{species}^2\)	0.16	0.05	0.04	0.26
\(\sigma_{study}^2\)	0.22	0.06	0.11	0.34

---

Trait sensitivity analysis

Calculating trait sensitivity

The following workflow uses the raw data set to estimate response-specific sensitivity to Bd infection. All response values were standardised to relative change in response.

# Extract mean, 10th and 90th quantile control values by species and responses
control_mean <- as.data.frame(trait_dat %>%
  dplyr::filter(Bd == 0) %>% 
  dplyr::group_by(species, response) %>% 
  dplyr::summarise(mean            = mean(trait_value, na.rm = TRUE),
                   quantile_10     = quantile(trait_value, .10),
                   quantile_90     = quantile(trait_value, .90)) %>%
  tidyr::unite("specific_response", c(species, response), sep = "_",  remove = FALSE)) # unique identifier by species and response

# Combine mean controls and calculate relative change of response
traits_exposed <- trait_dat %>% 
  tidyr::unite("specific_response", c(species, response), sep = "_",  remove = FALSE) %>%
  dplyr::mutate(control           = control_mean$mean[match(specific_response, control_mean$specific_response)], 
                change            = ifelse(trait_unit == "relative", trait_value, (trait_value - control) / control), # calculate relative change
                specific_response = as.factor(specific_response),
                quantile_90       = control_mean$quantile_90[match(specific_response, control_mean$specific_response)],
                quantile_10       = control_mean$quantile_10[match(specific_response, control_mean$specific_response)],
                control_90        = ifelse(trait_unit == "relative", trait_value, (quantile_90 - control) / control), # calculate relative 90th quantile
                control_10        = ifelse(trait_unit == "relative", trait_value, (quantile_10 - control) / control)) # calculate relative 10th quantile

# Run lm models to find ones with good R2 
R_sq_model <- as.data.frame(traits_exposed %>% 
                              dplyr::filter(!is.na(change)) %>% # remove NA's
                              dplyr::group_by(specific_response) %>%
                              dplyr::do(model = broom::glance(lm(change ~ lnBd, data = .))) %>%
                              tidyr::unnest(model)) %>%
  dplyr::select(specific_response, r.squared, adj.r.squared, p.value)

# extract equation (intercept and slope)
lm_eq <- as.data.frame(traits_exposed %>% 
                         dplyr::filter(!is.na(change)) %>%
                         dplyr::group_by(specific_response) %>%
                         dplyr::do(model = broom::tidy(lm(change ~ lnBd, data = .))) %>% 
                         tidyr::unnest(model)) %>%
  dplyr::select(specific_response, term, estimate) %>%
  reshape(idvar = "specific_response", timevar = "term", direction = "wide") %>% # from long to wide for estimate column
  dplyr::rename(intercept = "estimate.(Intercept)",
                slope     = "estimate.lnBd")

# Deal with relative units that were not converted previously (ifelse(trait_unit == "relative", trait_value, (quantile_90 - control) / control))
control_range <- as.data.frame(traits_exposed %>%
                                 dplyr::filter(Bd == 0) %>%
                                 dplyr::group_by(specific_response) %>%
                                 dplyr::summarise(control_90 = max(control_90),
                                                  control_10 = min(control_10)))

# Generate Bd load sequence from min and max and predict response
xseq       <- seq(0, max(traits_exposed$lnBd[!is.na(traits_exposed$lnBd)]), length.out = 50)
eq_bd_pred <- data.frame(lm_eq %>%
                           dplyr::group_by(specific_response) %>%
                           dplyr::summarise(pred = intercept + slope * xseq) %>%
                           dplyr::mutate(pred_bd    = xseq,
                                         trait      = traits_exposed$trait[match(specific_response, traits_exposed$specific_response)],
                                         control    = traits_exposed$control[match(specific_response, traits_exposed$specific_response)],
                                         control_90 = control_range$control_90[match(specific_response, control_range$specific_response)],
                                         control_10 = control_range$control_10[match(specific_response, control_range$specific_response)]))

# Keep specific responses with R2 >= 0.1 and p value <= 0.05
responses_kept <- R_sq_model %>% dplyr::filter(adj.r.squared >= 0.1 | p.value <= 0.05) %>% droplevels()
kept_pred      <- eq_bd_pred %>% dplyr::filter(specific_response %in% c(levels(responses_kept$specific_response))) %>% droplevels()

# Filter Bd load lower than 10th and higher than 90th quantile control values 
sensitivity_pos <- data.frame(kept_pred %>%
                                dplyr::group_by(specific_response) %>%
                                dplyr::filter(pred_bd != 0 & pred > 0) %>% # in positive direction
                                dplyr::filter(pred > control_90) %>%
                                dplyr::summarize(minBd = min(pred_bd)))

sensitivity_neg <- data.frame(kept_pred %>%
                               dplyr::group_by(specific_response) %>%
                               dplyr::filter(pred_bd != 0 & pred < 0) %>%  # in negative direction
                               dplyr::filter(case_when(control < 0 ~ pred < control_90, T ~ pred < control_10)) %>%
                               dplyr::summarize(minBd = min(pred_bd)))

# Combine sensitivity_pos and sensitivity_neg
sensitivity <- rbind(sensitivity_pos, sensitivity_neg) %>%
  dplyr::mutate(lnminBd     = minBd,
                response    = traits_exposed$response[match(specific_response, traits_exposed$specific_response)],
                trait       = traits_exposed$trait[match(specific_response, traits_exposed$specific_response)],
                bio_hier    = traits_exposed$bio_hier[match(specific_response, traits_exposed$specific_response)],
                species     = traits_exposed$species[match(specific_response, traits_exposed$specific_response)],
                species_OTL = traits_exposed$species_OTL[match(specific_response, traits_exposed$specific_response)],
                life_stage  = traits_exposed$life_stage[match(specific_response, traits_exposed$specific_response)],
                lnDose      = traits_exposed$lnDose[match(specific_response, traits_exposed$specific_response)],
                inv_temp    = traits_exposed$inv_temp[match(specific_response, traits_exposed$specific_response)])

sensitivity$species <- sub(" ", "_", sensitivity$species_OTL)
sensitivity$species <- factor(sensitivity$species, levels = tree_tip_label) # relevel order by tree

Sensitivity analysis

Run phylogenetic corrected, multi-level model to test the sensitivity of traits to Bd load.

set.seed(10)
sensitive_model <- brms::brm(lnminBd ~ -1 + trait + lnDose + inv_temp + (1 | species_OTL) + (1 | gr(species, cov = phylo)),
                             data    = sensitivity,
                             family  = gaussian,
                             data2   = list(phylo = phylo_cor),
                             prior   = c(prior(normal(0, 5), lb = 0), # mean of 0 and SD of 5
                                       prior(student_t(3, 0, 5), "sd"), # class of random effect deviation
                                       prior(student_t(3, 0, 5), "sigma")), # residual SD parameter
                             iter    = 1e4, warmup = 5e3, cores = 4, chains = 4,
                             control = list(adapt_delta = 0.99, max_treedepth = 18))

Fig. S6. Scatterplot of the observed data (y) vs the average simulated data (yrep) from the posterior predictive distribution. Dashed line represents a slope of 1.

Model output

Table S7. Mean parameter estimates, estimate error, and 95% Bayesian credible intervals for the sensitivity model, which includes the functional traits, natural logarithm of inoculation dose (lnDose), and the inverse of temperature as fixed effects. Group-level effects include the standard deviations (\(\sigma\)) for species identity (\(\sigma_{species}^2\)), and phylogenetic relatedness (\(\sigma_{phylogeny}^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
Behaviour	4.87	2.34	0.63	9.67
Body condition	5.06	1.95	1.19	8.80
Cardiovascular	3.79	2.38	0.23	9.14
Energy metabolism	4.23	2.67	0.26	10.11
Hormone level	2.51	1.86	0.11	6.87
Immune	2.82	2.20	0.10	8.15
Osmoregulation	2.52	1.47	0.19	5.71
Reproduction	10.60	3.26	4.06	16.93
Skin integrity	2.46	1.70	0.11	6.43
lnDose	0.08	0.07	0.00	0.26
Inverse temperature	1.34	1.23	0.03	4.53
Group-level effects
\(\sigma_{phylogeny}^2\)	2.05	1.52	0.08	5.66
\(\sigma_{species}^2\)	4.08	1.74	0.57	7.49

---

Survival analysis

The survival data set comprises of 33 species.

# Rename species based on species_OTL
surv_dat$species <- sub(" ", "_", surv_dat$species_OTL)
surv_dat$species <- factor(surv_dat$species, levels = tree_tip_label)  # relevel order by tree
surv_dat <- surv_dat %>%
    dplyr::mutate(species = factor(species))

Survival and Bd and time-dependent mortality analysis

Run phylogenetic corrected, multi-level models to test for the influence of Bd load on survival.

# Binary logistic regression model
set.seed(10)
surv_model <- brms::brm(alive ~ 1 + lnBd + time_post_days + lnDose + life_stage +
    inv_temp + species_OTL + (1 | study_ID:species_OTL) + (1 | gr(species, cov = phylo)),
    data = surv_dat, family = bernoulli(link = "logit"), data2 = list(phylo = phylo_cor),
    prior = c(prior(normal(0, 2), class = Intercept), prior(normal(0, 2), class = b)),
    iter = 5000, warmup = 2000, chains = 4, cores = 4, control = list(adapt_delta = 0.99,
        max_treedepth = 18))

time_bd_model <- brms::brm(lnBd ~ lnTime * life_stage + lnDose + inv_temp + (1 |
    study_ID:species_OTL), data = surv_dat %>%
    dplyr::filter(alive == 0), family = gaussian(), prior = c(prior(normal(0, 2),
    class = Intercept), prior(normal(0, 2), class = b)), iter = 5000, warmup = 2000,
    chains = 4, cores = 4, control = list(adapt_delta = 0.99, max_treedepth = 18))

Fig. S7. Posterior predictive check for (a) the survival model, and (b) the Bd-time model. The Bd-time model is represented by a scatterplot of the observed data (y) vs the average simulated data (yrep) from the posterior predictive distribution. Dashed line represents a slope of 1.

Model output

Table S8 - Survival model

Table S8. Mean parameter estimates, estimate error, and 95% Bayesian credible intervals for the probability of survival, which includes the natural logarithm of Bd load (lnBd), days post exposure, natural logarithm of inoculation dose (lnDose), life stage, and the inverse of temperature as fixed effects. Group-level effects include the standard deviations (\(\sigma\)) for phylogenetic relatedness (\(\sigma_{phylogeny}^2\)) and species-study interaction (\(\sigma_{species:study}^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	17.16	3.29	11.58	24.59
lnBd	-0.96	0.05	-1.06	-0.85
Days post-exposure	-0.03	0.01	-0.05	-0.02
lnDose	-0.82	0.27	-1.44	-0.41
Life stage - Juvenile	2.15	0.48	1.23	3.11
Inverse temperature	0.94	0.60	-0.20	2.17
Group-level effects
\(\sigma_{phylogeny}^2\)	2.39	2.05	0.09	7.56
\(\sigma_{species:study}^2\)	7.41	1.42	5.08	10.61

Table S9 - Threshold

Table S9. Predicted 50% and 90% survival threshold by species presented as zoospore equivalent (ZE) and the natural logarithm of ZE (ln ZE + 1).

Species	Family	50% (ln ZE + 1)	50% (ZE + 1)	90% (ln ZE + 1)	90% (ZE + 1)
Acris crepitans	Hylidae	8.1	3134	10.3	30486
Anaxyrus americanus	Bufonidae	7.8	2412	10.1	25431
Anaxyrus boreas	Bufonidae	7.8	2412	10.1	24788
Anaxyrus terrestris	Bufonidae	7.5	1886	9.9	20875
Atelopus zeteki	Bufonidae	8.7	6142	11.1	63494
Dendropsophus_elegans	Hylidae	8.6	5415	10.9	55199
Dendropsophus_minutus	Hylidae	8.7	5908	11.0	59648
Dryophytes versicolor	Hylidae	7.0	1055	9.4	11672
Ischnocnema parva	Brachycephalidae	8.7	5964	11.0	60506
Lechriodus fletcheri	Limnodynastidae	8.6	5274	10.9	51985
Limnodynastes ornatus	Limnodynastidae	8.3	4073	10.6	41840
Limnodynastes peronii	Limnodynastidae	8.1	3167	10.3	31350
Limnodynastes tasmaniensis	Limnodynastidae	8.4	4476	10.7	46123
Litoria aurea	Pelodryadidae	7.8	2458	10.2	25966
Litoria caerulea	Pelodryadidae	8.3	4024	10.6	40368
Litoria ewingii	Pelodryadidae	8.1	3167	10.4	32290
Litoria spenceri	Pelodryadidae	8.2	3545	10.5	36467
Litoria verreauxii	Pelodryadidae	7.5	1724	9.7	15618
Pseudacris feriarum	Hylidae	7.2	1341	9.5	13080
Pseudacris regilla	Hylidae	7.6	2077	10.1	24100
Pseudacris triseriata	Hylidae	7.2	1386	9.5	13506
Pseudophryne corroboree	Myobatrachidae	8.1	3134	10.3	30486
Rana aurora	Ranidae	7.5	1724	9.7	15618
Rana muscosa	Ranidae	8.7	6068	11.0	62587
Rana pipiens	Ranidae	7.2	1310	9.5	12883
Rana sierrae	Ranidae	8.7	6240	11.1	66762
Rana sphenocephala	Ranidae	8.1	3429	10.5	35145
Rana sylvatica	Ranidae	7.2	1379	9.5	13506
Average		8.0	3369	10.3	34349

Table S10 - Bd-time model

Table S10. Mean parameter estimates, estimate error, and 95% Bayesian credible intervals for the Bd-time model, which includes the natural logarithm of days post-exposure (lnTime), life stage, resistance, the inverse of temperature and the interaction of lnTime and life stage as fixed effects. Group-level effects include the species-study interaction (\(\sigma_{species:study}^2\)).

Parameter	Estimate	Est.Error	Q2.5	Q97.5
Fixed effects
\(\beta_0\)	7.24	3.96	-1.08	14.68
lnTime	-0.22	0.26	-0.74	0.27
Life stage - Juvenile	3.45	0.84	1.76	5.10
lnDose	0.26	0.34	-0.39	0.96
Inverse temperature	0.25	0.84	-1.39	1.92
lnTime : Life stage	-1.90	0.36	-2.60	-1.19
Group-level effects
\(\sigma_{species:study}^2\)	3.61	0.67	2.55	5.15

---

Figures for main text

Figures produced here were modified in Adobe Illustrator for publication.

Figure 2 - Trait model

# Extract marginal effects from model
model_me <- brms::conditional_effects(trait_model, c("trait", "life_stage"))
trait_me <- as.data.frame(model_me[[1]]) %>%
    dplyr::rename(estimate = estimate__, ci.lb = lower__, ci.ub = upper__)

# Plot model output
trait_me %>%
    ggplot(aes(x = trait, y = estimate)) + geom_hline(yintercept = 0, linetype = "dashed",
    alpha = 0.5) + ggforce::geom_sina(data = es_all_dat, aes(x = trait, y = Zr, size = Zr_inv),
    colour = "#cfcfcf") + geom_point(aes(colour = trait), size = 4, show.legend = FALSE) +
    geom_errorbar(aes(ymin = ci.lb, ymax = ci.ub, colour = trait), size = 0.8, width = 0.1,
        show.legend = FALSE) + colorspace::scale_colour_discrete_sequential(palette = "RedOr",
    nmax = 14, order = 4:14) + xlab(NULL) + ylab(expression("Effect size" ~ (italic(Z)[r]))) +
    coord_flip() + mytheme() + theme(legend.position = "bottom")

Fig. 2. Effect of Bd infection on all biological (cardiovascular, immune, hormone level, energy metabolism, osmoregulation, skin integrity) and functional (behaviour, body condition, locomotor capacity, reproduction, thermoregulation) traits measured. Only traits with more than five effect size were analysed. Numbers on the right side of the plot indicate the number of effect sizes analysed. Effect sizes (Zr) presented as estimates ± 95% CI. Grey points represent individual effect size, and the size of each point indicates study precision (inverse of standard error).

Figure 3 - Sensitivity model

# Extract marginal effects
sensitive_me <- brms::conditional_effects(sensitive_model)
sensitive_me <- as.data.frame(sensitive_me[[1]]) %>%
    dplyr::rename(estimate = estimate__, ci.lb = lower__, ci.ub = upper__) %>%
    dplyr::mutate(estimate = exp(estimate), ci.lb = exp(ci.lb), ci.ub = exp(ci.ub))

sensitive_post <- brms::posterior_samples(sensitive_model) %>%
    dplyr::sample_n(4000) %>%
    dplyr::select(b_traitBehaviour:b_traitSkinintegrity) %>%
    reshape2::melt(value.name = "sample")

sensitive_post$variable <- substring(sensitive_post$variable, 8)
source("https://raw.githubusercontent.com/datavizpyr/data/master/half_flat_violinplot.R")

sensitive_post %>%
    dplyr::group_by(variable) %>%
    dplyr::summarise(mean = mean(sample)) %>%
    dplyr::mutate(variable = forcats::fct_reorder(variable, mean)) %>%
    ggplot() + geom_point(aes(x = variable, y = mean), size = 2) + geom_flat_violin(data = sensitive_post,
    aes(x = variable, y = sample, fill = variable), colour = NA, show.legend = FALSE) +
    geom_point(aes(x = variable, y = mean), size = 3) + geom_point(aes(x = variable,
    y = mean, colour = variable), size = 2, show.legend = FALSE) + colorspace::scale_fill_discrete_sequential(palette = "RedOr",
    nmax = 12, order = 4:12) + colorspace::scale_colour_discrete_sequential(palette = "RedOr",
    nmax = 12, order = 4:12) + xlab(NULL) + ylab("Infection intensity (log ZE + 1)") +
    coord_flip() + mytheme()

Fig. 3. Posterior distribution density (4,000 subset iterations) of the minimum Bd load (ln ZE +1) to change trait response above or below the normal response range (5th or 95th quantile). Enclosed circles represent the mean estimate response.

Figure 4 - Survival model

# Plot
surv_plot <- surv_marg_eff %>%
    ggplot(aes(x = x, y = predicted, colour = life_stage)) + geom_vline(xintercept = log(10000),
    linetype = "dashed", alpha = 0.5) + geom_line(aes(group = group, linetype = life_stage),
    size = 0.5) + geom_point(data = surv_dat, aes(x = lnBd, y = alive, colour = life_stage),
    size = 1.5, alpha = 0.5) + colorspace::scale_colour_discrete_sequential(palette = "Blues 3",
    nmax = 5, order = c(3, 5)) + ylab("Survival probabilty") + xlab("Infection intensity (ln ZE + 1)") +
    mytheme() + theme(legend.position = "top") + facet_grid(rows = vars(life_stage))

inflict_plot <- data.frame(surv_marg_eff %>%
    dplyr::group_by(group, life_stage) %>%
    dplyr::filter(predicted < 0.52 & predicted > 0.48) %>%
    dplyr::summarise(threshold = mean(x)) %>%
    dplyr::mutate(bd = exp(threshold))) %>%
    dplyr::mutate(group = forcats::fct_reorder(group, -threshold)) %>%
    ggplot(aes(x = group, y = threshold, colour = life_stage)) + geom_point(aes(shape = life_stage),
    size = 2) + colorspace::scale_colour_discrete_sequential(palette = "Blues 3",
    nmax = 5, order = c(3, 5)) + ylab(expression(atop(italic("Bd") ~ "load (ln ZE + 1)",
    "with 50% mortality probabilty"))) + xlab(NULL) + scale_x_discrete(position = "top") +
    coord_flip() + mytheme() + theme(legend.position = "top", axis.text.y = element_text(size = rel(0.8)))

cowplot::plot_grid(surv_plot, inflict_plot, labels = c("a", "b"), ncol = 2, align = "h",
    axis = "tb", rel_widths = c(0.9, 1))

Fig. 4. Relationship between infection intensity (ln ZE + 1) and the probability of survival of 27 species. (a) Species-specific survival probability were grouped either as adult or juveniles. Solid lines represent adult frogs and dashed lines represent juvenile frogs. Dashed vertical line represents the Vredenburg’s “10,000 zoospore rule” (Kinney et al., 2011). (b) The predicted inflection points for probability of infection that causes 50% mortality. Circles represent adults and triangles represent juveniles.

Figure 5 - Bd-time dependent mortality model

# Extract marginal effects
time_marg_eff <- as.data.frame(ggeffects::ggpredict(time_bd_model, terms = c("lnTime[sample = 50]",
    "life_stage")))

# Create matrix of min and max values per group
time_range <- surv_dat %>%
    dplyr::filter(alive == 0) %>%
    dplyr::group_by(life_stage) %>%
    dplyr::summarise(min = min(lnTime[!is.na(lnTime)]), max = max(lnTime[!is.na(lnTime)])) %>%
    as.data.frame()

# Add min and max values to model df and keep predictions within data range
time_marg_eff$min <- time_range$min[match(time_marg_eff$group, time_range$life_stage)]
time_marg_eff$max <- time_range$max[match(time_marg_eff$group, time_range$life_stage)]
time_marg_eff <- time_marg_eff %>%
    dplyr::group_by(group) %>%
    dplyr::filter(x >= min & x < max) %>%
    dplyr::mutate(life_stage = group, days = exp(x), Bd = exp(predicted))

# Plot
time_bd_plot <- ggplot(data = time_marg_eff, aes(x = x, y = predicted, group = life_stage)) +
    geom_point(data = surv_dat %>%
        dplyr::filter(alive == 0), aes(x = lnTime, y = lnBd, colour = life_stage,
        shape = life_stage), size = 2, alpha = 0.4) + geom_ribbon(aes(x = x, ymin = conf.low,
    ymax = conf.high, fill = life_stage), alpha = 0.1, show.legend = FALSE) + geom_line(aes(colour = life_stage,
    linetype = life_stage), size = 1) + colorspace::scale_colour_discrete_sequential(palette = "Blues 3",
    nmax = 5, order = c(3, 5)) + colorspace::scale_fill_discrete_sequential(palette = "Blues 3",
    nmax = 5, order = c(3, 5)) + ylab("Infection intensity (ln ZE + 1)") + xlab("Time post exposure (ln days)") +
    mytheme() + theme(legend.position = "top")

time_plot <- surv_dat %>%
    dplyr::filter(alive == 0) %>%
    ggplot(aes(x = life_stage, y = lnTime)) + geom_flat_violin(aes(fill = life_stage),
    colour = NA, alpha = 0.5, show.legend = FALSE) + geom_point(aes(shape = life_stage,
    colour = life_stage), size = 1.5) + colorspace::scale_colour_discrete_sequential(palette = "Blues 3",
    nmax = 5, order = c(3, 5)) + colorspace::scale_fill_discrete_sequential(palette = "Blues 3",
    nmax = 5, order = c(3, 5)) + ylab("Time post exposure (ln days)") + xlab(NULL) +
    mytheme() + theme(legend.position = "top")

cowplot::plot_grid(time_bd_plot, time_plot, labels = c("a", "b"), ncol = 2, align = "h",
    axis = "bt", rel_widths = c(1, 0.5))

Fig. 5. (a) Relationship between days post exposure to Bd (ln days) and the infection intensity (ln ZE + 1) at death. Regression lines (solid = adults, dashed = juveniles) represent model estimates and shaded area represents 95% CI. (b) Difference in time to death post-exposure between juvenile and adults.

Supplmentary figures

Fig S8 - Bd-time dependent mortality model

es_all_dat %>%
    ggplot(aes(x = response, y = Zr, colour = trait)) + geom_hline(yintercept = 0,
    linetype = "dashed", alpha = 0.5) + geom_point(aes(shape = life_stage), show.legend = FALSE) +
    viridis::scale_colour_viridis(discrete = TRUE) + xlab(NULL) + ylab(expression("Effect size" ~
    (italic(Z)[r]))) + coord_flip() + facet_wrap(~trait, ncol = 3, scales = "free_y") +
    mytheme() + theme(axis.text = element_text(size = 5), legend.position = "bottom")

Fig. S8. Effect of Bd infection on specific responses across different biological and functional traits. Data presented as individual effect sizes (\(Z_r\)). Circles represent adults and triangles represent juveniles.

---

References

Berger, L., Hyatt, A. D., Speare, R. and Longcore, J. E. (2005). Life cycle stages of the amphibian chytrid Batrachochytrium dendrobatidis. Diseases of Aquatic Organisms 68, 51–63.

Borenstein, M., Hedges, L. V., Higgins, J. P. and Rothstein, H. R. (2009). Converting among effect sizes. In Introduction to meta-analysis (ed. Borenstein, M.), Hedges, L. V.), Higgins, J. P.), and Rothstein, H. R.), West Sussex, Uk: John Wiley & Sons.

Brannelly, L. A., Berger, L., Marrantelli, G. and Skerratt, L. F. (2015). Low humidity is a failed treatment option for chytridiomycosis in the critically endangered southern corroboree frog. Wildlife Research 42, 44–49.

Brutyn, M., D’Herde, K., Dhaenens, M., Rooij, P. V., Verbrugghe, E., Hyatt, A. D., Croubels, S., Deforce, D., Ducatelle, R. and Haesebrouck, F. (2012). Batrachochytrium dendrobatidis zoospore secretions rapidly disturb intercellular junctions in frog skin. Fungal Genetics and Biology 49, 830–837.

Campbell, C. R., Voyles, J., Cook, D. I. and Dinudom, A. (2012). Frog skin epithelium: Electrolyte transport and chytridiomycosis. International Journal of Biochemistry & Cell Biology 44, 431–434.

Grames, E. M., Stillman, A. N., Tingley, M. W. and Elphick, C. S. (2019). An automated approach to identifying search terms for systematic reviews using keyword co‐occurrence networks. Methods in Ecology and Evolution 10, 1645–1654.

Grogan, L. F., Robert, J., Berger, L., Skerratt, L. F., Scheele, B. C., Castley, J. G., Newell, D. A. and McCallum, H. I. (2018a). Review of the amphibian immune response to chytridiomycosis, and future directions. Frontiers in Immunology 9, 2536.

Grogan, L. F., Skerratt, L. F., Berger, L., Cashins, S. D., Trengove, R. D. and Gummer, J. P. A. (2018b). Chytridiomycosis causes catastrophic organism-wide metabolic dysregulation including profound failure of cellular energy pathways. Scientific reports 8, 1–15.

Hedges, L. V. and Olkin, I. (2014). Statistical methods for meta-analysis. London, UK: Academic Press Inc.

Hillyard, S. D., Møbjerg, N., Tanaka, S. and Larsen, E. H. (2008). Osmotic and ionic regulation in amphibians. In Osmotic and ionic regulation: Cells and animals (ed. Evans, D. H.), FL, USA: CRC Press.

Kindermann, C., Narayan, E. J. and Hero, J.-M. (2017). Does physiological response to disease incur cost to reproductive ecology in a sexually dichromatic amphibian species? Comparative Biochemistry and Physiology Part A: Molecular & Integrative Physiology 203, 220–226.

Lipsey, M. W. and Wilson, D. B. (2001). Practical meta-analysis. SAGE publications, Inc.

Martel, A., Spitzen-van der Sluijs, A., Blooi, M., Bert, W., Ducatelle, R., Fisher, M. C., Woeltjes, A., Bosman, W., Chiers, K. and Bossuyt, F. (2013). Batrachochytrium salamandrivorans sp. nov. causes lethal chytridiomycosis in amphibians. Proceedings of the National Academy of Sciences 110, 15325–15329.

Nakagawa, S. and Cuthill, I. C. (2007). Effect size, confidence interval and statistical significance: A practical guide for biologists. Biological Reviews 82, 591–605.

Nakagawa, S., Lagisz, M., Jennions, M., Koricheva, J., Noble, D. W., Parker, T. H., Sánchez-Tójar, A., Yang, Y. and O’dea, R. E. (2022). Methods for testing publication bias in ecological and evolutionary meta-analyses. Methods in Ecology and Evolution 13, 4–21.

Noble, D. W., Lagisz, M., O’dea, R. E. and Nakagawa, S. (2017). Nonindependence and sensitivity analyses in ecological and evolutionary meta‐analyses. Molecular Ecology 26, 2410–2425.

Ouzzani, M., Hammady, H., Fedorowicz, Z. and Elmagarmid, A. (2016). Rayyan—a web and mobile app for systematic reviews. Systematic reviews 5, 1–10.

Peterson, J. D., Steffen, J. E., Reinert, L. K., Cobine, P. A., Appel, A., Rollins-Smith, L. and Mendonça, M. T. (2013). Host stress response is important for the pathogenesis of the deadly amphibian disease, chytridiomycosis, in Litoria caerulea. PLoS ONE 8, e62146.

Richards-Zawacki, C. L. (2010). Thermoregulatory behaviour affects prevalence of chytrid fungal infection in a wild population of panamanian golden frogs. Proceedings of the Royal Society B: Biological Sciences 277, 519–528.

Rollins-Smith, L. A. and Le Sage, E. H. (2021). Batrachochytrium fungi: Stealth invaders in amphibian skin. Current Opinion in Microbiology 61, 124–132.

Rollins-Smith, L. A., Ramsey, J. P., Pask, J. D., Reinert, L. K. and Woodhams, D. C. (2011). Amphibian immune defenses against chytridiomycosis: Impacts of changing environments. Integrative and Comparative Biology 51, 552–562.

Van Rooij, P., Martel, A., Haesebrouck, F. and Pasmans, F. (2015). Amphibian chytridiomycosis: A review with focus on fungus-host interactions. Veterinary Research 46, 137.

Voyles, J., Berger, L., Young, S., Speare, R., Webb, R., Warner, J., Rudd, D., Campbell, R. and Skerratt, L. F. (2007). Electrolyte depletion and osmotic imbalance in amphibians with chytridiomycosis. Diseases of Aquatic Organisms 77, 113–118.

Voyles, J., Young, S., Berger, L., Campbell, C., Voyles, W. F., Dinudom, A., Cook, D., Webb, R., Alford, R. A. and Skerratt, L. F. (2009). Pathogenesis of chytridiomycosis, a cause of catastrophic amphibian declines. Science 326, 582–585.

Wu, N. C., Cramp, R. L. and Franklin, C. E. (2018). Body size influences energetic and osmoregulatory costs in frogs infected with Batrachochytrium dendrobatidis. Scientific Reports 8, 3739.

Wu, N. C., Cramp, R. L., Ohmer, M. E. and Franklin, C. E. (2019). Epidermal epidemic: Unravelling the pathogenesis of chytridiomycosis. Journal of Experimental Biology 222, jeb191817.

---

Session Information

R version 4.2.1 (2022-06-23)

Platform: aarch64-apple-darwin20 (64-bit)

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: pander(v.0.6.5), ggraph(v.2.1.0), litsearchr(v.1.0.0), ggnewscale(v.0.4.8), ggtree(v.3.4.1), phytools(v.1.2-0), maps(v.3.4.1), ape(v.5.6-2), rotl(v.3.0.14), ggeffects(v.1.1.4), performance(v.0.10.1), bayestestR(v.0.13.0), rstan(v.2.21.7), StanHeaders(v.2.21.0-7), brms(v.2.18.0), Rcpp(v.1.0.9), cowplot(v.1.1.1), colorspace(v.2.0-3), forcats(v.0.5.2), stringr(v.1.4.1), dplyr(v.1.0.10), purrr(v.0.3.5), readr(v.2.1.3), tidyr(v.1.2.1), tibble(v.3.1.8), ggplot2(v.3.4.0), tidyverse(v.1.3.2), Matrix(v.1.5-3) and knitr(v.1.41)

loaded via a namespace (and not attached): pacman(v.0.5.1), utf8(v.1.2.2), synthesisr(v.0.3.0), tidyselect(v.1.2.0), htmlwidgets(v.1.5.4), grid(v.4.2.1), combinat(v.0.0-8), munsell(v.0.5.0), codetools(v.0.2-18), DT(v.0.26), rentrez(v.1.2.3), miniUI(v.0.1.1.1), withr(v.2.5.0), Brobdingnag(v.1.2-9), highr(v.0.9), rstudioapi(v.0.14), stats4(v.4.2.1), bayesplot(v.1.10.0), labeling(v.0.4.2), emmeans(v.1.8.2), polyclip(v.1.10-4), mnormt(v.2.1.1), optimParallel(v.1.0-2), farver(v.2.1.1), datawizard(v.0.6.4), bridgesampling(v.1.1-2), coda(v.0.19-4), vctrs(v.0.5.1), treeio(v.1.20.1), generics(v.0.1.3), metafor(v.3.8-1), clusterGeneration(v.1.3.7), xfun(v.0.35), timechange(v.0.1.1), R6(v.2.5.1), markdown(v.1.4), graphlayouts(v.0.8.3), cachem(v.1.0.6), gridGraphics(v.0.5-1), assertthat(v.0.2.1), promises(v.1.2.0.1), scales(v.1.2.1), googlesheets4(v.1.0.1), gtable(v.0.3.1), processx(v.3.8.0), phangorn(v.2.10.0), tidygraph(v.1.2.2), rlang(v.1.0.6), scatterplot3d(v.0.3-42), lazyeval(v.0.2.2), stopwords(v.2.3), gargle(v.1.2.1), broom(v.1.0.1), checkmate(v.2.1.0), inline(v.0.3.19), yaml(v.2.3.6), reshape2(v.1.4.4), abind(v.1.4-5), modelr(v.0.1.10), threejs(v.0.3.3), crosstalk(v.1.2.0), backports(v.1.4.1), httpuv(v.1.6.6), tensorA(v.0.36.2), tools(v.4.2.1), bookdown(v.0.30), ggplotify(v.0.1.0), cubature(v.2.0.4.5), ellipsis(v.0.3.2), jquerylib(v.0.1.4), posterior(v.1.3.1), MCMCglmm(v.2.34), plyr(v.1.8.8), base64enc(v.0.1-3), progress(v.1.2.2), ps(v.1.7.2), prettyunits(v.1.1.1), viridis(v.0.6.2), zoo(v.1.8-11), ggrepel(v.0.9.2), haven(v.2.5.1), fs(v.1.5.2), magrittr(v.2.0.3), colourpicker(v.1.2.0), reprex(v.2.0.2), googledrive(v.2.0.0), mvtnorm(v.1.1-3), matrixStats(v.0.63.0), hms(v.1.1.2), patchwork(v.1.1.2), shinyjs(v.2.1.0), mime(v.0.12), evaluate(v.0.18), xtable(v.1.8-4), shinystan(v.2.6.0), XML(v.3.99-0.12), readxl(v.1.4.1), gridExtra(v.2.3), rstantools(v.2.2.0), compiler(v.4.2.1), crayon(v.1.5.2), htmltools(v.0.5.3), corpcor(v.1.6.10), ggfun(v.0.0.9), later(v.1.3.0), tzdb(v.0.3.0), aplot(v.0.1.9), expm(v.0.999-6), RcppParallel(v.5.1.5), lubridate(v.1.9.0), DBI(v.1.1.3), sjlabelled(v.1.2.0), formatR(v.1.12), tweenr(v.2.0.2), dbplyr(v.2.2.1), MASS(v.7.3-58.1), cli(v.3.4.1), quadprog(v.1.5-8), parallel(v.4.2.1), insight(v.0.18.8), igraph(v.1.3.5), pkgconfig(v.2.0.3), rncl(v.0.8.6), numDeriv(v.2016.8-1.1), xml2(v.1.3.3), ngram(v.3.2.2), dygraphs(v.1.1.1.6), bslib(v.0.4.1), stringdist(v.0.9.10), estimability(v.1.4.1), rvest(v.1.0.3), yulab.utils(v.0.0.5), distributional(v.0.3.1), callr(v.3.7.3), digest(v.0.6.30), rmarkdown(v.2.18), cellranger(v.1.1.0), fastmatch(v.1.1-3), tidytree(v.0.4.1), curl(v.4.3.3), shiny(v.1.7.3), gtools(v.3.9.3), lifecycle(v.1.0.3), nlme(v.3.1-160), jsonlite(v.1.8.3), viridisLite(v.0.4.1), fansi(v.1.0.3), pillar(v.1.8.1), lattice(v.0.20-45), loo(v.2.5.1), fastmap(v.1.1.0), httr(v.1.4.4), plotrix(v.3.8-2), pkgbuild(v.1.4.0), glue(v.1.6.2), xts(v.0.12.2), metadat(v.1.2-0), shinythemes(v.1.2.0), ggforce(v.0.4.1), stringi(v.1.7.8), sass(v.0.4.4) and mathjaxr(v.1.6-0)

---

1. Hawkesbury Institute for the Environment, Western Sydney University, NSW 2753, Australia, nicholas.wu.nz@gmail.com↩︎