Abstract
To estimate the basic reproduction number (R 0) of Borrelia lusitaniae and Borrelia afzelii, we formulated a mathematical model considering the interactions among the tick vector, vertebrate hosts, and pathogens in a 500-ha enclosed natural reserve on Le Cerbaie hills, Tuscany, central Italy. In the study area, Ixodes ricinus were abundant and were found infected by B. lusitaniae and B. afzelii. Lizards (Podarcis spp.) and mice (Apodemus spp.), respectively, are the reservoir hosts of these two Borrelia burgdorferi sensu lato (s.l.) genospecies and compete for immature ticks. B. lusitaniae R 0 estimation is in agreement with field observations, indicating the maintenance and diffusion of this genospecies in the study area, where lizards are abundant and highly infested by I. ricinus immature stages. In fact, B. lusitaniae shows a focal distribution in areas where the tick vector and the vertebrate reservoir coexist. Mouse population dynamics and their relatively low suitability as hosts for nymphs seem to determine, on the other hand, a less efficient transmission of B. afzelii, whose R 0 differs between scenarios in the study area. Considering host population dynamics, the proposed model suggests that, given a certain combination of the two host population sizes, both spirochete genospecies can coexist in our study area. Additional incompetent hosts for B. burgdorferi s.l. have a negative effect on B. afzelii maintenance, whose R 0 results > 1 only with high mouse population densities and/or low lizards abundance, but they do not seem to influence B. lusitaniae transmission cycle on Le Cerbaie. Secondly, our model confirms the importance of nymphs' infestation, of host population density and diversity, and spirochetes host association for the maintenance of the transmission cycle of B. burgdorferi s.l.
Introduction
The vector of LB in Europe is Ixodes ricinus, a generalist tick feeding on a wide variety of mammalian, avian, and reptilian species. The tick life cycle is regulated by biotic and abiotic factors that determine its presence, abundance, and phenology in different habitats (Randolph 1997, 2002, Bisanzio et al. 2008). Immature ticks generally parasitize small vertebrates and are the key stages for the maintenance of LB enzootic cycle (Van Buskirk and Ostfeld 1995). As the transovarial route is considered a negligible way of transmission (Patrican 1997, Matuschka et al. 1998), larvae acquire B. burgdorferi by feeding on infected hosts. The spirochete may also pass from nymphs to larvae that are active in the same period and feed in proximity on the skin of the host (cofeeding), even in the absence of systemic infection of the vertebrate (Randolph and Rogers 2006).
Hosts differ in their capacity of harboring ticks and being infected by B. burgdorferi s.l., whose genospecies have a peculiar association to various vertebrate species (Kurtenbach et al. 1998, 2002). Consequently, at a specific geographic location, the composition of the wildlife host population determines the presence and the relative frequency of each genospecies.
In Europe, small mammals are the major reservoir for Borrelia afzelii and B. burgdorferi sensu stricto, whereas birds are reservoirs for Borrelia garinii and Borrelia valaisiana (Hanincova et al. 2003b). Lizards are the most likely reservoirs for Borrelia lusitaniae (Younsi et al. 2005, Amore et al. 2007, Majlathova et al. 2008). This genospecies has been mainly found in host-seeking I. ricinus in foci in Southern Europe and Northern Africa, corresponding to the southern limit of the tick vector geographic range. This focal geographical distribution might be explained by the necessity of the simultaneous presence of I. ricinus and lizards. Indeed, on Le Cerbaie Hills, Tuscany, central Italy, B. lusitaniae was the dominant B. burgdorferi s.l. genospecies in I. ricinus in 2004, when it accounted for 82.9% of infected host-seeking ticks. Other genospecies were found in the same study location at lower frequency levels. In particular, B. afzelii was relatively infrequent in host-seeking ticks (2.4% of infection; Bertolotti et al. 2006). In 2005, wild rodents and lizards were trapped at the same location: B. lusitaniae was the only genospecies detected in feeding larvae and tissues from lizards, while B. afzelii was not found in tissues and feeding larvae from Apodemus spp., which is known to be the main reservoir host for this genospecies. Prevalence of infestation by nymphs was significantly greater on lizards than on mice, whereas levels of infestation by larvae were similar (Amore et al. 2007). The observed scenario was consistent with the hypothesis, generated in previous studies (Richter and Matuschka 2006), that relatively abundant lizard populations may act as a limiting factor for the transmission of genospecies other than B. lusitaniae, since they feed large fractions of immature ticks, thus subtracting these from other vertebrates such as small rodents.
A different scenario was observed in Baden-Wurttemberg, Germany, where both B. lusitaniae and B. afzelii were identified in vertebrate hosts with relatively high levels of infection (Richter and Matuschka 2006). Such a difference might be explained by a generally greater level of B. burgdorferi s.l. transmission in Germany than in Tuscany. Moreover, relative abundance of vertebrate hosts may differ across geographic locations, affecting intensity of B. burgdorferi s.l. transmission and accounting for the diffusion of different genospecies.
Models can be used to point out key features of a complex system, such as the transmission dynamics of B. burgdorferi s.l. genospecies. They can also guide further field research by generating hypotheses and identifying specific lack of knowledge.
It is well known that the basic reproduction number (R 0) is a useful parameter to study the transmission and maintenance of infectious agents. Several mathematical models focusing on R 0 have been generated to study tick-borne pathogens in the past few years, some referring to LB in particular. Randolph (1998) described a model that is easily understandable but not directly applicable to our scenario, since it does not to consider the simultaneous presence of different B. burgdorferi genospecies. Other models (Rosa et al. 2003, Hartemink et al. 2008, Pugliese and Rosa 2008) capture all factors influencing the pathogen transmission dynamics, but include complex formulations that are hardly applicable to fieldwork data.
In 2006, we estimated infestation by immature I. ricinus on mice and lizards in two different habitat types on Le Cerbaie, and we tested larval ticks for B. burgdoferi s.l. Data were included in a simple mathematical model to identify specific factors associated with host population dynamics affecting R 0 and thresholds in host abundance that determine the maintenance (R 0 > 1) or extinction (R 0 < 1) of B. lusitaniae and B. afzelii on Le Cerbaie.
Materials and Methods
Field work
The study was conducted in a natural reserve located on Le Cerbaie Hills, Pisa Province, central Italy. It is characterized by bottomlands, relatively humid and covered by deciduous trees, and uplands, which are drier and whose vegetation is typically Mediterranean (Bisanzio et al. 2008).
Small rodents were trapped as described in Amore et al. (2007) during monthly trapping sessions (March–August 2006) for two consecutive trap-nights. We used 180 live traps, set 10 m apart, in two 9-by-10 grids (9000 m2); one was located in a typical upland site and one in a bottomland site. The bottomland site was the same sampled in 2005 (Amore et al. 2007). Captured mice were anesthetized and examined; before releasing, they were marked by tattooing to permit individual identification. Population density was estimated from capture data using CAPTURE software (Otis et al. 1978, White et al. 1982). Lizards were captured by a noose affixed to a stick, in the same sites of mouse captures, and examined as described in Amore et al. (2007).
Prevalence and 95% exact binomial confidence intervals (CIs) of infestation by immature ticks were calculated by vertebrate species [BINOMIAL option, PROC FREQ; SAS Institute, (Cary, North Carolina) 1999]. Recaptures of the same rodents in the second trap-night of the same session were excluded from the analysis. Prevalence of infestation in lizards and mice was compared by the Fisher exact test. Mean numbers of ticks per host, 95% CIs, and negative binomial dispersion parameters were obtained using intercept-only generalized linear models (GLMs) in the SAS system. Negative binomial error was used to take into account the aggregated distribution of ticks among hosts. GLMs were used to compare mean tick infestation among host species. Model checking was accomplished by goodness-of-fit statistics (Littell et al. 2002). The effect of trapping site on the probability of infestation of mice by I. ricinus nymphs was tested by logistic regression analysis (PROC LOGISTIC; SAS Institute, 1999). To adjust for the effect of the month of capture, the seasonal pattern of infestation was included as a sinusoidal fluctuation, with amplitude of 1 (peak in April) and period of 1 year (Bisanzio et al. 2008).
The degree of concurrent infestation by at least one I. ricinus larva and nymph on the same individual host was tested by Kappa statistics (Fleiss 1981). Kappa is commonly used as a measure of agreement between categorical classification criteria, and a value not significantly different from zero indicates no agreement beyond chance.
Infection of B. burgdorferi s.l. ticks was investigated using polymerase chain reaction, as previously described (Mannelli et al. 1999, Amore et al. 2007). At least one tick from every captured animal was tested, to investigate the largest number of individuals. Infesting ticks were chosen considering their engorgement index, which is proportional to their host attachment duration and consequently to the probability of being infected (Yeh et al. 1995). Prevalence of polymerase chain reaction–positive ticks was calculated by host species, vector stage, and capture grid. The proportions of lizards and mice harboring B. burgdorferi s.l.–infected ticks were compared by the Fisher exact test. Analyses were performed using the R software (R Development Core Team 2007).
Mathematical model
A mathematical model was derived to qualitatively estimate R 0 for the different genospecies of B. burgdorferi s.l. from field data, capturing their interaction with I. ricinus immature stages and vertebrate hosts.
In the case of infections that spread via host and vector agents, R 0 may be defined as the expected number of infected hosts directly infected by a single host inserted in a susceptible population of hosts and vectors.
This refers to the number of infected hosts after a single generation of transmission. In fact, in such a scenario, some host-seeking uninfected larvae will feed on the infectious individual. A portion of them will get infected with the host's specific B. burgdorferi genospecies, and will drop off the host and molt, thus transmitting the infection to their next stage. Infected nymphs will then need to infest and feed on a susceptible competent host to transmit the infection. The expected number of resulting infected hosts represents, in our framework, a qualitative estimation of R 0 for the B. burgdorferi s.l. genospecies under consideration.
According to this scheme, R
0 can be expressed as the expected number of infected larvae from a single infected host times the reduction factor from infected larvae to infected nymphs, multiplied by the reduction factor from infected nymphs to infected hosts:
where Φ(L,Hc,αLH,DH,τL
) is the number of susceptible larvae that feed on the infected host during the host's infectious time period. It is a function of the number of larvae (L), the number of susceptible hosts (Hc
), the probability that a larva feeds on a host (α
LH
), the duration of the host infectious period (DH
), and the duration of the larva's feeding period (τ
L
);
βHL, βNH
, and βLN
are, respectively, the transmission probabilities from an infected host to a susceptible larva, from an infected nymph to a susceptible competent host, and from an infected larva to its nymph stage;
sN
is the survival probability from feeding larva to feeding nymph; Ψ(N,Hc,Hnc
) is the proportion of competent hosts on which the infected nymphs feed. It is a function of the expected number of nymphs (N) feeding on competent (Hc
) and incompetent (Hnc
) hosts.
The competition between hosts competent for different B. burgdorferi genospecies is captured by Ψ(N,Hc,Hnc ). In fact, hosts only compete for those nymphs that carry the infection they are competent for: host-seeking infected nymphs can feed on any host, but only those feeding on the competent ones lead to the pathogen transmission.
This model mainly resides on two assumptions: Every infected nymph mounts and feeds on a different host. Only systemic infection transmission is considered.
While the first assumption can be easily justified, since we consider a framework where a large number of potential uninfected hosts for the ticks exist, the other one should be carefully considered and discussed when applying this model to a real scenario. In fact, different host species can lead the system to strongly violate the latter assumption, reducing the predictive power of the model. In particular, for those ecological systems where cofeeding transmission plays an important role in the epidemiological dynamics (Ogden et al. 1997), the proposed model should be seen as a lower bound estimation of basic reproductive number.
Results
Mouse host populations were estimated in both the 9000 m2 trapping grid locations. A total of 90 and 54 Apodemus spp. were trapped in the upland and the bottomland sites, respectively. In both scenarios, collected data positively passed CAPTURE software's closure test, allowing the estimation of mouse population density. Under the appropriate estimator determined by CAPTURE (Chao for the M(th) model in the upland site, and Jackknife for the M(h) model in the bottomland grid), mouse population densities have been calculated. These resulted in 116.12 individuals per hectare
During 19 trapping sessions, 86 Podarcis spp. have been captured, 40 of them in the upland and 46 in the bottomland. Due to the subjectivity of the collection method, lizard population abundance could not be consistently estimated and we referred to bibliographic data in Mediterranean habitats, ranging from 100
We collected 638 larvae and 26 nymphs on mice and 343 larvae and 211 nymphs on lizards. Host infestation results are summarized in Table 1. Prevalence of larvae infestation on lizards and that on mice were not significantly different in both capture grids, while this difference was significant as regard nymphs infestation (p < 0.001). Logistic regression analysis showed that the probability of infestation by nymphs on mice was greater in the upland than in the bottomland trapping grid (adjusted odds ratio = 7.8), and such an effect was at the threshold of statistical significance (p = 0.05).
The mean numbers of nymphs feeding on lizards and mice in each capture grid are illustrated in Table 2. Based on goodness-of-fit statistics and residuals plots, the GLMs with negative binomial error were appropriate for the analysis of counts of immature ticks on vertebrate hosts. No significant difference was found between the mean numbers of larvae on lizards and mice (p = 0.6), while nymphs were significantly more abundant on lizards (p < 0.0001).
Confidence interval cannot be calculated.
Kappa, as a measure of coinfestation of lizards by immature I. ricinus, was −0.03 for lizards (95% CI: −0.23, 0.17) and 0.07 for mice (95% CI: 0.03, 0.12), indicating no evidence of cofeeding of larvae and nymphs on the same individual hosts.
B. lusitaniae and B. afzelii were identified, respectively, in Podarcis spp. and Apodemus spp. infesting ticks. In particular, B. lusitaniae was detected in ticks collected on lizards in both sampled grids, while B. afzelii only in ticks from mice in the upland site (Table 3). No significant differences in B. burgdorferi s.l. infection prevalence in ticks were found by comparing host species (Fisher Exact test, p > 0.05).
To apply the proposed formulation for R 0, the different parameters in the formula were estimated. When they could be considered independent from the peculiar ecological conditions, we decided to refer to published results. On the other hand, parameters influenced by habitat factors typical of the study area were derived from the field study data.
Transmission probabilities are commonly agreed to be: βHL × βLN = 0.5 for B. lusitaniae (Majlathova et al. 2006) and 0.4 for B. afzelii (Hanincova et al. 2003a), and βNH = 0.67 (Donahue et al. 1987). Another parameter that was not directly estimable for our field study data is sN , the expected fraction of feeding larvae that will feed in their nymphal stage on a host. For it we have used the value 0.1, which has been previously proposed (Randolph 1993, 1994, Hartemink et al. 2008).
The number of susceptible larvae that feed on the infected host during its infectious period can be estimated from field data as the expected number of larvae on captured hosts multiplied by the ratio between the duration of the infectious time period of a host DH
= 122 days (Randolph et al. 1996) and the duration of the nutrition of a larva on a host τ
L
= 2.5 (Hartemink et al. 2008):
where
As regard parameters taken from literature, unfortunately they have been reported only as point estimations by the authors, who did not discuss any variability of their values. It should be noted that the model estimation is directly proportional to βHL, βLN, βNH, sN , and DH , but inversely proportional to τ L . Therefore, any change in the evaluation of these parameters should directly reflect on the R 0 calculation. In such a case, the R 0 estimation must be multiplied by the ratio between the new and the initial values of the parameter, except for τ L , where the ratio shall be between the initial and the new values.
Finally, the proportion of competent hosts on which the infected nymphs feed Ψ(N,Hc,Hnc
) can be seen as the average number of nymphs feeding on captured competent hosts times the estimated number of competent hosts, over the estimated number of nymphs on all hosts in the ecosystem:
where
NH
is the number of nymphs on each host.
While the numerator is straightforwardly drawn from field data and population estimation of competent hosts, the denominator can be calculated as the sum of the estimated number of nymphs feeding on competent and on incompetent hosts. This takes into account the role of the population of incompetent vertebrates in the transmission dynamics of different B. burgdorferi s.l. genospecies. To highlight the role of lizards in B. afzelii transmission (and, on the other side, of mice for B. lusitaniae), the factor concerning incompetent hosts (in the denominator) can be seen as the sum of the estimated number of nymphs feeding on lizards (respectively, mice) and the estimated number of nymphs feeding on all other incompetent hosts (i.e., all those not being lizards, or mice). Thus, we propose this denominator to be composed of the sum of the number of nymphs feeding on known competent hosts (lizards or mice), on known incompetent hosts (mice or lizards, respectively), and on all other unknown incompetent hosts (as a function of their population density):
for the parameter relative to the R
0 estimation of B. afzelii (Ψ
afz
) and B. lusitaniae (Ψ
lus
), where the subscripts As, Pm, and nc refer, respectively, to mice, lizards, and other host noncompetent for both considered genospecies. In the Ψ
afz
and Ψ
lus
estimations in the numerator
The numerical framework that describes B. burgdorferi transmission dynamics is drawn in Figure 1 for the upland (a–c) and the bottomland (d–f) grids. It depicts the epidemic thresholds for B. afzelii and B. lusitaniae (i.e., R
0 = 1) in the two scenarios on Le Cerbaie, according to our model. Being not possible to report the point that corresponds to the combination of the exact estimations of both mouse and lizard populations, in all the figures only nymph-infested mouse population density is reported, together with its s.e. estimation, and nymph-infested lizard population is shown in the y-axis, with a maximum density value of 100

