Do pathogens become more virulent as they spread? Evidence from the amphibian declines in Central America

1. Introduction

Pathogen virulence can vary dramatically over space and time [1]. This variation in virulence can have profound implications for disease management and risk assessment [2], as well as for the success of biological control programmes [3,4]. Much of the variation in virulence that we observe is cyclical, being driven by density dependence and/or seasonal environmental shifts [5–8].

Less commonly, however, we observe directional shifts in virulence. These directional shifts are almost always patterns of declining virulence in biological control agents, a pattern attributable to rapid coevolutionary shifts in both host and pathogen [3,9]. In theory, virulence is usually thought not to increase over time, and when it does so this pattern is usually attributed to environmental shifts [6,10–12] rather than evolution.

It is, however, increasingly appreciated that spatial effects can cause very different evolutionary trajectories from those occurring in equilibrium populations [13,14]. Range-shifting populations, for example, experience selection pressures fundamentally different to those in stable populations. In particular, populations on an expanding range edge experience low conspecific density and are spatially assorted by dispersal ability [15]. These conditions produce strong evolutionary pressures favouring increased dispersal and population growth rates [16]. Although most systems in which these processes have been examined to date have involved range expansion in non-parasitic species [17–20], the same evolutionary forces should apply equally to parasitic species as they spread into a new environment (which for a parasite is a naive population of hosts) [21].

The evolutionary influence of invasive spread on parasite life history is unlikely to be simple, and we currently have little specific theory on the issue to guide us. Nonetheless, it seems likely that as a parasite invades a naive host population, parasites in the vanguard of that invasion will continuously experience a high density of susceptible hosts relative to that found behind the invasion front. When the parasite in question is a pathogen, these conditions of continuously high host density should select for increased virulence on the invasion front [22–24], particularly in species where dispersal is decoupled from local transmission. Thus, as a pathogen invades a naive host population, we might expect to see rapid evolution of increased virulence in the invasion vanguard. To our knowledge, the possibility that pathogens increase in virulence as a consequence of their spatial spread has never been empirically assessed. Indeed, most studies examining changes in virulence are sampled through time rather than through space, and most models of host–parasite coevolution are aspatial [14]. Here, we examine the possibility of increasing virulence during the spread of an invading pathogen: the amphibian chytrid fungus.

Since the late 1970s, many amphibian populations around the world have suffered major declines and extinctions, even inside protected areas [25]. The emerging infectious disease chytridiomycosis (caused by what is likely a virulent strain of the chytrid fungus, Batrachochytrium dendrobatidis, Bd) is now generally considered to be the main driver of these declines [26–28]. The evidence implicating Bd and chytridiomycosis in the amphibian declines is compelling, and consists of three lines: (i) Bd can be pathogenic to amphibians in both the field and laboratory [26,29]; (ii) genetic evidence suggests the emergence of a hypervirulent strain of Bd that shows genetic signal consistent with range expansion [30–32]; and (iii) amphibian population declines appear to have followed a broad wave-like pattern consistent with the spread of a novel pathogen (for example, from north to south in Central America and from south to north in eastern Australia) [27,33].

Given the likelihood of recent spread of a virulent strain of Bd in many areas of the world, we examined this system to assess the possibility of a directional increase in virulence during spread. To do this, we examined the relationship between the timing of Bd arrival and the timing of amphibian decline across 13 sites in Central America. Central America is one of the best-studied regions in the world with regard to amphibian declines and chytridiomycosis [27,34–41], so key data on the prevalence of Bd before and after population declines are often available. If Bd has evolved increased virulence over the course of its invasion through Central America, then we would expect to see a lag between pathogen arrival and population decline that diminishes over time.

2. Material and methods

(a) Amphibian decline and prevalence data for Central America

