Seasonal dynamics of tick burden and associated Borrelia burgdorferi s.l. and Borrelia miyamotoi infections in rodents in a Dutch forest ecosystem

Ixodes ricinus ticks transmit Borrelia burgdorferi sensu lato (s.l.) as well as Borrelia miyamotoi. Larvae become infected when feeding on infected rodents, with horizontal transmission of B. burgdorferi and horizontal and vertical transmission of B. miyamotoi. We studied seasonal dynamics of infection rates of I. ricinus and their rodent hosts, and hence transmission risk of these two distinctly different Borrelia species. Rodents were live-trapped and inspected for ticks from May to November in 2013 and 2014 in a forest in The Netherlands. Trapped rodents were temporarily housed in the laboratory and detached ticks were collected. Borrelia infections were determined from the trapped rodents and collected ticks. Borrelia burgdorferi s.l. and B. miyamotoi were found in ticks as well as in rodents. Rodent density was higher in 2014, whereas tick burden as well as the Borrelia infection rates in rodents were higher in 2013. The density of B. miyamotoi-infected nymphs did not differ between the years. Tick burdens were higher on Apodemus sylvaticus than on Myodes glareolus, and higher on males than on females. Borrelia-infection rate of rodents varied strongly seasonally, peaking in summer. As the larval tick burden also peaked in summer, the generation of infected nymphs was highest in summer. We conclude that the heterogeneity of environmental and host-specific factors affects the seasonal transmission of Borrelia spp., and that these effects act more strongly on horizontally transmitted B. burgdorferi spp. than on the vertically transmitted B. miyamotoi. Supplementary Information The online version contains supplementary material available at 10.1007/s10493-022-00720-z.


Introduction
Lyme borreliosis is caused by an infection with Borrelia burgdorferi sensu lato (s.l.), which is, in Europe, transmitted by the sheep tick Ixodes ricinus. The same tick species also transmits Borrelia miyamotoi, which can cause Hard Tick-borne Relapsing 1 3 Fever (HTRF). After hatching, I. ricinus larvae search for a vertebrate host. Rodents are commonly used as blood hosts by larvae and B. burgdorferi s.l. spirochaetes can be transmitted from a persistently infected rodent to feeding larvae (Hofmeester et al. 2016). Infected engorged larvae moult into infected nymphs, which can transmit the spirochaetes to rodents, thus closing the B. burgdorferi infection cycle among rodents. In contrast to B. burgdorferi s.l., B. miyamotoi is not only transmitted from transiently infectious rodents to larvae (Burri et al. 2014), but also transovarially in the ticks (Hauck et al. 2020;Scoles et al. 2001). The lifecycle of B. miyamotoi is, therefore, less dependent on rodents (Barbour et al. 2009). Ixodes ricinus can, at a low level, be co-infected with both pathogens (Kiewra et al. 2014;Kjelland et al. 2015). Ticks are highly aggregated among rodents and a minority of rodents feeds the majority of the ticks (Perkins et al. 2003), following the Pareto principle or the 20/80 rule (Woolhouse et al. 1997). A high burden of infected ticks increases the likelihood of rodents becoming infected with B. burgdorferi s.l., and a high larval tick burden on infected rodents increases the infection prevalence in questing nymphs. It is expected that the contribution of rodents to the density of infected nymphs is, therefore, also heterogeneous (Harrison and Bennett 2012;Krawczyk et al. 2020). The contribution of an individual rodent to the density of B. burgdorferi-infected nymphs is affected by (1) the number of larvae that feed on the rodent (larval tick burden), (2) the chance that the rodent is infected (rodent infection rate), (3) the chance that the spirochaetes are transmitted to feeding larvae (rodent infectivity), and (4) the chance that the engorged larvae moult into nymphs (moulting success). As the transmission route of B. miyamotoi goes via two routes, from adult tick to larva (via the egg) or from rodent to larva, the density of B. miyamotoi-infected nymphs is less dependent on rodent density. Transmission, therefore, may be affected by various biotic and abiotic factors, which may vary temporally and geographically (Li et al. 2012;Perez et al. 2012).
Rodent species differ in their suitability as hosts for ticks and for B. burgdorferi s.l. (Kurtenbach et al. 1994;Matuschka et al. 1992;Nilsson and Lundqvist 1978). Tick burden also varies within rodent species, which may be due to differences in rodent characteristics (age, sex, immune status), microclimate or vegetation composition (Gassner et al. 2011;Randolph and Storey 1999). Climate and season are known to influence the behaviour of both ticks and rodents, and can also affect tick burdens on rodents (Estrada-Peña et al. 2013;Hubalek et al. 1994). The spatial aggregation of questing larvae in the environment may also contribute to the heterogeneity in larval tick burdens. The tick burden on white-footed mice (Peromyscus leucopus) was, however, consistent over time per individual mouse and was only slightly affected by this spatial aggregation (Devevey and Brisson 2012), suggesting a large effect of the rodent characteristics on rodent tick burdens. The interactions between ticks, rodents and B. burgdorferi s.l. are complex (Van Duijvendijk et al. 2015), and little is known about the ecology of B. miyamotoi. Most studies investigated only one or a few variables that potentially explain tick burdens on rodents, and hence the transmission dynamics of Borrelia spp. These variables may, however, also interact with each other. The relative contribution of each variable to Borrelia transmission can, therefore, not be compared, but is nevertheless important for understanding the ecology of the parasite communities.
The aim of this study was to quantify the infection rates of two Borrelia species, B. burgdorferi s.l. and B. miyamotoi, as well as their transmission rates in I. ricinus and in their rodent hosts. The study was conducted over several successive months during 2 years, so as to include the effects of spatiotemporal variation of the ecosystem on the factors being examined.

