Spatial and seasonal determinants of Lyme borreliosis incidence in France, 2016 to 2021

Background Lyme borreliosis (LB) is the most widespread hard tick-borne zoonosis in the northern hemisphere. Existing studies in Europe have focused mainly on acarological risk assessment, with few investigations exploring human LB occurrence. Aim We explored the determinants of spatial and seasonal LB variations in France from 2016 to 2021 by integrating environmental, animal, meteorological and anthropogenic factors, and then mapped seasonal LB risk predictions. Methods We fitted 2016–19 LB national surveillance data to a two-part spatio-temporal statistical model. Spatial and temporal random effects were specified using a Besag-York-Mollie model and a seasonal model, respectively. Coefficients were estimated in a Bayesian framework using integrated nested Laplace approximation. Data from 2020–21 were used for model validation. Results A high vegetation index (≥ 0.6) was positively associated with seasonal LB presence, while the index of deer presence (> 60%), mild soil temperature (15–22 °C), moderate air saturation deficit (1.5–5 mmHg) and higher tick bite frequency were associated with increased incidence. Prediction maps show a higher risk of LB in spring and summer (April–September), with higher incidence in parts of eastern, midwestern and south-western France. Conclusion We present a national level spatial assessment of seasonal LB occurrence in Europe, disentangling factors associated with the presence and increased incidence of LB. Our findings yield quantitative evidence for national public health agencies to plan targeted prevention campaigns to reduce LB burden, enhance surveillance and identify further data needs. This approach can be tested in other LB endemic areas.


Table of contents
Supplementary Material S1 Calculation of the frequency of tick bite reports per department and per quarter ..

Supplementary Material S1 Calculation of the frequency of tick bite reports per department and per quarter
A total of 43,915 human tick bite reports were included in the analysis which was collected during 2017-2021 by the national CiTIQUE project (www.citique.fr) [1,2]. For each report, we extracted the information on reporting dates and GPS position (WGS84). Since the tick species in each report was not available, we used all tick bites reports. However, 96% of a random sample of 2009 of those human-attached ticks (at least 150 randomly selected from each region of France) were identified morphologically and genetically as Ixodes ricinus. This gives weight to our assumption that Ixodes ricinus was the main tick species responsible for the large majority of bites in the dataset.
We calculated the departmental and quarterly proportion of tick bite reports. For the period 2016-2019, we produced averaged quarterly values using the data available (i.e., from July 2017), and used these for each year of the period. For the year 2020 and 2021, we used corresponding 2020 and 2021 CiTiQUE data. Considering that human recreational activities are closely related to the probability of tick bite exposure, we assumed that the population at risk is proportional to the fraction of leisure-related green space in the department.
where =1, 2, …, for each department, and =1, 2, …, for each quarter. Here = 4, we used the available data during 2017-19 to calculate the quarterly average for four quarters (in order January to March (winter), April to June (spring), July to September (summer), October to December (autumn)). is the number of inhabitants in department , Gd represents the cumulative percentages of recreational areas in that department where people are likely to encounter tick bites. They are derived from the Corine land cover dataset ( (1), we adjusted the population weight of department to obtain the number of quarterly reported tick bites per 100,000 inhabitants, denoted as ℎ ( , ). In equation (2) the population-adjusted tick bite reports ℎ was aggregated for all department and all quarters. We then calculated the quarterly proportion of tick bite reports ( , ) for each department. To be consistent with the size of the grid cells analysed, we rasterized the department polygons and assigned the same value ( , ) to grid cells belonging to the same department. All operations were performed using the raster library in R version 4.0.5 [4].
Population data were obtained from the 2017 national census by the National Institute of Statistics and Economic Studies (INSEE) [5]. French administrative boundaries shapefiles were downloaded from the Global Administrative Areas (GADM) Data [6]. Geographical information was retrieved from the National Institute of Geography and Forestry (IGN) [7].

Supplementary Material S2: Spatial interpolation by ordinary kriging method
To apply the ordinary kriging method, we assumed that the difference of LB incidences between two locations depends only on their distance [8]. We first estimated the semivariance function, which is defined by the following equation: where Z (xi, yi) represent the observed LB incidence value at coordinates (xi, yi) and N(h) is the number of pairs of observations at a distance h from each other. The computation of ̂(ℎ) is repeated sequentially for 2h, 3h, ... , kh, By default, k is set up to 10 and all operations were performed using the kriging library in R version 4.0.5. [4] We then fitted the observed semivariogram with a spherical model and estimated the parameters of the theoretical model.
0 is the value of nugget parameter, representing the variability observed at a smaller scale than the minimum sampling interval or the measure errors, while the sill parameter is given by 0 + 1 that is, the maximum semivariance (ℎ). 2 suggests that beyond this distance, the observations are independent of each other and are no longer spatially correlated [8]. We plotted experimental and fitted models using the average LB incidence from 2016 to 2019 and observed spatial autocorrelation in the distribution of LB within 110 km ( Figure S1).      Supplementary Material S6: Cross-validation of the gamma model Figure S8. Histogram of the probability integral transform (PIT), 2020. Figure S9. Histogram of the probability integral transform (PIT), 2021. Figure S10. The overall national predicted and observed annual Lyme borreliosis incidence rate per 100,000 inhabitants, France, 2016-2021.
The predicted mean incidence values and their confidence interval (95% CI) were calculated from all grid cells and quarters of the gamma model and for each year (in green). The observed mean incidence values and their 95% CI were those calculated by the Réseau Sentinelles (in blue) [9].