Does the early frog catch the worm? Disentangling potential drivers of a parasite age-intensity relationship in tadpoles.

The manner in which parasite intensity and aggregation varies with host age can provide insights into parasite dynamics and help identify potential means of controlling infections in humans and wildlife. A significant challenge is to distinguish among competing mechanistic hypotheses for the relationship between age and parasite intensity or aggregation. Because different mechanisms can generate similar relationships, testing among competing hypotheses can be difficult, particularly in wildlife hosts, and often requires a combination of experimental and model fitting approaches. We used field data, experiments, and model fitting to distinguish among ten plausible drivers of a curvilinear age-intensity relationship and increasing aggregation with host age for echinostome trematode infections of green frogs. We found little support for most of these proposed drivers but did find that the parsimonious explanation for the observed age-intensity relationship was seasonal exposure to echinostomes. The parsimonious explanation for the aggregated distribution of parasites in this host population was heterogeneity in exposure. A predictive model incorporating seasonal exposure indicated that tadpoles hatching early or late in the breeding season should have lower trematode burdens at metamorphosis, particularly with simulated warmer climates. Application of this multi-pronged approach (field surveys, lab experiments, and modeling) to additional parasite-host systems could lead to discovery of general patterns in the drivers of parasite age-intensity and age-distribution relationships.

Introduction

Understanding the factors that drive patterns of changing parasite intensity and aggregation with host age is central to parasitology and disease ecology and has repercussions for understanding the temporal and spatial dynamics of other trophic interactions (Hassell 2000; Raffel et al. 2008; Wilson et al. 2002). In particular, the pattern of changing parasite intensity (the average number of parasites infecting individual hosts, after Hudson and Dobson 1995) with host age has become a fundamental tool of epidemiology, allowing estimation of the parasite reproductive number, R 0, and providing insights into dynamical processes, such as parasite-induced mortality and the ability of hosts to develop acquired immune memory (Anderson and Gordon 1982; Anderson and May 1985; Dietz 1993; Duerr et al. 2003; Farrington et al. 2001; Ferguson et al. 1999; Hudson and Dobson 1995; Wilson 2002). Age–intensity relationships can take a variety of forms, including Type I, a constant increase in mean intensity with age as a consequence of constant infection and no parasite mortality; Type II, an asymptotic increase to a constant mean intensity in older hosts (generally assumed to reflect constant infection and parasite mortality rates); and Type III, a convex relationship in which older hosts have lower mean intensity than hosts in intermediate age classes (Hudson and Dobson 1995). A second important component of the infection process is the age–distribution relationship: the pattern of changing parasite aggregation with host age. Most parasites exhibit aggregated distributions (i.e., most parasites are in a few hosts; Shaw and Dobson 1995; Shaw et al. 1998), and the degree of aggregation can vary with age and infection history, with important implications for designing disease control measures (Barger 1985; Hassell 1980). Determining the mechanisms that drive such ecological patterns are crucial for (1) designing public and veterinary health measures to combat parasitic diseases, (2) deriving basic insights about the biology of host–parasite systems, such as selective pressures on parasite and/or host traits (e.g., timing of breeding; Vandegrift et al. 2008), and (3) predicting impacts of perturbations on host–parasite dynamics (e.g., climate change; Lafferty 2009).

Although age–intensity and age–distribution patterns represent useful concepts and tools in disease ecology, the large number of individual and combined processes that can generate similar types of age–intensity and age–distribution relationships also pose major challenges for identifying causal factors (Duerr et al. 2003; Fulford et al. 1992). For example, mechanisms that can influence the curvature of age–intensity relationships include (1) changes in the rate of exposure to infectious stages through time (i.e., seasonal or interannual variation in parasite abundance), (2) age-related changes in behavioral avoidance of parasites, (3) innate changes in susceptibility with age, (4) acquired immunity following repeated exposure, (5) density-dependent parasite establishment due to competition (e.g., for space in the infected tissue/organ of the host), (6) density-dependent facilitation of infection (i.e., current infection facilitates establishment of future infections), and (7) parasite-induced host mortality. The degree of parasite aggregation can be influenced by three related processes, namely, (8) heterogeneity in exposure to parasites (i.e., individual variation in behaviors or habitat use leading to higher exposure in some individuals than in others), (9) heterogeneity in susceptibility (i.e., individual variation in susceptibility to parasite establishment), and (10) clustering of parasites into discrete groups in the environment.

Each of these proposed mechanisms is predicted to have specific effects on the age–intensity profile and the degree of parasite aggregation, as measured by the variance-to-mean ratio of parasite intensity (Table 1). Nevertheless, researchers seldom test among competing hypotheses to explain the shape of age–intensity curves, especially in wild host populations where experimental manipulations can be challenging (but see Hudson and Dobson 1997; Quinnell 1992). By combining information about the age–intensity and age–distribution relationships, it is sometimes possible to narrow down which processes are the primary contributors to the infection dynamics of a given host–parasite system (Duerr et al. 2003), especially with the benefit of manipulative experiments and quantitative models. Here, we present a case study of how results from field studies, experiments, and models can be integrated to tease apart mechanisms generating a given age–intensity or age–distribution relationship, using echinostome infections of green frog tadpoles as a model system.