Two datasets were used for the analysis. The first contains the decline dates for 13 localities from southern Mexico to central Panama, which were sourced from the literature as follows. Southern Mexico: Parque Nacional El Chico, Hidalgo [41]; Cerro Chicahuaxtla, Veracruz [41]; Cerro San Felipe, Oaxaca [41]. Honduras: Pico Bonito National Park [42,43]; El Cusuco [43]. Guatemala: San Marcos [41]. Costa Rica: Monteverde [38,44]; Las Tablas [34]. Panama: Fortuna [35]; Santa Fe [45]; El Copé [27]; El Valle [46]; Campana [46]. For most localities, the date of decline is not very precise. To achieve some consistency, we defined the earliest date of decline (^ED) across all localities as the last time a population was sampled before a decline was observed and the latest date of decline (^LD) as the first sampling period in which it was obvious that the population had gone through a population decline. Where there was ambiguity (e.g. one species declines but another does not decline for a few years), we have opted for choices that maximize our bounds around decline date. Thus, we have opted for bounds on decline date that capture our true uncertainty. These bounds represent the earliest and latest possible dates at which decline could have occurred.

The second dataset we use documents localities at which Bd prevalence has been surveyed in the region. This latter dataset on Bd prevalence spans a much larger range of localities than the decline data [27,36,39,41,45–48]. To deal with this, spatial data were collapsed into spatial bins (herafter, ‘sites’) for the purposes of estimating Bd arrival time at each site. Bins were defined around localities where declines are known to have occurred and for which estimates of earliest and latest dates of decline were available. We placed Bd sampling localities into the site representing the closest decline locality. This ensured that each sampling locality was always associated with the closest decline locality (for map of binned data and the full dataset, see the electronic supplementary material).

(b) Estimating the timing of Batrachochytrium dendrobatidis arrival and amphibian decline across the region

The basic premise of our analysis was to first estimate date of arrival and post-arrival prevalence of Bd, and then regress the decline date against these estimated arrival dates and post-arrival prevalences. A regression slope of 1 for the partial coefficient of decline date on arrival date would indicate a consistent time lag between Bd arrival and population decline across multiple sites. Importantly, a partial coefficient different from 1 in this analysis indicates that the time lag between Bd arrival and population shows directional shifts through time.

In our analysis, we assume that both decline date and the infection status of frogs depends on an underlying process: the arrival of Bd in the system. The process of Bd arrival at each site, i, is represented by a simple threshold model in which, for each site at some point in time (T_i), Bd switches from being absent to being present with some mean prevalence (P_i) [49]. We treat the number of infected frogs per sampling time within site (^infN_i,t) as a draw from a binomial distribution with a sample size equal to the number of frogs sampled at that site.time (^totN_i,t)

where p_i,t describes the thresholdWe used a uniform prior (spanning dates between 1950 and 2010) for each T_i and a flat beta prior (alpha = beta = 1) for each P_i.

In our model, the underlying process of Bd arrival also affects the date of decline. We treat the date of decline at each site as an unobserved event d_i, occurring between the earliest and latest possible decline dates. We treat this unobserved date of decline as a draw from a normal distribution, censored at the earliest and latest possible decline dates:

where I enforces the censoring interval (by setting the likelihood to zero outside the interval defined by earliest and latest decline dates), φ is the inverse of the error variance (the precision) andis the multiple regression model linking arrival time and prevalence to decline date at each site. In this last equation, α is the intercept and the β terms are the partial regression coefficients. A graphical representation of our model is provided in the electronic supplementary material.

For each site, then, we estimated the arrival time of Bd, and simultaneously estimated the overall linear regression parameters linking decline date to arrival date and post-arrival prevalence.

The model was fitted in a Bayesian framework using the JAGS Gibb's sampler [50], with flat priors for arrival time, T (between −30 and 30) and P (between 0 and 1 using a flat Beta distribution); and minimally informative Gaussian priors (mean = 0, variance = 10⁵ for the regression coefficients; table 1). To speed convergence, we roughly zero-centred our time data by subtracting 1980 from all values of year. Therefore, the priors for T actually represent the interval between 1950 and 2010. The uniform prior on T was chosen for logistical reasons (the posterior of T is also uniform) rather than representing an a priori assumption that Bd was not present in the region before 1950. Because of this, 95% credible interval bounds that approach either 1950 or 2010 should be treated as conveying no real information on actual arrival time at that bound.

