Distribution of clones among hosts for the lizard malaria parasite Plasmodium mexicanum

Distributions

Two distributions are commonly used to model count data like the number of infectious bites per host: the Poisson distribution and the NB distribution (Hilbe, 2014). These distributions model the number of times an event occurs in a given interval of time or space. When applied to the distribution of malaria clones, the ‘events’ are the clones themselves (or, equivalently, bites by an infected insect vector) and the units of space are the hosts. In practice, the Poisson and NB distributions used are usually conditional (they are zero-truncated, meaning that they do not include not infected hosts, Schneider, 2021), in part to account for the fact that many data sets do not include robust sampling of not infected hosts (Hill & Babiker, 1995).

The simpler of the two is the Poisson distribution. The Poisson distribution is described by a single parameter, µ, which is both the mean and variance of the distribution. The Poisson distribution has a number of assumptions, including the independence of clones (i.e., having a certain clone does not affect the chances of having another) and the equivalence of the mean and variance (i.e., as the average number of clones per infection increases, so does the variance in clones per infection).

In reality, most parasite distributions are overdispersed relative to the Poisson distribution, meaning that the variance-to-mean ratio is greater than 1. For microparasites like malaria, overdispersion often means that parasite clones show significant clumping or aggregation, with a smaller number of hosts containing a larger number of clones than would be expected if clones were distributed at random. This clumping can occur for a variety of reasons, including if individual hosts vary in their susceptibility or exposure or if parasite propagules are distributed unequally in the environment (e.g., due to variation in vector abundance or vectors carrying multiple clones; Boulinier, Ives & Danchin, 1996). The negative binomial (NB) distribution is often a better fit in such cases (Anderson & May, 1992). The NB is described by two parameters: the mean, µ, and the dispersion parameter, α. When α is 0, the NB is equivalent to the Poisson. Higher values of α indicate greater aggregation (i.e., parasite clones are concentrated in a smaller proportion of the host population). Note that the statistical software R, which is used in this paper, uses θ (or ‘size’) in place of α.θ is inversely related to α, thus small values of θ indicate substantial clumping of clones and large values indicate similarity to a Poisson distribution (Hilbe, 2014). Hilbe (2014) provides useful coverage of both the Poisson and NB distributions.

Deviation from these basic distributions could occur for a number of reasons. For example, clones may not be independent of one another, such as if multiple malaria clones are ingested by a vector and transmitted simultaneously to a new host. Hill & Babiker (1995) dismissed this as unlikely, but experimental data suggest not only that multiple clones can transfer to the vector and develop there (e.g., Huber et al., 1998; Vardo-Zalik, 2009), but also that clones that were underrepresented in the vertebrate may actually tend to increase their relative representation in the vector (Nwakanma et al., 2008; Taylor, Walliker & Read, 1997; Vardo-Zalik, 2009). Hill & Babiker (1995) outline an approach to model the distribution of clones in such a case, but note that we do not have sufficient information about the transmission biology of Plasmodium to build a meaningful model at this time. Nonetheless, it seems that this situation would lead to an overabundance of multi-clone infections relative to the Poisson or negative binomial distribution.

Transmission may also be extremely localized such that some portion of an apparently continuous population has little to no chance of contracting the parasite, while another portion is at much higher risk than the overall prevalence of the parasite would otherwise suggest. Although some variation in exposure or susceptibility can generate the clumping and overdispersion characteristic of the negative binomial distribution (Boulinier, Ives & Danchin, 1996), major differences would cause deviations. Data suggest that such heterogeneity in exposure is common. Transmission of human malaria is known to show hotspots, where transmission in small areas is much higher than the surrounding region (Gaudart et al., 2006; reviewed in Bousema et al., 2012). Similar patterns are seen in wildlife parasites. For example, Eisen & Wright (2001) measured multiple landscape features at different spatial scales and found that ground cover within 10m of a host’s point of capture was the best predictor of infection status with the lizard malaria parasite Plasmodium mexicanum.

An extreme case of heterogeneity in parasite exposure in which a portion of the population is not exposed at all and a portion has ordinary variation in exposure can be modeled by a zero-inflated distribution. This distribution assumes some proportion of the population has zero probability of becoming infected while the remainder of the population has clones distributed in accordance with the Poisson or NB distribution. Zero-inflated models are mixed models that include a parameter (the ‘zero-probability’) describing the proportion of observed 0s (uninfected hosts) that are not a result of the processes that generate the other distribution (Poisson or NB; Hilbe, 2014); in other words, some portion of uninfected hosts are not infected because they are not exposed to infectious vectors at all and other uninfected hosts just happen to have not been bitten (but had the same chance of being bitten as the infected hosts). The number of not exposed hosts can be estimated as the zero probability by fitting this zero-inflated distribution. Note that fitting a zero-inflated model should predict a very similar distribution of clones among infected hosts to fitting a zero-truncated (conditional) model, but the former may provide greater insight on whether or not infected hosts are subject to the same underlying forces (e.g., exposure to infected vectors) as infected hosts.