Table 1

Predictions of mechanisms proposed to influence the shape of age–intensity (1–7) and age–distribution (4–10) relationships in host–parasite systems

Mechanism no.	Process	Mean intensity	Variance-to-mean ratio
1.	Seasonal or inter-annual variability in exposure	System-specific	No effectb
2.	Behavioral avoidance improves with age	Asymptotic (Type II)	No effect
3.	Innate changes in immunity with age	System-specific	No effect
4.	Acquired immunity following repeated exposure	Asymptotic (Type II)a	Reduced aggregation (<1), decreasing with agea
5.	Density-dependent parasite competition	Asymptotic (Type II)b	Reduced aggregation (<1), decreasing with ageb
6.	Density-dependent parasite facilitation	Exponential increaseb	Aggregated (>1), increasing with ageb
7.	Host mortality due to infection	Asymptotic (Type II) or convex (Type III)b,c	Reduced aggregation, decreasing with ageb,c
8.	Heterogeneity in exposure	No effectb,c	Aggregated (>1), increasing with ageb,c
9.	Heterogeneity in susceptibility	No effectb,c	Aggregated (>1), increasing with ageb,c
10.	Clumped infection	No effectb	Aggregated (>1), decreasing with ageb

aAnderson and May (1985)

bDuerr et al. (2003)

cAnderson and Gordon (1982)

Echinostoma spp. are common trematode parasites whose larvae encyst in the kidneys of larval and adult amphibians. They can cause generalized edema and 30–40% mortality in young (Gosner stage 25) tadpoles due to kidney damage (Holland et al. 2007; Schotthoefer et al. 2003). This parasite has a complex life cycle. It utilizes snails as first intermediate hosts, where the parasite reproduces asexually, releasing the cercarial stage that is infectious to tadpoles (Huffman and Fried 1990). Cercariae enter tadpoles through the urinary ducts and infect the kidneys where they form durable metacercarial cysts that must be ingested by a mammalian or avian predator to mature (Huffman and Fried 1990). These cysts remain visible after the parasite itself dies and can be distinguished from live metacercariae (Martin and Conn 1990), allowing quantification of parasite mortality. These infections can also cause mortality in young tadpoles (Holland et al. 2007; Rohr et al. 2008a; Schotthoefer et al. 2003). Due to the increased abundance of echinostomes and other trematodes in response to urbanization and agrochemicals, echinostomes might play a significant role in limiting amphibian conservation efforts (Rohr et al. 2008b; Skelly et al. 2006).

We determined the shapes of the age–intensity and age–distribution relationships for a naturally occurring echinostome species infecting green frog tadpoles in a Pennsylvania pond. We then used results from experimental infection studies with another echinostome species as well as a model parameterized from laboratory and field experiments to test among the suite of competing mechanisms postulated to drive the observed age–intensity and age–distribution patterns (Table 1). Finally, we used predictive modeling to assess how hatching date and seasonal parasite exposure might influence frog parasite burdens at metamorphosis.

Methods

Animal collection

Green frog tadpoles between stages 25 and 39 (n = 96) were collected on 14 October 2005 from Beaver 1, a landlocked permanent pond in the PA State Game Lands #176 (40°45′52.6″N, 78°0′43.6″W). In an effort to minimize any bias in tadpole collection, we dip-netted parallel to the shoreline while zig-zagging in and out of deeper waters. Tadpoles were held in dechlorinated tap water and fed ground alfalfa pellets (rabbit chow) and fish flakes ad libitum. Echinostoma trivolvis cercariae for the infection experiment (see below) were obtained from naturally infected Planorbella trivolvis snails collected from a pond near Harrisburg, Pennsylvania (40°4′56.0″N, 76°46′2.0″W), which were maintained in bubbled 35-L aquaria containing artificial spring water (ASW, as described by Cohen et al. 1980) and fed boiled lettuce and fish flakes ad libitum. All animals were maintained at room temperature under a constant 12/12-h light/dark cycle.

Experimental infections

Cercarial collection occurred on 24 October 2005 following the procedures described by Kiesecker (2002). Approximately 4 h after emergence from the snails, 20 cercariae were transferred into 84 urine sampling cups (100 mL) containing individual tadpoles and ASW. Twenty hours later, the tadpoles were transferred into larger deli cups containing 600 mL of dechlorinated tap water and fed ground alfalfa and fish flakes ad libitum. Since no visible cercariae remained at this time point, all were assumed to have entered the tadpoles, so that variation in the number of encysted cercariae should reflect variation in physiological resistance to infection by individual tadpoles.

Tadpoles were euthanized 7 days post-exposure by immersion in 1% benzocaine and fixed in 10% buffered formalin. Mass and developmental stage (Gosner 1960) were recorded for each tadpole following fixation in 10% buffered formalin. Twelve additional control tadpoles were examined for naturally occurring trematodes.

Trematode identification and enumeration