Estimation of the epidemic outbreak threshold (i.e., basic reproduction number [R
0] = 1) on Le Cerbaie. The solid line corresponds to R
0 = 1 for Borrelia lusitaniae, while the dashed line corresponds to R
0 = 1 for Borrelia afzelii. Light gray areas around the two lines represent the effects of the negative binomial 95% confidence intervals in mean nymph's infestations. The area above the solid line corresponds to R
0 > 1 for B. lusitaniae, while the area below the dashed line corresponds to R
0 > 1 for B. afzelii. The area between the two lines corresponds to R
0 > 1 for both B. lusitaniae and B. afzelii. Densities of nymph-infested mouse population estimated in the study area are reported, together with their standard error estimations (dash-dotted vertical lines and dashed vertical areas). The upland scenario is depicted in the upper row (
Discussion
The perpetuation of the complex LB system in the environment needs an equilibrium between its main components: tick infesting and questing populations, host ecological dynamics, and B. burgdorferi s.l. infection. Our model was built to identify the key points for the maintenance of two B. burgdorferi genospecies in a focus of B. lusitaniae in central Italy. Our aim was exploring the mechanisms underlying field observations in two real scenarios on Le Cerbaie, and pointing out priorities for future investigations on LB eco-epidemiology.
In previous studies on Le Cerbaie, extensive host-seeking tick collection highlighted the differences in the acarologic risk between the two main habitat types in the area and the predominance of B. lusitaniae (Bisanzio et al. 2008). This genospecies was constantly found in the study area in host-seeking ticks and in lizards' ticks, both in 2005 in the bottomland site and in 2006 in the bottomland and upland sites. On the other hand, B. afzelii was only detected with low prevalence values in host-seeking ticks and it was not present in the ticks from mice collected in bottomland capture grid in 2005 (Amore et al. 2007), neither in 2006 in the same location. This genospecies was only identified in immature ticks feeding on mice captured in the upland site during the 2006 field work. Moreover, in the upland we detected higher mouse densities and a significantly higher infestation by nymphs on mice compared to bottomland. These observations suggest that the maintenance of B. burgdorferi s.l. is influenced by factors such as reservoir host population densities and their infestation by nymphs. Indeed, fluctuations in vertebrate host populations might change the relative importance of different genospecies. Apodemus spp. mice, the known reservoirs for B. afzelii, are widespread in different habitats at all latitudes (Ouin et al. 2000, Michelat and Giraudoux 2006), but their population dynamics are characterized by within- and among-year fluctuations (Jamon 1986, Marcstrom et al. 1990, Butet 1994) influencing B. burgdorferi s.l. maintenance (Mather et al. 1989, Brisson and Dykhuizen 2006). Lizards are a major component of wildlife in Mediterranean areas (Avery 1978), where they allow B. lusitaniae maintenance in habitats that are suitable to its vector (Majlathova et al. 2008).
Another essential element to be considered in the LB cycle is the vertebrate's suitability to host nymphs. On Le Cerbaie, we found out that lizards feed more nymphal ticks than rodents (Tables 1 and 2), and this result is supported by the 2005 study, when only 3 out of 53 mice were infested with 1 nymph each (Amore et al. 2007). Lizards have a zooprophylactic effect on genospecies other than B. lusitaniae (Lane 1990, Matuschka et al. 1998, Kuo et al. 2000) and, when abundant, can reduce the transmission of other genospecies (Richter and Matuschka 2006). Anyway, the lizards' zooprophylactic effect seems not sufficient to explain the low and inconstant levels of B. afzelii in our study area. As evaluated by dragging, host-seeking nymphs are more abundant in the bottomland, which appears to be the most suitable habitat for ticks on Le Cerbaie (Bisanzio et al. 2008), but nymph infestation on captured mice is higher in the upland. This leads to assume a different host-seeking behavior of nymphs according to different habitats. For example, the lower humidity in the upland could determine the displacement of nymphs to the litter zone to avoid desiccation, or a greater activity in colder moments of the day (Perret et al. 2003) so that nymphs encounter and feed more easily on mice in the upland than in the bottomland.
Moreover, the presence of other vertebrate hosts that are not competent reservoirs for B. burgdorferi s.l. can dilute the infection rate in ticks. In our study area, a scarce ungulate population is present, including <15 roe deer (Capreolus capreolus L.) and a few wild boars (Sus scrofa L.). As we could not quantify their role as hosts for ticks, we considered scenarios where populations of incompetent hosts harbor up to 1/3 and 1/2 of all feeding nymphs.
Our results show that B. lusitaniae transmission cycle is robust on Le Cerbaie, where it seems to be maintained even in the case of low lizard densities and presence of noncompetent vertebrate hosts for the vector, especially in the upland (Fig. 1a–c). On the other hand, B. afzelii shows to be more sensitive to host population fluctuations, and cannot be maintained in the presence of low mouse densities, few contacts between nymphs and mice, and a large population of incompetent hosts. Indeed, the model shows that the presence of B. afzelii in the bottomland is not likely according to our mouse population estimations (Fig. 1d–f). On the contrary, B. afzelii seems to be maintained in the upland, especially if mouse densities are closer to the upper bound of our estimation and there is a weak competition of other hosts for nymphs. As the presence of these competing vertebrate species cannot be excluded and we detected B. afzelii in the upland, we might presume that mouse density is higher than what we estimated by captures, or that lizard population is not so abundant.
All these considerations lead us to think that B. afzelii can be maintained together with B. lusitaniae only in areas where high transmission burdens of the B. burgdorferi s.l. complex exist.
In conclusion, our model helped to gain insight into the ecological conditions influencing the transmission cycles of two different B. burgdorferi genospecies on Le Cerbaie, where B. lusitaniae exists in a stable focus. It confirmed that the habitat suitable to lizards and I. ricinus ticks is essential for the perpetuation of B. lusitaniae, but, on the other side, limits the transmission cycle of other B. burgdorferi genospecies, such as B. afzelii. Indeed, the maintenance of B. afzelii differs between different scenarios on Le Cerbaie, being strongly influenced by the densities of vertebrate host population and their pattern of infestation by nymphs.
The model is in agreement with field observations and explains the mechanisms at the basis of the presence of LB in the two studied sites. However, further investigations are needed to determine the robustness of the model to variations in the parameter values that we drawn from bibliographic data. For example, it would be important to assess the infectious time period of lizards and mice, as it could influence the probability of the two genospecies' transmission and maintenance. Moreover, studies on an extensive area are needed to evaluate the habitat effect on the biological cycles of B. burgdorferi on Le Cerbaie. In particular, it would be interesting to determine if the observed data refer to the single upland and bottomland studied sites, or if they reflect a pattern in the reserve and can be generalized to other habitats. More information on the lizard population density, specific spatial and temporal mouse population fluctuations, host-seeking behavior, and the role of other vertebrate species as hosts for nymphs would be also useful to better understand the pathogen ecology and epidemiology.
Footnotes
Disclosure Statement
No competing financial interests exist.