One final example of why standard distributions may not accurately predict clonal distributions relates to intraspecific dynamics within infected hosts; negative interactions among Plasmodium clones could lead to a paucity of multi-clone infections. For example, there is ample evidence that established infections can prevent the entry of additional clones, a phenomenon known as premunition (e.g., Vardo, Kaufhold & Schall, 2007; De Roode et al., 2005; reviewed for P. falciparum in Smith et al. 1999). Blocking the entrance of additional clones would reduce the number of multi-clone infections because a host might harbor only a single parasite clone even if it had received multiple infectious bites. This might alter the distribution of clones in various ways depending on the transmission intensity and strength of premunition, but it seems likely that it would lead to data that are underdispersed relative to the Poisson. At the extreme, with very high transmission and perfect premunition (and ignoring the possibility that vectors may transmit multiple clones at once), most hosts would be infected by a single clone; the mean number of clones per infection would be low (near 1 clone per host), but the variance would be even lower.

Research question

Obtaining reliable data for testing the distribution of clones among hosts can be challenging. One reason is that some parasites (e.g., Plasmodium falciparum) are known to sequester in capillaries, making it necessary to take multiple samples from a patient to ensure that sequestered clones are not missed (Tusting et al., 2014). Additionally, control efforts like anti-malarial prophylactics (Tusting et al., 2014) and vaccine candidates (reviewed in Arnot, 2002) can alter MOI. Third, reliable counts of uninfected hosts are not always available due to study design (e.g., uninfected hosts may not be sampled) or risk of false negatives (Hill & Babiker, 1995).

I therefore pursue a study of the distribution of clones among hosts using the lizard malaria parasite Plasmodium mexicanum in its host the western fence lizard Sceloporus occidentalis. This natural host-parasite system has been studied for decades (Schall & St. Denis, 2013) and is free from many of the challenges outlined above: the parasite is not known to sequester (Ford & Schall, 2011), the hosts are not exposed to chemotherapy or vaccination, and false-negatives are known to be quite rare (Perkins, Osgood & Schall, 1998). Additionally, despite differences in biology between P. mexicanum and human parasites (e.g., P. mexicanum is transmitted by sand flies rather than mosquitoes), many of the biological complexities potentially impacting the distribution of clones among hosts have shown similar patterns between P. mexicanum and other malaria parasites (see papers referenced under ‘Distributions’). My aim was to determine whether the distribution of clones among hosts differs significantly from the distributions outlined in Table 1.

Study sites and sampling

The malaria parasite Plasmodium mexicanum infects Sceloporus occidentalis lizards in western North America and has been studied since 1978 at the University of California Hopland Research and Extension Center (HREC) near the town of Hopland in Mendocino County California (Schall & St. Denis, 2013; Vardo & Schall, 2007). Lizards are captured using a fishing pole with attached noose, and a small sample of blood is taken before the lizard is released at its point of capture. Blood samples are used to make thin blood smears for processing with Giemsa stain and a small amount of blood is stored frozen on filter paper (Whatman^® qualitative circles or similar) for genetic analysis.

Sampling at HREC is conducted at a number of well-established sites. The parasite shows genetic differentiation among these sites, some of which are only a few hundred meters apart (Fricke, Vardo-Zalik & Schall, 2010). Because combining data from multiple sites or years might cause deviations from standard statistical distributions, I tested the fit of the selected distributions using data taken from the four individual sites sampled in 2009 or 2010 that had the most infections: Water Tank (WT, 2009), Parson’s Creek (PC, 2009), Middle Lower Horse (MLH, 2010) and Gorge (GOR, 2010). The largest of these sites, MLH, is less than 4.5 ha; Eisen & Wright (2001) describe this site in detail, though the area sampled in 2009–2010 is a bit smaller than their study area. The other sites are each closer to 1ha (Fig. 1).

Plasmodium infections in this system are known to be chronic; as a result, variation in the duration of exposure may also cause deviations from standard distributions (i.e., older lizards might harbor more parasite clones because they have been exposed to infectious bites over a longer period). To limit the impact of variation in lizard age on clonal distributions, I retained only lizards in a limited size range (40–64 mm) in the final analysis. For lizards collected in early summer (the vast majority of lizards in this data set), this range of sizes corresponds to lizards that hatched in late summer of the previous year with a high level of confidence (Davis, 1967).

Genetic analysis