Encysted echinostome metacercariae were counted using the clearing and staining procedure described by Hanken and Wassersug (1981). After whole tadpoles had been cleared and stained, the kidneys were dissected and examined more closely to ensure accurate cyst counts. Two distinct morphospecies of encysted echinostome larvae were present in the kidneys of many experimentally infected tadpoles, and these cysts were counted separately for all tadpoles. We have not observed any obvious change in cyst shape or cyst wall thickness of E. trivolvis with time since infection (T. Raffel and J. Rohr, personal observation), and thus we are confident that these two morphospecies are indeed distinct species. In the article, we refer to the naturally occurring echinostome species as “environmental cysts” and to the experimentally induced infections as “experimental cysts”. The life cycle of an experimental cyst was completed as described previously to obtain an adult specimen for specific identification as E. trivolvis (Raffel et al. 2009b). All identifications of trematode metacercariae were performed by researchers blind to the groups to which the tadpoles belonged, and voucher specimens were deposited in the US National Parasite Collection (E. trivolvis adult: USNPC 101920; “Experimental” metacercariae: USNPC 103086, USNPC 103087; “Environmental” metacercariae: USNPC 103088, USNPC 103089). We take advantage of the relatedness and similar location and modes of infection of the two morphospecies and the fact that the environmental cysts were established before the experimental infections to evaluate the level of support for acquired immunity/shared immunity to these infections and negative or positive density-dependence of these infections.

Modeling methods

Tadpole ages were estimated from individual masses using the temperature-dependent tadpole growth model of Berven et al. (1979) and the temperature profile of Beaver 1 pond, as described in detail in Appendix S-1 of the Electronic Supplementary Material (ESM). We then determined that the trematode age–intensity relationship was non-linear (Type II). We next developed a probabilistic model of infection intensity incorporating two potential drivers of the Type II age–intensity profile, i.e., seasonal exposure and stage-dependent susceptibility, as described in detail in ESM Appendix S-2, because most other potential mechanisms had little support based on the experimental results and the shape of the age–intensity relationship (see “Discussion”). Briefly, models with and without seasonal exposure and stage-dependent susceptibility were fitted to the observed age–intensity relationship. Relative seasonal exposure levels (Fig. 2a) were estimated from literature values for seasonal P. trivolvis population dynamics and echinostome infection prevalence (Sapp and Esch 1994; Wetzel and Esch 1996). Lymnaea and Physa spp. snails, the only other potential intermediate hosts in this region (Kostadinova and Gibson 2005), show similar seasonal patterns of abundance and echinostome prevalence, with peak abundances in July and August (Clampitt 1974; Eisenberg 1966; Snyder and Esch 1993). Across multiple years, the seasonality and actual densities of trematode-infected and uninfected P. trivolvis snails have been documented to be strikingly similar (Lemly and Esch 1984), and thus we assume similar seasonality and snail and parasite densities across years in our model. Relative susceptibility of tadpoles to infection at different developmental stages was estimated from the observed relationship between the Gosner (1960) stage and experimental echinostome encystment rates (Fig. 2b). Dates at which tadpoles passed through various stages were estimated based on a modification of Berven et al.’s (1979) temperature-dependent tadpole development model and the Beaver 1 temperature profile. Models were fitted to the observed age–intensity relationship using maximum likelihood and selected using Akaike’s information criterion (AIC). Their relative support was assessed with bootstrapping using Matlab v6.1 (Mathworks, Cambridge, MA; ESM Appendix S-2).

Fig. 2

Seasonal exposure and stage-dependent susceptibility. a Seasonal changes in echinostome exposure rates of tadpoles (i.e., the estimated number of infected snails; filled circles, left axis), calculated as the product of total Planorbella trivolvis snail abundance (open diamonds, left axis) and prevalence of snails producing cercariae (open circles, right axis) (data from Sapp and Esch 1994). b Stage-dependent susceptibility to echinostome encystment as determined by the experimental infection study (mean ± SE). The number of cercariae-producing snails was assumed to be zero from December until March, during which time snails move into deep water or into the substratum (Sapp and Esch 1994), so that any cercariae produced during these months are unlikely to infect tadpoles. The exposure rate for November was estimated by averaging the estimated numbers of infected snails for October and December

To determine how hatching date influences parasite burden at metamorphosis, we constructed a predictive simulation model in R statistical software (R Development Core Team 2006), examining a range of potential temperature profiles that span the latitudinal range of green frogs and predicted climate change scenarios (ESM Appendix S-3). To simulate warmer or cooler climates, temperature was added or subtracted uniformly to the measured seasonal temperature profile of Beaver 1 pond. This model incorporated (1) a degree-day model of time to metamorphosis and (2) seasonal exposure (corresponding to the preferred model based on AIC model selection), making the simplifying assumptions that seasonal exposure rates and the number of degree-days required for metamorphosis are constant across a range of temperatures (ESM Appendix S-3).

Statistical methods