Posterior densities for the model coefficients were estimated from three chains, each of 100 000 samples following a burn-in of 10 000 iterations. Convergence was assessed both by trace plots and with the Gelman–Rubin convergence statistic [51]. The ability of the model to estimate parameters was confirmed using simulated data. All data manipulation and analysis was conducted in R [52], and the scripts and data are available in the electronic supplementary material.

(c) Assessment of parameter bias

Our prevalence dataset exhibited large variations in the intensity and frequency of sampling over time. Some populations were sampled regularly, others occasionally, some with very large samples and other with very small samples (see the electronic supplementary material, figure S7). This raises the possibility that uneven sampling through time or space was affecting our parameter estimates. To control as much as possible for this possibility, we ran our model on two subsets of the data. In the first, we excluded the extremely intensively studied site (El Copé), to even out sampling effort, and in the second, we excluded the sites at which spatial binning occurred (Monteverde and Las Tablas), to remove a possible binning effect.

Additionally, because we were particularly interested in the partial regression coefficients of decline date on arrival date, it was important to understand the extent of possible bias in these parameter estimates driven by uncertainty in arrival time (which might cause us to underestimate the slope of the true relationship through regression to the mean). To assess this, we took two approaches. First, we simulated data in which 13 hypothetical ‘true’ arrival and decline dates were defined, but in which our uncertainty around the arrival dates varies. By examining the way in which estimation bias scales with uncertainty in arrival date in our model, we can assess the probable impact of regression to the mean on our final inference. These ‘true’ times were considered equal (i.e. decline time equals arrival time) and were placed regularly over a 40-year interval (a period of time similar to that spanned by the data). We simulated the interval censoring in our decline data by adding and subtracting random draws from a uniform distribution (U(0,7), in years) to the ‘true’ date to create an interval. To control observation error in arrival time, we then defined a sampling interval, I, which was centred over the true arrival time. At each site, a sample of 50 ‘frogs’ with zero prevalence was created at the beginning of this interval, and a sample again of size 50, but of prevalence 0.3, was created at the end of this interval. The resulting uncertainty in arrival time is then approximately the size of the defined interval. By changing the value of I (from 1 to 20 years), we can introduce increasing levels of uncertainty into our arrival times and assess the effect of regression to the mean on our slope estimate.

To do this, for each integer value of I between 1 and 20, we generated 1000 datasets and fitted our model to these simulated datasets. Because this was a computationally intensive exercise, we simplified the model by dropping the effect of post-arrival prevalence from the model and only sampled the posterior distribution of the partial coefficient of arrival time. Convergence was not checked on these 20 000 model fits.

Our second approach to assessing the bias in our estimate of the true relationship between decline date and arrival time was to make use of the SIMEX algorithm [53]. In this algorithm, we stepwise inflate the observation error in our real data, calculating the parameters at each level of inflation. By examining the relationship between level of inflation and parameter estimates, we can extrapolate back to estimate the parameter estimate we could expect with zero observation error. To do this, we inflated the time intervals at which prevalence samples were undertaken; inflating the interval by a factor of {1, 1.4, 1.8, … , 3} times the actual interval. To keep the model coherent, we centred this inflation for each site on the point in time between the first observation of Bd at each site and the previous sampling period (in which Bd was not detected). We then estimated regression slopes under the full model and observed how these scaled with our inflation of sampling interval. We then fitted a linear regression to these data to estimate the parameter value we could expect under zero observation error.

3. Results

Our analysis consists of two parts, operating simultaneously. First, we estimate arrival date of Bd within each site using a simple threshold model. Second, we fit a multiple regression model across sites between decline date (the response) and arrival time and post-arrival prevalence (the independent variables). Information flows between the levels of this hierarchy, and in our case, this had the effect of providing bounds on arrival date at each site that were more precise than estimates made purely at the site level. Indeed, the model was able to estimate arrival dates at one site, EV, which has only a very small sample of Bd prevalence (n = 4) and at which Bd has yet to be reported in the literature (figure 1). Nonetheless, estimates of arrival date were still typically broad at each site (see figure 1; electronic supplementary material, table SI.1).

• Download figure
• Open in new tab
• Download PowerPoint

