Introduction
Haemosporidian parasites (Plasmodium and related genera) exploit thousands of vertebrate host species including mammals, squamate reptiles and birds on all warm continents (Galen et al., Reference Galen, Borner, Martinsen, Schaer, Austin, West and Perkins2018). A vexing issue in the study of human malaria focuses on prevalence patterns and their variation over time and space (Hackett, Reference Hackett1937). A recent study has examined >50 000 prevalence surveys of Plasmodium falciparum in sub-Saharan Africa since 1900 (Snow et al., Reference Snow, Sartorius, Kyalo, Maina, Amratia, Mundia, Bejon and Noor2017), and found no comprehensive explanation for changes in parasite abundance. Understanding the interplay between the parasites, their hosts and the abiotic conditions of malaria infection has implications for human public health, but also for wildlife conservation. For example, a single strain of P. relictum caused sharp population declines of endemic birds in the Hawaiian islands (Beadell et al., Reference Beadell, Ishtiaq, Covas, Melo, Warren, Atkinson, Bensch, Graves, Jhala, Peirce, Rahmani, Fonseca and Fleischer2006). Also, climate change has now brought the parasites to northern latitudes and to bird species previously not infected by Plasmodium (Martinsen et al., Reference Martinsen, Sidor, Flint, Cooley and Pokras2017). One pointed question focuses on the stability in prevalence patterns; that is, are the dynamics of these parasite–host systems stable through time and robust to environmental fluctuations (MacDonald, Reference Macdonald1952)? Classic theoretical models on malaria predict temporally stable dynamics (review in Aron and May, Reference Aron and May1982). Large environmental changes, though, may well disrupt stable endemic patterns to unpredictable outcomes (Chiyaka et al., Reference Chiyaka, Tatem, Cohen, Gething, Johnston, Gosling, Laxminarayan, Hay and Smith2013).
With ongoing climate changes in temperature and rainfall patterns, could the prevalence of malaria parasites be driven into unstable patterns? Malaria parasites are vector-borne, and thus sensitive to environmental fluctuations (Harvell et al., Reference Harvell, Mitchell, Ward, Altizer, Dobson, Ostfeld and Samuel2002; Altizer et al., Reference Altizer, Ostfeld, Johnson, Kutz and Harvell2013), particularly through the effect of environmental variation on various stages of the vector and parasite life cycle (LaPointe et al., Reference LaPointe, Goff and Atkinson2010; Paaijmans et al., Reference Paaijmans, Blanford, Bell, Blanford, Read and Thomas2010; Mordecai et al., Reference Mordecai, Paaijmans, Johnson, Balzer, Ben-Horin, de Moor, McNally, Pawar, Ryan, Smith and Lafferty2013). There is little consensus on the expected influence of climate change on the dynamics of vector-borne parasites including Plasmodium (Rohr et al., Reference Rohr, Dobson, Johnson, Kilpatrick, Paull, Raffel, Ruiz-Moreno and Thomas2011). Expected temperature increases of 1–5 °C may result in expanded malaria geographic ranges (e.g. Martinsen et al., Reference Martinsen, Sidor, Flint, Cooley and Pokras2017), increased vector biting and parasite replication rates promoting higher transmission rates (Pascual et al., Reference Pascual, Ahumada, Chaves, Rodó and Bouma2006). However, an increase in temperature may also increase vector mortality. Therefore, the relationship between abiotic factors, disease transmission and its ecological consequences remains an important open question in disease ecology and parasitology (Lively et al., Reference Lively, de Roode, Duffy, Graham and Koskella2014).
Temporal instability could push a cascade of subtle consequences, including changes in the composition of coexisting Plasmodium species in the same or in closely related sympatric hosts, possibly affecting parasite virulence (Elliot et al., Reference Elliot, Blanford and Thomas2002; Wolinska and King, Reference Wolinska and King2009). These issues are complex and are best probed with a long-term view, analysing long data series that could allow an understanding of climate-driven changes in parasite dynamics and their consequences (Clutton-Brock and Sheldon, Reference Clutton-Brock and Sheldon2010). Most long-term studies on malaria parasites focus on those of human medical importance and are influenced by public health efforts that result in changing environments for the parasites, superimposed on natural fluctuations. Long-term studies of malaria parasites in wildlife hosts present a cleaner alternative. For example, a study characterizing bird haemosporidian parasite lineages at a millennial time frame found indications of genetic lineage switches at this time scale on Caribbean islands (Soares et al., Reference Soares, Latta and Ricklefs2017). At a shorter time frame, a 17-year study of haemosporidian parasites infecting great reed warblers (Acrocephalus arundinaceus) found temporally stable composition of parasite lineages and no significant changes to host fitness (Bensch et al., Reference Bensch, Waldenström, Jonzén, Westerdahl, Hansson, Sejberg and Hasselquist2017). At northern latitudes, a 10-year study showed that Plasmodium infection of black-capped chickadees is temporally variable and associated with keratin disorder (Wilkinson et al., Reference Wilkinson, Handel, Van Hemert, Loiseau and Sehgal2016). Thus, these long-term studies focusing on avian haemosporidian systems provide invaluable insights into the dynamics of host–parasite systems, but are scarce in other hosts (but see Schall and Denis, Reference Schall and Denis2013).
Plasmodium infection in lizards offers a useful system to understand the interplay between environmental change and the stability of host–parasite interactions. Lizards are ectothermic hosts with life spans that allow studying multiple generations in a short period of time. Here we present a long-term study of three lizard malaria parasites (P. azurophilum, P. floridense and P. leucocytica) infecting Anolis gundlachi in the tropical rainforest of Puerto Rico. We leverage data assembled from 1990 to 1999 (Schall et al., Reference Schall, Pearson and Perkins2000), add extensive recent sampling and reanalyse the data using statistical methods appropriate to a longer-term approach. There are several notable features for this study: the lizard life span is typically 1 year or shorter, and thus many host generations are followed in the study. The system has not been disturbed by deliberate human intervention, such as parasite control efforts, logging or human-induced land-cover change during the study period. The parasite assemblage includes three species exploiting the same host, and the dynamics of all three were followed over time. Changes in relative proportion of the Plasmodium species would be particularly interesting because P. leucocytica infects several classes of white blood cells and the other two infect erythrocytes, with possible complex competitive interactions (Perkins, Reference Perkins2000). Last, a well-defined measure of body condition allows a measure of changes in host health condition over time. We seek resolution of four important questions: (1) was the parasite prevalence and risk of infection stable during the long-term study period? (2) If not, was temporal variability driven by abiotic factors (e.g. temperature and precipitation) or were there secular trends apparently independent of environmental changes [resembling the ‘autonomous’ pattern of human malaria suggested by Hackett (Reference Hackett1937)]? (3) Was parasite species composition stable during this long-term period, or were there signs of shifts on probability of infection among species? (4) Do abiotic factors mediate the potential negative consequences (virulence) of infection to the host? This system thus offers a rare view of a malaria system in a nonhuman host – over >15 generations of the vertebrate host.
Materials and methods
Study system, field sampling and diagnostics
We studied the lizard A. gundlachi (Fig. 1a), and three species of Plasmodium parasites: P. azurophilum (Fig. 1b), P. floridense (Fig. 1c) and P. leucocytica (Fig. 1d) at the El Verde field station located at the Luquillo Experimental Forest, in Puerto Rico (central point 18°19.263′N–65°49.146′W). Even though P. azurophilum infects red blood cells and P. leucocytica white blood cells, these species are closely related (sister taxa) in the most recent published phylogeny of the Plasmodium group (Galen et al., Reference Galen, Borner, Martinsen, Schaer, Austin, West and Perkins2018). In contrast, P. floridense is more related to a group including P. mexicanum and P. chiricahue, which to date have not been found in Puerto Rico. The host, A. gundlachi, is a medium-sized lizard (mean snout–vent length (SVL) of 53.5 mm, mean mass 3.89 g) and is the most common anole in the forest understory, reaching population sizes of 2000 ha−1 (Reagan, Reference Reagan1992). This anole is among seven other species living at the site, but the other anoles are rarely infected (Schall et al., Reference Schall, Pearson and Perkins2000). We sampled lizards during 15 periods from 1990 to 2017: summers (May–August) 1990, 1996, 1997, 1998, 2015, 2016, 2017, and winters (January–March) 1991, 1997, 1998, 1999, 2001, 2002, 2016 and 2017.
(A) The host Anolis gundlachi is infected by three malaria parasites. The figure shows examples of Giemsa-stained blood smears of (B) Plasmodium azurophilum, (C) P. floridense and (D) P. leucocytica.
To ensure consistency over time we replicated rigorously the protocols of field sampling and diagnostics used in the original study by Schall et al. (Reference Schall, Pearson and Perkins2000). Specifically, we consulted sketch maps of study sites made during the early years to sample the same sites within the forest (between 4 and 10 trails each year) and collected similar sample sizes (N > 650 per survey). The experimental forest is well mapped, including a 16-ha plot with every tree identified and numbered, facilitating return to the earlier collection sites. To ensure uniform scoring of blood films, LO and MAA worked with the initiator of the early samples (JJS) with a dual-viewing microscope to become proficient in the scoring of infected lizards and identification of the parasites.
We searched for anoles on tree trunks, branches, rocks, ground and other perches. We captured lizards by hand or noose and kept them in individual bags to transport them to the laboratory where we determined sex, measured SVL, mass and sampled blood using toe clipping (Schall et al., Reference Schall, Pearson and Perkins2000). Toe clipping also ensures that individuals are not sampled more than once in a season. Captured lizards were released within 24 h after capture in the same areas where they were collected. We made thin blood smears for each individual and kept dried blood samples on filter papers for future molecular analyses. We fixed the smears using methanol (100%) and stained the slides using Giemsa stain at pH 7.5 for 50 min following the protocol by Schall et al. (Reference Schall, Pearson and Perkins2000). To determine if a lizard was infected we diagnosed each sample by scanning stained blood smears at 1000× for 6–10 min and identified parasite species based on morphological traits and cell class infected (Telford, Reference Telford2016). To estimate host body condition, we determined the SVL and mass of each lizard.
The 15 surveys included in these analyses comprise the longest lizard–malaria study in the tropics. These surveys took place from 1990 to 2002 and additional surveys in 2015–2017, including a gap from 2002 to 2015. The data, though, are appropriate to answer our stated questions. First, each survey is an independent sample because the average life span of these lizards is ~1 year. Second, the 15 surveys sample the full range of environmental variables including years under El Niño Southern Oscillation events (extreme droughts on one extreme, and high rainfall due to hurricanes on the other) and mild years that represent average conditions. Third, the methods were rigorously consistent among samples, which prevents any potential bias due to observational variability. Therefore, the surveys are a thorough representation of the environmental conditions experienced by this host–parasite system in the long-term period.
Environmental variables
We compiled temperature and rainfall data from El Verde weather station records and estimated daily mean minimum, maximum, cumulative (rain) and variance of registered temperatures and rainfall 30 and 120 days prior to the sampling month (Schall et al., Reference Schall, Pearson and Perkins2000). Hurricane events occurred in September 1995 (two), July 1996 and September 1998. While Hurricane Maria affected Puerto Rico in September 2017, surveys took place in June and July. Therefore, the effect of the storm is not captured in this dataset. We did not include hurricanes per se in the analysis, but high rainfall marked those periods in the precipitation data.
Analyses
We followed Chamberlin's (Reference Chamberlin1965) approach of comparing multiple working hypotheses for each one of our four research questions. This approach contrasts with the more common Popperian method in which a single alternative hypothesis is compared against a null one. The key difference between the two approaches is that the multiple working hypotheses framework allows for the possibility that more than one hypothesis may be simultaneously true (Elliott and Brook, Reference Elliott and Brook2007; Betini et al., Reference Betini, Avgar and Fryxell2017). Here we operationalize this approach developing a list of a priori hypotheses for each of the questions, which are represented by a model. We compared each model – including a null (intercept only) model – using Akaike Information Criterion adjusted for small sample sizes (AICc). The most parsimonious model has the lowest AICc score. In the case of multiple models having similar low AICc scores (<2 AICc units), we use AIC weights to assess the relative contribution of each hypothesis to explain the observed relationship. If two or more models had similar AICc scores we averaged the top models by weighting each shared parameter estimate by the model's AICc weight (Burnham and Anderson, Reference Burnham and Anderson2004). We fitted all models using maximum likelihood and conducted model selection and averaging using the AICcmodavg package in R 3.0.3 (R Core Team, 2013).
To determine if the individual probability of infection was stable through the long-term study period, we compared eight binomial models. These models (with exception of the null model) include season, sex and/or SVL as controlling covariates because Schall et al. (Reference Schall, Pearson and Perkins2000) found strong evidence showing that the probability of infection was higher in the summer than in the winter season. They also found that bigger males had a higher probability of infection than smaller males or females. We compared models predicting the probability of infection as a function of (1) SVL, year and sex separately, (2) the additive effect of sex, SVL and season, (3) a similar model, but including an additive effect of year, (4) a model including sex, SVL and season with an interactive effect between sex and SVL, (5) a similar model including an additive effect of year and (6) an intercept-only model that represented the null hypothesis that none of the variables explains the observed temporal variability in probability of infection (Supplementary Table S1). If there was significant temporal variation in the dynamics, the most parsimonious model (lowest AICc score) would include the effect of sampling year.
We followed a similar model selection approach to test which environmental variables (i.e. rainfall and/or temperature) better predict the proportion of infected individuals through time. We compared 35 models that included various combinations of mean rainfall 30 days and 120 days before the field sampling (Schall et al., Reference Schall, Pearson and Perkins2000). These also included models with the variance of rainfall or temperature 30 or 120 days before the field sampling (Vasseur et al., Reference Vasseur, DeLong, Gilbert, Greig, Harley, McCann, Savage, Tunney and O'Connor2014). To allow for a possible non-linear effect we also fitted individual and additive models with a second-degree polynomial. Last, we fitted a null (intercept only) model to describe the case where none of the tested variables better explains the observed patterns in prevalence (Supplementary Table S3).
To assess the relationship between environmental factors, infection status and body condition we estimated body condition using the residual index R i (Cox and Calsbeek, Reference Cox and Calsbeek2015). We calculated this index using the residuals of the linear regression of log10 mass on log10 SVL. Lizards with positive residuals are heavier than average (better body condition), while lizards with negative residuals are skinnier relative to their SVL than average (Schall and Pearson, Reference Schall and Pearson2000). We made a separate analysis for each sex, including data from the years for which body mass and SVL data were available (1996, 1997, 1998, 2015, 2016 and 2017). Previous studies show that A. gundlachi lizards are particularly sensitive to maximum environmental temperatures (Huey and Webster, Reference Huey and Webster1976; Hertz et al., Reference Hertz, Huey and Stevenson1993) and cumulative rain (Schall and Pearson, Reference Schall and Pearson2000). Therefore, we compared 14 models that predicted variability in body condition as a function of maximum temperature 30 days of the census, cumulative rain 6 months before the survey, infection state and their additive combinations. To account for potential non-linear effects, we also compared models that incorporated a second-degree polynomial effect of temperature and rainfall. We also compared a null (intercept-only) model that represents the case where none of these variables explains variability in body condition (Supplementary Tables S5 and S7).
To assess changes in composition of the three Plasmodium species through time we compared 10 multinomial logit models in their ability to predict the individual probability of being non-infected or infected by one of three Plasmodium parasites (four categories; co-infections were not considered because they were infrequent). This modelling approach is an extension of a logistic regression for multinomial response variables. We included sex and SVL in the models as controlling variables because these may influence the probability of infection by different Plasmodium species. We compared models including (1) the single effect of sex, SVL or year; (2) models considering the additive and interactive effect of year and SVL; (3) a model considering the additive effect of sex and SVL and (4) the additive effect of these three variables, with either interactions of year with SVL, or sex. We also fitted a null (intercept-only) model that represents the case where none of these variables explains the probability of getting infected by any of the three Plasmodium parasites (Supplementary Table S10). The models were fitted using maximum likelihood applying the nnet package in R (Venables and Ripley, Reference Venables and Ripley2013).
Results
Long-term dynamics in probability of infection
A total of 8991 A. gundlachi lizards were sampled in 15 surveys from 1990 to 2017. The most parsimonious model explaining variability in the individual probability of infection included the additive effect of sex, SVL, season, year and the interaction between sex and SVL (ΔAICc = 10.31 between this model and the next best, AICc weight = 0.99; Supplementary Table S1). In the summer, the probability of infection and prevalence were relatively stable from 1990 to 1998 (Supplementary Fig. S1) with an estimated individual probability of infection fluctuating between 0.28 and 0.38 for males and between 0.17 and 0.27 for females. The probability of infection decreased during the 2015–2017 period to 0.10–0.19 for males, and 0.06–0.12 for females (Fig. 2A, Supplementary Table S2). In the winter, the probability of infection and prevalence was low in 1991 (Supplementary Fig. S1; average probability of infection of 0.15 in males, and 0.09 in females). Then it increased and remained relatively stable from 1997 to 2002 with an individual probability of infection fluctuating between 0.28 and 0.38 in males and between 0.19 and 0.26 in females. This apparent stability was disrupted in 2016–2017 where the probability of infection decreased (average probability of infection of 0.11–0.14 in males, and 0.07–0.08 in females; Fig. 2B, Supplementary Fig. S1 and Table S2).
The predicted partial relationship between probability of infection of A. gundlachi by three malarial parasites and time. This probability of infection is relatively constant in the early years but decreases in 2015–2017. Points represent mean estimates and the bars represent 95% CIs.
Environmental drivers of temporal variability
The most parsimonious model explaining temporal variations in the proportion of infected individuals through time included the additive quadratic effect of mean maximum daily temperature and mean daily rainfall through 120 days before the sampling (ΔAICc = 25.72 between this model and the next best, AICc weight = 1; Supplementary Table S3). The probability of infection followed a nonlinear response with temperature and rainfall. The proportion of infected individuals increased with increasing temperature and maximized at ~26 °C after which the relationship decelerates (Fig. 3A, Supplementary Table S4). The proportion of infected individuals maximized at an intermediate measures of rainfall (~10 mm; Fig. 3B, Supplementary Table S4).
Predictions of the most parsimonious model explaining temporal variation in the proportion of infected A. gundlachi by malaria parasites. This model predicted the proportion of infected individuals by the additive effect of squared (A) mean maximum temperature and (B) mean daily rainfall 120 days before sampling. Lines represent estimated means and shaded areas 95% CIs.
Environment-mediated effects of infection on body condition
The most parsimonious model explaining the relationship between body condition (R i) of females and abiotic variables included a quadratic additive effect of mean daily maximum temperatures and cumulative rain 1 month before the survey (ΔAICc = 0.49 between this model and the next best, AICc weight = 0.56; Supplementary Table S5). Yet, the following best model, which included the same variables in addition to the infection state, had also similar explanatory power (AICc weight = 0.44). Therefore, we conducted model averaging between these two models (Supplementary Table S6). This average model predicts a non-linear partial relationship between body condition and temperature where body condition is maximized at ~26 °C (Fig. 4A). This model also predicts a slow decrease in body condition with increasing cumulative rainfall (Fig. 4B, Supplementary Table S6). While average body condition was higher in non-infected females, the confidence intervals (CIs) overlapped suggesting lack of difference between the two groups. The most parsimonious model explaining these relationships for males included the effect of cumulative rain a month before the survey (ΔAICc = 0.34 between this model and the next best, AICc weight = 0.47; Supplementary Table S7). The following best model had a similar explanatory power included the same variable in addition to the infection state (AICc weight = 0.39). Therefore we conducted model averaging between these two models (Supplementary Table S8). This average model predicted that male's body condition maximized at intermediate levels of cumulative rain (~152 mm of rain in a month; Fig. 5). While average body condition was predicted higher for infected individuals, the CIs overlapped with those predicted for non-infected individuals suggesting lack of difference in body condition between the two groups (Fig. 5).
Partial predictions from the most parsimonious model describing the relationship between A. gundlachi females body condition and (A) mean daily maximum temperature, and (B) cumulative rain during the month previous to the survey. Red shaded area represents 95% CI for infected individuals, the blue shaded area represents 95% CI for the non-infected individuals and the grey shaded area represents overlap between them. This area of overlap indicates that the predicted relationship was similar between infected and non-infected females.
Partial predictions from the most parsimonious model describing the relationship between A. gundlachi males’ body condition and cumulative rain during the month previous to the survey. Red shaded area represents 95% CI for infected individuals, the blue shaded area represents 95% CI for the non-infected individuals and the grey shaded area represents overlap between them. This area of overlap indicates that the predicted relationship was similar between infected and non-infected females.
Long-term changes in parasite species composition
The most parsimonious model explaining the probability of an individual either being infected or non-infected by P. azurophilum, P. leucocytica or P. floridense was best described by the additive effect of SVL, sex and year (ΔAICc = 9.02 between this model and the next best, AICc weight = 0.99; Supplementary Table S9). This model predicts a decrease in the probability of infection of all Plasmodium parasites after 2001. The probability of infection by P. azurophilum increases to the levels before 2001 in 2017 (Fig. 6, Supplementary Table S10). The most common species was consistently P. azurophilum. The remaining parasite species, P. floridense and P. leucocytica changed their relative dominance, but this apparently was due to a secular decline in P. floridense over the entire study period (Fig. 6, Supplementary Table S10). This model predicted no differences in the partial relationship between SVL and probability of infection by the different Plasmodium parasites (Supplementary Fig. S3). Whereas the model predicts no difference between sexes in the probability of getting infected by P. azurophilum or P. floridense, the probability of becoming infected by P. leucocytica was greater in females (Supplementary Fig. S4).
Partial predictions from the most parsimonious multinomial model showing a general decrease through time in the probability of infection of A. gundlachi by three Plasmodium parasites. The model predicts a secular decrease in the probability of infection by P. floridense. Dots represent the infection probability estimates for each Plasmodium species and the bars represent 95% CIs.
Discussion
We assembled a long-term dataset including 15 surveys conducted from 1990 to 2017 to identify the environmental drivers of Anolis–Plasmodium temporal dynamics and their consequences for host body condition and parasite species composition. Our results show that the probability of infection for the anole by Plasmodium parasites varied little early in the study and declined substantially in 2015–2016. This drop in infection probability coincided with a severe drought in Puerto Rico. However, the long-term approach allowed detection of subtler, but important patterns. Overall, temporal variation in probability of infection was associated with fluctuations in temperature and rainfall 4 months before each sampling period. The relationship with both environmental variables was nonlinear predicting maximum probability of infection at intermediate rainfall and a decelerating relationship with increasing maximum temperature. Similarly, this long-term analysis showed that host body condition maximized at the intermediate levels of maximum temperature and/or rainfall. Moreover, the optimal temperature for female's body condition was very similar to the optimal temperature for transmission (26 °C), suggesting that the environmental conditions that maximize host health may also catalyse parasite transmission. During the study period, P. azurophilum was consistently the most common parasite while P. floridense declined to almost local extinction. The study by Schall et al. (Reference Schall, Pearson and Perkins2000), which analysed the first decade of data, concluded that the dynamics were relatively stable over time because these patterns were not apparent during the early years of the study. Therefore, a longer-term approach (>10 years) was needed to reveal the response of the parasite dynamics to environmental changes and its consequences.
The delayed link between rainfall and temperature with parasite prevalence (120 days) most likely is a consequence of shifting vector abundance. During the severe 2015 drought, we noted an overall dry aspect to the forest, with little standing water in puddles, Heliconia flowers, or epiphytic bromeliads that could have restricted breeding locations for mosquitoes. Most strikingly, despite the general effect of variation in temperature and humidity over the study period, no dramatic effect on prevalence was seen after short environmental disruptions from dry periods or severe hurricanes in the 1990s. Hurricanes caused severe physical damage to the forest including massive defoliation (Reagan and Waide, Reference Reagan and Waide1996; Schall et al., Reference Schall, Pearson and Perkins2000). Thus, the drop in the probability of infection in the later samples could indicate that the 2015 drought was of sufficient duration and strength to shift the system to a new dynamic state. Substantial theory and empirical evidence support the conclusion that natural ecological systems can experience such alternative stable states (Petraitis, Reference Petraitis2013).
This new stable state may be the consequence of shifts in the dominance of polymorphic forms of each species as a result of co-evolution of the host and its parasites. Strain theory has a long history in human malaria research recognizing that these parasite species are highly polymorphic, often resulting in a collection of almost discrete forms within species (Gupta and Day, Reference Gupta and Day1994; McKenzie et al., Reference McKenzie, Smith, O'Meara and Riley2008). The severe drought of 2015 may have shifted dominance patterns favouring strains that can still spread under these generally unfavourable conditions at the expense of having low transmission rates. A non-mutually exclusive alternative hypothesis includes a decrease in favourable environmental conditions for vector transmission.
The identity of the vectors that transmit Plasmodium among Anolis lizards in Puerto Rico is not yet known. Two of our findings may suggest that there may be more than one dipteran species transmitting the parasite. First, if there were a unique vector species transmitting the three Plasmodium parasites we may find similar infection patterns among species. However, our results show that the probability of infection of female lizards is greater for P. leucocytica and P. floridense, but there is no difference between sexes for P. azurophilum (Supplementary Fig. S3). This difference is not likely to be the result of similar natural history traits between the parasites species. The latest Plasmodium phylogeny shows P. leucocytica and P. azurophilum as sister taxa and P. floridense more related to a different group (Galen et al., Reference Galen, Borner, Martinsen, Schaer, Austin, West and Perkins2018). Second, we found different patterns of temporal variation in the probability of infection among parasite species (Fig. 6). This may be the result of different responses to climate variability by multiple vector species. Our research group is currently testing mosquitoes for the presence of the three Plasmodium species infecting the lizards in El Verde. To date, in the field station we have only detected these parasites in Culex mosquitoes (Acevedo, M, unpublished), which are also the most common mosquito genus at the site (Yee, D, unpublished). Moreover, P. floridense is known to infect Culex mosquitoes in Florida (Klein et al., Reference Klein, Young and Telford1987). Yet, amplifying Plasmodium from a mosquito shows that the vector fed on an infected lizard, but does not necessarily imply vector competence. Therefore, more thorough work is still needed to identify the vector(s).
Why are males more likely to get infected than females? Higher rates of infection in males are common in many species (Roberts et al., Reference Roberts, Walker and Alexander2001). Ecological and physiological mechanisms have been proposed to explain this pattern. While the ecological mechanisms explaining this pattern in Puerto Rican anoles are not yet understood, we know that males can defend territories that are two times larger than females (Losos, Reference Losos2011). To defend their territories males move larger distances increasing their exposure to infectious vector bites by experiencing a larger diversity of vector habitat. Physiology may also play an important role. Males are characterized by having higher testosterone levels that have been shown to suppress immune responses against protozoan infections in many species (Wunderlich et al., Reference Wunderlich, Mossmann, Helwig and Schillinger1988; Zuk and McKean, Reference Zuk and McKean1996).
What caused the differences in parasite prevalence for the summer vs winter seasons? Rainfall is greater in the summer, which increases potential habitats for vector breeding, increasing the ratio of vectors to hosts. Also, summer is the mating season for A. gundlachi when lizards are more active and defending territories (Gorman and Licht, Reference Gorman and Licht1974; Reagan, Reference Reagan1992). If the vectors are daytime active, the lizards could be more exposed to vector bites. Alternatively, if the parasite reduces reproduction in the blood when vector activity ceases, the lower prevalence in winter may only be apparent. There is no evidence that lizards can clear an infection. Also, using a polymerase chain reaction-based method to detect very low-parasitemia infections showed cryptic infections were more common in the winter season (Vardo et al., Reference Vardo, Wargo and Schall2005).
Our study adds to the growing evidence that the relationship between malaria transmission and climatic measures is complex and often nonlinear (Garamszegi, Reference Garamszegi2011; Mordecai et al., Reference Mordecai, Paaijmans, Johnson, Balzer, Ben-Horin, de Moor, McNally, Pawar, Ryan, Smith and Lafferty2013). The observed patterns could be the result of complex interactions between variation in environmental conditions and its effects on various stages in the parasite life cycle. Variations in temperature and rainfall can affect host susceptibility, which is particularly important in ectothermic hosts such as anoles. Our results showing a non-linear effect of environmental fluctuation and anole body condition partially supports this hypothesis. Environmental fluctuations can also affect other variables predicted by the Ross-Macdonald theoretical model that influence transmission such as vector abundances, biting rates, the probability of infection from the host to vector and from the vector to host, vector death rates, incubation periods and recovery rates (Smith et al., Reference Smith, Battle, Hay, Barker, Scott and McKenzie2012). Adding to this complexity is the role of spatial scale. For example, in the Sceloporus occidentalis–P. mexicanum system in California, when studying the relationship between landscape features and infection prevalence, the type of ground cover (grass, rocks and leaf litter) influenced the probability of infection (Eisen and Wright, Reference Eisen and Wright2001). Note that changes in temperature and rainfall could alter microhabitat quality, including the production of leaf litter and ground cover. Therefore, suggestions that rising temperatures alone will increase the prevalence of malaria parasites ignore the more complex story.
Temperature and rainfall also mediated Anolis body condition. Measures of body condition are commonly used as a proxy of the energetic state of lizards and their fitness (Peig and Green, Reference Peig and Green2009; Cox and Calsbeek, Reference Cox and Calsbeek2015). The host A. gundlachi is a thermoconformer-shaded forest specialist whose metabolic performance depends closely on temperature (Huey and Webster, Reference Huey and Webster1976). We found that females had an optimal body condition at an environmental temperature of ~26 °C, which is consistent with the experimental voluntary maximum temperature tolerated by the species (Hertz et al., Reference Hertz, Huey and Stevenson1993). The non-linear relationship between A. gundlachi’s body condition and rainfall likely reflects the conflict between the levels of rainfall that promote abundance of food resources and periods of extremely heavy precipitation that result in physiological stress for the lizards. The optimal temperature for female body condition was similar to optimal temperature for transmission. This result suggests support for the theoretical prediction that environmental conditions mediate the interaction and potential co-evolution between the host and parasites (Elliot et al., Reference Elliot, Blanford and Thomas2002; Wolinska and King, Reference Wolinska and King2009). Because the parasite reproduces at the expense of the host resources, having ambient temperatures that promote optimal host body condition will likely result in optimal resources for the parasite to exploit for within-host reproduction (Seppälä et al., Reference Seppälä, Liljeroos, Karvonen and Jokela2008). Optimal temperatures may also mediate other mechanisms that promote transmission. For example, increasing temperatures may result in higher gametocyte production, but if production increases beyond a certain threshold it may ultimately decrease vector survival (assuming virulence to the vector). Also, temperature may mediate transmission-blocking immune responses or processes related to competition among co-infecting parasites (Mideo and Day, Reference Mideo and Day2008). Field experiments will be a natural follow-up to confirm these observational results.
We found no significant difference in body condition between the infected and non-infected individuals. This is consistent with previous studies that found little evidence of a relationship between male lizard health and parasite infection in the Caribbean (Schall and Pearson, Reference Schall and Pearson2000; Schall et al., Reference Schall, Staats and Gatten2002), but contrasts with studies of lizard malaria in California or Africa (Schall, Reference Schall2002). This pattern in lizards may be explained by three non-mutually exclusive explanations. First, virulence to male hosts can be expressed in other traits that may not be captured by a body-condition index such as changes in aggressive behaviour or stamina. Second, parasitemia of Caribbean Plasmodium species is lower than those of California or Africa (Schall pers obs.), which may result in a lower rate of blood cell turnover. Third, low virulence could be the result of lower transmission intensity as predicted by the theory of virulence evolution in the host–parasite systems that have co-evolved for many generations (Ewald, Reference Ewald1995). The pattern in avian malaria is similar in which some host species Plasmodium infection has negligible effects on body condition – mostly attributed to low parasitemia (e.g. Granthon and Williams, Reference Granthon and Williams2017), but there are significant negative effects in other species leading to dramatic increases in mortality (e.g. van Riper et al., Reference van Riper, van Riper, Goff and Laird1986).
The El Verde lizard malaria system has an added interest because we could follow through time the relative proportions of three species of Plasmodium that infect a single lizard host. Despite the major drop in the probability of infection during the later sample periods, and the major disruption of the forest by hurricanes in the early years, P. azurophilum remained at a constant higher proportion of infections. P. floridense gradually decreased over the years almost approaching local extinction. What could account for this secular change? P. floridense and P. leucocytica infect different cell classes (erythrocytes for P. floridense and several classes of white blood cells for P. leucocytica), and therefore are not likely to be competitors for within-host habitats. Yet, P. floridense would be competing with P. azurophilum – the most dominant species – for erythrocytes. The competitive interactions of malaria parasite species that exploit the same host are poorly studied (but see Schall and Bromwich, Reference Schall and Bromwich1994). Here, the gradual drop over in the prevalence of P. floridense over both hurricane and drought conditions is particularly perplexing. We suspect the explanation falls to the biology of the vectors, either different species or how parasites may compete within a single vector. Similar parasite species turnovers have been observed in long-term studies of bird malaria. For instance, Plasmodium and Haemoproteus lineages in Caribbean birds show significant temporal variation at the scale of decades, which suggests frequent parasite's local colonization and extinction events (Fallon et al., Reference Fallon, Ricklefs, Latta and Bermingham2004).
The most significant finding of this study is the value, even the necessity, of a long-term approach to understand the ecology of host–parasite systems. In our study, the cycle of parasites between hosts likely follows an annual pattern because the lizards live ~1 year. Studying the long-term dynamics of the A. gundlachi–Plasmodium system allowed us to gain a better picture of the patterns and possible mechanisms behind prevalence changes, shifts in the blood parasite assemblage and the effects on the body condition of the hosts. The influence of environmental variables became apparent only after taking this long-term approach.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182018001750.
Acknowledgements
We are grateful to the Editor and two anonymous reviewers who provided invaluable comments and suggestions that greatly improved the presentation of our work. We thank Stephanie Aguila, Larry Alejandro, David Clark, Génesis Colón, Laura Davo, Ashleigh Farmer, Paola Iglesias, Adriana Ortiz, Víctor Ramos, Judith Reyes, Gabriela Silva and Norberto Torres for their help in the field. We thank the NSF-LTER group at El Verde for providing lodging and logistics.
Financial support
MAA and LO were supported by the University of Puerto Rico Fondo Institucional para la Investigación (FIPI) from the University of Puerto Rico and a grant from the Puerto Rico Science and Technology Trust. JJS was supported by the USA National Science Foundation.
Conflict of interest
None.
Ethical standards
This research was conducted under the permits of Institutional Animal Care and Use Committee protocols from the University of Puerto Rico.