The null model typically utilized in studies of age–intensity curves (in the presence of constant exposure to infectious stages and in the absence of additional factors) is an asymptotic Type II curve based on the assumption of a constant parasite mortality rate which eventually comes to balance the rate of infection (Duerr et al. 2003; Hudson and Dobson 1995). However, parasite mortality appears to have been negligible in this system since we found no dead cysts of either echinostome (but see Holland 2009), and none of our tadpoles were developed enough to have lost cysts via pronephros absorption (Belden 2006); therefore, the null model for the age–intensity curve in this study was a linear Type I relationship for all statistical analyses.

Because smaller sample sizes lead to underestimation of the mean parasite intensity and the degree of aggregation (Gregory and Woolhouse 1993; Lloyd-Smith 2007), effects of age on the intensity of naturally occurring echinostomes were analyzed using a generalized linear model with a negative binomial error distribution and a log link, which fully utilizes the individual-level data (Wilson and Grenfell 1997). Since individual-level analysis was not possible for the analysis of changing aggregation with age, we created age classes for the estimation of variance-to-mean ratios according to sample size, instead of equal age ranges, in order to avoid a decrease in sample size with age. To construct these age classes, we ranked individuals by estimated age and lumped them into groups of approximately four individuals, with minor modifications in this number to accommodate blocks of tadpoles with identical ages. Because the number of experimental cysts was bounded by the constant exposure level of 20 cercariae, these data were converted to proportions and arcsine-square-root transformed, and significant predictors of the experimental echinostome encystment rate were determined with a linear model. For both naturally occurring and experimental echinostomes, curvature in the effect of mass, stage, or age was analyzed by testing for significant quadratic effects of these variables. Significant predictors were determined by assessing their contribution to the final models using likelihood ratio tests (change in deviance) for naturally occurring cysts or analysis of variance for experimental cysts (McCullagh and Nelder 1989). To confirm that curvature in the age–intensity profile was not due to the small sample size of second-year tadpoles, we also conducted a simulation analysis to confirm that the parasite intensities in second-year tadpoles were truly lower than would be expected if parasite intensity increased linearly with age (ESM Appendix S-4). Briefly, this analysis involved extrapolating the age–intensity relationship for the first-year animals out to the second year, assuming a Type I age–intensity relationship, re-sampling from a negative binomial distribution based on this mean and aggregation, and evaluating how many simulations out of 10,000 had quadratic coefficients as negative as in the observed data. Finally, to determine whether the probabilistic models (ESM Appendix S-2) fully explained the curvature of the age–intensity relationship and the unexpectedly low intensity of infection in older tadpoles, we tested for curvature in the residuals from these models and used a bootstrap analysis to test for significant linear trends (ESM Appendix S-5). Curvature in the residuals would indicate that a given model did not fully explain the curvature of the age–intensity relationship.

Results

Wild-caught tadpoles were naturally infected with metacercariae of an echinostome species infecting the mesonephros, which were morphologically distinct from the E. trivolvis metacercariae resulting from experimental infections. These environmental cysts were larger and more spherical, with a thinner cyst membrane, and they had narrower, longer, and lighter staining immature reproductive organs (ESM Fig. S-2). Double-blind examination of encysted metacercariae showed that the 12 control tadpoles were only infected with environmental cysts, verifying our ability to distinguish the two species. This finding also indicates that E. trivolvis was absent or at least at very low abundance in this population of tadpoles, almost certainly less than a mean intensity of 0.32 cysts per tadpole (P < 0.05), based on 10,000 random datasets generated from a negative binomial distribution with the same degree of aggregation as the environmental cysts (k = 0.57). This procedure was repeated over a range of postulated mean intensities to find a value at which there was greater than 95% probability of detecting at least one cyst with a sample size of 12. There was no statistical relationship between the intensities of the two trematode species in individual tadpoles, regardless of which was used as the response variable or what additional covariates were added to the models (all P > 0.5). All metacercariae had fully intact internal structures, revealing that they were alive upon fixation and indicating that there was no detectable mortality of established metacercariae for either echinostome species.

When estimated age was used as a predictor of trematode burden instead of mass, the quadratic term was highly significant (age: coefficient 1.3 × 10−2, X 12 = 8.7, P = 0.003; age2: coefficient −4.4 × 10−5, X 12 = 7.2, P = 0.007). A conservative simulation analysis, based on extrapolation and resampling of the negative binomial distribution for the first-year tadpoles, showed that this degree of curvature in the age–intensity relationship was highly unlikely to arise from the small sample size of second-year tadpoles causing underestimation of their trematode burdens. Indeed, only 2.0% of 10,000 simulations had quadratic coefficients as negative as that observed, despite a conservative estimate of the slope for the age–intensity relationship of first-year tadpoles (ESM Appendix S-4). All other predictors of environmental cyst numbers and their quadratic terms, including mass, Gosner stage, and the number of successfully encysted experimental echinostomes, were non-significant in the final age model (all P > 0.1). The distribution of environmental cysts was significantly different from a Poisson distribution (X 132 = 536.3, P < 0.001) but not from the negative binomial distribution (X 142 = 9.5, P = 0.796; ESM Fig. S-3), with a variance-to-mean ratio greater than one, indicating aggregation (mean 18.9, variance-to-mean 35.2). The variance-to-mean ratio of the environmental cysts increased linearly with age (coefficient 0.083, F 1,24 = 5.0, P = 0.034; Fig. 1a), with no significant quadratic relationship (F 1,23 = 1.2, P = 0.281).