Rodent trapping
Rodents were trapped with Heslinga live traps provided with hay, grain, carrot, and a few mealworms. Traps were placed at 3-week intervals at 15:30 h and inspected the next day at 08.30 h, from May until November in 2013 and 2014 in both plots (10 trapping weeks each year). In 2013, plot B was started in July (7 trapping weeks). Traps were placed at 5 m inter-trap distance in grids of 12 × 12 traps in 2013 and 6 × 12 traps in 2014. Trapped rodents were identified to species level, sexed, weighed, and counted for ticks by searching the ears, snout, head, belly, legs, armpits, throat and tail. A small ear biopsy was collected from each captured male rodent with sterile scissors and tweezers during the first capture in 2013 and from males and females during each capture in 2014 (Sinsky and Piesman 1989). Ear biopsies were stored in 70% alcohol at -20 ºC until further analysis. Rodents were individually marked by trimming their fur in a unique pattern. The fur was cut at nine locations on the body of a rodent (sides and back, at the front, middle or rear), allowing for 256 unique combinations (2 8 ) per rodent species. A reference number and trap number were recorded for each rodent. A maximum of 15 male rodents per species was transported to the laboratory for the collection of feeding ticks (see also below). Captured shrews were released immediately. All experiments were approved by the ethical committee of Wageningen University (permit numbers 2013017 and 2014064).

Rodent housing and tick collection
Trapped rodent males that were transported to the laboratory were housed individually in Makrolon type II cages placed over water-filled pans, with petroleum jelly on the edge to prevent ticks escaping. Cages were provided with standard rodent bedding and a tissue and toilet roll for shelter. Food and water were provided ad libitum. After 3 and 6 days, engorged ticks were collected from the water, dried on filter paper for 2 h and housed in ventilated 0.2-ml Eppendorf tubes at 20 ºC, 90-95% relative humidity and L14:D10 photoperiod. Engorged ticks were checked weekly and were stored at -20 ºC after they moulted to the next developmental stage until further analysis. In 2013, all emerged ticks were identified to species level. After 6 days, each rodent was released at the place of capture.

