INTRODUCTION
Plants and animals are not evenly distributed across their geographical ranges: the abundance of a species varies in space in response to variations in the suitability of the habitat. Abundance is usually highest in parts of the range where conditions are near optimum, resulting in an uneven distribution across the species' geographical range (Gaston, 2003). The habitat of a parasite differs from that of most free-living organisms; it is not spatially continuous but consists of similar patches represented by its hosts, whereas the environment between these patches is absolutely unsuitable for parasites (Poulin, 1998a). The distribution of a parasite population across a host population is characterized by its aggregation, or overdispersion. In other words, most parasite individuals occur in a few host individuals, while most host individuals have only a few, if any, parasites (Anderson and May, 1978; Poulin, 1993; Shaw and Dobson, 1995; Wilson et al. 2001). This particular distribution of parasite individuals among hosts is caused by a variety of factors (Poulin, 1998a) and has important consequences for different aspects of the evolutionary ecology of parasites (e.g. Morand et al. 1993). Overdispersion of parasites is an almost universal phenomenon (Anderson and May, 1978; May and Anderson, 1978; Anderson and Gordon, 1982; Shaw and Dobson, 1995; Shaw, Grenfell and Dobson, 1998). This ubiquity suggests that similar processes may be involved in the generation of the same pattern in different host-parasite systems. In particular, overdispersion is commonly thought to arise as a result of heterogeneities in host populations and/or infection pressure (Anderson and May, 1979; Shaw and Dobson, 1995; Wilson et al. 2001). However, while all parasites appear to be aggregated, strong differences exist in the degree of aggregation and the causes of this are not well understood.
In spite of sharp differences in the pattern of spatial distribution between free-living organisms and parasites, the relationship between abundance and distribution is surprisingly similar in both groups. Indeed, a positive relationship between local abundance and occupancy is one of the most pervasive macroecological patterns (Gaston, 2003). The positive correlation between local abundance and occupancy has been documented at a variety of scales, across various habitats, in different biogeographical regions and for various taxa (see Gaston, 2003 and references therein). In particular, in the application of this relationship to host-parasite systems, a positive correlation between the mean number of parasite individuals per host (mean abundance) and the percentage of infected hosts (host occupancy by parasites, i.e. prevalence) has been supported in many studies (Shaw and Dobson, 1995; Morand and Guégan, 2000; Krasnov, Khokhlova and Shenbrot, 2002a; Simkova et al. 2002). This pattern can be explained by a number of mechanisms that are not all mutually exclusive (see review in Gaston, Blackburn and Lawton, 1997; Gaston, 2003).
Morand and Guégan (2000) tested 4 of 8 existing hypotheses (see reviews in Gaston et al. 1997; Gaston, 2003) explaining the positive relationship between abundance and prevalence of nematodes, and found that this pattern could be explained by demographic and stochastic mechanisms revealed by simple epidemiological models without invoking complex mechanisms such as a degree of ecological specialization (host specificity) (Brown, 1984). In contrast, studies of abundance and prevalence of monogeneans of freshwater fish of Slovakia (Simkova et al. 2002) and fleas parasitic on small mammals of Holarctic (Krasnov et al. 2004) demonstrated that the degree of specialization of a parasite was negatively correlated with its mean abundance and, thus, can explain positive abundance-prevalence relationship. These studies, as well as a comparative analysis of distribution patterns in macroparasite infections of wildlife (Shaw and Dobson, 1995), were carried out using interspecific comparisons. However, a trend of positive abundance-prevalence relationships that arises from interspecific comparisons can represent a net result of different intraspecific patterns and, thus, mask true relationships in particular host-parasite systems.
Another pattern of abundance and distribution that is astonishingly similar in both free-living and parasitic organisms is the power relationship between mean abundance and its variance. This empirical relationship known as Taylor's power law (Taylor, 1961) is supported by numerous data on both free-living and parasitic species (Taylor and Taylor, 1977; Taylor and Woiwod, 1980; Anderson and Gordon, 1982; Perry and Taylor, 1986; Shaw and Dobson, 1995; Morand and Guégan, 2000; Simkova et al. 2002). Exponent (parameter b or slope of Taylor's relationship) of this power function usually varies among species as 1<b<2, but the causes of the variation of this parameter among species (e.g. Kilpatrick and Ives, 2003) as well as within species (at spatial or temporal scale) are far from being understood. For parasites, it has been shown to be an inverse indicator of parasite-induced host mortality (Anderson and Gordon, 1982) as an increase in b suggests that at least some of the hosts are infected with heavy burdens of parasites. Kilpatrick and Ives (2003), using simulation models, demonstrated that the between-species variation in this parameter can be caused by negative interactions among species in a community. However, variation of the slope of Taylor's law has been rarely tested by real data.
Here, we report on a study of abundance and distribution of fleas parasitic on rodents in the desert region of Israel. Fleas (Siphonaptera) are characteristic ectoparasites of mammals and are most diverse and abundant on small and medium-sized burrowing species. Fleas usually alternate between periods when they occur on the host body and periods when they occur in its burrow or nest. In most cases, pre-imaginal development is entirely off-host. The larvae are usually not parasitic and feed on organic debris in the burrow and/or nest of the host. The cycle of activity and parasitism of imago fleas varies from being continuously active all year to short-term strictly summer or winter periods of activity. In the study region, flea species range from being highly host specific to being host opportunistic (Krasnov et al. 2004) and are active either all year round or only during the winter (Krasnov et al. 1997, 2002b).
The goals of this study were to investigate possible causes of variability in aggregation and abundance-prevalence relationships of fleas among different host-flea systems as well as among different host and flea species. Specifically, we aimed to (a) confirm aggregated distribution of fleas within populations of host species as well as across whole host communities and to understand how the degree of aggregation varies among different flea-host systems; (b) confirm the existence of a positive relationship between local parasite abundance (mean number of fleas per individual host) and their prevalence and to test the hypothesis of ecological specialization (Brown, 1984; Gaston et al. 1997), predicting that mean abundance and prevalence will be higher in host-opportunistic fleas than in host-specific fleas and in fleas that are active all year round than in fleas that are active as adults only in winter; (c) understand if and how the degree of flea specialization and flea community richness affects the pattern of the relationship between mean abundance and its variance testing the hypothesis of Kilpatrick and Ives (2003) that the slope of Taylor's power law is negatively correlated with flea community size; and (d) test whether a simple epidemiological model based on mean flea abundance and degree of aggregation can predict the observed pattern of prevalence. We also tested 2 different methods of the aggregation parameter estimation and fitted models based on these methods to the observed patterns of prevalence.
MATERIALS AND METHODS
Study area
Rodent trapping and flea collection were carried out in 3 regions of southern Israel (Ramon erosion cirque, Hatzeva Nature Reserve and Ashalim sands). The Ramon erosion cirque is situated at the southern boundary of the Negev Highlands (30 °35′N, 34 °45′E). The altitude of the north rim of the cirque ranges between 900 and 1200 m above sea level, while the south rim is about 510 m above sea level. The climate is characterized by hot, dry summers (mean daily air temperature of July is 26·2 °C) and relatively cool winters (mean daily temperature of January is 12·5 °C). There is a sharp decrease in annual rainfall from 100 mm on the north rim to 56 mm in the bottom of the cirque. The dominant vegetation consists of Retama raetam, Zygophyllum dumosum, Artemisia monosperma, Anabasis articulata, Atriplex halimus, and Noaea mucronata.
Hatzeva Nature Reserve (30 °48′N, 35 °15′E) is situated in the central part of the Arava valley (60 m above sea level). The climate in this area is similar to that of the Ramon cirque, although mean daily temperatures are slightly higher (29·8 °C in July and 14 °C in January) and annual rainfall is half (25·8 mm). Dominant vegetation is represented by Haloxylon persicum and Calligonum comosum.
Ashalim sand dunes are situated in the western Negev desert (30 °59′N, 34 °41′E). The landscape is represented by linear stabilized dunes, covered by A. monosperma, Echiochilon fruticosum and R. raetam. Average annual rainfall is 108 mm, and the mean daily air temperature is 25 °C in July and 11·1 °C in January.
Rodent trapping and flea collection
In 1992–2004 we carried out 226 short-term (1–5 days) rodent trapping surveys on 1 ha plots (179 surveys in the Ramon cirque, 18 surveys in Hazeva Nature Reserve and 28 surveys in Ashalim sands). Sampling was evenly balanced across seasons. Each plot was sampled with 25–50 Sherman live-traps placed in a 5×5 grid. We collected fleas from each individual rodent only when captured the first time. The animal's fur was combed thoroughly, using a toothbrush, over a white plastic pan and fleas were carefully collected. The fur was brushed several times until no fleas were recovered. Fleas were stored in 70% alcohol. After this, each rodent was marked by either toe-clipping or by implanting a micro-transponder subcutaneously with a Trovan® ID-100, and released. For subsequent individual identification, we used Allflex® ISO compatible RFID portable reader for passive electronic identification of transponders. Data only on common host species and common flea species were included in the analyses (in total, 3127 individual adult hosts and 16551 individual fleas; see Results section for details).
Analysis of general patterns of abundance and prevalence
We calculated mean abundance, variance of abundance and prevalence for each trapping survey for each flea species and each host species; for each host species across all flea species (except for Acomys cahirinus, Gerbillus andersoni allenbyi and Gerbillus pyramidum that were infested by 1 flea species only; see below); and for each flea species across all host species on which this flea was found (except for Parapulex chephrenis and Nosopsyllus pumulionis which infested 1 host species only).
We tested for negative binomial distribution of each flea species on each host species, as well as within host species across all fleas, and within flea species across all hosts using the maximum likelihood method (see Wilson et al. 2001 for details). The distribution of the prevalences of a community of parasites among host populations has been reported to be either bimodal (Morand and Guégan, 2000) or negative binomial (Simkova et al. 2002). We looked at the distributions of the prevalences of fleas within host species across all flea species as well as within flea species across all host species. There was no hint of the bimodality in the data. Therefore, frequencies of mean flea abundances (calculated for each survey) were tested against log-normal distributions for each flea-host association, each host across all fleas and each flea across all hosts using the Kolmogorov-Smirnov test for normality on log-transformed mean abundances. Analyses within hosts across fleas and within fleas across hosts were not carried out for hosts A. cahirinus, G. a. allenbyi, and G. pyramidum, and for fleas P. chephrenis and N. pumilionis (see above). Analyses for seasonal fleas (see Results section) were carried out only for seasons when these species occurred.
Variance/mean analysis
Mean abundance (M) and variance of abundance [V(M)] of an organism's distribution are related by a power law (Taylor, 1961) V(M)=aMb. We regressed log-transformed variance of flea abundance against log-transformed mean of flea abundance (both calculated within a trapping survey) for each flea-host association, each host across all fleas and each flea across all hosts and compared among slopes of the resulting relationships. We used t-tests to test whether the slope differs significantly from 1. To test for the relationships between the slope of Taylor's relationship and species richness of flea assemblages, we calculated the number of co-occuring flea species for each flea species across all hosts and in the entire sampling period (assemblage richness) and regressed values of b against log-transformed values of the assemblage richness using both conventional regression and method of independent contrasts (Felsenstein, 1985) that controls for the confounding effect of phylogeny. A phylogenetic tree for fleas was constructed using morphological taxonomy (see Krasnov et al. 2004 for details). The same flea species occurring on different hosts were considered as sister branches. To compute independent contrasts, we used the PDAP:PDTREE program (Midford, Garland and Maddison, 2004) implemented in Mesquite Modular System for Evolutionary Analysis (Maddison and Maddison, 2004). The procedure of the method of independent contrasts followed Garland, Harvey and Ives (1992). Regressions were forced through the origin as is standard practice with contrasts, because of the arbitrariness associated with the sign of independent contrasts (see Garland et al. 1992 for details). To test for the effect of the degree of ecological specialization of fleas on the slope of Taylor's relationship, we compared this slope between host specific or seasonal (‘specialist’) and host opportunistic (‘generalist’) or year-round-active fleas (see below), respectively, using Mann-Whitney U-test.
Analysis of abundance/prevalence relationships
To estimate the proportion of variance in flea prevalence explained by their mean abundances, we fitted the observed relationship between prevalence and log mean abundance to the logistic curve. This is because a logistic curve, in general, fits well to an empirical relationship between the fraction of sites where a given species occurs (prevalence) and its mean abundance (Gaston, 1999). This was carried out for each flea on each host, each flea across all hosts and all fleas within each host. We applied least squares estimation procedures via the Levenberg-Marquardt algorithm (Moré, 1977).
Test of ecological specialization hypothesis
We classified fleas as specialists and generalists with respect to the level of their host specificity and as seasonal and year-round-active according to their seasonality. A flea species was considered as a host specific (specialist) if it parasitized either only 1 host species or 2 consubgeneric host species and as a host opportunistic (generalist) if it parasitized 2 or more host species (belonging to the different subgenera, in the case of 2 hosts) (see Poulin and Mouillot, 2003). We compared prevalence, abundance, variance of abundance, slopes of Taylor's power relationship, and the indicator of aggregation k calculated by 2 methods (see below) between host specific and host opportunistic fleas and between fleas with year-round and winter periods of activity using Mann-Whitney U-tests.
Test of the epidemiological model
Epidemiological models (Anderson and May, 1985) predict that the probability distribution of parasite numbers per host individual, being negative binomial, determines the relationship between the prevalence of infection P(t) (proportion of infected hosts) at any given time and the mean number of parasites per host individual M(t) at time t as:
where k is the parameter of the negative binomial distribution inversely indicating degree of aggregation. When k is large (>20), the distribution converges on the Poisson.
There are several methods for estimation of k (Southwood, 1966; Elliott, 1977; Wilson et al. 2001). For example, k can be estimated using parameters a and b of Taylor's power law (Taylor, Woiwod and Perry, 1979) as
Another method to estimate k is to use the moment estimate of Elliot (1977), corrected for sample size:
where M is mean abundance, V(M) is variance of abundance and n is host sample size. We calculated k using both methods. Then we calculated the expected prevalence (Pexp) for each flea-host association across surveys, each host across all flea species and each flea across all hosts based on the two estimates of k (Pexp1 and Pexp2, respectively), and compared the estimated prevalence with the observed prevalence using linear regression. We used t-tests to test whether the slopes of the resulted regression differed significantly from 1.
We avoided an inflated Type I error by performing Bonferroni adjustments of alpha across all analyses. Significance is only recorded at the adjusted level.
RESULTS
Hosts and fleas
In total, we recorded 14 rodent species, namely Psammomys obesus, Meriones crassus, Gerbillus dasyurus, Gerbillus henleyi, Gerbillus gerbillus, Gerbillus nanus, Gerbillus andersoni allenbyi, Gerbillus pyramidum, Sekeetamys calurus, Acomys cahirinus, Acomys russatus, Mus musculus, Jaculus jaculus, and Eliomys melanurus. Species that either were uncommon (E. melanurus) or required special methods of sampling (P. obesus, J. jaculus) were excluded from the analysis. The same was true for rodents with extremely low flea infestation (G. henleyi, A. russatus, and M. musculus).
Fourteen flea species were collected including Xenopsylla conformis mycerini, Xenopsylla dipodilli, Xenopsylla ramesis, Xenopsylla cheopis, Parapulex chephrenis, Synosternus cleopatrae pyramidis, Coptopsylla africana, Nosopsyllus iranus theodori, Nosopsyllus pumilionis, Myoxopsylla laverani traubi, Stenoponia tripectinata medialis, Rhadinopsylla masculana, Ophthalmopsylla volgensis and Leptopsylla algira costai. Details of the relationships between particular rodent species and particular flea species have been published previously (Krasnov et al. 1997, 1999). Of these, the following common species were included in the analyses: X. c. mycerini, X. dipodilli, X. ramesis, P. chephrenis, S. c. pyramidis, N. i. theodori, N. pumilionis, R. masculana and S. t. medialis. In total, there were 17 flea species/host species associations (Table 1). Of these flea species, P. chephrenis, S. c. pyramidis, N. pumilionis, and R. masculana were considered host specific, whereas X. c. mycerini, X. dipodilli, X. ramesis, N. i. theodori and S. t. medialis were considered host opportunistic. Imagoes of both Nosopsyllus species, R. masculana and S. t. medialis were active in 3–4 months of the cool season (winter fleas), whereas all other fleas were active as adults all year round (Krasnov et al. 1997, 2002a).
General patterns of flea abundance and prevalence
The negative binomial models provided a statistically satisfactory fit to the observed frequency distribution of each flea species on each host species (Kolmogorov-Smirnov tests, d=0·0005–0·08, P>0·4 for all) as well as those of all fleas within host species (Kolmogorov-Smirnov tests, d=0·007–0·08, P>0·1 for all; see Fig. 1A for the illustrative example with G. dasyurus) and of each flea across all hosts (Kolmogorov-Smirnov tests, d=0·0004–0·01, P>0·7 for all).
Fig. 1. (A) Frequency distribution of number of fleas per individual Gerbillus dasyurus (n=1802 and k=0·48). (B) Frequency distribution of log-transformed mean abundances of fleas among G. dasyurus across trapping surveys.
Frequency distributions of mean flea abundances fitted well to a log-normal pattern among flea-host associations (Kolmogorov-Smirnov tests, d=0·09–0·29, P>0·2 for all; see Fig. 1B for the illustrative example with G. dasyurus). The same was true for distributions of flea abundances both within host species across all fleas (Kolmogorov-Smirnov tests, d=0·09–0·22, P>0·2 for all) and within flea species across all hosts (Kolmogorov-Smirnov tests, d=0·04–0·17, P>0·2 for all).
Frequency distribution of prevalences of fleas among samples was found to be normal in most cases (Kolmogorov-Smirnov tests, d=0·08–0·23, P>0·1 for all; see Fig. 2 for the illustrative example of N. i. theodori on G. dasyurus). However, frequency distribution of X. ramesis on M. crassus and of X. dipodilli across all hosts fitted to normal distribution after log-transformation only (Kolmogorov-Smirnov tests, d=0·21 and d=0·09, respectively; P>0·2 for both). In addition, S. c. pyramidis attained more that 90% prevalence on G. a. allenbyi in 22 of 25 surveys and 100% prevalence on G. pyramidum in 11 of 13 surveys. Consequently, no characteristic frequency distribution was found for this flea on either host.
Fig. 2. Frequency distribution of prevalences of Nosopsyllus iranus theodori on Gerbillis dasyurus across surveys.
Relationships between variance and mean abundance of flea infestation
The slope of the relationship between log-transformed variance and mean abundance of fleas was significantly greater than unity in all flea-host associations as well as within host species across all flea species and within flea species across all host species (Table 2; see Fig. 3 for illustrative example of G. dasyurus and X. dipodilli). In 1 case (S. c. pyramidis on G. pyramidum), a slope of log variance against log mean abundance exceeded, albeit non-significantly, a value of 2. Slopes of other regressions of log-transformed variances against log-transformed mean abundances varied from a low of 1·19 in N. pumilionis on G. nanus to a high of 1·93 in S. c. pyramidis on G. a. allenbyi.
Fig. 3. Relationship between log variance of abundance and log mean abundance of Xenopsylla dipodilli on Gerbillus dasyurus.
Table 2. Summary of regression analyses of relationships between log variance and log mean abundance of fleas on their rodent hosts (P<0·001 for all cases) (Slopes are significantly higher than unity (t-tests, P<0·01 for all cases). Abbreviations of host names are Ac (A. cahirinus), Ga (G. a. allenbyi), Gd (G. dasyurus), Gn (G. nanus), Gp (G. pyramidum), Mc (M. crassus) and Sc (S. calurus). Abbreviations of flea names are Ni (N. i. theodori), Np (N. pumilionis), Rm (R. masculana), Sc (S. c. pyramidis), St (S. t. medialis), Xd (X. dipodilli), Xc (X. c. mycerini) and Xr (X. ramesis).)
The slope of Taylor's relationship decreased with the increase of species richness of flea assemblages (number of additional co-occurring flea species across all hosts, ASR) as b=1·71–0·49log(ASR) (r2=0·39, F1,15=9·4, P<0·007). Comparative analysis of these two variables using independent contrasts produced a non-significant relationship (r=−0·28, P=0·26). However, when an outlier, N. pumilionis, was omitted from the analysis, a negative correlation between contrasts in b and contrasts in species richness of flea assemblages became significant (r=−0·53, P<0·03). The slope of Taylor's power relationship between mean abundance and its variance was significantly higher in host specific than host opportunistic fleas as well as in seasonal than in year-round-active species (Mann-Whitney test, Z=3·79 and Z=9·47, respectively; P<0·0001 for both).
Relationships between flea abundance and prevalence
A logistic function described the relationship between prevalence of fleas and their mean abundance except for (a) S. c. pyramidis on G. pyramidum and on both G. pyramidum and G. a. allenbyi and (b) for N. i. theodori on S. calurus (Table 3; see Fig. 4 for illustrative example of fleas on G. dasyurus). As in previous analyses, causes of these exceptions can be the absence of variance in prevalence (G. pyramidum and G. a. allenbyi) and small sample size (S. calurus). In general, prevalences were low at low flea abundances (as in Shaw and Dobson (1995)) and rose rapidly to high asymptotes (although varied markedly) with an increase in abundances (Fig. 4).
Fig. 4. Logistic function fit (black line) for prevalence of fleas parasitizing Gerbillus dasyurus plotted against log mean flea abundance (see details of fit in Table 3).
Table 3. Summary of fitting logistic function to the relationships between prevalence and log mean abundance of fleas (All regressions, except those marked by *, are significant, P<0·001. See Table 2 for the abbreviations of flea and host names.)
Ecological specialization and abundance/prevalence
Prevalence, mean abundance and variance of abundance were significantly higher in host-specific than host-opportunistic flea species (Mann-Whitney tests, Z=5·80, Z=5·64 and Z=5·06, respectively; P<0·0001 for all). Values of k calculated using moment estimation corrected for host sample size were significantly lower in host-specific than host-opportunistic fleas (Mann-Whitney tests, Z=9·47 and Z=3·95, respectively; P<0·0001 for both). In contrast, k values calculated using a and b coefficients of Taylor's power law did not differ between these two flea categories (Mann-Whitney test, Z=0·96, P>0·3).
If, however, ecological specialization was considered based on the seasonal pattern of activity, the comparison of abundance and distribution parameters provided different results. Prevalence, mean abundance and its variance were significantly higher in year-round-active than seasonal species (Mann-Whitney tests, Z=−2·56, Z=−2·23 and Z=−2·16, respectively; P<0·01 for all). Values of k calculated by either method did not differ between seasonal and year-round-active fleas (Mann-Whitney tests, Z=−0·93 and Z=−1·77, respectively; P<0·1 for both).
Flea distribution as predicted from a simple epidemiological model
In general, regressions of flea prevalences expected from epidemiological models with k values calculated from Taylor's power law with observed flea prevalences produced slopes significantly less than 1 (Table 4). This indicated that the models, in general, under-estimated prevalences and, in most cases, expected prevalences, Pexp1, were lower than observed prevalences, Pobs. However, if moment estimates of k corrected for host number were used instead, expected prevalences, Pexp2, did not differ significantly from observed prevalences (Table 4), except for P. chephrenis on A. cahirinus, S. c. pyramidis on G. a. allenbyi, X. dipodilli on S. calurus and pooled fleas on M. crassus (in these cases, expected prevalences were significantly lower than observed prevalences). Furthermore, the model with k calculated using Taylor's power law tended to over-estimate relatively low prevalences and under-estimate relatively high prevalences (see Fig. 5 for illustrative example of N. i. theodori). No significant relationship between expected and predicted prevalences was found for S. c. pyramidis on G. pyramidum and for N. i. theodori on S. calurus. No relationship in the former case can be due to the absence of variance in observed prevalences (100% in almost all surveys), whereas no relationship in the latter case can be explained by extremely small sample size (N. i. theodori was recorded on S. calurus in 4 surveys only). In addition, r2 values for the relationships between Pobs and Pexp1 were, in general, lower than those for Pobs and Pexp2.
Fig. 5. Relationship between observed and expected (from the epidemiological models, with different k estimation) prevalence of Nosopsyllus iranus theodori across all rodent hosts. Pexp1 – open circles, solid line; Pexp2 – closed circles, dashed line. See text for explanations.
Table 4. Summary of regressions of expected (Pexp1 and Pexp2) from the epidemiological model against observed prevalences of fleas infesting rodent hosts (All, except those marked by *, are significant, P<0·001). k values for Pexp1 and Pexp2 were calculated using either Taylor's power law or corrected for host number moment estimate, respectively. See Table 2 for abbreviations of host and flea names. * – see text for explanations for S. c. pyramidis on G. pyramidum and N. i. theodori on S. calurus. ** – slope does not significantly differ from 1 (t-tests, P>0·05). See Table 2 for the abbreviations of flea and host names.)
If, however, mean flea abundance, its variance and flea prevalence were calculated for pooled data for each host/flea species association across all surveys, models with k values estimated using both methods predicted flea prevalence equally well (slope of regression of expected against observed prevalence was 0·93±0·11 for Pexp1 and 0·96±0·04 for Pexp2, P<0·001 for both). Nevertheless, r2 value was higher for the model with k estimation corrected for host sample size (r2=0·81 versus r2=0·96).
DISCUSSION
General patterns of abundance and distribution
The results of this study confirm the aggregated distribution of fleas across their host individuals (Anderson and May, 1978; Shaw and Dobson, 1995; Shaw et al. 1998). This aggregation was well described by a negative binomial distribution as is the case with a variety of other parasites (see Shaw and Dobson, 1995 for reviews). It should be noted that the negative binomial distribution of fleas across host individuals holds on whatever scale of consideration we selected, namely within host within flea species, within host across flea species and within flea species across hosts. The causes and consequences of aggregated distribution in parasites will not be discussed further in this paper because they have been discussed extensively in previous publications (see Poulin, 1998a; Combes, 2001 for reviews).
The log-normal frequency distribution of flea abundances found in this study also supports previous findings (Shaw and Dobson, 1995; Morand and Guégan, 2000). Again, this pattern held at various scales; within and across flea species and within and across host species. This pattern was observed previously in communities of different parasites of a single host species (Dobson, 1990; Simkova et al. 2002) and among host species (Morand and Guégan, 2000). This confirms the ubiquity of the log-normal distribution of abundance in both endo- and ectoparasites and in both water and terrestrial ecosystems. The log-normal distribution of abundances is also the case for free-living organisms (May, 1975; Williamson and Gaston, 1999). This distribution arises because of interaction of a variety of factors affecting species abundances (see Shaw and Dobson, 1995 for the description of these factors in terms of host-parasite systems). Furthermore, clear lower and upper limits of flea abundance can be envisaged. Indeed, mean minimum abundance among flea species varied from 0·04 to 0·1 fleas per individual host. This suggests that there is some species-specific minimum abundance level below which a flea population cannot persist over time (Keymer, 1982). In addition, most flea species attained maximum mean abundance of 1·0–8·6 fleas per individual hosts, whereas mean abundances over 20 fleas per individual host were observed for 2 fleas only (S. c. pyramidis on G. pyramidum and X. conformis on M. crassus). The maximum mean abundances can reflect either flea-specific maximum transmission rate possible (Shaw and Dobson, 1995) and/or host-specific carrying capacity. The host-focused explanation seems to be more feasible in our case because maximum mean abundances of these two fleas were lower on other hosts (16·0 S. c. pyramidis per individual G. a. allenbyi and 3·3 X. conformis per individual G. dasyurus).
The distribution of flea prevalences both within host species among flea species and within host and flea species among surveys, however, did not support previous findings on this issue in other host and parasite taxa. Indeed, Poulin (1999) and Simkova et al. (2002) found negative binomial distributions of prevalences in helminth parasites of birds and fish, respectively. Morand and Guégan (2000) reported bimodal distribution of prevalences in interspecific comparative study of mammalian nematodes. Furthermore, bimodality was reported to be a common type of distribution of occupancy (equivalent of prevalence) for free-living animals (Hanski, 1991) [although this commonness has been questioned for parasites (Nee, Gregory and May, 1991)] and has been used for testing core-satellite species hypothesis (see below). Our findings, thus, do not support core-satellite hypothesis.
Taylor's power relationship
Relationships between flea mean abundance and its variance fit Taylor's power law for all scales we used. Taylor's law has been found to describe accurately the observed spatial distribution of organisms (Taylor and Taylor, 1977; Taylor and Woiwod, 1980; Perry and Taylor, 1986). The value of the exponent b is suggested as an indicator of a tendency of organisms to be mutually attracted (Perry, 1988) which seemed to be the case in our study (b values were significantly greater than unity). Furthermore, high r2 values of the regressions (0·82–0·98 in most cases) suggest that fleas seemed to be constrained in the degree of variation observed for any given mean abundance. This variation is lower than that found in free-living animals (e.g., Taylor, Woiwod and Perry, 1978) but close to that found for a variety of parasites (Shaw and Dobson, 1995). This tightness of the relationship in parasites has been suggested to result from a trade-off between being too aggregated (and hence risking extinction because of high mortality of heavily infected hosts) and being too random (and hence lowering mating opportunities) (Anderson and Gordon, 1982; Shaw and Dobson, 1995).
In addition, the b value of Taylor's relationship varied in different flea species on different hosts and even in the same flea species parasitizing different hosts (e.g., X. ramesis on M. crassus on G. dasyurus; 1·45±0·01 versus 1·27±0·10, respectively). The value of b was suggested to reflect regulation processes in host-parasite systems and to be inversely indicative of parasite-induced host mortality (Anderson and Gordon, 1982; Madhavi and Anderson, 1985). Indeed, parasitism by X. ramesis leads to higher energy expenditure of G. dasyurus (Khokhlova et al. 2002), but not of M. crassus (I. S. Khokhlova, B. R. Krasnov and A. A. Degen, unpublished data). In addition, the highest b values were observed for host specific fleas, P. chephrenis and S. c. pyramidis. Moreover, mean values of b were significantly higher in host-specific than host- opportunistic fleas. This can be explained by relatively low virulence of a host-specific parasite because a specific host and its descendants are necessary to maintain a parasite population. However, the common view that parasite-induced host mortality can decline during long evolutionary history of the relationship between a particular host species and a particular parasite species has been strongly criticized in the context of optimal virulence (Poulin, 1998a; Combes, 2001). The selection of a strategy of host exploitation by a parasite should depend on a concern for maximization of fitness rather than the welfare of the host (Poulin, 1998a). Indeed, the b value for a host-specific N. pumilionis was relatively low suggesting a highly negative effect on the host. This flea is active in 3–4 winter months only. Perhaps, it faces a trade-off between (a) the risk of being lost together with over-parasitized hosts and (b) the necessity of finishing the life-cycle during a very short time. Nevertheless, in general, mean values of b were higher in seasonal than in year-round-active fleas. This contradiction between high mean b in seasonal fleas and low b in N. pumilionis can be related to this flea being both highly host specific and strictly seasonal. Possibly, these two traits invoke different constraints on the pattern of aggregation of N. pumilionis with the effect of seasonality being stronger. In other words, variation in b values among different flea-host systems can depend on a variety of factors related to particular life-history traits that can affect demographic patterns of both parasites and hosts (Downing, 1986).
The value of b for Taylor's relationships for fleas was also negatively affected by the number of co-occurring species, supporting the results of the simulation study of Kilpatrick and Ives (2003). They demonstrated that the slopes of Taylor's power relationship average 2 where there is a lack of interspecific competition; direct or apparent (negative interactions mediated by a shared competitor) competition leads to slopes of between 1 and 2, whereas the increase of competition between species decreases the slopes from 2 to 1 (see Kilpatrick and Ives, 2003). However, the slope of the power relationship seemed to be insensitive to community size exceeding 5 species (Kilpatrick and Ives, 2003). In addition, according to Kilpatrick and Ives (2003) the slope might be the product not only of species interactions but also of how species respond to the environment relative to each other.
Relationship between abundance and prevalence
The relationship between parasite abundance and prevalence is analogous to the abundance-occupancy relationship of free-living organisms which is overwhelmingly positive (see Gaston, 1999, 2003 and references therein), although contrary cases also have been reported (Boeken and Shachak, 1998). In general, an empirical relationship between the fraction of sites (or samples) where a given species occurs and its mean density can be modelled as a logistic curve (Gaston, 1999). This relationship is supported for various taxa of parasites including both endo- (Poulin, 1999; Morand and Guégan, 2000) and ectoparasites (Krasnov et al. 2002a; Simkova et al. 2002; this study). Positive abundance-occupancy relationships were explained using various hypotheses (listed in Gaston, 2003), the most common being the sampling artifact hypothesis, core-satellite hypothesis and ecological specialization hypothesis (see Hanski, Kouki and Halkka, 1993).
The sampling hypothesis states that positive abundance-occupancy relationship can be a result of sampling artifacts simply because locally rare species are more difficult to detect than locally abundant species. Although the sampling model has been successfully applied to various sets of data on free-living organisms (Hanski et al. 1993), it was strongly criticized (Brown, 1995) and has been rejected in studies of various parasite-host systems (Morand and Guégan, 2000; Krasnov et al. 2002a).
The core-satellite hypothesis predicts bimodal distribution of organisms in their environment (Hanski, 1982) with most species occurring either in a large proportion or only in a small proportion of patches. Although we did not test the core-satellite hypothesis in this study, the absence of bimodality in prevalence distribution suggested rejection of this hypothesis. Studies of other parasite-host systems also rejected core-satellite hypothesis as a possible explanation of positive prevalence-abundance relationships (Poulin, 1998b; Simkova et al. 2002).
Ecological specialization
Abundance and prevalence of ‘specialist’ fleas were either higher or lower than those of ‘generalist’ fleas depending on which parameter of flea biology was selected as an indicator of specialization. This can explain, at least partly, the discrepancy in the results of the studies that compared abundance between ‘specialists’ and ‘generalists’ in different taxa (e.g., Morand and Guégan, 2000 vs Simkova et al. 2002) or at a different scale within the same taxon (this study vs Krasnov et al. 2004).
Moreover, studies of different taxa that used the same parameter to describe specialization provided sharply different results. For example, Poulin (1998b) observed a negative relationship between the number of fish host species used by 188 species of metazoan parasites and their average abundance in hosts. This apparent trade-off between the number of host species exploited and the abundance achieved by parasites in these hosts was explained by the presumably high cost of parasite adaptations to multiple host defence mechanisms (trade-off hypothesis – Poulin, 1998b). In contrast, neither prevalence nor abundance were affected by host specificity for intestinal nematodes of mammals (Morand and Guégan, 2000), indicating that the diversity of host defence mechanisms does not seem to affect abundance of nematodes. Yet, studies of Barger and Esch (2002) and Simkova et al. (2002) demonstrated positive relationships between host specificity and abundance, thus supporting the resource range hypothesis of Brown (1984).
The divergence in the abundance-specialization relationship in dependence on parameter of specialization again suggests that a variety of life-history traits can affect the pattern of parasite abundance and distribution.
Simple epidemiological model predicts flea distribution
Our results demonstrated that a simple epidemiological model can, in general, successfully predict the occurrence of fleas in populations of their hosts. However, this appeared to be true only after controlling for host sample size. The epidemiological model we used (e.g. Anderson and May, 1978) takes into account demographic parameters interacting with spatial heterogeneity (Anderson et al. 1982). Parameters of this model appear to be, generally, the most parsimonious set of factors that explain much of the variance in flea prevalence without involving other mechanisms such as the degree of flea host specificity or the level of host resistance against flea parasitism. Nevertheless, in some cases (S. c. pyramidis on G. al. allenbyi and G. pyramidum, P. chephrenis on A. cahirinus, X. dipodilli and N. i. theodori on S. calurus), observed flea prevalence was higher than that predicted from the epidemiological model. In other words, there were more infested animals than expected from the model. In most of these cases, the discrepancy of the observed prevalence from that predicted by the model can be explained by either a relatively low negative effect of flea parasitism on a host and/or strong resistance of a host to flea parasitism. For example, parasitism of S. c. pyramidis has no negative effect of metabolic rate on G. a. allenbyi (H. Hawlena, Z. Abramsky, B. R. Krasnov and I. S. Khokhlova, unpublished data), whereas it does in the case of P. chephrenis and A. cahirinus (I. S. Khokhlova, M. Sarfati, B. R. Krasnov and A. A. Degen, unpublished data). G. a. allenbyi generates a strong immune response to the antigens from S. c. pyramidis (Khokhlova et al. 2004). Higher than expected prevalence can be also be related to a close tie between a particular flea and a particular host. G. a. allenbyi and G. pyramidum are very closely related species (belonging to the same subgenus) and are the only hosts of S. c. pyramidis, whereas A. cahirinus is a specific host of P. chephrenis. S. calurus is the main (and sometimes the only) host for X. dipodilli in rocky habitats (this host is absent in other habitats inhabited by X dipodilli). In addition, the unexplained portion of variance by the epidemiological model (which can attain as much as 39%) can be explained by specific relationships in a particular flea-host system as well as some, still unknown, biological properties of flea and/or host.
No deviation from the model was found when the data were pooled across hosts for each flea species. This supports the notion that a trend that arises from interspecific comparison can mask true relationships, in particular host-parasite systems. In addition, the epidemiological model predicted accurately the observed prevalence only when it was corrected for host sample size. This confirms the important effect of host density on flea distribution (e.g. Krasnov et al. 2002a; Stanko et al. 2002). Models that used k values calculated using moment estimation predicted prevalence better than those that used k values calculated from Taylor's power relationship. This suggests that calculation of k from Taylor's relationship should be used cautiously, although it was used successfully in other studies (Morand and Guégan, 2000; Simkova et al. 2002).
In conclusion, the results of this study demonstrated that (a) flea distribution is aggregated within populations of host species as well as across whole host communities; (b) flea abundance is positively correlated with flea prevalence; both these parameters were higher in host-specific fleas than in host-opportunistic fleas and in fleas that are active all year round than in seasonal fleas; (c) the slopes of the Taylor's relationships and, consequently, aggregation levels tend to be higher in specialized flea species and tend to decrease with increasing richness of flea assemblages; and (d) accounting for the host sample size in the estimation of the aggregation parameter improves the estimations of the prevalence given by the simple epidemiological model.
We thank Allan Degen, A. Marm Kilpatrick and two anonymous referees for helpful comments on the earlier draft of the manuscript. This study was partly supported by the Israel Science Foundation (Grant no. 249/04 to B.R.K., G.I.S. and I.S.K.). This is publication no. 183 of the Ramon Science Center and no. 487 of the Mitrani Department of Desert Ecology.