Fig. 1

Age–distribution and age–intensity relationships for naturally occurring echinostomes in a population of green frog (Rana clamitans) tadpoles based on days since hatching. Raw datapoints (small circles) and averages for age classes [large squares, mean ± standard error (SE), n = 4 for all except the largest age class] are shown. a Linear change in the degree of aggregation with age, as estimated by variance-to-mean ratios. b Infection intensity as a function of age, showing the predicted age–intensity relationship for the null model of constant exposure compared to models incorporating stage-dependent susceptibility (thin lines) and seasonal exposure (solid lines)

There was strong evidence for stage-dependent susceptibility to echinostome infection, via a significant quadratic effect of Gosner stage (constant 0.648; stage: coefficient 0.007, F 1,77 = 0.7, P = 0.394; stage2: coefficient −0.007, F 1,77 = 8.1, P = 0.006) on the encystment rate for experimental echinostomes (Fig. 2b). Both early- and late-stage tadpoles had lower encystment rates than middle-stage tadpoles (Fig. 2b). All other predictor variables and their quadratic terms, including mass, age, and number of environmental cysts, were not significant in the final model (all P > 0.1).

We tested for a significant quadratic effect of age on residuals from each of the predictive models to determine how much curvature remained after fitting a given model to the observed age–intensity relationship. These residuals were log-transformed to improve normality as described in ESM Appendix S-5. There was a significant quadratic effect of age on the log-residuals of the null model of constant exposure (F 1,89 = 7.1, P = 0.009; ESM Fig. S-4A), indicating that this model failed to fully explain the Type II age–intensity curve (Fig. 1b). Stage-dependent susceptibility alone provided a slightly better fit to the data than constant exposure (Fig. 1b) and reduced the curvature of the log-residuals, reducing the statistical significance of the quadratic term (F 1,89 = 3.6, P = 0.060; ESM Fig. S-4B). In contrast, both of the models including seasonality of exposure explained most of the curvature of the age–intensity curve, leading to straight-line relationships between age and the log-residuals with no significant quadratic effects (both P > 0.4; ESM Fig. S-4C, D). We used a bootstrap analysis to test for significant trends in the residuals for each of these models (ESM Appendix S-5), revealing significant trends in the residuals for the constant exposure or stage-dependent susceptibility models but no trend for either model incorporating seasonal exposure, indicating that there was no significant signal in the age–intensity pattern beyond that explained by seasonal exposure (ESM Appendix S-5).

Both models incorporating seasonal exposure had AIC values at least 3 units lower than those of the other models considered (Table 2). It is conventional to consider one model clearly preferable to another if its AIC value is lower by 2 or more units (Burnham and Anderson 2002), indicating a strong degree of support for the seasonal exposure models on information–theoretic grounds. Furthermore, the bootstrap analysis selected a model including seasonal exposure in 97.8% of 10,000 randomizations (Table 2). In contrast, models including stage-dependent susceptibility were selected for in fewer than half of the bootstrapped datasets (Table 2), indicating that stage-dependent susceptibility seems to be less influential than seasonality. However, we cannot distinguish between models with seasonality alone and seasonality combined with stage-dependent susceptibility (Table 2). These results provide support for the hypothesis that seasonality of exposure is an important determinant of the observed age–intensity pattern. They also suggest that stage-dependent susceptibility has little effect on the age–intensity relationship, although it probably has other important effects on the ecology of this host–parasite relationship. Additional research would be necessary to more conclusively determine its importance.

Table 2

Model outputs for the null model of constant exposure as compared to models incorporating stage-dependent susceptibility and seasonality in exposure

Model	Parameter c	Parameter K	AIC	∆AIC	Akaike weight	Bootstrap supporta
Null (constant exposure)	0.8248	0.598	711.7	5.1	0.038	0.016
Susceptibility stage-dependent	0.7989	0.606	710.2	3.6	0.080	0.006
Seasonality in exposure	0.4315	0.629	706.6	0	0.483	0.501
Seasonality and susceptibility	0.4208	0.627	707.0	0.4	0.399	0.477

AIC (Akaike information criterion) values are compared to the best-fit model, that incorporating only seasonality in exposure

aProportion of times selected using 10,000 bootstrap datasets

Additional information from experiments described in the published literature is considered along with these results in the following “Discussion” section to address the full range of proposed hypotheses outlined in Table 1.

Predictive modeling incorporating seasonal exposure and estimated time to metamorphosis indicated that tadpoles hatching early (May) or late (August–September) in the green frog breeding season should have lower trematode burdens at metamorphosis than tadpoles hatching in June or July (Fig. 3). Trematode burden at metamorphosis was generally higher for tadpoles hatching near the time of peak exposure (Fig. 2a) or for tadpoles hatching late enough in the season to require an extra year to metamorphose (Fig. 3, ESM Fig. S-5). Assuming no change in the seasonality of exposure with latitude, the proportional difference in cumulative exposure between early/late tadpoles and those hatching in mid-season was more pronounced with simulated warming/decreased latitude (Fig. 3), particularly with uniform addition of at least 2°C so that tadpoles could metamorphose in less than 1 year (ESM Fig. S-5). In Fig. 3, the +2°C line is nearly flat because with this temperature profile tadpoles took almost exactly 1 year to metamorphose (ESM Fig. S-5), thus all tadpoles received the same total exposure regardless of hatching date.