Identification of Borrelia burgdorferi s.l. and B. miyamotoi infections
Ear biopsies and emerged nymphs were analysed individually as described by Sinsky and Piesman (1989). Briefly, nymphs were heated for 20 min at 99 °C in 100 μl 1 M ammonium hydroxide (NH 4 OH), after which they were centrifuged and heated for 20 min at 99 °C with open lids to evaporate the ammonia. The lysates were stored at 4 °C. DNA was extracted from the ear biopsies using the Qiagen DNeasy Blood & Tissue Kit (Jahfari et al. 2014). Borrelia burgdorferi s.l. DNA was detected using qPCR in the IQ Multiplex Powermix with a volume of 20 μl, containing iTaq DNA polymerase (Bio-Rad Laboratories, USA), 200 nM of each primer and 3 μl of template DNA (Heylen et al. 2013). Outer surface protein A gene (OspA) (forward primer: 5′-AAT ATT TAT TGG GAA TAG GTC TAA-3′; reverse primer: 5′-CTT TGT CTT TTT CTT TRC TTA CA-3′ and probe: 5′-Atto520-AAG CAA AAT GTT AGC AGC CTT GA-BHQ1-3′) and the B. burgdorferi s.l. flagellin gene (flaB) (forward primer: 5′-CAG AIA GAG GTT CTA TAC AIA TTG AIA TAG A-3′; reverse primer: 5′-GTG CAT TTG GTT AIA TTG YGC-3′ and probe: 5′-Atto425-CAA CTI ACA GAI GAA AXT AAI AGA ATT GCT GAI CA-Pho-3′, where X stands for an internal BHQ-1 quencher attached to thymine) were used as targets. The qPCR cycling program (using a light cycler 480 real-time PCR system; Hoffmann-La Roche, Switzerland) was performed using a two-step PCR program: Taq activation for 5 min at 95 °C followed by 60 cycles of 5 s at 94 °C and 35 s at 60 °C involving a single point measurement at 60 °C with corresponding filters, finishing with one cycle of 20 s at 37 °C. In the same qPCR, B. miyamotoi was detected with primers and probe based on the flagellin gene (flaB) for detection of the bacteria.

Data analysis
Data were analysed with a generalized linear mixed model with a binomial distribution. Rodent identification number and trap number were used as random factors in each model. First, the effects of the main factors (year, rodent species, rodent sex, rodent weight, week, week 2 , rodent infection and plot) and all possible two-factor interactions between the main factors were analysed in four full models (Appendix).
The number of larvae counted in the field was used as dependent variable in the model to analyse larval tick burden on rodents. This model included all main factors and interactions. Data from ear biopsies were used as dependent variable in the model to analyse rodent infection rate. Rodent infection and the year × rodent sex interaction were excluded from this full model. The number of B. burgdorferi s.l. infected ticks out of the number of analysed ticks that were collected from the male rodents which were found B. burgdorferi s.l.-positive by PCR was used as dependent variable in the model to analyse rodent infectivity. Rodent sex and rodent infection were excluded from this full model. The number of moulted ticks out of the number of analysed ticks that were collected from the male rodents was used as dependent variable in the model to analyse larval moulting success. Rodent sex was excluded from this full model. Second, each model was reduced by backwards elimination of the fixed factors, starting with the highest P-values until all factors were either significant, or were included in an interaction factor that was significant. Significance threshold was α = 0.05. All analyses were performed using SAS v.9.3 statistical software (SAS Institute, Cary, NC, USA).

Tick burden
In total, 3157 larvae were counted in the field on the trapped rodents (Table S1). Mean larval tick burden was 7.02 ± 0.87 in 2013 and 3.08 ± 0.22 in 2014. Larvae were highly aggregated among rodents; in 2013, 20% of the rodents fed 62.3% of the larvae, and in 2014, 19% of the rodents fed 73.4% of the larvae. In 2013, 1128 larvae were collected in the laboratory from male rodents, from which 1 086 were identified to species level after they moulted into nymphs. One of these was I. trianguliceps, all others were I. ricinus. In 2014, 1410 larvae were collected in the laboratory from male rodents, which were not identified to species level. In 2013, 33 nymphs were collected in the laboratory from male rodents, from which 17 were identified to species level after they had moulted into adults. Two of these were I. trianguliceps, all others were I. ricinus. In 2014, 32 nymphs were collected in the laboratory from male rodents, which were not identified to species level. No adult ticks were collected from the water basins, even though four (unidentified) adult ticks were observed feeding on the rodents during trapping.
Nymphal tick burden was higher in 2013 (0.11, 95% CI 0.05-0.25) than in 2014 (0.03, 95% CI 0.02-0.05; P = 0.0015; Figure S2). Rodents having a higher body weight Estimated mean larval tick burden on rodents in week 32 per A year, B rodent species, C rodent sex, and D plot. Error bars represent 95% confidence intervals. Significant differences between means are indicated with ***P < 0.001 (generalized linear mixed model) had larger nymphal tick burdens than rodents weighing less (P = 0.0004; Figure S3). With each gram of body weight, larval tick burden was multiplied by 1.19.
Rodent infection rate with B. burgdorferi s.l. was higher in 2013 (0.68, 95% CI 0.45-0.85) than in 2014 (0.16, 95% CI 0.12-0.22; P < 0.0001; Fig. 6). The difference in rodent infection rate between years varied over time (P = 0.0190; Fig. 7). There was a shorter, but higher peak in 2013 than in 2014. Rodent infection rate also varied over time (P = 0.0034; Fig. 7), peaking in summer. Rodents weighing more had a higher infection rate than those weighing less (P = 0.0011; Figure S4). For each increasing gram of body weight, rodent infection rate increased by a factor 1.09. The effect of rodent sex on rodent infection rate varied over time (P = 0.022; Figure S5). Rodent infection rate was highest in week 32-36 for both sexes.

