Predicting the growth of the amphibian chytrid fungus in varying temperature environments

Abstract Environmental temperature is a crucial abiotic factor that influences the success of ectothermic organisms, including hosts and pathogens in disease systems. One example is the amphibian chytrid fungus, Batrachochytrium dendrobatidis (Bd), which has led to widespread amphibian population declines. Understanding its thermal ecology is essential to effectively predict outbreaks. Studies that examine the impact of temperature on hosts and pathogens often do so in controlled constant temperatures. Although varying temperature experiments are becoming increasingly common, it is unrealistic to test every temperature scenario. Thus, reliable methods that use constant temperature data to predict performance in varying temperatures are needed. In this study, we tested whether we could accurately predict Bd growth in three varying temperature regimes, using a Bayesian hierarchical model fit with constant temperature Bd growth data. We fit the Bayesian hierarchical model five times, each time changing the thermal performance curve (TPC) used to constrain the logistic growth rate to determine how TPCs influence the predictions. We then validated the model predictions using Bd growth data collected from the three tested varying temperature regimes. Although all TPCs overpredicted Bd growth in the varying temperature regimes, some functional forms performed better than others. Varying temperature impacts on disease systems are still not well understood and improving our understanding and methodologies to predict these effects could provide insights into disease systems and help conservation efforts.