Fig. 3

Predicted trematode burden at metamorphosis for a green frog (Rana clamitans) tadpole hatched on a given date, based on a model incorporating seasonal exposure and using a degree-day model to estimate time to metamorphosis (ESM Fig. S-5). The model output using the Beaver 1 pond temperature profile (0°C) is indicated by a bold solid line, and narrower (solid, dashed or dotted) lines indicate model outputs given simulated warming (+1 to +6°C) or cooling (−1 to −3°C). Gray shading indicates the period outside the known green frog breeding season, when green frog tadpoles are unlikely to hatch because hatching occurs 3–6 days after breeding (Schalk et al. 2002)

Discussion

Age–intensity relationship

Despite there being several mechanisms that could drive the observed asymptotic age–intensity relationship between green frogs and echinostomes (see Table 1), our study revealed that a single mechanism, seasonal exposure, was the parsimonious driver of this age–intensity relationship. Changes in susceptibility with tadpole development were detected in this study, with intermediate-stage tadpoles being more susceptible to experimental echinostome infections than early- or late-stage tadpoles. This stage-dependent susceptibility to infection might be important to other aspects of the host–parasite relationship (e.g., Raffel et al. 2010), but it added little explanatory power to our mechanistic models of the age–intensity relationship. The observed convex relationship between the stage of tadpoles and their susceptibility to echinostomes was also apparent in a previous study on echinostomes and green frogs (Holland et al. 2007), although not explicitly discussed by these authors. Lower encystment rates in late-stage tadpoles were also observed in studies of Rana pipiens and Bufo americanus (Raffel et al. 2010; Schotthoefer et al. 2003), suggesting that this might be a general pattern for tadpole susceptibility to E. trivolvis. This pattern might be due to improvements in immune function or changes in kidney morphology associated with later tadpole development (Flajnik et al. 1987; Raffel et al. 2010; Rollins-Smith 1998), whereas the lower susceptibility of early-stage green frog tadpoles might be due to the absence of well-developed mesonephral kidneys. The mesonephral kidneys begin to develop late in Gosner stage 25 and provide more area for echinostome encystment than the smaller pronephral kidneys that tadpoles have upon hatching (Holland et al. 2007).

Although the observed asymptotic age–intensity relationship can be explained by seasonal exposure alone, this does not necessarily demonstrate that seasonality was the only causal factor. However, based on experimental results and the observed age–distribution relationship, there is little support for alternative hypotheses, such as density-dependent parasite competition, acquired immunity, and parasite-induced host mortality (Table 1). Density-dependent parasite competition and acquired immunity following repeated parasite exposure both act to reduce the degree of parasite aggregation, so both hypotheses predict that the variance-to-mean ratio will decrease with age, which is in direct contrast to what we observed. Additional evidence against the hypothesis of density-dependent competition is provided by experimental infection studies showing that the echinostome encystment rate was independent of the number of cercariae to which ranid tadpoles were exposed (Holland et al. 2007; Schotthoefer et al. 2003), except for very young tadpoles (Gosner stage 25–26) prior to development of the mesonephros (Schotthoefer et al. 2003). Both studies included treatments that exposed tadpoles to higher numbers of echinostomes than were observed in any of our wild-caught tadpoles. This apparent absence of density-dependent infection in older tadpoles seems to rule out this potential driver of the asymptotic age–intensity relationship.

The ability of tadpoles to acquire enhanced resistance to helminth infection following repeated exposure has not, to our knowledge, been tested experimentally. However, given the incomplete allograft responsiveness and restricted T and B lymphocyte diversities of tadpoles compared to adult amphibians (Flajnik et al. 1987; Rollins-Smith 1998), not to mention the complete lack of significant cross-immunity between the two echinostome species in this study (as indicated by the lack of a negative correlation between the two cyst types), it seems reasonable to conclude that acquired immunity plays at most a relatively minor role in tadpole parasite dynamics. This conclusion is consistent with a study of red-spotted newt parasites that found no signature of acquired immunity in age–intensity relationships of encysted helminths (Raffel et al. 2009a).