Moulting success of larvae
In 2013, 1128 engorged larvae were collected from the male rodents in the laboratory, of which 1086 (96%) moulted into nymphs (Table S2). In 2014, 1410 engorged larvae were collected from the male rodents in the laboratory, of which 1349 (96%) moulted into nymphs. Moulting success was lower for plot A (0.91, 95% CI 0.88-0.94) than for plot B (0.95, 95% CI 0.93-0.96; P = 0.044; Figure S6). Moulting success did not differ between rodent species (P = 0.057). The effect of rodent species was affected by year (P = 0.018); in 2013, moulting success was not affected by rodent species (P = 0.74) but in 2014, moulting success was higher for larvae that fed on bank voles compared to larvae that fed on wood mice (P = 0.0016; Figure S7). Moulting success of larvae varied over time (P < 0.0001).

Contribution of rodents to the indirect estimation of the density of infected nymphs
The calculated relative numbers of infected questing nymphs per year, rodent species, plot and rodent sex are shown in Table S2. The relative number of B. burgdorferi s.l.infected questing nymphs that emerged from the larvae that fed on the rodents was 0.058 for 2013 and 0.081 for 2014. The relative number of B. miyamotoi-infected questing nymphs that emerged from the larvae that fed on the rodents was the same for 2013 and 2014 (0.00006). year. Error bars represent 95% confidence intervals. Significant differences between means are indicated with *P < 0.05, **P < 0.01; N.S.: P > 0.05 (generalized linear mixed model)