All infections included in this study were genotyped previously; these data were published by Neal & Schall (2014). Genotyping was performed for 4 microsatellite loci- Pmx306, Pmx732, Pmx747 and Pmx839- using fluorescently labeled PCR primers and conditions described by Schall & Vardo (2007). The PCR product was diluted and underwent fragment analysis at Cornell University Biotechnology Resource Center.

The parasite is haploid in vertebrate blood, so each clone will have a single allele at each locus and produce a single peak on an electropherogram. Clones that are at a low relative frequency in the blood may be hard to detect in infections due to a weak signal on the electropherograms. No peak-size cut-off was applied to correct for this, but a previous study demonstrated that even alleles with 10-fold differences in concentration could usually be detected (Vardo-Zalik, Ford & Schall, 2009). Additionally, alleles were classified by relative electropherogram peak height (highest, >1/3 highest, <1/3 highest [small]); these designations were not used in my analysis, but are available in the data set at osf.io/6yvh5.

I estimated the number of clones per host as the maximum number of alleles seen at any one locus (a common measure also used for human malaria (Zhong et al., 2018)). Loci with fewer alleles may result if clones share alleles at that locus. This measure provides a minimum estimate of the number of clones per infection because some clones may share alleles, and will be referred to hereafter as the observed minimum number of clones, or observed MC. Anderson et al. (2010) showed that it was possible to distinguish about 90% of Plasmodium falciparum clones using only four microsatellites, and the loci in that paper had lower average heterozygosity than the loci I used. This suggests that observed MC, while a minimum estimate of MOI, may not greatly underestimate the number of distinct clones in most cases.

Fitting distributions

All analyses were performed in R (R Development Core Team, Vienna, Austria). The following analyses were performed on observed MC data for each of the four single sites individually (GOR [2010], MLH [2010], PC [2009] and WT [2009]). I fit the Poisson distribution using the glm() function and the NB using the glm.nb() function. The zero-inflated models were fit using the zeroinfl() from the R package “pscl”. I assessed dispersion using the Pearson dispersion statistic as suggested by Hilbe (2014) and compared each zero-inflated model to the corresponding standard model with the Vuong test using the function vuong(). I performed a chi-squared test comparing the observed vs. expected values of each distribution after grouping multi-clone infections to reduce the number of categories with fewer than five observed infections following (generally) the suggestion of Bliss & Fisher (1953; two sites unfortunately had fewer than five single-clone infections: PC had four and WT had three). Due to the relatively small number of infections at each site considered individually, I used Fisher’s combined probability test to obtain an overall assessment of whether the distribution of observed MC at the four sites collectively differed significantly from each distribution.

As mentioned above, observed MC is known to provide a minimum (and therefore biased) estimate of the true MOI. True MOI is almost certainly higher because some distinct clones will share alleles and some hosts may have been infected by the same clone multiple times (Schneider, 2018). If the bias is slight, it may not greatly affect the fit of distributions, but if the bias is greater, it could cause MOI estimates to deviate from expected distributions even if true MOI does not. To account for this possibility, I used the likelihood approach of Schneider & Escalante (2014) to obtain an unbiased (Schneider, 2018) estimate of average MOI at each site, simulated the observable number of alleles per locus using a Poisson distribution of clones per infection and allele frequencies observed at each site, and compared this expected distribution to the distribution of observed MC. True MOI and true allele frequencies were estimated from genetic data for each site using the MLE() R function provided by Schneider (2018; supplemental materials). A separate MOI estimate was obtained for each of the four microsatellite loci, and these values were averaged to obtain a single overall estimate of MOI for each site. I then simulated the distribution of clones among n_s hosts, where n_s is the number of lizards collected for the observed data from each site s. Simulations were performed by drawing n_s random numbers from a Poisson distribution with lambda equal to the estimated MOI. These random numbers represent the number of infectious bites received by each of the n_s simulated hosts. I then used the function mnom(), also from the R code in Schneider (2018; supplemental materials), to randomly draw x_i alleles from a multinomial distribution with allele frequencies estimated from the observed data (described above), where x_i is the random number of infectious bites for each host drawn from the Poisson distribution in the previous step. Sampling of the multinomial distribution was performed for each of the n_s hosts at each of the four microsatellite loci. After allele sampling, a count of the number of distinct alleles per locus was tallied (note that some alleles may be sampled more than once, but would only be counted once) and the maximum number of distinct alleles at any locus was recorded as simulated MC. Taken together, these steps were meant to simulate the expected underlying biology (Poisson distribution of infectious bites) and known observational limitations (overlap in alleles among clones, repeated infection with the same clones) of the observed data. The simulation was repeated 1,000 times to determine the mean, minimum and maximum number of infections observed with each possible number of clones to determine if observed data fell in this range (estimation of a 95% confidence interval did not seem necessary; see Results).