Another possible mechanism for the asymptotic age–intensity curve is host mortality due to infection. Parasite-induced host mortality has the potential to generate asymptotic age–intensity relationships for highly aggregated persistent parasites, such as the echinostomes in this study, because individuals with high parasite intensities would be expected to die at a higher rate, leading to the absence of heavily infected individuals in older age classes (Duerr et al. 2003; Hudson and Dobson 1995). This is unlikely to be a relevant mechanism for our system because echinostomes would need to cause significant mortality in older, more heavily infected tadpoles in order to influence the overall age–intensity curve, and experiments so far have found no evidence of mortality due to echinostome infection in any but the youngest (Gosner stage 25) ranid tadpoles, even with exposure levels comparable to the highest parasite burdens observed in this study (Belden 2006; Holland et al. 2007; Schotthoefer et al. 2003). Although these studies were conducted in the laboratory where many mortality factors that can interact with infections are limited (e.g., predation), field patterns also do not support parasite-induced host mortality as a considerable driver of the age–intensity relationship. There was no decrease in the variance-to-mean ratio for older tadpoles collected from the field, as predicted for over-dispersed parasites that induce host mortality (Anderson and Gordon 1982; Duerr et al. 2003). It is theoretically possible that an unobserved increase in parasite intensity or aggregation occurred during the gap between the first- and second-year tadpoles (Fig. 1), followed by a decrease to the levels observed in the second-year tadpoles. However, this seems unlikely since the gap in sampled tadpole ages is due to the gap in amphibian breeding during winter (i.e., no tadpoles of this age range existed to be sampled), which coincided with a period of virtually no exposure to cercariae (Fig. 2a). Hence, while field data on echinostome-induced host mortality are lacking, there is presently limited support for this hypothesis as an explanation of the observed age–intensity relationship.

One final alternative explanation for the asymptotic age–intensity curve, that tadpoles improve their ability to avoid exposure with age by changing their behavior patterns, cannot be dismissed because this hypothesis was not tested experimentally. This hypothesis is feasible because ranid tadpoles are known to change microhabitat use and to increase activity levels as they get older (Golden et al. 2001), tadpoles change activity levels in response to cercariae (Rohr et al. 2009), and activity at least appears to reduce tadpole probability of infection by echinostome cercariae (Koprivnikar et al. 2006; Thiemann and Wassersug 2000).

Although seasonal exposure appears to be the parsimonious driver of this age–intensity relationship, this conclusion depends on a several important assumptions. First, in estimating tadpole ages, we assumed that rates of temperature-dependent growth and development were similar to those observed by Berven et al. (1979) for this frog species (ESM Appendix S-1). Relaxing this assumption would have required more information about factors influencing the rates of tadpole growth and development in this pond, such as predation or competition (Raffel et al. 2010). Even if relaxing this assumption had been possible, however, we doubt it would have changed our conclusions regarding the age–intensity relationship because the cut-off between first- and second-year tadpoles, which was the strongest determinant of the shape of the age–intensity relationship, was robust to changes in model parameters (ESM Appendix S-1). Furthermore, there was a gap in masses—but not in developmental stages—between tadpoles identified as first- and second-year individuals (ESM Appendix S-1), which is best explained by the cessation of tadpole development, but not growth, during the winter (Crawshaw et al. 1992). For these reasons, we feel confident in our identification of the larger individuals as second-year tadpoles and in our characterization of the shape of the age–intensity relationship. Second, in estimating seasonal exposure, we assumed that seasonal changes in snail abundance and infection prevalence were similar to those observed for P. trivolvis by previous authors (Sapp and Esch 1994; Wetzel and Esch 1996). We think this assumption is reasonable because observed patterns of seasonal snail abundance and echinostome infection prevalence are broadly similar among multiple populations and species of snails in the northeastern USA (Clampitt 1974; Eisenberg 1966; Lemly and Esch 1984; Sapp and Esch 1994; Snyder and Esch 1993; Wetzel and Esch 1996). Furthermore, small changes in the pattern of assumed seasonal exposure, such as a shift in peak exposure by a month forward or backward, would not have greatly influenced the fit of the seasonal exposure model.

Age–distribution relationship

Analysis of the age–distribution relationship provides further opportunity to discriminate among possible mechanisms in our system. We proposed seven mechanisms that might influence the age–distribution relationship (Table 1), of which only four are capable of generating the observed pattern of increasing aggregation with age. Of these, heterogeneity in exposure, meaning consistent differences in exposure rates for individual tadpoles, is the only candidate mechanism that is consistent with our data and other observations, predicting both the aggregated parasite distribution and the increase in aggregation with age (Duerr et al. 2003; Quinnell et al. 1995). Heterogeneity in exposure seems probable in the tadpole–echinostome system given the highly aggregated distributions of snails, and thus cercarial production in ponds (Sapp and Esch 1994, T. Raffel and J. Rohr, personal observation), and the high degree of site fidelity exhibited by amphibians within ponds (Bellis 1968). Hence, tadpoles living and foraging in one part of the pond will have a different exposure history from tadpoles living in other areas. Furthermore, amphibians exhibit significant individual-level variation in activity (Rohr and Madison 2003; Smith and Doupnik 2005), which has been found to influence echinostome infection rates in laboratory studies (Koprivnikar et al. 2006; Thiemann and Wassersug 2000). Aggregated spatial patterns of infectious stages and heterogeneity in host behavior have both been shown to be capable of generating high degrees of aggregation in laboratory host–parasite systems (Anderson et al. 1978; Keymer and Anderson 1979).