When we examined the relationships between dates of decline, pathogen prevalence and arrival dates, two interesting patterns arise. First, post-arrival prevalence was typically low (with only four sites having post-arrival prevalences greater than 0.4; see electronic supplementary material), and had no substantial effect on date of decline (95% credible interval for effect on decline date: −4.9 to 4.5).

Second, we find a significant positive relationship between decline dates and Bd arrival dates across southern Mexico and Central America. This is the expected result if the arrival of Bd is associated with population declines, and is consistent with the notion that Bd has spread through the region. Our analysis also revealed that the slope of this positive relationship is clearly less than one (95% credible interval for effect of arrival date on decline date: 0.65–0.89). This latter result was qualitatively unchanged when we fitted the model to the data subsets with reduced variance in sampling intensity and bin size (see electronic supplementary material).

Importantly, our data simulations suggest that the hierarchical model performs well in the face of observation error, consistently producing parameter estimates that are unbiased at reasonable levels of uncertainty in arrival date. When arrival date becomes sufficiently uncertain (greater than ±10 years), however, parameter estimation bias does still become a problem. At levels of observation error around arrival time similar to that in our data (less than ±5 years), estimation bias is likely to be very small (slope estimates biased downwards by less than 0.02; see electronic supplementary material). Moreover, the SIMEX algorithm suggests that inflating the observation error around arrival time in our data makes little difference to our estimate of the regression slope (see electronic supplementary material). In total, these results suggest that our estimated slope values of less than one are not statistical artefacts, and that the lag between Bd arrival and naive host population decline has, in fact, decreased over time. Thus, we find strong evidence that amphibian decline is linked temporally to the arrival of Bd across Central America, but also that Bd appears to have become more virulent as it spread southwards (figure 2).

• Download figure
• Open in new tab
• Download PowerPoint

4. Discussion

The estimated arrival dates we observe for Bd across Central America are consistent with a general pattern of spread of Bd moving from north to south (figure 1). Furthermore, when we examined the correlation between Bd arrival and population decline across sites there was an unequivocal positive correlation between the two: Bd appears to have been spreading, and frog populations clearly decline after Bd becomes prevalent in the system (figure 2). This link between Bd arrival and population decline across an entire region has proved difficult to make in other areas due to insufficient time series of pathogen sampling [49]. Having made this link between pathogen arrival and population declines, however, what can we infer about changes in virulence during pathogen spread?

First, our analysis shows that post-arrival prevalence has little effect on decline date. Although this seems counterintuitive, it probably just reflects the highly variable nature of Bd prevalence and the vagaries of the timing of sampling. Bd prevalence varies strongly with season [54], and probably also varies with host density [55]. If the only post-arrival samples occurred during a moment of peak prevalence and population die-off, for example, then this population will have a high estimated post-arrival prevalence. If, on the other hand, this same population was sampled some time later when host densities were greatly reduced and the season was poor for Bd, then the same population would have a low estimated post-arrival prevalence. Clearly, the population's decline date is the same in both cases. It seems likely, then, that the non-effect of prevalence on decline date in our analysis probably just reflects ad hoc sampling of a dynamic system.

More interestingly, our analysis shows that, while positive, the slope of the relationship between date of population decline and Bd arrival is almost certainly less than one (figure 2 and table 1). We initially thought sampling biases inherent in the data might drive this low slope value (e.g. changes in the frequency and intensity of sampling over time, or an uneven binning range through space). Subsetting our data to reduce or remove these possible biases, however, failed to qualitatively change our results (see electronic supplementary material). We then were concerned that the low slope values might represent an instance of regression to the mean rather than a true pattern. Analysis of simulated data and variance-inflated real data, however, led us to conclude that the bias in the estimated slope values is negligible in our case. In the light of these results, we feel confident that the low slope values are not driven by regression to the mean, nor by uneven sampling or binning across sites. The unexpectedly low slope value between dates of population decline and Bd arrival probably represents a real pattern: one of a diminishing lag period between pathogen arrival and population decline. Virulence appears to have increased during the invasion of Bd through Central America.

