Climate and tree seed production predict the abundance of the European Lyme disease vector over a 15-year period
Parasites & Vectors volume 13, Article number: 408 (2020)
To predict the risk of tick-borne disease, it is critical to understand the ecological factors that determine the abundance of ticks. In Europe, the sheep tick (Ixodes ricinus) transmits a number of important diseases including Lyme borreliosis. The aim of this long-term study was to determine the abiotic and biotic factors driving the annual abundance of I. ricinus at a location in Switzerland where Lyme borreliosis is endemic.
Over a 15-year period (2004 to 2018), we monitored the abundance of I. ricinus ticks on a monthly basis at three different elevations on Chaumont Mountain in Neuchâtel, Switzerland. We collected climate variables in the field and from nearby weather stations. We obtained data on beech tree seed production from the literature, as the abundance of Ixodes nymphs can increase dramatically two years after a masting event. We used AIC-based model selection to determine which ecological variables drive annual variation in tick density.
We found that elevation site, year, seed production by beech trees two years prior, and mean annual relative humidity together explained 73.2% of the variation in our annual estimates of nymph density. According to the parameter estimates of our models, (i) the annual density of nymphs almost doubled over the 15-year study period, (ii) changing the beech tree seed production index from very poor mast (1) to full mast (5) increased the abundance of nymphs by 86.2% two years later, and (iii) increasing the field-collected mean annual relative humidity from 50.0 to 75.0% decreased the abundance of nymphs by 46.4% in the same year. Climate variables collected in the field were better predictors of tick abundance than those from nearby weather stations indicating the importance of the microhabitat.
From a public health perspective, the increase in nymph abundance is likely to have increased the risk of tick-borne disease in this region of Switzerland. Public health officials in Europe should be aware that seed production by deciduous trees is a critical driver of the abundance of I. ricinus, and hence the risk of tick-borne disease.
Climate change has resulted in dramatic changes in the distribution and abundance of medically important arthropod vectors, such as mosquitoes, ticks and biting flies [1, 2]. These range expansions are expected to increase the incidence of vector-borne diseases that are transmitted by these arthropod vectors with serious consequences for human health . For example, the range expansion of Ixodes ticks into northern latitudes has resulted in a dramatic increase in the incidence of Lyme borrelosis and other tick-borne diseases in Canada and Scandinavia [4,5,6,7,8,9]. In contrast, whether climate change is affecting the abundance of Ixodes ticks in areas where they are endemic is not clear [10,11,12,13,14,15,16,17].
Climate change will influence the population ecology of Ixodes ticks via its effects on both abiotic factors (e.g. climate variables) and biotic factors (e.g. vegetation or host community). Ixodes ticks spend ~ 98% of their time off the host and have to cope with seasonal changes in temperature and precipitation. Their life history traits (development, survival and reproduction) are highly sensitive to climate variables that are changing due to global warming . For example, tick development rates and survival increase with temperature and relative humidity, respectively [19,20,21,22]. However, we currently do not have a good understanding of how climate change will influence the abundance of Ixodes ticks in endemic areas [2, 20].
With respect to biotic factors, climate change is expected to drive changes in plant and animal communities that will feed back on the population ecology of Ixodes ticks. For example, climate influences seed production in beech and oak trees [23,24,25,26], which affects the population dynamics of animal species that use these seeds as a food resource, such as rodents and birds [27, 28]. In turn, fluctuations in rodent abundance can drive variation in the abundance of Ixodes ticks because ticks feed on rodents [10, 27, 29,30,31]. Previous work in North America has shown that the density of Ixodes ticks and the risk of Lyme disease depend on the abundance of rodents in the previous year and on the abundance of acorns two years previously [10, 31,32,33]. In contrast, in Europe there is no direct evidence demonstrating that tree seed production is driving the population dynamics of Ixodes ticks [34,35,36].
In Europe, the sheep tick (Ixodes ricinus) transmits a number of important diseases including Lyme borreliosis . Ixodes ricinus is a generalist tick species that feeds on a variety of vertebrate hosts, including mammals, birds and reptiles . Ixodes ticks have three active stages: larva, nymph and adult; each stage requires a single blood meal to develop to the next stage (female ticks take a blood meal to produce eggs). Blood meals often occur in different years, and the total life-cycle can range from 3 to 5 years. Thus, Ixodes tick population ecology is characterized by cohorts of larvae, nymphs and adults that were born in different years, but that search for hosts at the same time of year. In addition, Ixodes ticks have highly seasonal activity patterns; they search for hosts during the spring, summer, and fall and enter diapause during the winter. As most tick sampling methods only capture questing ticks (i.e. ticks searching for a blood meal), within-year field studies inevitably capture the dramatic seasonal changes in the abundance of questing ticks.
The aim of this study was to better understand how climate and altitude influence the ecology of I. ricinus ticks. We used a long-term study to test whether the density of I. ricinus has increased along an altitudinal gradient in Switzerland, and whether any such increase in tick abundance was correlated with climate change. Specifically, we used data from a 15-year study that monitored the monthly abundance of I. ricinus nymphs and adult ticks at four different elevations on Chaumont Mountain, in the canton of Neuchâtel, Switzerland. We collected climate variables that are known to be important for tick ecology (e.g. temperature, relative humidity, saturation deficit and precipitation) and we obtained data on seed production by beech trees from the literature. We expected tick abundance to be higher at the lower elevations compared to the higher elevations. We predicted tick abundance to increase over time at all elevations, and that this increase would be more pronounced at the higher elevations compared to the lower elevations. Finally, we predicted that the intensity of seed production by beech trees would drive the abundance of nymphs two years later and the abundance of adult ticks three years later.
The study was conducted on the south-facing slope of Chaumont Mountain, which is part of the Jura mountains, and is located in the canton of Neuchâtel, in western Switzerland. The forest on Chaumont Mountain is mainly composed of European beech (Fagus sylvatica; 28.6%), Norway spruce (Picea abies; 28.5%), European silver fir (Abies alba; 20.4%), sycamore maple (Acer pseudoplatanus; 5.9%), European ash (Fraxinus excelsior; 3.7%), Scots pine (Pinus sylvestris; 2.3%), sessile oak (Quercus petraea; 2.3%), willows (Salix spp.; 2.1%), common whitebeam (Sorbus aria; 1.6%), and European hornbeam (Carpinus betulus; 1.0%). There is logging activity in the area, and there are also hiking trails and recreation areas for the public.
Four tick sampling sites, referred to as low, medium, high and top, were established at elevations of 620, 740, 900 and 1073 m above sea level (ASL), respectively. The four sites have been previously described [39, 40]. During the study, the canton of Neuchâtel approved the construction of a network of mountain bike trails and an adventure park, which were located at a distance of 25 m from the top site. The construction of the network of mountain bike trails occurred from 2006 to 2010, and maintenance work on the trail is done every year. The construction of the adventure park started in 2011 and each year it is heavily used by the public from April to October. We believe that the construction of these recreation facilities reduced the quality of the forest habitat at the top elevation and compromised its integrity as a long-term field site.
Questing I. ricinus nymphal and adult ticks were collected on a monthly basis over a period of 15 years from January 2004 to December 2018 at each of the four elevations. The sampling protocol has been previously described . Briefly, a 1-m2 cotton flag was dragged across low vegetation over a distance of 100 m (low elevation) or 120 m (other elevations). The flag was inspected every 20 m and nymphs and adults were counted separately. The same person (Olivier Rais) conducted all of the 3427 drags (15 years × 12 months × 5 or 6 drags × 4 elevations = 4140 drags). Drags were not performed when there was snow on the ground (713 drags).
Field-collected climate data
Temperature (T; units are °C) and relative humidity (RH; units are %) were recorded at 60 cm above ground at one moment in time on the day of tick collection (usually between 10:00 am and 2:00 pm) at each sampling site using a thermohygrometer (Model 615, Testo SA, Lonay, Switzerland). Thus, for each combination of elevation and year, we had a total of 12 field-collected measurements of temperature and relative humidity. The saturation deficit (SD) is a measure of the drying power of the atmosphere (units are mm of mercury) and is calculated using temperature (T) and relative humidity (RH) as follows: SD = (1 − RH/100) × 4.9463 × e0.0621T [42, 43]. We confirmed that our field-collected climate data were accurate by comparing them to the Climap-net data (Additional file 1: Sections 1 and 2).
In addition to the field-collected data, we obtained climate data from the Climap-net database of the Swiss meteorological center (MeteoSuisse). Two weather stations close to our four elevation sites are located in Neuchâtel at 485 m ASL and in Chaumont at 1136 m ASL. These weather stations sample at 200 cm above ground the temperature and relative humidity every hour, and the total precipitation each day. We used the daily mean temperature (average of the 24-hourly measurements), the daily mean relative humidity, and the daily total precipitation. Thus, for each year, we had a total of 365 Climap-net measurements of the daily mean temperature, the daily mean relative humidity, and the daily total precipitation. The saturation deficit was calculated as previously described. For each of the four elevation sites, we calculated site-specific climate variables by interpolating the data between the two weather stations (Additional file 1: Section 3).
Tree masting data
Mast refers to the fruit of forest trees, such as the acorns of oak trees and the beech nuts of beech trees. Data on masting (or seed production) were obtained for Neuchâtel from the MASTREE database for two important tree species: European beech (Fagus sylvatica) and Norway spruce (Picea abies) . These two species accounted for 57.1% of the trees at our study location. In this database, the mast intensity is classified into five classes: 1, 2, 3, 4, and 5, which refer to very poor mast, poor mast, moderate mast, good mast, and full mast, respectively.
Statistical analyses were restricted to the three lower elevation sites (low, medium and high). The top elevation was excluded because tick abundance at this site declined dramatically following the construction of the network of mountain bike trails and the adventure park. In Additional file 1: Sections 4 and 5, the data are analyzed with all four elevation sites together.
Annual cumulative tick density
The cumulative nymph density (CND) is a measure of the total annual abundance of questing nymphs per 100 m2 and was estimated by integrating the area under the curve (AUC) of the 12 monthly questing nymph densities for each year [43, 45]. We used this AUC approach because it is less likely to be biased by missing data compared to calculating a simple average for each year. We assume that the CND represents a small unknown fraction of the nymphs that were actually present in the area. No dragging was performed on days when there was snow on the ground (hereafter referred to as snow days). Over the study period (15 years × 12 months = 180 days), a total of 34 snow days occurred, and these were coded as missing data. In summary, tick abundance data from 2581 individual drags were collapsed into 45 estimates of CND (15 years × 3 elevations = 45 annual estimates of abundance). The same approach was used to calculate the cumulative adult tick density (CAD).
Annual mean climate variables
To investigate the relationship between climate and our annual estimates of tick abundance, we collapsed our monthly or daily weather data into a single annual value. For the field-collected data, the annual means were calculated over the 12 monthly measurements (i.e. a single measurement for each month). For the Climap-net data, the annual means were calculated over 365 daily means (i.e. a total of 365 days × 24 measurements/day = 8760 hourly measurements). Thus, the Climap-net annual means were based on 730 times more data than the field-collected annual means. However, an important advantage of the field-collected data was that they were specific for each of the four elevation sites. In contrast, the Climap-net data came from two weather stations that were located at some distance from the four elevation sites. To facilitate comparison between the slopes of the climate variables, we standardized the climate variables to z-scores (mean of 0, standard deviation of 1).
Annual tree masting variables
The model by Ostfeld et al.  predicts a 2-year time lag between masting events and the CND and a 3-year time lag between masting events and the CAD. To validate the model by Ostfeld et al. , we compared five different time lags where the CND is predicted by mast events that occurred 0, 1, 2, 3 and 4 years into the past (Additional file 1: Section 9). This analysis confirmed that the 2-year time lag is an order of magnitude better than all of the other time lags considered. The 2-year time lag has a partial r2 value of 26.3% and a P-value of 0.0002. In contrast, the partial r2 values of the other four time lags (0, 1, 3 and 4 years in the past) are much lower (range: 0.0–5.1%) and none of the P-values are statistically significant. This preliminary analysis validates our decision to model the CND and CAD as a function of the tree mast scores (of European beech and Norway spruce) two years previously (year y-2) and three years previously (year y-3), respectively. For example, tree mast scores from the year 2001 predicted the CND in year 2003 (2 years later) and the CAD in year 2004 (three years later).
Analysis of the annual variation in nymph and adult tick abundance using linear models
Count data follow a Poisson distribution or a negative binomial distribution. However, our estimates of the CND (or CAD) are summary statistics (sums or integrals) that are based on the counts of ~ 60 drags (12 dates × 5 drags per date). According to the central limit theorem of statistics, summary statistics will follow a normal distribution even if the observations on which they are based are drawn from a non-normal distribution. We therefore assumed that the residuals of our CND values follow a normal distribution. The CND values were log10-transformed to further improve their fit to the normal distribution. The log10-transformed CND values (n = 45) were analyzed using linear models (LM) with normal errors. In section 11 of Additional file 1, we show that the analysis of the CND using generalized linear models with negative binomial errors (which are appropriate for overdispersed count data) gives the same results.
The explanatory variables included the elevation site (3 levels: low, medium, high), the covariate year (rescaled as 1, 2, 3, … 15), the covariate beech mast score (range: 1–5), the covariate spruce mast score (range: 1–5), mean annual temperature, mean annual RH, mean annual SD, and mean annual precipitation (all annual climate variables were transformed to z-scores). We used mean annual climate variables based on the field-collected data and on the Climap-net data to determine which climate variables were better for explaining variation in tick abundance. As time lags are known to be important in tick ecology, we modelled the CND as a function of the mean climate variable in the present year or the previous year. The same approach was used to model the CAD (n = 45).
Model selection approach
We used a model selection approach based on the Akaike information criterion (AIC) to find the most parsimonious model. Models were ranked according to their AIC values and the Akaike weights were calculated for each model. We used the Akaike weights to calculate the model-averaged parameter estimates and their 95% confidence intervals (CIs) (Additional file 1: Section 6). For the best model from the model selection table for the log10-transformed CND (and CAD), we confirmed that the residuals met the assumptions of normality and equal variances (Additional file 1: Section 7). This result shows that we were justified in using linear models to analyze the CND and CAD.
We used R version 3.6.1 for all statistical analyses . We used the auc() function (integrates the area under the curve) in the flux package to estimate the CND and the CAD for each year. We used the lm() function in the basic package to model the log10-transformed CND as a linear model of the explanatory variables. We used the mod.sel() function and the model.av() function in the MuMIn package. The mod.sel() function creates an automated model selection table using previously generated models (linear models with normal errors in our case), and the mod.av() function calculates model-averaged parameter estimates and their standard errors. The raw data used for these statistical analyses can be found in Additional file 2: Table S1.
Differences in climate between the elevation sites
The three elevation sites differed in climate. Across the 15 years of the study, the mean annual temperature was warmest at the low elevation site (15.4 °C), intermediate at the medium elevation site (14.3 °C), and coldest at the high elevation site (13.2 °C) and these differences were significant (ANOVA: F(2, 42) = 9.431, P < 0.0001). The mean annual relative humidity was lowest at the low elevation site (61.9%), intermediate at the medium elevation site (63.7%), and highest at the high elevation site (66.5%), but these differences were not significant (ANOVA: F(2, 42) = 1.948, P = 0.155). The mean annual saturation deficit was highest at the low elevation site (6.3 mmHg), intermediate at the medium elevation site (5.5 mmHg), and lowest at the high elevation site (4.6 mmHg), and these differences were significant (ANOVA: F(2, 42) = 7.199, P = 0.002). These analyses show that the three elevation sites differed with respect to temperature, humidity and saturation deficit.
Total number of ticks collected
Over the 15 years of the study, we collected a total of 39,255 I. ricinus ticks: 31,067 nymphs and 8188 adult ticks (4168 males and 4020 females) at the three lower elevation sites.
Variation in nymph abundance among elevation sites
We used a one-way ANOVA to compare the mean CND between the three elevations (i.e. this analysis ignores the other explanatory variables). The interpretation of the CND is the theoretical number of questing nymphs that would have been collected in an area of 100 m2 over the course of the year if we had sampled the tick population each day (i.e. 365 sampling occasions per year). The mean CND was significantly different between the three elevation sites (ANOVA: F(2, 42) = 6661.4, P < 0.0001). The mean CND (and 95% CI) for the low, medium, and high elevations were as follows: 21,532 (95% CI: 16,936–27,374), 17,793 (95% CI: 13,996–22,621), and 11,535 (95% CI: 9073–14,665). The CND can be converted to a daily mean number of nymphs per 100 m2 by dividing by 365 days. Thus, the mean number of nymphs collected in an area of 100 m2 (and 95% CI) for the low, medium, and high elevations were as follows: 59.0 (95% CI: 46.4–75.0), 48.7 (95% CI: 38.3–62.0), and 31.6 (95% CI: 24.9–40.2). In summary, the mean nymph density was inversely related to the altitudinal gradient and it was highest at the low elevation site and lowest at the high elevation site.
Effect of the climate variables and mast years on variation in nymph abundance among years
For the log10-transformed CND, the best six models had a combined support of 95.0% (Table 1). The best model had 76.0% of the support (Table 1), explained 73.2% of the variation in the log10-transformed CND, and contained the explanatory variables of elevation site (partial r2 = 26.6%), year (partial r2 = 14.8%), beech mast score 2 years prior (partial r2 = 26.9%), and the field-collected relative humidity from the same year (i.e. no time lag; partial r2 = 7.6%). The five next-best models all contained the same explanatory variables of elevation site, year, and beech mast score, but differed with respect to the fourth explanatory variable, which included the site:year interaction, field-collected SD from the same year, field-collected temperature from the same year, field-collected RH from the previous year, and field-collected SD from the previous year (Table 1). The full model selection analysis with 52 models is presented in Additional file 1: Section 6.
Using the parameter estimates from the top models in Table 1, the effect sizes were calculated on the original scale of the CND with respect to the following baseline: the site was low elevation, the year was 2004, the beech tree mast index 2 years prior was set to 1, and the field-collected relative humidity from the same year was set to 50.0%. For simplicity, we present the parameter estimates from the top model in Table 1 (Additional file 1: Section 8), but we note that these are qualitatively similar to the model-averaged parameter estimates for the 52 models in the full model selection analysis (Additional file 1: Section 6).
The log10-transformed CND was significantly different between the three elevation sites (Figs. 1 and 2; Additional file 1: Section 8). The CND (on the original scale) at the low elevation was 11.4% higher than the medium elevation (Medium-Low contrast in Y-intercept = -0.053, SE = 0.045, P = 0.253). The CDN (on the original scale) at the low elevation was 43.1% higher than the high elevation (High-Low contrast in Y-intercept = -0.245, SE = 0.005, P < 0.0001). In summary, the effect of elevation site in this multiple regression approach was the same as the results from the one-way ANOVA. The CND was inversely related to the altitudinal gradient; it was highest at the low elevation and lowest at the high elevation (Figs. 1 and 2; Additional file 1: Section 8).
The slope of the covariate year was positive (slope = 0.020 per year, SE = 0.005, P < 0.0001; Additional file 1: Section 8), indicating that the log10-transformed CND was increasing over time at the three elevation sites. According to the estimate of the general slope for the three elevation sites on Chaumont Mountain, the CND (on the original scale) increased by 88.4% over the 15-year period of the study (2004–2018) (Fig. 3).
The beech mast score had a strong positive effect on the log10-transformed CND two years later (slope = 0.067 per class, SE = 0.011, P < 0.0001; Fig. 1 and Additional file 1: Section 8). According to the estimate of the general slope for the three elevation sites on Chaumont Mountain, increasing the beech mast score from 1 (poor mast) to 5 (full mast) increased the CND (on the original scale) by 86.2% (Fig. 4). Over the 15-year study period, beech trees produced 6 years of good or full mast scores (2004, 2006, 2009, 2011, 2014 and 2016), where high CND years were expected to occur two years later (2006, 2008, 2011, 2013, 2016 and 2018, respectively). At the low, medium, and high site, these 6 good mast years produced 3 (2016, 2018 and 2006), 4 (2006, 2013, 2018 and 2016), and 4 (2018, 2006, 2016 and 2011) of the six-highest CND values two years later (Fig. 1). For the adult ticks, the beech tree mast score had a strong positive effect on the log10-transformed CAD three years later (slope = 0.101, SE = 0.019, P < 0.001, partial r2 = 21.0%; Additional file 1: Section 5).
The field-collected mean annual relative humidity at 60 cm above ground had a negative effect on the log10-transformed CND in the same year (slope = -0.074 per standard deviation, SE = 0.021, P = 0.001; Additional file 1: Section 8). According to the estimate of the general slope for the three elevation sites on Chaumont Mountain, increasing the field-collected mean annual relative humidity from 50.0% to 75.0% decreased the CND from the same year (on the original scale) by 46.4% (Fig. 5).
The models with other climate variables had much less support (Table 1): field-collected mean annual SD from the same year (5%), field-collected mean annual temperature from the same year (5%), field-collected mean annual RH from the previous year (2%), and field-collected mean annual SD from the previous year (1%). We note here that field-collected mean annual RH (regardless of the time lag) always had a negative effect on nymph abundance, whereas field-collected mean annual temperature and field-collected mean annual SD always had a positive effect on tick abundance (Additional file 1: Section 6).
The analysis of the log10-transformed CND for the four elevation sites gave very similar results (Additional file 1: Section 4). The only difference was that the best model contained an interaction between elevation and year. This interaction was driven by the dramatic decline in CND over time at the top site. The best model for the four-site analysis contained the effects of elevation site, year, beech mast score, and mean annual relative humidity in the same year (i.e. no time lag), which is identical to the three-site analysis. The parameter estimates were also very similar between the 3-site and 4-site analyses. We believe that the dramatic decline in CND at the top site was caused by the construction of the recreation facilities during the study. For this reason, we presented the 3-site analysis of the CND in the main text. The analysis of the log10-transformed cumulative adult density (CAD) for all four elevation sites is also presented in Additional file 1: Section 5.
In summary, the CND increased significantly over time and with tree seed production two years earlier while it decreased significantly with elevation and with the field-collected relative humidity in the same year.
Climate change is driving the range expansion of arthropod vectors and the increased incidence of vector-borne diseases [3, 47]. In Europe and North America, there is currently much interest in determining which abiotic and biotic factors are influencing the abundance of Ixodes ticks and their associated pathogens. To address this question, we quantified the abundance of I. ricinus ticks over a period of 15 years in an area of Switzerland that is endemic for Lyme borreliosis. Variation in our annual estimates of cumulative nymph density was best explained by elevation site, year, abundance of beech nuts 2 years prior, and mean annual relative humidity in the same year. Nymph density increased by 88.4% (i.e. almost doubled) over the 15-year study period at the three elevation sites.
One of the most important results of this study is our demonstration of a strong positive association between seed production by beech trees and the density of nymphs two years later. The effect of beech seed production two years prior on the log10-transformed CND was highly statistically significant (P < 0.000001). According to our parameter estimates, the effect size was large; changing the beech mast score two years prior from 1 to 5 increased the CND by 86.2%. Beech nut abundance 2-years prior explained at least 26.9% of the variation in our annual estimates of nymph abundance. By accounting for large inter-annual fluctuations in the CND, the inclusion of the beech mast index in our models revealed the effects of more subtle ecological factors, such as mean annual relative humidity and year. In keeping with the slow life-cycle of Ixodes ticks, the beech seed production increased the abundance of adult ticks three years later (Additional file 1: Section 5). The synchronized seed production in beech trees is driven by climate variables such as temperature and precipitation [23,24,25,26]. Some recent studies have suggested that climate change is increasing the frequency of masting events [23, 25, 48]. If true, this suggests that climate change could be increasing tick abundance in endemic areas via indirect effects on tree seed production; more frequent masting events would increase the long-term abundance of reservoir hosts for ticks and tick-borne pathogens. Future long-term studies investigating whether climate change is affecting the abundance of Ixodes ticks and the incidence of tick-borne diseases should include the proper time-lagged estimates of tree seed production into their models.
The 2-year time lag between masting events and the abundance of nymphs was first discovered by Ostfeld and colleagues in the eastern USA where the blacklegged tick (Ixodes scapularis) is the vector of Lyme borreliosis [10, 31, 32]. The chain of causality is as follows. Mast seeding in year y increases the abundance of small mammals, deer, and larval ticks in year y + 1 [10, 27, 29, 31, 49, 50]. Higher densities of larval ticks coincide with and feed on higher densities of small mammals, which in turn increases the abundance of nymphs in year y + 2 [10, 31,32,33, 51,52,53]. In summary, there is good evidence in the eastern USA that mast events drive the abundance of I. scapularis nymphs with a 2-year time lag. In contrast, there are few studies in Europe that have demonstrated this phenomenon for I. ricinus. A 7-year study in Poland found the expected 2-year time lag between oak mast seeding and the incidence of Lyme borreliosis, but the causal link between tree seed production and nymphal abundance was not demonstrated . A 10-year study in the Netherlands found that the mast scores of common oak trees and beech trees in year y had a negative effect on larvae, nymphs and adults in year y + 1 . This result is expected because once ticks obtain a blood meal, they are no longer available to be sampled via dragging. An 18-year study in central Europe found a strong correlation between rodent densities in year y and tick-borne diseases in year y+1, but did not demonstrate a direct link to ticks or tree seed production . Thus, our study is the first demonstration that mast seeding events are strongly associated with the inter-annual abundance of I. ricinus ticks in Europe, which in turn influences the risk of tick-borne diseases .
An important contribution of our study is the demonstration that only the 2-year time lag was important for explaining variation in the CND (Additional file 1; Section 9). Masting with a 2-year time lag had a P-value of < 0.000001. When the time lag between masting and the CND was changed to 0, 1, 3 or 4 years, the relationship was no longer statistically significant. Masting with a 2-year time lag had a partial r2 value of 26.3%. In contrast, for the other time lags (0, 1, 3 and 4 years in the past), the partial r2 values of masting ranged from 0.0 to 5.1%. This analysis clearly shows that variation in the CND at our study location was strongly associated with masting with a 2-year time lag and not with any other time lag.
Our long-term study is a rare demonstration that the abundance of Ixodes ticks has increased in a Lyme disease-endemic area. We found that the abundance of nymphs at the three lower elevation sites increased by 88.4% from 2004 to 2018. A 10-year study in the Netherlands found an increase in the abundance of I. ricinus larvae and adults . In contrast, a 15-year study at a site in Neuchâtel, Switzerland, that is very close to the location of the present study found a significant decrease in the abundance of I. ricinus nymphs . Another 8-year study in central New Jersey, USA, found no significant directional change in the abundance of I. scapularis nymphs . Other long-term studies have found large inter-annual fluctuations in density, but none of them found that Ixodes tick density was increasing [10, 13, 14]. Increased tick abundance could have important consequences for the incidence of tick-borne diseases and public health. All else being equal, our study suggests that the risk of tick-borne disease has doubled at this location in Switzerland.
Theoretical models have predicted that Ixodes ticks will increase their distribution and abundance under global warming [6, 54,55,56]. Numerous empirical studies have shown that Ixodes ticks and tick-borne diseases are emerging public health problems at the northern limit of their geographic distribution [8, 57,58,59]. In contrast, fewer studies have investigated whether global warming has influenced Ixodes tick abundance in areas where Lyme disease is endemic [10,11,12,13,14,15,16,17]. The present study did not find direct links between climate change and the observed doubling in tick abundance. Mean annual relative humidity in the same year was the climate variable that had the greatest impact on the CND, but there was no evidence that it had undergone directional change over the 15-year study period (Additional file 1: Section 10). Conversely, the mean annual temperature increased significantly over the 15-year duration of the study (i.e. evidence of directional climate change; Additional file 1: Section 10), and while temperature had a significant and positive effect on the CND (Additional file 1: Section 6), models that contained temperature had low support (5.7%; Additional file 1: Section 6). In summary, despite finding a significant increase in temperature and ticks over the 15-year study period, the evidence for a causal link between the two is weak.
An interesting result is our demonstration that the mean annual relative humidity had a negative effect on the abundance of nymphs in the same year. This result was unexpected, as numerous studies have shown that survival of immature Ixodes ticks increases with relative humidity [18, 19, 60]. However, this result is not without precedent and other studies in Europe have found a negative relationship between relative humidity and the abundance of I. ricinus nymphs [61,62,63,64,65]. One explanation is that humid environments are favorable for the development of entomopathogenic fungi, which have been shown to cause high mortality in Ixodes ticks [66, 67]. A striking result was that the field-collected climate variables (measured at 60 cm above ground) were much more important for predicting tick abundance than climate variables from nearby weather stations (measured at 200 cm above ground). This result is even more remarkable when one considers that the annual means from the Climap-net data are based on 730 times more data (365 days × 24 hours = 8760 measurements per year) compared to the field-collected data (12 measurements per year). A recent long-term study at a nearby site in Neuchâtel  found the same result: field-collected climate variables were better than weather station climate variables at explaining inter-annual variation in tick abundance. The explanation is that the field-collected data measure the local conditions, whereas the weather stations are some distance away from the elevation sites. We suspect that microclimate factors like soil and vegetation are influencing the field-collected annual climate means to such an extent that they have only a weak resemblance to the weather station annual climate means (Additional file 1: Section 10). Our study suggests that researchers should make the effort to measure local climate variables rather than relying on data from nearby weather stations.
Ticks were found at all three elevations (620, 740 and 900 m ASL), but there were important differences in their abundance. Nymphal tick abundance was negatively correlated with elevation, which is in agreement with previous studies [39, 40]. A mechanistic explanation for this phenomenon is the relationship between temperature and tick development rates [18, 21]. At higher and colder elevations, eggs and larvae have much slower development rates, which ultimately reduces the number of larvae that reach the nymphal stage [20, 68, 69].
The nymphal tick abundance was lowest at the top site (1073 m ASL). At the top site, the CND decreased by 81.0% over the 15 years of the study, whereas at the other three elevation sites, the CND increased by 88.4% over the same time period (Additional file 1: Section 4). Thus, the population dynamics at the top site were very different from the low, medium, and high elevation sites. One explanation is the construction of recreation facilities at a distance of ~ 25 meters from our top site. These recreation facilities, which include a network of mountain bike trails (constructed from 2006 to 2010) and an adventure park with a zip line and outdoor laser games (constructed in 2011), have greatly increased the number of human visitors to the top of Chaumont Mountain. Thus, the destruction of the forest habitat and subsequent human disturbance may have caused the decline of the I. ricinus tick populations over time at our top site. An interesting alternative explanation is that the monthly tick sampling over a period of 15 years decreased the CND at the top site. Field studies typically assume that dragging removes a small fraction of the available tick population, but this assumption may not be true in habitats where tick density is already low.
We found that the abundance of I. ricinus nymphs almost doubled over 15 years at our study location in Switzerland. From a public health perspective, this increase in nymph abundance is likely to have increased the risk of tick-borne disease in this region. Beech mast years at 2 and 3 years prior were strongly and positively associated with the abundance of nymphs and adult ticks, respectively. We found that relative humidity had a negative effect on nymph abundance, which was surprising because immature Ixodes ticks are known to be sensitive to desiccation. We found no direct association between climate change and the doubling of nymphal tick abundance. However, there could be indirect effects if climate change is increasing the frequency of masting events and thereby increasing the abundance of reservoir hosts and ticks. Public health officials should be aware that mast years are an important time-lagged predictor of tick abundance and the incidence of tick-borne diseases.
Availability of data and materials
The raw data for this study are stored in the Additional file 2. The climate data are available from the Climap-net database of the Federal Office for Meteorology and Climatology (http://www.meteosuisse.admin.ch/home/service-et-publications/conseil-et-service/portail-de-donnees-dedie-aux-specialistes.html. Accessed 4 Mar 2019). The tree seed production dataset is available in the Ecology–Ecological Society of America repository (http://onlinelibrary.wiley.com/doi/10.1002/ecy.1785/suppinfo).
Akaike information criterion
above sea level
cumulative adult density
cumulative nymph density
Medlock J, Leach S. Impact of climate change on vector-borne disease in the UK. Lancet Infect Dis. 2015;15:159–99.
Gage KL, Burkot TR, Eisen RJ, Hayes EB. Climate and vectorborne diseases. Am J Prev Med. 2008;35:436–50.
Mills JN, Gage KL, Khan AS. Potential influence of climate change on vector-borne and zoonotic diseases: a review and proposed research plan. Environ Health Perspect. 2010;118:1507–14.
Lindgren E, Tälleklint L, Polfeldt T. Impact of climatic change on the northern latitude limit and population density of the disease-transmitting European tick Ixodes ricinus. Environ Health Perspect. 2000;108:119–23.
Tälleklint L, Jaenson TGT. Increasing geographical distribution and density of Ixodes ricinus (Acari: Ixodidae) in central and northern Sweden. J Med Entomol. 1998;35:521–6.
Ogden NH, Maarouf A, Barker IK, Bigras-Poulin M, Lindsay LR, Morshed MG, et al. Climate change and the potential for range expansion of the Lyme disease vector Ixodes scapularis in Canada. Int J Parasitol. 2006;36:63–70.
Brownstein JS, Holford TR, Fish D. Effect of climate change on Lyme disease risk in North America. EcoHealth. 2005;2:38–46.
Jaenson TGT, Jaenson DG, Eisen L, Petersson E, Lindgren E. Changes in the geographical distribution and abundance of the tick Ixodes ricinus during the past 30 years in Sweden. Parasit Vectors. 2012;5:8.
Tokarevich NK, Tronin AA, Blinova OV, Buzinov RV, Boltenkov VP, Yurasova ED, et al. The impact of climate change on the expansion of Ixodes persulcatus habitat and the incidence of tick-borne encephalitis in the north of European Russia. Global Health Action. 2011;4:8448.
Ostfeld RS, Canham CD, Oggenfuss K, Winchcombe RJ, Keesing F. Climate, deer, rodents, and acorns as determinants of variation in Lyme-disease risk. PLoS Biol. 2006;4:e145.
Hauser G, Rais O, Morán Cadenas F, Gonseth Y, Bouzelboudjen M, Gern L. Influence of climatic factors on Ixodes ricinus nymph abundance and phenology over a long-term monthly observation in Switzerland (2000–2014). Parasit Vectors. 2018;11:289.
Schulze TL, Jordan RA, Schulze CJ, Hung RW. Precipitation and temperature as predictors of the local abundance of Ixodes scapularis (Acari: Ixodidae) nymphs. J Med Entomol. 2009;46:1025–9.
Brugger K, Walter M, Chitimia-Dobler L, Dobler G, Rubel F. Forecasting next season’s Ixodes ricinus nymphal density: the example of southern Germany 2018. Exp Appl Acarol. 2018;75:281–8.
Daniel M, Malý M, Danielová V, Kříž B, Nuttall P. Abiotic predictors and annual seasonal dynamics of Ixodes ricinus, the major disease vector of central Europe. Parasit Vectors. 2015;8:478.
Berger KA, Ginsberg HS, Dugas KD, Hamel LH, Mather TN. Adverse moisture events predict seasonal abundance of Lyme disease vector ticks (Ixodes scapularis). Parasit Vectors. 2014;7:181.
Hayes LE, Scott JA, Stafford KC. Influences of weather on Ixodes scapularis nymphal densities at long-term study sites in Connecticut. Ticks Tick Borne Dis. 2015;6:258–66.
Burtis JC, Sullivan P, Levi T, Oggenfuss K, Fahey TJ, Ostfeld RS. The impact of temperature and precipitation on blacklegged tick activity and Lyme disease incidence in endemic and emerging regions. Parasit Vectors. 2016;9:606.
Kilpatrick A, Dobson A, Levi T, Salkeld D, Swei A, Ginsberg H, et al. Lyme disease ecology in a changing world: consensus, uncertainty and critical gaps for improving control. Philos Trans R Soc Lond B. 2017;372:20160117.
Gray JS. Review The ecology of ticks transmitting Lyme borreliosis. Exp Appl Acarol. 1998;22:249–58.
Randolph SE. Tick ecology: processes and patterns behind the epidemiological risk posed by ixodid ticks as vectors. Parasitology. 2004;129:37–65.
Ogden NH, Lindsay LR, Beauchamp G, Charron D, Maarouf A, OʼCallaghan CJ, et al. Investigation of relationships between temperature and developmental rates of tick Ixodes scapularis (Acari: Ixodidae) in the laboratory and field. J Med Entomol. 2004;41:622–33.
Gray JS, Dautel H, Estrada-Peña A, Kahl O, Lindgren E. Effects of climate change on ticks and tick-borne diseases in Europe. Interdiscip Perspect Infect Dis. 2009;2009:593232.
Drobyshev I, Niklasson M, Mazerolle MJ, Bergeron Y. Reconstruction of a 253-year long mast record of European beech reveals its association with large scale temperature variability and no long-term trend in mast frequencies. Agric For Meteorol. 2014;192:9–17.
Drobyshev I, Övergaard R, Saygin I, Niklasson M, Hickler T, Karlsson M, et al. Masting behaviour and dendrochronology of European beech (Fagus sylvatica L.) in southern Sweden. For Ecol Manage. 2010;259:2160–71.
Övergaard R, Gemmel P, Karlsson M. Effects of weather conditions on mast year frequency in beech (Fagus sylvatica L.) in Sweden. Forestry. 2007;80:555–65.
Piovesan G, Adams JM. Masting behaviour in beech: linking reproduction and climatic variation. Can J Bot. 2001;79:1039–47.
Clotfelter E, Pedersen A, Cranford J, Ram N, Snajdr E, Nolan V, et al. Acorn mast drives long-term dynamics of rodent and songbird populations. Oecologia. 2008;154:493–503.
Schnurr JL, Ostfeld RS, Canham CD. Direct and indirect effects of masting on rodent populations and tree seed survival. Oikos. 2002;96:402–10.
Ostfeld RS, Jones CG, Wolff JO. Of mice and mast. Bioscience. 1996;46:323–30.
Clement J, Vercauteren J, Verstraeten WW, Ducoffre G, Barrios JM, Vandamme AM, et al. Relating increasing hantavirus incidences to the changing climate: the mast connection. Int J Health Geogr. 2009;8:1.
Ostfeld RS, Schauber EM, Canham CD, Keesing F, Jones CG, Wolff JO. Effects of acorn production and mouse abundance on abundance and Borrelia burgdorferi infection prevalence of nymphal Ixodes scapularis ticks. Vector Borne Zoonotic Dis. 2001;1:55–63.
Ostfeld RS, Levi T, Keesing F, Oggenfuss K, Canham CD. Tick-borne disease risk in a forest food web. Ecology. 2018;99:1562–73.
Schauber EM, Ostfeld RS, Evans J, Andrew S. What is the best predictor of annual Lyme disease incidence: Weather, mice, or acorns? Ecol Appl. 2005;15:575–86.
Hartemink N, van Vliet A, Sprong H, Jacobs F, Garcia-Marti I, Zurita-Milla R, et al. Temporal-spatial variation in questing tick activity in the Netherlands: The effect of climatic and habitat factors. Vector Borne Zoonotic Dis. 2019;19:494–505.
Bogdziewicz M, Szymkowiak J. Oak acorn crop and Google search volume predict Lyme disease risk in temperate Europe. Basic Appl Ecol. 2016;17:300–7.
Tkadlec E, Václavík T, Široký P. Rodent host abundance and climate variability as predictors of tickborne disease risk 1 year in advance. Emerging Infect Dis. 2019;25:1738.
Jongejan F, Uilenberg G. The global importance of ticks. Parasitology. 2005;129:3–14.
Gern L, Humair PF. Ecology of Borrelia burgdorferi sensu lato in Europe. In: Gray JS, Kahl O, Lane RS, Stanek G, editors. Lyme Borreliosis: biology, epidemiology and control. Wallinford: CABI International; 2002. p. 149–74.
Jouda F, Perret JL, Gern L. Ixodes ricinus density, and distribution and prevalence of Borrelia burgdorferi sensu lato infection along an altitudinal gradient. J Med Entomol. 2004;41:162–9.
Morán Cadenas F, Rais O, Jouda F, Douet V, Humair PF, Moret J, et al. Phenology of Ixodes ricinus and infection with Borrelia burgdorferi sensu lato along a North- and South-facing altitudinal gradient on Chaumont Mountain, Switzerland. J Med Entomol. 2007;44:683–93.
Morán Cadenas F, Rais O, Humair PF, Douet V, Moret J, Gern L. Identification of host bloodmeal source and Borrelia burgdorferi sensu lato in field-collected Ixodes ricinus ticks in Chaumont (Switzerland). J Med Entomol. 2007;44:1109–17.
Randolph SE, Storey K. Impact of microclimate on immature tick-rodent host interactions (Acari: Ixodidae): implications for parasite transmission. J Med Entomol. 1999;36:741–8.
Perret JL, Guigoz E, Rais O, Gern L. Influence of saturation deficit and temperature on Ixodes ricinus tick questing activity in a Lyme borreliosis-endemic area (Switzerland). Parasitol Res. 2000;86:554–7.
Ascoli D, Maringer J, Hacket-Pain A, Conedera M, Drobyshev I, Motta R, et al. Two centuries of masting data for European beech and Norway spruce across the European continent. Ecology. 2017;98:1473.
Eisen RJ, Eisen L, Castro MB, Lane RS. Environmentally related variability in risk of exposure to Lyme disease spirochetes in northern California: effect of climatic conditions and habitat type. Environ Entomol. 2003;32:1010–8.
R Development Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2016. https://www.R-project.org/
Caminade C, McIntyre KM, Jones AE. Impact of recent and future climate change on vector-borne diseases. Ann N Y Acad Sci. 2019;1436:157–73.
Schmidt W. Temporal variation in beech masting (Fagus sylvatica L.) in a limestone beech forest (1981–2004). Allgemeine Forst und Jagdzeitung. 2006;177:9–19.
Wolff JO. Population fluctuations of mast-eating rodents are correlated with production of acorns. J Mammal. 1996;77:850–6.
Jones CG, Ostfeld RS, Richard MP, Schauber EM, Wolff JO. Chain reactions linking acorns to gypsy moth outbreaks and Lyme disease risk. Science. 1998;279:1023–6.
Gray JS, Kahl O, Janetzki C, Stein J. Studies on the ecology of Lyme disease in a deer forest in County Galway. Ireland. J Med Entomol. 1992;29:915–20.
Gilbert L, Maffey GL, Ramsay SL, Hester AJ. The effect of deer management on the abundance of Ixodes ricinus in Scotland. Ecol Appl. 2012;22:658–67.
Perez G, Bastian S, Agoulon A, Bouju A, Durand A, Faille F, et al. Effect of landscape features on the relationship between Ixodes ricinus ticks and their small mammal hosts. Parasit Vectors. 2016;9:20.
Ogden NH, St-Onge L, Barker IK, Brazeau S, Bigras-Poulin M, Charron DF, et al. Risk maps for range expansion of the Lyme disease vector, Ixodes scapularis, in Canada now and with climate change. Int J Health Geogr. 2008;7:24.
Li S, Gilbert L, Harrison PA, Rounsevell MDA. Modelling the seasonality of Lyme disease risk and the potential impacts of a warming climate within the heterogeneous landscapes of Scotland. J R Soc Interface. 2016;13:20160140.
Alkishe AA, Peterson AT, Samy AM. Climate change influences on the potential geographic distribution of the disease vector tick Ixodes ricinus. PLoS ONE. 2017;12:e0189092.
Ogden NH, Lindsay LR, Morshed M, Sockett PN, Artsob H. The emergence of Lyme disease in Canada. CMAJ. 2009;180:1221–4.
Leighton PA, Koffi JK, Pelcat Y, Lindsay LR, Ogden NH. Predicting the speed of tick invasion: an empirical model of range expansion for the Lyme disease vector Ixodes scapularis in Canada. J Appl Ecol. 2012;49:457–64.
Korotkov Y, Kozlova T, Kozlovskaya L. Observations on changes in abundance of questing Ixodes ricinus, castor bean tick, over a 35-year period in the eastern part of its range (Russia, Tula region). Med Vet Entomol. 2015;29:129–36.
Rodgers SE, Zolnik CP, Mather TN. Duration of exposure to suboptimal atmospheric moisture affects nymphal blacklegged tick survival. J Med Entomol. 2007;44:372–5.
Hubálek Z, Halouzka J, Juricova Z. Host-seeking activity of ixodid ticks in relation to weather variables. J Vector Ecol. 2003;28:159–65.
Schwarz A, Maier WA, Kistemann T, Kampen H. Analysis of the distribution of the tick Ixodes ricinus L. (Acari: Ixodidae) in a nature reserve of western Germany using Geographic Information Systems. Int J Hyg Environ Health. 2009;212:87–96.
Li S, Heyman P, Cochez C, Simons L, Vanwambeke SO. A multi-level analysis of the relationship between environmental factors and questing Ixodes ricinus dynamics in Belgium. Parasit Vectors. 2012;5:149.
James M, Bowman A, Forbes K, Lewis F, McLeod J, Gilbert L. Environmental determinants of Ixodes ricinus ticks and the incidence of Borrelia burgdorferi sensu lato, the agent of Lyme borreliosis, in Scotland. Parasitology. 2012;140:1–10.
Kiewra D, Kryza M, Szymanowski M. Influence of selected meteorological variables on the questing activity of Ixodes ricinus ticks in Lower Silesia, SW Poland. J Vector Ecol. 2014;39:138–45.
Benjamin MA, Zhioua E, Ostfeld RS. Laboratory and field evaluation of the entomopathogenic fungus Metarhizium anisopliae (Deuteromycetes) for controlling questing adult Ixodes scapularis (Acari: Ixodidae). J Med Entomol. 2002;39:723–8.
Hartelt K, Wurst E, Collatz J, Zimmermann G, Kleespies RG, Oehme RM, et al. Biological control of the tick Ixodes ricinus with entomopathogenic fungi and nematodes: Preliminary results from laboratory experiments. Int J Med Microbiol. 2008;298:314–20.
Ogden NH, Bigras-Poulin M, O’Callaghan CJ, Barker IK, Lindsay LR, Maarouf A, et al. A dynamic population model to investigate effects of climate on geographic range and seasonality of the tick Ixodes scapularis. Int J Parasitol. 2005;35:375–89.
Eisen RJ, Eisen L, Ogden NH, Beard CB. Linkages of weather and climate with Ixodes scapularis and Ixodes pacificus (Acari: Ixodidae), enzootic transmission of Borrelia burgdorferi, and Lyme disease in North America. J Med Entomol. 2016;53:250–61.
We would like to thank Lise Gern for her financial support and for generously giving us these data. We are grateful to Jan Boni, forest engineer of the city and district of Neuchâtel, for providing useful information on the construction and maintenance of the recreational areas on Chaumont Mountain. This manuscript was improved by comments from Dolores Genné, Christopher Zinck, Jessica Thoroughgood, and Matthew Yunik. This study was part of the Ph.D. thesis of Cindy Bregnard.
This study was supported by grants obtained by Lise Gern from the Swiss National Science Foundation: FN 32-57098.99, FN 3200B0-100657, FN 320000-113936 and FN 310030-127064 and by grants obtained by Lise Gern from the Federal Office of Public Health National Reference Center: 2009/10 (Projekt (911) 316) and 2011/13 (11.006911/ 304.0001-707). The doctoral salary of CB was supported by the University of Neuchâtel. This research was also supported by two grants awarded to MJV: a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada (RGPIN-2019-04483) and an Establishment Grant from the Saskatchewan Health Research Foundation (4583).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Validation of the field-collected temperature data. Section 2. Validation of the field-collected relative humidity data. Section 3. Interpolation of the Climap-net climate data. Section 4. Full statistical analysis of the cumulative nymph abundance (CND) for the four elevation sites. Section 5. Full statistical analysis of the cumulative adult abundance (CAD) for the four elevation sites. Section 6. Full model selection table, support of each individual explanatory variable, and model-averaged parameter estimates of the nymph abundance for the three lowest elevation sites. Section 7. Assumptions of the linear models for the best models from the AIC-based model selection approach of the nymph and adult abundance. Section 8. Parameter estimates of the top model in the model selection table of the nymph abundance for the three lowest elevation sites. Section 9. Analysis of different time lags between the cumulative nymph abundance and beech masting. Section 10. Climate change over the 15-year study period. Section 11. Analysis of CND using Generalized Linear Models with negative binomial errors.
Raw data used for all statistical analyses.
About this article
Cite this article
Bregnard, C., Rais, O. & Voordouw, M.J. Climate and tree seed production predict the abundance of the European Lyme disease vector over a 15-year period. Parasites Vectors 13, 408 (2020). https://doi.org/10.1186/s13071-020-04291-z