. Therefore, it is not surprising that temperature effects on organisms have been studied in many different systems.
Studies that examine temperature effects on organisms often do so in the laboratory under highly controlled conditions where temperature is set to one or multiple constant levels for the duration of an experiment. Many of these studies examine an organism's performance at several temperatures and compare the results at different conditions (Adamo & Lovett, 2011;Bieri et al., 1983;Damos & Savopoulou-Soultani, 2008;Fielding & Ruesink, 1988;Stevenson et al., 2013;Voyles et al., 2017). Such studies have been used to answer questions about temperatures that optimize performance/ traits, the thermal ranges of traits, and thermal adaptation of traits.
One way that the results of these thermal experiments are compared or summarized is by fitting a thermal performance curve (TPC) to the performance data across temperature treatments. TPCs are mathematical functions that are used to describe how an organism performs over a range of temperatures (Huey & Stevenson, 1979). Most of these functions are assumed to be unimodal, meaning that they exhibit an optimal temperature (T opt ) where the trait being described is maximized. These functions also typically include maximum and minimum temperature thresholds denoted as T max and T min , respectively, where the trait being measured approaches, or is equal to, zero (Angilletta Jr, 2006;Huey & Kingsolver, 1989). TPCs have been fit for numerous species across a variety of traits using constant temperature data (Deutsch et al., 2008;Niehaus et al., 2012;Stevenson et al., 2013;Voyles et al., 2017). These functions describe how performance changes over constant temperatures and can be used to infer performance at temperature not directly measured (Huey & Stevenson, 1979).
Thermal performance curves can be applied to microscopic organisms, including pathogens (Ratkowsky et al., 1983;Verant et al., 2012;Voyles et al., 2017). Thus, TPCs are often used in disease studies to quantify how environmental temperatures may regulate disease dynamics. Temperatures can influence both pathogens and hosts and their relationship (Cohen et al., 2017;Gehman et al., 2018;Kirk et al., 2005;Linder et al., 2008). Pathogenic characteristics, such as their distributions, growth rates, and survival rates, can be affected by environmental temperatures (Harvell et al., 2002;Lafferty, 2009;Stevenson et al., 2013).
These data can be used to predict how an organism will perform in more natural thermal conditions because experiments using more complex, fluctuating temperatures are often rare (Greenspan et al., 2017).
However, extrapolating from these constant temperature experiments and TPCs to varying temperature regimes (as are experienced in almost all natural thermal conditions) is difficult and the accuracy of doing so is still debated (Bernhardt et al., 2018;Liu et al., 1995;Niehaus et al., 2012). It is important to improve methods to generalize from constant temperature experiments to more natural varying temperature environments due to the numerous host and pathogen characteristics that temperature can influence (Adamo & Lovett, 2011;Bailey et al., 2017;Stevenson et al., 2013;Ward et al., 2007).
The amphibian chytrid fungus, Batrachochytrium dendrobatidis (Bd), is a fungal pathogen that causes chytridiomycosis (Berger et al., 1998;Fisher et al., 2009). This disease has been linked to amphibian population declines worldwide (Scheele et al., 2019). Field studies have found that Bd infections are strongly affected by temperatures.
For example, infection prevalence and intensity tend to be lower in warmer seasons, lower elevations, and in warmer aquatic and terrestrial habitats (Berger et al., 2004;Forrest & Schlaepfer, 2011;. Laboratory studies have shown that Bd growth and survival are limited at temperatures over 30°C and optimal growth occurs at temperatures between 17 and 25°C (Piotrowski et al., 2004;Stevenson et al., 2013;Woodhams et al., 2008). Bd can grow at temperatures as low as 4°C, although growth at this temperature is slow (Piotrowski et al., 2004;Voyles et al., 2012). Recently, research has also shown how temperature fluctuations can impact Bd's growth rate (Lindauer et al., 2020;Stevenson et al., 2020). Fluctuating temperatures and heat pulses can reduce mortality and morbidity in amphibian hosts (Greenspan et al., 2017;Raffel et al., 2015;Woodhams et al., 2003).
In this study, we assessed the capability of a hierarchical model fit to constant temperature data to predict measurements taken under known fluctuating temperature regimes, using Bd grown in vitro as our study system. We used data collected by Stevenson et al. (2013) and Stevenson et al. (2020) and previously modeled in a different manner by Greenspan et al. (2017). Specifically, we used Bd optical density growth data collected at 10 constant temperature data by Stevenson et al. (2013) to fit a hierarchical logistic growth model. In this model, we constrained the logistic growth rate by a TPC. To test how the choice in TPC influenced the model, we fit the model five times with the same constant temperature data and altered the TPC used in the model. We then used these five fitted models to make predictions about growth in three varying temperature scenarios.
We then validated these predictions by using Bd data grown in the three varying temperature scenarios.
Given previous work (Kingsolver & Woods, 2016;Ma et al., 2015; Tomanek, 2010), we expected that our simple method of generalizing from constant to time-varying temperatures would perform least well for varying temperatures near the TPC peak, due to the nonlinearity. We also expected that the generalization would perform best under the varying temperature regime with the smallest daily fluctuations. Lastly, we expected that the TPCs that had a typical left-skewed shape would make more accurate predictions due to their common use and model comparisons in the literature (Huey & Kingsolver, 1989;Rebaudo et al., 2018;Shi & Ge, 2010).

| Bd cultures
The Bd strain Paluma-Lgenimaculata #2-2011-CO was used for both constant and varying temperature trials in this study. This strain was isolated from a green-eyed tree frog tadpole (Litoria serrata) near Paluma, Queensland, Australia, using the protocol described in Stevenson et al. (2013). The strain was cultured in TGhL broth (8 g tryptone, 1 g gelatin hydrolysate, and 2 g lactose in 1 L of distilled water) and passaged 24 times before the constant temperature experiments and 12 times for the varying temperature experiments.

| Constant temperature experiments
Here, we briefly review the protocols used by Stevenson et al. (2013) to collect constant temperature optical density growth measurements. Bd was grown on TGhL agar plates prior to the study. They allowed the Bd to grow for three days, on these plates, before flooding the plates with TGhL broth. The TGhL broth used to flood the plates was recollected and filtered to remove zoosporangia. Zoospores in the filtered solution were quantified with a hemocytometer. The zoospore concentration in the filtered solution was found to be 5.75 × 10 7 zoospores per ml. The zoospores solution was then used to inoculate 96-well plates used in the constant temperature growth experiments.
Each 96-well plate consisted of 18 positive wells (100 μl of the filtered zoospore solution), six negative wells (100 μl of heat-killed Bd, 60°C for 45 min), and 24 wells containing 100 μl of TGhL. The plates were randomly assigned to 1 of 10 constant temperature treatment (13, 15, 17, 19, 21, 23, 25, 26, 27, and 28°C). were checked daily for contamination, usually seen by high optical density readings and/or discoloration in the well. Contaminated wells were not included in the final analysis. All optical density readings were adjusted by subtracting out the mean of the negative control wells (heat-killed Bd wells).

| Bd optical density logistic growth model
We assumed that the Bd optical density growth pattern could be described by a logistic growth model that had a delay period before Bd's exponential growth ( Figure 1). Specifically, we used the following equation: to fit Bd optical density growth patterns, where Y(t i ) is the initial optical density measurement, d is the delay period, K is the maximum optical density, Y 0 is the initial optical density, r T is the logistic growth rate, and t i is the ith time. This logistic growth model was first described in Voyles et al. (2017) and used to fit Bd optical growth patterns ( Figure   1). The logistic growth equation (Equation 1) was modified to add a delay period due to the time it takes zoospores to settle in the well and develop into zoosporangia. This phase is known as the delay phase and lasts for d amount of time. Next, wells enter the exponential phase where zoosporangia release new zoospores and those new zoospores settle and develop into zoosporangia. The rate at which this happens is controlled by r T ; the higher the value of r T , the faster the growth and development of Bd. Lastly, there is the stationary phase when optical density is not increasing anymore due to new zoosporangia not releasing new zoospores into the wells.

| Thermal performance curve
We found several thermal performance curves (TPCs) in the literature.
The TPCs found varied in functional form and were developed for a wide range of organisms. Most of the TPCs found were phenomenological models, but some were mechanistic models based on enzyme kinetics. We chose five TPCs from the literature that varied in function and purpose for development. The models we selected were as Figure shows a graphical representation of Equation (1) used to fit the Bd optical density data and also shows the three phases of the model. The model starts out at time 0 and Y 0 and enters the delay phase. The length of the delay phase is controlled by d. The steepness of the exponential growth phase is controlled by r T and the stationary phase is reached once the optical density value is equal to K. The logistic growth rate, r T , is temperature sensitive is regulated by a thermal performance curve (shown in the bottom right corner) follows: Briere 2 (Briere et al., 1999), Ratkowsky (Ratkowsky et al., 1983), Ikemoto (Ikemoto, 2005), Logan 10 (Logan et al., 1976), and Stinner (Stinner et al., 1974). The TPCs selected were used to explore how different TPC shapes and mathematical functions could affect predictions made about Bd growth in a varying temperature environment. All TPCs chosen had the logistic growth rate (r T ) go to zero at temperatures below T min and above T max . This choice was made because Bd growth was measured with optical density, which measures everything in the sample, including both live and dead cells. Therefore, we do not expect Bd optical density measures to decrease. All parameter definitions for the chosen TPCs can be found in Appendix S1.
The five TPCs chosen ranged in the number of parameters (4-8).
The five functions also varied in whether they included critical thermal parameters, such as T min , T opt , and T max . Briere 2, Logan 10, and Ratkowsky each have two critical thermal parameters while Ikemoto and Stinner each have only one (Briere et al., 1999;Ikemoto, 2005;Logan et al., 1976;Ratkowsky et al., 1983;Stinner et al., 1974). Most functions came from the entomology literature, with Logan 10 being one of the most popular. Ratkowsky is an outlier in the TPCs we consider, coming from the microbial literature. The Stinner model, the oldest model that we included, has a distinct tabletop shape that differed from the other TPCs. Lastly, some of the TPCs chosen do not have true T min or T max values because they asymptotically approach zero. Thus, for these functions, the logistic growth rate (r T ) will never be exactly zero. In order to allow comparisons between all models, we define effective T min and T max for each TPC, which are the lower and upper temperatures where the logistic growth rate (defined by the TPC) becomes ≤ 0.01.

| Fitting the logistic growth model
We fit Equation (1) to the constant temperature Bd optical density growth data using a hierarchical Bayesian approach. That is, we fit all the constant temperature optical density data simultaneously using one logistic growth model (Equation 1). However, instead of placing a simple prior on the logistic growth rate (r T ), we constrained the shape of the growth rate to conform to one of the five TPC functional forms as presented in the previous section. Each of the five was fit in turn, and the predictions of the fitted models were compared. Mathematically, the hierarchical model is given by: This is the general formulation of the model where f(T, ) denotes the functional form of a chosen TPC with parameter set and the priors distributions of the parameters determining the shape of the TPC are denoted by g(ϴ). Logistic growth parameters priors were based on Voyles et al. (2017) priors. Priors for the parameters in each of the TPC were based on values found in previous publications, which are given in devRate (Rebaudo et al., 2018). TPC priors are shown in Appendix S1.
We chose to fit the model hierarchically instead of fitting 10 separate logistic growth models (1 per constant temperature treatment) and then using posterior distribution samples to fit the thermal performance curves. Fitting this model hierarchically allows the logistic growth parameters and thermal performance curve parameters to be fit together and influence each other.
This also gives the model extra information about the expected value for the logistic growth rate (r T Parameters from the model were estimated via MCMC using the rjags package (Plummer, 2003) in R (Team & R. C., 2013). For each fitted model, five chains were run for 20,000 iterations and the first 10,000 samples discarded for burn-in for a total of 50,000 MCMC samples from the posterior distribution for each model. Trace plots of the MCMC chains were checked visually for convergence. We also checked chains for autocorrelation. Two of the hierarchical models exhibited high autocorrelation. In these cases, the number of iterations was increased from 10,000 to 100,000 then thinned (keeping every 100th sample) to maintain the same number of posterior samples.
To compare the different hierarchical model fits to the constant temperature optical density data, we used deviance information criterion (DIC), which is defined by: where D is the expected deviance and is added to (pD), effective number of parameters (Spiegelhalter et al., 2002). The pD was defined as the one proposed by (Plummer, 2002). We used DIC for model comparison due to its use for Bayesian model comparison by taking into account the fit (expected deviance) and complexity (effective number of parameters) of the model (Spiegelhalter et al., 2002). The DIC calculations for each hierarchical model were done in rjags using the dic.sample function (Plummer, 2003). The DIC values were compared, and the model with the lowest DIC value was considered the best fit to the constant temperature data (Table 1).

| Varying temperature experiments
We validated our predictions made from the constant temperature hierarchical models with varying temperature Bd optical density data (Stevenson et al., 2020). The varying temperature regimes were realistic simulations of regimes experienced by amphibians based on temperature loggers placed in the environment. Thermal regimes were recorded during the wet summer season and dry winter season, at low and high elevations, and in the water, air, or in a frog model (made from agar, using models with both perfect and zero resistance to evaporative water loss; Roznik & Alford, 2014. The agar frog models were placed in daytime and nighttime locations used by tracked frogs . The incubator temperature regimes were created by taking the logged temperature information and averaging it (separately for daytime and nighttime locations) and splitting it into 4-h intervals. A temperature logger was placed in each incubator to record the temperature during the experiments (Figure 2).
We only used data from 3 out of the 12 varying temperature regimes (all reported in Stevenson et al. (2020)) because of natural weather fluctuations in the other nine datasets made it difficult to evaluate consistency in thermal regimes across days. This constraint highlights that our method and other methods using integration (Worner, 1992) are essentially unable to evaluate 75% of our data, which was collected under natural weather conditions. We  Notes: The T min , T max , and r T columns show the posterior mean and the 95% highest posterior density interval in parentheses. Adjusted T min and T max are defined at the temperatures at which the posterior medians of the logistic growth rate reach 0.01. Maximum r T is the temperature at which the logistic growth rate is maximized. Penalized deviance column shows the values calculated by rjags using the dic.sample function. Lastly, the ΔDIC is defined by ΔDIC = DIC i − DIC min . The DIC indicates that the Stinner model is the best fit to the constant temperature data.

F I G U R E 2
Plotted are the varying temperature regimes that were used to grow Bd in incubators. The three regimes are based on temperatures recorded in either low or high altitude, wet or dry seasons, and air or on a frog. Black lines are temperatures recorded in the incubator over time, while the dashed gray lines represent the piece-wise functions fit to the incubator temperatures. There were some inconsistencies in the thermal data, as seen in the spike in the wet high frog regime and when the wet low frog cycle went out of sync. The dry low air also had tiny spikes that the piece-wise function did not take into account due to the short time frame of the spikes

| Predicting Bd performance in varying temperature regimes
To make predictions about how Bd grew in the three varying temperature regimes, we used the following logistic growth function: where x is optical density, r T is the temperature-dependent logistic growth rate, and K is the maximum optical density. With this model, we We wanted to compare predictions about the parameters in the logistic growth model so we refit Equation (1), once per varying temperature regime (n = 3), without constraining the logistic growth rate by a TPC. This model was fit with rjags and used similar priors to the hierarchical model and r T had a prior this time (Gamma(1, 1)). We also assigned Y i a normal distribution due to negative optical density values in the varying temperature treatments. We ran the model for the same number of iterations and removed the same burn-in. We then compared our predicted parameter values to these parameter values, found by just fitting a non-temperature-dependent logistic growth model.

| Thermal performance curve models
The relationships between the logistic growth rate (r T ) and temperature in the hierarchical models were constrained by the thermal performance curve (TPC) used. The structure of the Briere 2, Logan

10, and Ratkowsky models results in similar relationships between
Bd growth rate and temperature-that is, these three models have a left-skewed shape. In contrast, the Ikemoto model is more symmetrical ( Figure 3). The Stinner model is also symmetrical, but is nearly constant at intermediate temperatures.
We compare predictions of key thermal parameters such as T min , T opt , and T max which differed due to using different TPCs to constrain the shape of the logistic growth rate over temperature. At lower temperatures, the predicted logistic growth rates (r T ) patterns were similar in the Logan 10, Ikemoto, and Ratkowsky models. These three models have a less drastic reduction in r T as temperature decreases compared to the other two models. The Briere 2 and Stinner models F I G U R E 3 Thermal performance curves were created from 1000 samples from the posterior distribution of each of the hierarchical models. The output from five thermal performance curves used to constrain the logistic growth rate. Dashed lines represent 95% credible intervals and solid lines represent medians had higher adjusted T min and a more distinctive decrease in r T as temperatures become cooler. The Stinner model had the highest median adjusted T min at 10.95°C (Table 1). There were similar patterns at the upper end of the TPCs. All models except the Ikemoto model predicted a sharp drop in r T at temperatures over the predicted thermal optimums (Figure 3). The Ikemoto model also had a much higher T max than the other four models (Table 1).

| Hierarchical model parameters
Along with the TPC parameters, the hierarchical models also had differences in the logistic growth model parameters, specifically the carrying capacity, K, and the length of the delay period, d. Stinner models predicted a significantly higher K value (Table 2).
In contrast, the estimated initial optical density, Y 0 , was the same across all models.

| Predicting optical density under varying temperature conditions
We attempted to use the fitted models from the previous sections to predict optical density measurements under three time-varying temperature regimes (predictions shown in Figure 4)  Notes: For each model, we report the posterior median values and 95% highest posterior density intervals values calculated from samples from the posterior distributions (N = 1000). Medians and highest posterior density intervals are also shown for the critical thermal values (T min , T opt , and T max ) if these were estimated parameters in the given model (adjusted T min and T max are shown in the supplemental material. All median values for Y 0 are 1 × 10 −3 and have all models have an interval that falls between (1 × 10 −3 , 1.04 × 10 −3 ).
Lastly, the WLF temperature regime had the highest mean and maximum temperature, with a range that spans most of the thermal maxima estimated by the hierarchical models. The Stinner model, which had the best Bd growth predictions in the other two varying temperatures, made the worst prediction at this temperature regime. The Ikemoto model, which had the worst predictions in the other two temperature regimes, made the best prediction in the WLF temperature regime. The Logan 10 and Ratkowsky models, which made similar predictions, did overpredict Bd growth but not by a significant amount. Briere 2 again overpredicted Bd growth at this temperature regime but made a more accurate prediction at this temperature regime when compared to the WHF and DLA temperature regime predictions.

| D ISCUSS I ON
In this paper, we explored the extent to which TPCs estimated under multiple constant temperatures can be generalized to make predictions of growth of a fungal pathogen, Bd, under varying temperature regimes. We also explored how changing the assumed form of the TPC used to constrain the logistic growth rate (r T ) in the hierarchical model affected predictions. Surprisingly, we found that our method did not make the most accurate predictions with the temperature regime that had the smallest fluctuations or with temperature regimes with more intermediate temperatures. We also found that which TPCs were used in the hierarchical model made a difference and that the more typical left-skewed TPCs, like Briere 2, did not make the best predictions. This reinforces the inadequacy of current methods for predicting thermal performance under time-varying temperature conditions, even for these small organisms, from measurements taken at constant temperatures.
All the hierarchical models generally overpredicted the growth  (Piotrowski et al., 2004;Stevenson et al., 2013;Voyles et al., 2017 The impact of temperature variability on the performance of organisms is being explored in the literature, and more studies are trying to predict performance in varying environments (Bernhardt et al., 2018;Denny, 2019;Ferguson & Sinclair, 2020;Khelifa et al., 2019;Kingsolver et al., 2015;Niehaus et al., 2012). The accuracy of methods, such as rate summation, which integrates over time-varying temperature to determine growth or trait performance at a given time, is still under debate (Bernhardt et al., 2018;Liu et al., 1995;Niehaus et al., 2012 Thermal regimes are changing around the world (Easterling et al., 2000). Many publications (Chen et al., 2011;Hinder et al., 2014;Marshall et al., 2010;Narum et al., 2013) have suggested that if the thermal envelopes of areas shift beyond the tolerances of the organisms inhabiting them, those organisms will either die out, emigrate, or perhaps adapt. Changing temperatures could accentuate mismatches between host and pathogen performance curves and produce large changes in interactions (Cohen et al., 2017). However, this may underestimate the affects of the changing climate by not considering variation in temperature. Previous research has shown that temperature variation can affect species' biology and interspecific interactions such as disease (Greenspan et al., 2017;Lindauer et al., 2020;Vasseur et al., 2014;Woodhams et al., 2003 We were unable to accurately predict an organism's performance in vitro in a varying temperature environment from constant temperature data, even using different TPCs and a hierarchical modeling approach that allows us to account for more uncertainty. This highlights a gap in our current models for understanding how organisms, even simple organisms in in vitro studies, react in varying temperature environments. Trying to make predictions with larger organisms we assumed that methods would need to be improved; however, we showed that even when trying to make predictions on smaller simple organisms, there is still a knowledge gap. Building off current methods, future studies could take some of these phenomenological functions and build more mechanistic functions explaining how an organism performance in varying temperature environments. With these models, we can start to determine what factors affect varying temperature performance and when current methods, like our hierarchical modeling approach, might be sufficient. For example, including an acclimation period for larger organisms might be needed instead of allowing the performance or a trait rate to immediately adjust based on the new temperature (Rohr et al., 2018). When temperatures exceed or drop below the thermal maximum and minimum, the organism could also accumulate damage and need to recover from that when temperatures are more optimal, further slowing the performance of a certain trait. Models will also have to make predictions with more irregular and sporadic temperature regimes, caused by weather events or thermoregulatory behavior. Building these new models will not only help improve our predictions of how organisms perform in varying temperature environments, but also improve our understanding of how temperature impacts organisms.

ACK N OWLED G M ENTS
We thank Lisa Belden, Meryl Mims, and Eric Smith for their feedback on the project and manuscript. ZG was partially supported by NIH

DATA AVA I L A B I L I T Y S TAT E M E N T
Constant temperature data used in this manuscript were from Stevenson et al. (2013) and varying temperature data were from Stevenson et al. (2020). Data and code for this manuscript will be publicly available on Zenodo at: https://doi.org/10.5281/zenodo.4525954.