The evolution of increased virulence in this case is a possibility for two reasons. To begin with, we expect it based on the simple theory outlined in §1. It is a routine observation that pathogens evolve very rapidly in response to host defences as well as to trade-offs between transmission and virulence [56–58]. These trade-offs are greatly weakened on the advancing front of an invasive pathogen population, because susceptible hosts are constantly at high density relative to places where the pathogen has been long-established. In long-established populations, susceptible hosts are likely to be at lower density, because some portion of the host population is infected already, and some portion of the population may have acquired resistance.

This differential in host density is even stronger in the case where the arrival of the pathogen causes population declines, as commonly happens with the arrival of Bd [27,59]. As a consequence, the discrepancy in the density of susceptible hosts between vanguard and older populations of Bd may be highly pronounced. Invasion-front Bd potentially never experiences the low density of susceptible hosts that endemic Bd experiences. Relative to endemic Bd, then, Bd on the invasion front consistently experiences easier transmission, and so Bd variants with high growth rates (and, as a consequence, greater virulence) would be favoured [22,24]. In such a situation, virulence would be under strong directional selection [15,21]. Recent work by Rosenblum et al. [60] suggests that Bd has a very dynamic genome, which could translate into rapid functional shifts over short periods of time. Moreover, the generation time of Bd is in the order of days [61], so actually there have been thousands of generations elapsed since it was first detected in Mexican samples collected in the 1970s. Thus, evolution would not be unexpected, nor have to be particularly rapid to cause significant changes in the virulence of Bd as it spread.

The second argument for the role of evolution in driving changes in virulence is that we can infer that normal coevolutionary processes are occurring in this system. If evolution is playing a role in this system, we would expect the initial virulence of Bd to be high at a newly invaded site, but to diminish at that site as time goes on and the usual coevolutionary forces come into play. This expectation is consistent with the observation that many populations of amphibians that survive the arrival of Bd continue to persist, and many recover, in the time following its arrival [38,62,63]. Indeed, even species that were thought to be extinct and have not been observed in the wild for decades are currently reappearing [64–67]. Some of these rediscoveries reflect populations in areas that are climatically unsuitable for Bd, but in others Bd is highly prevalent and causes ongoing mortality, suggesting the emergence of a stable host–pathogen system [63,68]. These population recoveries are expected in a system in which host–pathogen coevolution occurs, but are difficult to explain if evolution does not occur and virulence is driven purely by directional environmental shifts. Although we cannot say whether these population recoveries are driven by evolution in the host or pathogen, or (as seems most likely) both, they do point to an evolutionarily labile system in which virulence can change rapidly.

Nonetheless, shifts in the environment clearly play an important role in the virulence of Bd. It is clear that the virulence of Bd is strongly related to its abiotic environment [61,69] as well as its biotic environment [55,70]. Indeed, much of the work on Bd has focused on environmental drivers of virulence, and shows that systematic environmental changes such as climate change [44,65,71] and altitudinal gradients [39,45] also contribute to variation in Bd virulence in space and time. It remains possible that the pattern of increasing virulence we observe here is driven purely by Bd steadily moving into a more suitable environment in Central America (where that increased suitability is a spatio-temporal effect). More probable, however, is that both evolutionary and environmental factors are at play, with broad evolutionary trends underlying those driven by the environment. Nonetheless, more data and targeted experiments would be required to confirm our contention that Bd has evolved increased virulence during its invasion of Central America. In particular, examining virulence in standardized hosts across longitudinal time series and transects through Bd's invasion history will be of upmost importance.

If emergent pathogens experience important non-equilibrium selection forces as they spread, then we might expect systematic changes in pathogen virulence to be common. That these systematic increases in virulence are rarely observed may say more about our expectations from non-spatial theory (and our resultant sampling designs) than reality. Here, we have simply looked for a pattern of increasing virulence in a messy ad hoc dataset, which imposes severe limitations on our ability to separate environmental and evolutionary effects. Sampling designs that anticipate shifts in virulence in space and time are critical if we are to increase our understanding of the evolutionary dynamics of spreading pathogens. This understanding is important: if true, a general pattern of increasing virulence during pathogen spread could have profound consequences for our understanding of disease dynamics.

Funding statement

We thank the Australian Research Council for support through its fellowships programme. Wayne Mallet at James Cook University's high-performance computing facility provided support for the computationally demanding aspects of the analysis.

Footnotes