Heritability of Batrachochytrium dendrobatidis burden and its genetic correlation with development time in a population of Common toad (Bufo spinosus)

Despite the important threat that emerging pathogens pose for the conservation of biodiversity as well as human health, very little is known about the adaptive potential of host species to withstand infections. We studied the quantitative genetic architecture responsible for the burden of the fungal pathogen Batrachochytrium dendrobatidis in a population of common toads in conjunction with other life‐history traits (i.e., body size and development rate) that may be affected by common selective pressures. We found a significant heritable component that is associated with fungal burden, which may allow for local adaptation to this pathogen to proceed. In addition, the high genetic correlation found between fungal burden and development time suggests that both traits have to be taken into account in order to assess the adaptive response of host populations to this emerging pathogen.

have to be taken into account in order to assess the adaptive response of host populations to this emerging pathogen.

K E Y W O R D S : Batrachochytrium dendrobatidis, chytridiomycosis, fungal burden, genetic correlation, local adaptation, quantitative genetics.
Emerging diseases have roughly quadrupled over the last 50 years (Jones et al. 2008) threatening biodiversity and human health (Daszak et al. 2000). Pathogens have the potential to induce evolutionary changes in the hosts but the current rates of anthropogenic alteration of the environment may impede this evolutionary response (Longo et al. 2014). Furthermore, while genetic diversity is expected to be high in outbred populations, host organisms may not have enough adaptive potential to develop a rapid evolutionary response to a new selection pressure such as a novel pathogen (Falconer and Mackay 1996). To date, the adaptive potential of wild threatened taxa to respond to pathogen-induced selection is poorly understood.
Previous work estimating the quantitative genetic basis of susceptibility to pathogens has mostly focused either on animal strains for biomedical research (Flint et al. 1995;Rubattu et al. 1996;Råberg et al. 2009) or commercially important species (Price 1985;Roy and Kirchner 2000;Wilfert and Schmid-Hempel 2008). In turn, the complex life-history of many of the studied host species makes it difficult to acquire detailed information about pedigree structure or to undertake informative breeding designs. As a result, most evidence for a genetic basis in susceptibility to pathogens comes from observations of differences in pathogen load among strains/populations and quantitative genetic estimates of broad sense heritabilities (e.g., Mackintosh et al. 2000;McKinney et al. 2011). There is a clear gap in knowledge when it comes to the estimation of the causal quantitative genetic components (i.e., henceforth genetic architecture, Merilä and Sheldon 1999) for pathogen susceptibility in natural populations. This is in spite of the fact that the additive genetic variance is the ultimate determinants of short-term adaptive responses (Falconer and Mackay 1996;Bürger and Lynch 1997;Mousseau et al. 2000).
Fungal diseases have become a major concern for many taxa in the last decades (Fisher et al. 2012). One of the most important emerging fungal diseases is chytridiomycosis, caused by Batrachochytrium dendrobatidis (Bd), a pathogenic, virulent, and highly transmissible fungus. At a global scale, Bd has infected more than 500 species of amphibians (Aanensen 2007), driving many of them to extinction (e.g., La Marca et al. 2005;Schloegel et al. 2006) and the loss of vertebrate biodiversity associated with chytridiomycosis is the most severe in recorded history owing to a pathogen (Skerratt et al. 2007). Nevertheless, while the spatial epidemiology of the Bd panzootic has been the focus of much recent research (Morgan et al. 2007;Farrer et al. 2011;Rosenblum et al. 2013), to our knowledge, no study has yet estimated the quantitative genetic architecture of the fungal load in any amphibian host.
Because traits do not evolve in isolation, it is also necessary to estimate genetic covariances (i.e., the so-called G matrix) with other life-history traits that may be under similar selective pressures (Lande 1979;Rose and Charlesworth 1981;McGuigan 2006). In the case of amphibians, body size at metamorphosis and development rate are well-studied life-history traits (Berven and Gill 1983;Rowe and Ludwig 1991). They are under the selective influence of environmental factors that may also act on immune response, body condition, and exposure time to Bd, affecting the burden and susceptibility of the amphibians to the pathogen. For instance, environmental conditions both in high altitudes and latitudes have been shown to select for faster embryonic and larval development rates in amphibians (Martin and Miaud 1999;Merilä et al. 2000;Miaud and Merilä 2001;Laurila et al. 2002;Morrison and Hero 2003;Muir et al. 2014). Furthermore, amphibians might present trade-offs between immune system and development time or body size. Costs of accelerating larval development may entail a depression in immune response (Gervasi and Foufopoulos 2008) and the activation of immune defences may engender compromises in the condition of hosts (Garner et al. 2009).
In addition, alpine amphibians are likely to be more affected by chytridiomycosis due to the increased survival and persistence of the fungus and the reduced ability of the host to tolerate the infection in these cold regions (Daszak et al. 1999;Walker et al. 2010). In Central Spain, at the Peñalara Massif (altitude ß2000 m), Bd caused the near-extirpation of the common midwife toad, Alytes obstetricans (Bosch et al. 2001) and mass mortalities in the common toad, Bufo spinosus (formerly Bufo bufo (Recuero et al. 2012) (Bosch and Martínez-Solano 2006;Garner et al. 2009;Bosch et al. 2014). Negative population trends have been attributed to chytridiomycosis in a pond of the Peñalara Massif (Bosch et al. 2014). Thus, the Peñalara Massif, within the Guadarrama National Park (GNP), is a unique system that provides the opportunity to monitor the evolution of an amphibian population (i.e., Bufo spinosus) from the initial stages of an outbreak of chytridiomycosis.
The aim of our study is to assess whether a wild population has significant adaptive potential to face an emergent pathogen, and determine the phenotypic and genetic correlations with other relevant life-history traits that might influence the evolutionary response to the disease. To that end, we estimated the relative contribution of additive, maternal, and dominance effects on the overall phenotypic (co)variation of Bd-load, development time and body weight in an amphibian host (Bufo spinosus). We strove to estimate the genetic architecture in its most natural setting (in contrast with most quantitative genetic studies of amphibian evolution that estimate heritable parameters in artificial environments (e.g., Berven 1987;Laurila et al. 2002;Uller et al. 2002;Gomez-Mestre et al. 2004), by conducting our experiment under seminatural conditions in situ in ponds that are exposed to natural temperatures and fluctuations in photoperiod.

STUDY AREA
The study area is a set of three ponds (i.e., Laguna Grande, Laguna Chica, Laguna de Pájaros; detailed location and description in Table 1, Appendix 1) of glacial origin located in the Peñalara Massif in Central Spain (from 1800 to 2200 m; Fig. 1). In early spring, Bufo spinosus breeds in some of the largest and more permanent ponds in this area and each pond acts as a subpopulation.

SAMPLING AND EXPERIMENTAL DESIGN
To estimate additive genetic variance and other causal components of phenotypic variation, we collected 21 adults of Bufo spinosus in Laguna Grande from July 3rd to 9th 2013, and conducted controlled crosses. We performed 23 crosses with 14 males and seven females in the laboratory. When approximately a quarter of the eggs had been laid by a given female, the male fertilizing the clutch was replaced by a new one. Since the common toad exhibits external fertilization, this protocol assures that the offspring of each female results in half sib families. Each female was crossed with four different males, sharing two males with the next female (a scheme of the families that completed the experiment is presented in Table 2, Appendix 1). This design allows us to estimate the contribution of each parent, and their interaction, to the phenotypic variance in the offspring. Thus, providing the   means to separate dominance, additive and maternal effects in the animal model. This crossed design is the most appropriate in an outbreeding species with no parental care, large clutch size and long prereproduction period (sexual maturity at 4-6 years old) (Conner and Hartl 2004). After fertilization each full-sib family was divided in four replicates of ca. 150 eggs.
To estimate phenotypic divergence among ponds, two clutches from Laguna Chica and seven clutches from Laguna de Pájaros, fertilized in the wild, were collected during the same time period. They were mixed to obtain a pool of eggs for each pond. Every pool was divided in four replicates of ca. 150 eggs.
All offspring, both from the quantitative genetic artificial crosses and the clutches collected in the wild, were grown in the laboratory at 15°C until they reached Gosner stage 25 and then were transferred to Laguna Grande. Fifty randomly selected tadpoles from each replicate were raised in a meshed plastic container of 2L in volume placed randomly into groups of 14 containers (henceforth referred as "block"). A total of eight blocks formed a raft floating in Laguna Grande almost completely submerged. The raft was secured to the shore with a rope and tadpoles were fed ad libitum with ground fish food. A data-logger measured temperature in the pond at 10-minute intervals for duration of the experiment. Temperature differences among containers were slight (maximum difference 0.5ºC, thermometer error ± 0.1ºC). Mesh holes of the containers permitted water exchange but prevented the escape of the tadpoles. This impeded any interaction with other amphibian species. However, transmission of the fungus usually involves infected animals (Piotrowski et al. 2004;Medina et al. 2015) since Bd zoospores have low mobility, thus, the containers may significantly reduce exposure to Bd spores. Therefore, to mimic natural contact between species, one overwintered salamander larva from the same pond was introduced into each container for 15 days. As overwintered Salamandra salamandra has been shown to have an infection prevalence of 100% in spring in this pond system (Medina et al. 2015), it was a guaranteed, in this way, to expose experimental animals to infection. Four salamander larvae were sampled after their extraction from the containers and all tested positive for Bd as expected (loads ranging from 1.6 and 6.2 GE).
After 45 days, we adjusted the density inside the containers to avoid the potential effects of varying densities due to mortality. The tadpoles were collected when they reached Gosner stage 42 (four limbs and tail); they were weighed, photographed, and, finally, euthanized with an overdose of benzocaine and conserved in 96% ethanol.

PHENOTYPIC VARIATION
To quantify Bd-loads, DNA was extracted with the reagent Prep-Man Ultra (according to Boyle et al. 2004) from 1063 individual hind limb feet-clips. The extractions were diluted 1/10 before real-time polymerase chain reaction (qPCR) amplification, performed in duplicates, and with Bd genomic equivalent (GE) standards of 100, 10, 1, and 0.1, as well as a negative control. The qPCR were performed in a CFX96 TM Real-Time PCR Detection System, BIO-RAD. Bd-load was measured as the mean of two replicates from the same individual (maximum disparity between replicates = 69,500 GE; mean coefficient of variation per individual = 6,63%; 95% confidence interval = 6.20-7.02%). We considered detection as positive when both duplicates of an individual were positive, amplification curves presented the expected sigmoidal shapes, and the mean was above 0.1 GE.
We considered development time as the time inside the pond, thus, from Gosner stage 25-42. Some individuals were introduced in the experimental pond within few days of each other due to disparities in reaching Gosner stage 25. Accounting for potential effects of variation in thermal conditions caused by these days, individual development time was measured as accumulated degree-days. We calculated daily temperature averages.
Individual development time was the sum of these averages of the days spent inside the pond.
For the body mass, each tadpole was weighed (i.e., fresh weight) at Gosner stage 42 with a precision balance (±1 mg).

NEUTRAL GENETIC VARIATION
Toe-clips were obtained for 108 breeding adults from the three study ponds (40 from Laguna Grande, 27 from Laguna Chica, and 41 from Laguna de Pájaros) from 2011 to 2013. DNA was extracted with a salt protocol (modified from Aljanabi and Martinez (1997), Appendix 2).
Microsatellites BbU14, BbU54, BbU13, BbU49, BbU47, BbU24, BbU39, BbU62, and BbU23 (Brede et al. 2001) were amplified and divided in three groups to perform three multiplexes following a polymerase chain reaction (PCR) protocol of 2 μl of total volume (Kenta et al. 2008) (protocol and PCR conditions in Appendix 2). Microsatellite sequencing was performed using an ABI 3100 automatic DNA Sequencer at the DNA Sequencing Unit of the University of Oviedo.

DATA ANALYSIS
Distributions of the raw data for Bd-load, body weight, and development time were slightly log-normally distributed and normality was achieved by a base 10 logarithmic transformation. Differences in traits between ponds (i.e. phenotypic divergence) were tested with the ANOVA test implemented in the Stats package of the R software v. 3.0.2 (R Core Team 2013). To test for potential covariation, each pair of traits was fitted sequentially into a linear model in which one of the traits was the response variable and the other was a fixed factor; block was used as a random factor. The linear models were fitted with the lme4 package (Bates et al. 2014).
To investigate the genetic components of phenotypic variance, we fitted an animal model (Lynch and Walsh 1998) to the data from the offspring of Laguna Grande. This model enabled the separation of genetic (i.e., additive and dominance) from environmental (i.e., common environment/maternal and residual) variance components in Bd-load, body weight, and development time. Variances, narrow sense heritabilities, and the proportion of total variance explained by each component were estimated through linear-mixed models, using the Bayesian implementation in the R package MCMCglmm (Hadfield 2010) in R v. 3.0.2 (R Core Team 2013). Univariate models for each of the three traits were fitted as well as a trivariate model, analyzing all traits simultaneously and providing estimates of the variance components for each trait and the components of covariance between them. Block was used as random factor in all models to account for potential variation caused by the position in the raft. To ensure that the prior distributions did not affect the posterior estimates, we performed a sensitivity analysis. Different priors were tested in all models, that is Inverse Wishart priors, conservative priors, flat priors, and parameter expanded priors (Supporting Information). Estimates were obtained from 5,000,000 iterations, with a thinning of 100 and a burn in of 1,000,000 iterations. Convergence was estimated by visual inspection and the Heidelberg and Welch diagnostic. The same convergence criteria were applied for all Bayesian analyses presented in this article. The trivariate animal model was specified as follows: ⎡ Where y i a column vector containing the phenotypic values for Bd-load, y j a vector containing the phenotypic values for development time and y k a vector containing the phenotypic values for body weight. The three normally distributed after the normalization by logarithmic transformation. X is a design matrix linking fixed predictors to the data. This matrix has associated vector β, which is the vector of fixed effects, in our case, there were no fixed effects. Z 1 , Z 2 , Z 3, and Z 4 are the incidence matrices for random effects: additive genetic, maternal, genetic dominance, and block effects, respectively. These matrices have associated vectors of coefficients for each animal (individuals from 1 to 1063) representing each random effect contributing to the phenotype. Such vectors of normally distributed random effects are designated as: a for additive genetic, m for maternal effect, d for dominance effect, and b is the block effect. The vector e designates residual variation not explained by the factors in the model (Lynch and Walsh 1998). For instance, the model for individual 1 would be: (y i1 y j1 y k1 ) = X β + Z 1 (a 1 ) + Z 2 (m 1 ) + Z 3 (d 1 ) To ensure that univariate models could separate additive from dominance effects in the traits, we checked the correlation between the posterior distribution of the additive and the dominance variances.
For the microsatellite data, scoring errors, large allele dropout, and null alleles were checked using Micro-Checker (Van Oosterhout et al. 2004). We used the software Genepop v.4.2 (Raymond and Rousset 1995) to infer the effective number of migrants in the system and its Markov Chain Algorithm (1000 dememorisation steps, 1000 batches, 1000 iterations per batch) to test for deviations from Hardy-Weinberg Equilibrium (HWE) We used microsatellites, pedigree, and phenotypic data to study the relative roles of drift and selection on subpopulation divergence. Microsatellite data from adults were analyzed with the RAFM package (Karhunen and Ovaskainen 2012) using the R v. 3.0.2 software (R Core Team 2013). A coancestry coefficient matrix for each subpopulation was obtained from 1,000,000 iterations, with a thinning of 10 and a burn in of 1000 iterations. The coancestry matrices were used as the prior for the Driftsel package (Karhunen et al. 2013) also implemented in R v. 3.0.2 software. Driftsel detects whether divergence in trait means deviates from that expected under random drift and could be attributed to natural selection. Driftsel was run for 1,000,000 iterations, with a thinning of 10 and a burn-in of 1000 iterations. Experimental block was used as random factor and 1,000,000 tmpmax was fixed for sparse Cholesky decomposition.

PHENOTYPIC DIVERGENCE
The three traits revealed significant statistical differences among subpopulations (Bd-load P < 0.0001, F 2, 1153.8 = 14.217; Development time P < 0.01, F 2, 1152.7 = 6.181; Body mass P < 0.0001, F 2, 1160.3 = 10.035). The most conspicuous differences among ponds occurred for Bd-load, with the offspring from Laguna Chica having a much lower load ( Fig. 2 Analysis of phenotypic covariance between traits resulted in a significant positive covariance between Bd-load and development time (correlation = 0.38; P-value < 2.2 × 10 -16 ; Fig. 3) but not between any other pairs of traits.

QUANTITATIVE GENETIC ARCHITECTURE
Sequential fitting of the data to animal models with increasing complexity (i.e., additive + residual, additive + dominance + residual, additive + maternal + residual, and the complete model estimable: additive + dominance + maternal + residual) showed that the complete model (i.e., additive + dominance + maternal + residual) had the best fit both for the univariate and trivariate cases (i.e., lowest DIC values). All models reached good convergence.
Univariate animal models were insensitive to the prior used for all traits. These models separated additive from dominance effects well in all traits except for body weight. Body weight exhibited a strong negative correlation between the posterior distribution of the additive and the dominance variances, -0.82 (95% confidence interval: -0.8299 --0.8157). Thus, the estimates of additive and dominance variances for body weight are informative when it comes to the broad-sense heritability but our model could not separate properly the additive from the dominance effects, hence, separate estimates provided in Table 1 have to be taken with caution.
While the posterior distribution of the genetic correlation between Bd-load and development time was always accurate and consistent irrespectively of the prior used, the trivariate model yielded flat uninformative posterior distributions for individual trait estimates. This may be due to the difficulty of estimating the joint distribution of components of the G-matrix as compared to the distribution on its individual component separately, given the data available (Ovaskainen et al. 2008;Teplitsky et al. 2014). In this article, we present the narrow-sense estimates of additive, maternal, and dominance variances obtained from the univariate analyses and the genetic correlation from the trivariate model, both using a weak and unbiased prior with an Inverse Wishart distribution (see methods section). Variance estimates (modes) and the proportions of total variance explained are reported in Table 1. All three traits had a strong heritable component. In particular, the mode of the heritability of Bd-load (h 2 = 0.22) for additive effects indicated an important adaptive potential for this trait. Since our tested animals were exposed to Bd all the time and, to be sure that the additive estimates of Bd-load were not influenced by variation in exposure time we also estimated the heritability of Bd-load with development time as covariate. The adjusted heritability still resulted significant (h 2 ß0.15) and close to the previous value (h 2 = 0.22). The genetic dominance component of development time (Dß0.10) was of similar magnitude as the additive (h 2 ß0.11). Body weight has a strong broadsense heritable basis (i.e., additive + dominance) but we could not separate reliably additive from dominance effects. Development time and Bd-load presented a significant positive genetic correlation of 0.17 (95% HPDI, 0.0482 -0.2989). Thus, 0.17 of the overall phenotypic covariation between Bd-load and development time in Laguna Grande is due to additive genetic effects.

GENETIC DRIFT IN THE POPULATION
No evidence of scoring errors, allelic dropout, or stuttering was detected with Micro-Checker software except for possible null alleles at locus BbU14 in Laguna Grange and Laguna de Pájaros, and loci BbU13 and BbU49 in Laguna de Pájaros. Microsatellite loci conformed with HWE in most cases; only loci with possible null alleles (i.e., heterozygote deficit) had deviation from HWE. We also found slight linkage disequilibrium between BbU13 and BbU47. Analyses with and without these loci resulted in very similar conclusions; therefore, we present the results of the analyses including all microsatellite loci.
Microsatellite data results indicate low divergence between subpopulations, with low F ST = 0.021 (±0.007) and D est = 0.022. Both F ST and D est pairwise values were < 0.04 (Table 2). Both estimators found the lowest neutral genetic divergence between the nearest subpopulations: Laguna Chica and Laguna Grande.
Regarding overall migration among ponds inferred by microsatellite data, the effective number of migrants was six per generation.

POTENTIAL ROLE OF NATURAL SELECTION
Both RAFM and Driftsel Bayesian models reached convergence. No significant effect of natural selection on differentiation among subpopulation means was revealed by the Driftsel package. All the subpopulation trait means were located within the variation range expected under genetic drift. Values of the neutrality test close to 0.5 indicate drift whereas values close to 0 or 1 indicate stabilizing or diversifying selection, respectively. In our case neutrality test resulted in S = 0.76, a value not high enough to consider differentiation among subpopulations as caused by diversifying selection.

Discussion
Our study provides a detailed estimate of the quantitative genetic architecture of infection susceptibility to an emergent pathogen in a natural host population. To our knowledge, this is the first estimate of the adaptive potential to Bd-load in an amphibian host. Furthermore, our results revealed a significant genetic correlation between development time of the host and susceptibility to the infection. Further research is needed to better understand the potential joint evolution of these important life-history traits in other populations and taxa.

ADAPTIVE POTENTIAL AND GENETIC CORRELATION
We show that a host population can harbor significant evolutionary potential when facing a novel pathogen. Our results suggest that Bd-load might respond fast to strong selection since its heritable basis exhibited a predominantly additive component. The maternal effect and genetic dominance for Bd-load explained less than 0.01 of the overall phenotypic variation. Development time presented a quantitative genetic architecture with a lower additive component but larger dominance than Bd-load. This indicates that, under similar selection pressure, development time might display a slower evolutionary response (Roff 1997;Lynch and Walsh 1998). However, the response to selection of development time and Bd-load should not be considered independently. Bd-load and development time showed a significant positive phenotypic covariation in the offspring of Laguna Grande. A portion of this covariation could be due to Bd-exposure time, because the slower that development occurs, the greater is the time for keratin to accumulate Bd zoospores. In fact, it is well known that the Bd panzootic most aggressively affects those amphibian species with longer larval periods and, thus, greater exposure to the fungus (Blaustein et al. 2004;Carey et al. 2006). Our results demonstrated that around 0.17 of the phenotypic covariation is due to genetic correlation. Therefore, the selection on either Bd-load or development time may entail an evolutionary response on the other trait. For instance, the strong selection for shorter development time, as could be expected in situations of restricted growth opportunities as mountain environments, high latitudes, and temporal ponds (Miaud and Merilä 2001;Laurila et al. 2002;Muir et al. 2014), might produce genotypes with lower propensity to reach high Bd-loads.
Conversely, body weight appears to be able to evolve independently from the other traits since it did not show significant genetic correlations nor phenotypic covariation. Body weight had the largest broad sense heritability and may have a strong evolutionary response to selection. However, our data did not allow for a reliable separation of additive from dominance effects in this trait.
We would like to note that our detailed quantitative genetic estimates come from the population at Laguna Grande and, thus, could be specific to this population. Further population replication, including also lowland and Bd-free populations, would be needed to assess the generality of our findings.

EVOLUTIONARY PROCESSES IN THE STUDY SYSTEM
The host population in our system may be dealing with diverse selection pressures. In addition to the new pathogen, common toads are encountering novel environmental conditions in the recently colonized ponds. Laguna de Pájaros and Laguna Chica were occupied by Bufo spinosus after the Bd outbreak, while Laguna Grande was an established population long before (Bosch and Rincón 2008). Laguna Grande is the largest pond both in size and with respect to effective population size (28 breeding females captured in 2013 in contrast with 17 in Laguna de Pájaros and six in Laguna Chica; J. Bosch, unpubl. data), therefore, selection pressure being equal in all ponds, a faster rate of adaptation would be expected in Laguna Grande (Soulé 1976;Frankham 1996). Nevertheless, when it comes to phenotypic divergence, Laguna Chica is the most differentiated with the lowest Bd-load and development time (Fig. 2). This is the smallest pond, suffering strong thermal fluctuations and risk of desiccation at the end of the summer (J. Bosch, unpubl. data), which might impose a strong selection toward reducing development time. Due to the genetic correlation between development time and Bd-load, low fungal burden could be indirectly selected for. Furthermore, common toads in Laguna Chica might face less infection risk than the other study ponds because of the lack of overwintering salamander larvae, which form the main reservoir of Bd in the system after the disappearance of Alytes obstetricans (Medina et al. 2015).
Because adaptive divergence can occur even under moderate gene flow (Endler 1973;Rice and Hostert 1993;Garant et al. 2007) we assessed the possibility of adaptive divergence among the toad populations in the three study ponds. Comparison between quantitative genetic and neutral genetic differentiation may help us to assess whether the observed phenotypic divergence is likely to be caused by directional selection (Merilä and Crnokrak 2001;McKay and Latta 2002;Leinonen et al. 2008). In our study, the comprehensive tests implemented in Karhunen et al. (2013), combining multivariate coancestry and quantitative genetic coefficients, did not provide support for a role of divergent selection in the observed differences. This is not to say that selection is not acting on the studied traits but rather that we cannot discard genetic drift as a relevant cause of the current phenotypic differentiation.

DEVELOPMENT OR Bd-LOAD
Body size in ectotherms is expected to have a positive phenotypic relationship with development time (e.g., Laurila et al. 2002;Morey and Reznick 2004;Rudolf and Rödel 2007). However, occasionally, this relationship can also be negative or nonexistent (Travis 1981;Pfennig et al. 1991) that has been related to variation in food availability and mortality rates (Alford and Harris 1988;Rudolf and Rödel 2007). In our study, we did not find any significant phenotypic covariation between the two traits although food availability was not constrained. Due to an unusual heat peak, a high mortality episode was recorded soon after the transfer of the experimental animals to the pond. Dead larvae in the pond decomposed, or were ingested, too fast to measure body size. After this mortality episode, unrelated to Bd, no further significant mortality was detected until the termination of the experiment.
We also explored potential relationships between body weight and Bd-load. Previous studies found a negative effect of Bd-load (Voyles et al. 2012) and Bd-exposure (Garner et al. 2009) on body size in amphibian adults and juveniles, respectively. The lack of covariation of Bd-load with body weight in our results may imply that negative effects of high Bd-loads could become apparent later on in development, although differences in experimental design among studies could also play a role. For instance, since all our tadpoles were exposed to Bd, we cannot assess any potential effects of exposure per se on body weight. Further studies, linking directly Bd-load with mortality, are needed to better understand the defense mechanisms against the pathogen (e.g., resistance and/or tolerance) and its evolutionary relevance.

Conclusions
Based on the significant genetic correlation between Bd-load and development time found in this study, further work in other populations and taxa is required to better disentangle the evolutionary significance of the relationship between development rates and susceptibility to diseases. This is especially so for taxa with complex life cycles that may present a trade-off between immunity and development rate in time-constrained situations (Rolff et al. 2004) as suggested for amphibians (Gervasi and Foufopoulos 2008).
The heritable component associated with Bd-load may allow common toad populations to adapt to chytriodiomycosis. In this context, it is necessary to study the strength of Bd-induced selection (i.e., estimation of selection coefficients) and the relationship between Bd-load, host condition, development rate, and mortality in nature.
Given the detailed information since the outbreak of the infection and establishment of new common toad subpopulations, this study sets an initial time point for long-term evolutionary monitoring of the system. The significant genetic component of Bd-load, coupled with the current ease of deep sequencing technologies, makes this species a candidate suitable for mapping the molecular basis responsible for the fungal burden.

Ethics statement
Animal welfare, handling, and field permits to conduct this research were obtained from the Consejería de Medio Ambiente, Vivienda y Ordenación del Territorio de la Comunidad de Madrid (permit number: 10/071126.9/13).