Meteorological and climatic variables predict the phenology of Ixodes ricinus nymph activity in France, accounting for habitat heterogeneity
Sampling sitesLongitudinal observation campaigns for I. ricinus nymph activity were carried out at 11 sampling sites in forest areas from seven different tick observatories across France. Tick observatories are located at the following French municipal areas, where the coordinates of the centre of each municipal area and the climatic types29 are also provided as: (1) La Tour de Salvagny (45° 48′ 50.6″ N 4° 42′ 53.2″ E; Mixed climates); (2) Saint-Genès-Champanelle (45° 43′ 23.8″ N 3°01′ 08.0″ E; Mountain climate); (3) Etiolles (48° 37′ 59.9″ N 2° 28′ 00.1″ E; Degraded oceanic climate); (4) Carquefou (47° 17′ 58.5″ N 1° 29′ 26.0″ W; Oceanic climate); (5) Gardouch (43°23′ 25.7″ N 1° 41′ 02.1″ E; South-West Basin climate); (6) Velaine-en-Haye (48° 42′ 13.4″ N 6° 01′ 16.1″ E; Semi-continental climate); (7) Les Bordes (47° 48′ 47.3″ N 2° 24′ 01.3″ E; Degraded-oceanic climate) (Fig. 1). The observation campaigns were carried out from April/June 2014 to May/June 2021 in most observatories, except for Les Bordes, which began in April 2018.Figure 1The map was created using QGIS version 3.8, Zanzibar (https://www.qgis.org). The climatic region types were previously classified by Joly et al.29.The distribution of tick observatories according to the climatic region types of continental France: (1) Etiolles (degraded oceanic); (2) Velaine-en-Haye (semi-continental); (3) Les Bordes (degraded oceanic); (4) Carquefou (oceanic); (5) La Tour de Salvagny (mixed); (6) Saint-Genès-Champanelle (mountain); (7) Gardouch (south-west basin). Phenological patterns observed at each observatory were also indicated.Full size imageEach tick observatory corresponds to one sampling site except La Tour de Salvagny, Gardouch, and Les Bordes (Table S1). In La Tour de Salvagny, we had to withdraw the observations at the original site (La Tour de Salvagny A) in September 2016 because the site became no longer accessible. In April 2017, we continued our observations at a nearby site, approximately 2 km apart (La Tour de Salvagny B). In Gardouch, the activity of questing nymphs was observed both inside and outside the enclosed area of an experimental station on roe deer (Capreolus capreolus), referred to as Gardouch Inside and Gardouch Outside, respectively. The estimated population density of roe deer in Gardouch Inside (50 individuals per 100 ha) was higher than Gardouch Outside (less than 20 individuals per 100 ha) (H. Verheiden, personal communication, 15th October 2021). Furthermore, three sampling sites in Les Bordes, approximately 1.2 km apart, were referred to as Les Bordes A, B, and C, respectively. Additional sampling sites of these observatories were considered and reported as distinct sampling sites in further analyses, resulting in a total to 11 sampling sites from 7 observatories. Furthermore, due to their geographical proximity, meteorological/climatic factors of different sampling sites from the same observatories were considered identical in subsequent statistical analyses, whereas land cover and topography factors could be varied.Field observation campaigns were planned and carried out by local investigators who had been trained on the sampling protocol. The locations of forests, sampling sites, and passages were chosen where their biotopes are known to be suitable for I. ricinus tick populations around each observatory at the time the field observation campaigns started30. The observations were never carried out during the daytime when the weather was highly unfavourable to questing ticks, e.g., heavy rain, snow, or snow cover.Sampling protocol for questing Ixodes ricinus nymphsActivity of questing I. ricinus nymphs was observed by a cloth-dragging sampling technique31. Within a 1-km radius, a 1 m × 1 m white cloth was dragged over 10 observation units of 10 m short-grass vegetative forest floors, called transects. For each transect, a repeated removal sampling design was used27. The cloth-dragging sampling process was successively repeated three times per sampling. All nymphs found on white cloth in each campaign were removed and collected in a vial for subsequent morphological identification32 by the same acarologists at the corresponding laboratories. As a result, the questing nymph activity of each sampling site was monitored as a total number of confirmed I. ricinus nymphs collected from three repeated sampling on 10 transects, equivalent to a surface area of 100 m2. This measure was considered as an indicator for tick abundance on the day of sampling. The same transects were repeatedly sampled throughout the study period at approximately 1-month intervals.Environmental dataWe tested 28 environmental variables to explain the observed I. ricinus nymph activity (Table 1). These variables could be categorized as: (1) Daytime duration and meteorological variables (time-dependent, 9 variables); (2) Land cover, topography, and bioclimatic variables (time-independent, 19 variables).Table 1 Environmental variables (meteorological, land cover, topography, and bioclimatic variables) used to explain I. ricinus nymph counts per 100 m2 in regression analysis.Full size tableDaytime duration and meteorological variablesDaytime duration ((daytime)) from January 2013 to June 2021 at each sampling site was obtained from the corresponding latitude using geosphere package33. Hourly meteorological data (2-m temperature and relative humidity) were recorded locally at each forest. Subsequently, daily mean, minimum, and maximum values of temperature (({T}_{M}), ({T}_{N}), and ({T}_{X}); in °C) and relative humidity (({U}_{M}), ({U}_{N}), and ({U}_{X}); in %) were derived from these hourly records. The meteorological seasons of the temperate area in northern hemisphere are defined as: (1) Spring, 1st March to 31st May; (2) Summer, 1st June to 31st August; (3) Autumn, 1st September to 30th November; (4) Winter, 1st December to 28th or 29th February.Missing values found on these local daily-level variables were imputed by the random forest algorithm in mice package34. External daily meteorological data, i.e., daily average temperature and relative humidity, derived from neighbouring weather stations (Météo-France or INRAE), as well as month and year information, were used as auxiliary variables (Table S2). As a result, the imputation process creates a total of 500 iterated values for each variable. The median values of 500 imputations were used to replace the missing values.The imputed daily meteorological data were subsequently used to calculate the averaged values in different lagged time intervals for further analysis, called interval-average variables15. The interval-average variables were generated to reduce the uncertainty that might arise during the imputation process and to capture the cumulative effects of the meteorological variables, which were mean temperature ({T}_{M}) and minimum relative humidity ({U}_{N}). The interval-average variables were defined as the average values of a meteorological variable (Min) {({T}_{M}), ({U}_{N})} during a period between ({t}_{1}) to ({t}_{2}) month(s) before the sampling, denoted as ({M}^{{t}_{1}:{t}_{2}}), where 1 month consists of 28 days. As temperature conditions affect several ecological processes of tick populations, particularly developmental and questing rates3, the mean temperature ({T}_{M}) was selected for further analysis to reflect the overall temperature effects. While the minimum relative humidity ({U}_{N}) was chosen for the following reasons: (1) the survival of I. ricinus is highly sensitive to desiccation conditions6,7,8. As a result, when compared to mean or maximum relative humidity, minimum relative humidity is a relatively strong indicator of the effects of desiccation stress; (2) the variation of minimum relative humidity among all sites was higher than that of the mean and maximum relative humidity. This high variation allowed us to better describe meteorological characteristics of each sampling site.Here, we hypothesized that interval-average meteorological conditions influence the dynamics of observed nymph activity at different time lags in different manners. Short-term lags may have an impact on immediate responses, such as the probability of questing. At the same time, long-term lags may influence the dynamics of nymph abundance, which is associated with development and survival rates. Therefore, we explored the impact of each meteorological variable at following time lags on the observed nymphs activity in subsequent regression analysis: (1) 1-month moving average condition, ({M}^{0:1}); (2) previous 3-to-6-month moving average condition, ({M}^{3:6}); (3) 6-month moving average condition, ({M}^{0:6}); (4) 12-month moving average condition, ({M}^{0:12}). For instance, ({T}_{M}^{0:1}) denotes 1-month moving average temperature, representing an average of temperature between 0 and 1 months (0–28 days) before the day of sampling.In addition to the interval-average variables, monthly and seasonal average values of mean temperature and minimum relative humidity during the observation period were also calculated to describe the characteristics of meteorological conditions of each sampling site.Land cover, topography, and bioclimatic variablesWe obtained land cover, topography, and bioclimatic data from a 1-km radius buffer area around the center of each sampling site to capture habitat characteristics across all 10 transects. All the variables were handled and obtained by using QGIS version 3.8.035. The digital elevation model (DEM) data derived from the Shuttle Radar Topography Mission (SRTM) database36 was used to describe the topographic features of sampling sites, which included the mean (({mean}_{elv})) and standard deviation (({sd}_{elv})) of the elevation (in m above sea level), the proportion of flat area (({p}_{flat}); defined by the slope ≤ 2.5%37), the proportion of area facing north (({p}_{north})), east (({p}_{east})), west (({p}_{west})), and south (({p}_{south})), and the catchment area ((catchment)) as a proxy variable for moisture. Bioclimatic variables for each site (historical average conditions during 1970–2000) were derived from the WorldClim database38, including the annual mean temperature (({BIO1}_{Temp}); in °C), the mean diurnal range (({BIO2}_{Diur}); in °C), the maximum temperature of the warmest month (({BIO5}_{maxTemp}); in °C), and the annual precipitation (({BIO12}_{Prec}); in mm). The land cover features of each sampling site were described using the CORINE Land Cover (CLC) 201839, while the characteristics of forests were explained by the BD forêt version 2 data40. The forest fragmentation was characterized by the percentage of forest-covering area (({p}_{Forest})), the forest edge density (({ED}_{Forest}); in m/km2), and the number of forest patches (({n}_{Forest})). While the diversities of the land cover types (level-1 and level-2 CLC) and the forest types were calculated by using the Shannon’s diversity index41 ((H)) as (H=sum_{i=1}^{S}{p}_{i}mathrm{ln}{p}_{i}), where (S) is the total number of land cover/forest types and ({p}_{i}) is the proportion of land cover/forest type (i) within the 1-km radius buffer area. The Shannon’s diversity index for level-1 CLC, level-2 CLC, and forest types were denoted as ({H}_{CLC1}), ({H}_{CLC2}), and ({H}_{Forest}), respectively. Finally, the soil pH data (({pH}_{soil})) was retrieved from the European Soil Data Centre (ESDC) database42.Statistical analysisAll the statistical analyses were carried out using the programming language R version 3.6.043. The variations of questing nymph population of each site were described by using (1) baseline annual nymph counts (spatial variation); (2) phenological patterns (seasonal variation). A baseline annual nymph count of site (i) (({{N}_{base}}_{i})) was defined as a summation of monthly median nymph counts ({varvec{tilde{N}}}_{i}={{tilde{N }}_{i,t}}) across all 12 months (tin left{mathrm{1,2},dots ,12right}) and expressed as: ({{N}_{base}}_{i}=sum_{t=1}^{12}{tilde{N }}_{i,t}). Subsequently, the monthly median nymph counts of each site ({varvec{tilde{N}}}_{i}) were transformed into normalized monthly median nymph counts ({varvec{tilde{N}}}_{i}^{*}={{tilde{N }}_{i,t}^{*}}) following Eq. (1) to have a range value of 0 to 1, which allows us to compare phenological patterns among all sites that have different annual baseline nymph counts.$${tilde{N }}_{i,t}^{*}=frac{{tilde{N }}_{i,t}}{mathrm{max}({stackrel{sim }{{varvec{N}}}}_{i})}$$
(1)
The term (mathrm{max}({stackrel{sim }{{varvec{N}}}}_{i})) denoted the maximum monthly median nymph counts. The normalized median nymph count ({tilde{N }}_{i,t}^{*}) of 1 indicates the maximum nymph activity (peak), while the value ({tilde{N }}_{i,t}^{*}) of 0 designates the absence of nymph activity. Afterwards, the phenological patterns were descriptively classified using the following criteria: (1) the season which the peaks of activity arrive; (2) evidence of reduced activity during winter (November–January); (3) the number of activity waves in a year, whether the pattern is unimodal or bimodal. After assigning phenological patterns to each site, the overall trends of different patterns were derived from medians of the normalized monthly median nymph count ({tilde{N }}_{i,t}^{*}) from all sites that belonged to each pattern. Furthermore, the directional changes in the maximum nymph counts were tested using a Spearman’s rank correlation coefficient, a p-value More