The other three potential drivers of an increase in aggregation with age, namely, density-dependent facilitation, clumped infection, and heterogeneities in susceptibility (Table 1), are unlikely mechanisms for the observed age–distribution relationship. Density-dependent facilitation of parasite establishment, meaning that current infections make it easier for additional echinostomes to successfully infect the same tadpole, predicts an exponential increase in parasite intensity with age (Duerr et al. 2003). This prediction is contrary to our results; furthermore, Holland et al. (2007) found no evidence of density-dependence in experimental infections of green frogs with echinostomes when tadpoles were infected with varying doses of cercariae. Clumped infection, i.e., grouping of infectious stages into discrete clusters (Duerr et al. 2003), seems unlikely given that cercariae appear to swim independently of each other once released from snails (T. Raffel and J. Rohr, personal observation). Furthermore, clumped infection predicts that the degree of aggregation decreases with age (Duerr et al. 2003), which is opposite to the pattern we observed.

Heterogeneity in susceptibility did not appear to be an important driver of parasite aggregation in this study. If similar mechanisms are used to resist infection by both echinostome species in this study, which seems reasonable given their shared route of infection and tissue tropism, then individual heterogeneity in susceptibility would predict higher encystment rates of experimental echinostomes in tadpoles with more environmental cysts. The lack of any relationship between the two cyst types suggests that heterogeneity in susceptibility did not contribute substantially to the age–distribution pattern of the environmental cysts.

In summary, seasonal variation in population-wide exposure rates and individual heterogeneity in exposure appear to be the parsimonious drivers of the age–intensity and age–distribution relationships for echinostomes in this tadpole population. While a number of the other six proposed mechanisms might influence the age–intensity relationship to a small extent, their effects were not strong enough to be detected in this study. Amphibian life-cycles are often products of responses to seasonal changes in environmental conditions (e.g., Rohr et al. 2002, 2003), and parasites too are often highly seasonal (Altizer et al. 2006; Raffel 2006). As a result, it seems probable that many age–intensity relationships are driven largely by seasonality in parasite exposure, particularly for host species that are exposed to parasites for only one or a few years.

Parasites, host life-history, and climate

Parasites can have important selective effects on host reproductive decisions, including seasonal breeding strategies (Raffel 2006; Vandegrift et al. 2008). For example, gray tree frogs (Hyla versicolor) select oviposition sites according to the risk of trematode infection to their tadpole offspring (Kiesecker and Skelly 2000). Green frogs might have evolved similar responses to seasonal patterns of trematode infection. Our model predicted that tadpoles hatching either early (May) or late (August–September) in the green frog breeding season should have lower echinostome burdens at metamorphosis (Fig. 3). This conclusion was robust to changes in parameter values which, when varied, had similar effects to raising or lowering the temperature. This difference in trematode risk should provide a selective pressure for frogs to breed early or late in the breeding season, because offspring with lower echinostome exposure as tadpoles are more likely to survive and reproduce. The difference in trematode burdens between early/late and mid-season was greater with simulated warming, suggesting that this selection pressure for a dual breeding season might be stronger for frog populations at lower elevations or latitudes. Interestingly, Berven et al. (1979) observed dual peaks in breeding activity (mid-May and early August) in their lowland green frog populations but not in montane populations, not unlike dual peaks in egg-laying observed for the Beaver 1 green frog population by J.R. Rohr in 2007 (personal observation). Of course parasites are only one of many selection pressures that might influence green frog reproductive decisions, and warmer climatic conditions will probably influence the phenology of green frog breeding or the time to metamorphosis in ways unrelated to parasitism.

These results have potential implications for effects of climate change on amphibian parasite dynamics. The predictive model suggests that warming should lead to stronger selection pressure for a dual breeding season, along with lower trematode burdens at metamorphosis due to increased frog developmental rates. However, the latter effect might easily be offset by changes in snail or trematode population dynamics, such as increased cercarial shedding rates by snails in warmer conditions, as proposed by Poulin (2006). In addition, green frogs can adapt to cooler local climatic conditions by decreasing the number of degree-days it takes to metamorphose (parameter K in the model; Berven et al. 1979) and relaxing the assumption of an invariant K would reduce the magnitude (though not the seasonality) of the predicted temperature effects. Furthermore, our model does not account for potential changes in seasonal exposure, which could be influenced by phenological shifts in snail breeding, trematode egg inputs from vertebrate definitive hosts, or cercarial shedding. To accurately predict effects of climate change on overall tadpole exposure to trematodes, future researchers should incorporate responses of snails and trematodes to temperature in addition to effects of local adaptation on tadpole developmental rates at particular temperatures.

Conclusions

This study is among only a few to have assessed various drivers of age–intensity and age–distribution relationships in natural host–parasite systems. However, the combination of field, experimental, and modeling approaches described here should be applicable to other host–parasite systems, especially those in which researchers can estimate host age and conduct experimental infections under controlled conditions. Perhaps by using this powerful approach we will finally be able to elucidate general patterns in the drivers of parasite age–intensity and age–distribution relationships, improving our ability to design control measures for parasitic diseases and to predict effects of perturbations such as climate change.

Below is the link to the electronic supplementary material.

Supplementary material 1 (PDF 1984 kb)