Discussion
The most common rodent species in The Netherlands, bank voles and wood mice, were the only rodent species trapped in this study (de Boer et al. 1993;Gassner et al. 2008).
Rodent density was about one order of magnitude lower in 2013 than in 2014. The contribution of the 2013 rodent population to the density of B. burgdorferi s.l.-infected nymphs was, however, only a little less than that of the 2014 population. The contribution of these rodents to the density of B. miyamotoi-infected nymphs did not differ between years. The total number of feeding larvae was higher, whereas the estimated larval and nymphal tick burdens per rodent were lower in 2014 than in 2013. These lower tick burdens caused a lower rodent infection rate with B. burgdorferi s.l. in 2014. As a result, the contribution of the 2014 population to the density of B. burgdorferi s.l.-infected nymphs increased only slightly. The lower mean larval tick burden in 2014 was probably due to a dilution effect caused by the high rodent density compared to that of 2013 (Kiffner et al. 2011;Krasnov et al. 2007;Rosa et al. 2007). Rodent density is largely affected by food availability, which also affects home range size (Stradiotto et al. 2009). A large home range size of rodents can result in more tick encounters and therefore, a higher tick burden (Devevey and Brisson 2012;Sonenshine and Stout 1968). The contribution of the rodent populations to the density of B. miyamotoi-infected nymphs, however, did not differ between the years, supporting the importance of transovarial transmission for this spirochaete (Hauck et al. 2020;Krause et al. 2015).
Our data showed that larval and nymphal tick burdens were higher on wood mice than on bank voles, which is supported by previous studies in different ecological settings (Gassner et al. 2013;Mysterud et al. 2015;Tälleklint and Jaenson 1997). In addition, more larvae were found on male rodents than on female rodents and on rodents with a large body weight, which are in general older than smaller ones, which is also supported by previous studies from different ecological settings (Devevey and Brisson 2012;Nilsson and Lundqvist 1978;Sonenshine and Stout 1968;Tälleklint and Jaenson 1997). The higher tick burdens on male compared to female rodents may be the effect of larger home ranges of males (Stradiotto et al. 2009) and therefore tick encounter rates.
Larval and nymphal tick burdens were higher in 2013 compared to 2014. Considering Lyme borreliosis risk, only infected rodents can infect feeding larvae. Rodents can become infected with B. burgdorferi s.l. through the bite of an infected larva or nymph (Radolf et al. 2012;Van Duijvendijk et al. 2016). A high tick burden on rodents, therefore, results in a high exposure to the spirochaetes. As a result of this higher tick burden, rodent infection rate with B. burgdorferi s.l. was higher in 2013 compared to 2014. The lower rodent infection rate in 2014 contributed to a small increase in the density of B. burgdorferi s.l.infected nymphs fed by the 2014 rodent population compared to the 2013 rodent population. Variation in rodent infection rate between years was also found by Kurtenbach et al. (1995) and Krawczyk et al (2020) and is most likely caused by the annual variations in the ratio of rodent/questing tick density, with a dilution effect at high rodent densities, leading to lower B. burgdorferi s.l. infections in rodents. In our study, rodent infection rate peaked during summer, which would increase the infection rate in larvae that feed on the rodent population. An increase in rodent infection rate in the spring was also found by Kurtenbach et al. (1995). These authors did not, however, find a decrease in the infection rate of feeding larvae in the autumn and speculated that this difference may have been caused by a difference in longevity of the rodents. In our study, the rodent infection rate showed a similar peak over time to that of the larval tick burden. Larval tick burden is likely to be affected by larval questing activity. Likewise, tick questing activity affected the infection rate of white footed mice (Hofmeister et al. 1999). We did not find a difference in infection rate between rodent species, which is in contradiction to the results of Kybicova et al. (2008), Humair et al. (1999) and Krawczyk et al. (2020), who found a higher infection rate in bank voles than in wood mice. It is unclear what caused this difference. Gassner et al (2013), however, reported variable infection rates between bank voles and wood mice, collected from three geographically different areas. Our molecular analysis did not distinguish between the B. burgdorferi s.l. genospecies, but previous research showed that rodents and larvae collected from rodents are mainly infected with B. afzelii (Hanincova et al. 2003 Kybicova et al. 2008), which can cause skin problems in humans (Stanek et al. 2012;Strle and Stanek 2009).
For the infected rodents to have an effect on the density of infected nymphs, the spirochaetes should also be transmitted to feeding larvae. Rodent infectivity was 0.43 based on the raw data, and estimated rodent infectivity was 0.38 when corrected for trap number and rodent number. These values are in line with the findings of Burri et al. (2014). Rodent infectivity was not affected by any of the measured variables. This is in contrast with the results of Perez et al. (2012) and Humair et al. (1999), who found a higher infectivity of bank voles than of wood mice. Infected rodents do not lose their infection, but the infectivity varies between individual rodents and was enhanced by successive infestations with larval ticks (Gern et al. 1994).
Overall, moulting success of larvae was always > 90% and was not affected by rodent species. This is in line with the findings of Tälleklint and Jaenson (1997), but in contradiction to the results of Humair et al. (1999), who found a higher moulting rate of ticks that fed on wood mice (33.3%) compared to bank voles (6.2%). We found that the effect of rodent species varied between years. Dizij and Kurtenbach (1995) showed that moulting success of larvae that fed on bank voles was reduced with successive tick infestations, whereas this was not found for larvae feeding on wood mice. Years with high tick burdens could, therefore, result in reduced moulting success for ticks that fed on bank voles, possibly caused by immunity to tick antigens. We found, however, that moulting success was significantly decreased in bank voles in 2014, when tick burdens were low.
Understanding the factors that cause the aggregation of ticks on rodents helps to understand the transmission dynamics of B. afzelii and to identify the 20% of the rodents that contribute to 80% of the B. afzelii transmission cycle (Lloyd-Smith et al. 2005;Stein 2011;Woolhouse et al., 1997). Larval density, and not rodent density was found to be the best explanatory variable for nymphal density in the following year (Rosa et al. 2007).
Larval tick burdens on rodents show high temporal variation within and between years. A peak in rodent density (e.g., after a mast year) results in reduced tick burdens on rodents and a reduced rodent infection rate.
We conclude that years with high rodent densities do not lead to an increase in the contribution of the rodent population to the density of B. burgdorferi s.l.-infected nymphs and B. miyamotoi-infected nymphs. The risks for Lyme borreliosis and HTRF are, therefore, not directly affected by a temporal variation in rodent density. The strong seasonal variation in tick burden on rodents as well as the variation in rodent infection rates, however, cause large differences in Borrelia infection risk and should be considered in public health campaigns aimed at the prevention of tick-borne disease.

3
Main factors and two-factor interactions that were included in the generalized linear mixed models for analysing the effects on larval tick burden, rodent infection, rodent infectivity, and tick moulting success and their corresponding F-and P-values