# Meteorological and climatic variables predict the phenology of Ixodes ricinus nymphs activity in France, accounting for habitat heterogeneity

Phrutsamon Wongnak, Severine Bord, Maude Jacquot, Albert Agoulon, Frédéric Beugnet, Laure Bournez, Nicolas Cèbe, Adélie Chevalier, Jean-François Cosson, Naima N. Dambrine, et al.

► To cite this version:

Phrutsamon Wongnak, Severine Bord, Maude Jacquot, Albert Agoulon, Frédéric Beugnet, et al. Meteorological and climatic variables predict the phenology of Ixodes ricinus nymphs activity in France, accounting for habitat heterogeneity. Scientific Reports, 2022, 12 (1), 10.1038/s41598-022-11479-z. hal-03510954v3

HAL Id: hal-03510954

https://hal.inrae.fr/hal-03510954v3

Submitted on 13 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

# Meteorological and climatic variables predict the phenology of *Ixodes ricinus* nymph activity in France, accounting for habitat heterogeneity

Phrutsamon Wongnak<fnref n="1" />,², Séverine Bord³, Maude Jacquot<fnref n="1" />,²,¹³, Albert Agoulon⁴,
Frédéric Beugnet⁵, Laure Bournez⁶, Nicolas Cèbe⁷,⁸, Adélie Chevalier⁹,
Jean-François Cosson¹⁰,¹⁴, Naïma Dambrine<fnref n="1" />,², Thierry Hoch⁴, Frédéric Huard¹¹,
Nathalie Korboulewsky⁹, Isabelle Lebert<fnref n="1" />,², Aurélien Madouasse⁴, Anders Mårell⁹,
Sara Moutailler¹⁰, Olivier Plantard⁴, Thomas Pollet¹⁰,¹², Valérie Poux<fnref n="1" />,²,
Magalie René-Martellet<fnref n="1" />,², Muriel Vayssier-Taussat¹⁰, Hélène Verheyden⁷,⁸,
Gwenael Vourc'h<fnref n="1" />,² & Karine Chalvet-Monfray<fnref n="1" />,²✉

*Ixodes ricinus* ticks (Acari: Ixodidae) are the most important vector for Lyme borreliosis in Europe. As climate change might affect their distributions and activities, this study aimed to determine the effects of environmental factors, i.e., meteorological, bioclimatic, and habitat characteristics on host-seeking (questing) activity of *I. ricinus* nymphs, an important stage in disease transmissions, across diverse climatic types in France over 8 years. Questing activity was observed using a repeated removal sampling with a cloth-dragging technique in 11 sampling sites from 7 tick observatories from 2014 to 2021 at approximately 1-month intervals, involving 631 sampling campaigns. Three phenological patterns were observed, potentially following a climatic gradient. The mixed-effects negative binomial regression revealed that observed nymph counts were driven by different interval-average meteorological variables, including 1-month moving average temperature, previous 3-to-6-month moving average temperature, and 6-month moving average minimum relative humidity. The interaction effects indicated that the phenology in colder climates peaked differently from that of warmer climates. Also, land cover characteristics that support the highest baseline abundance were moderate forest fragmentation with transition borders with agricultural areas. Finally, our model could potentially be used to predict seasonal human-tick exposure risks in France that could contribute to mitigating Lyme borreliosis risk.

Ongoing anthropogenic climate change impacts redistributions of vector populations and vector-borne diseases have become increasingly concerned as an alarming public health threat<fnref n="1" />. The effects of climate change on the population of *Ixodes ricinus* ticks (Acari: Ixodidae), the most important arthropod vector species in Europe,

<footnote n="1">¹Université de Lyon, INRAE, VetAgro Sup, UMR EPIA, 69280 Marcy l'Etoile, France. ²Université Clermont Auvergne, INRAE, VetAgro Sup, UMR EPIA, 63122 Saint-Genès-Champanelle, France. ³Université Paris-Saclay, AgroParisTech, INRAE, UMR MIA-Paris, 75005 Paris, France. ⁴INRAE, Oniris, UMR BIOEPAR, 44300 Nantes, France. ⁵Global Technical Services, Boehringer-Ingelheim Animal Health, 69007 Lyon, France. ⁶Nancy Laboratory for Rabies and Wildlife, The French Agency for Food, Environmental and Occupational Health and Safety (ANSES), 54220 Malzéville, France. ⁷Université de Toulouse, INRAE, UR CEFS, 31326 Castanet-Tolosan, France. ⁸LTSER ZA PYRénées GARonne, 31326 Auzeville-Tolosane, France. ⁹INRAE, UR EFNO, 45290 Nogent-sur-Vernisson, France. ¹⁰ANSES, ENVA, INRAE, UMR 956 BIPAR, 94701 Maisons-Alfort, France. ¹¹INRAE, US 1116 AgroClim, 84914 Avignon, France. ¹²INRAE, CIRAD, UMR ASTRE, 34398 Montpellier, France. ¹³Present address: Ifremer, RBE-SGMM-LGPMM, 17390 La Tremblade, France. ¹⁴Jean-François Cosson is deceased. ✉email: karine.chalvet-monfray@inrae.fr</footnote>
have become the topic of interest in this region². As ixodid ticks are ectothermic vertebrates, their life cycles and key ecological processes, such as mortality, development rate, and host-seeking (questing) behaviour, are sensitive to weather conditions, particularly temperature and relative humidity³,⁴. Therefore, variations in these meteorological factors contribute to the regulation of the dynamics of tick abundance, phenology of questing activity, and eventually, human-tick exposure risks throughout the year. However, the direct impacts of climate change on tick population dynamics related to abiotic factors, such as temperature, which impacts the development rate⁵, and relative humidity, which impacts the mortality rate⁶⁻⁸, are not always straightforward. Rising temperatures could promote the probability of questing activity and interstadial developmental rate; however, ticks may also face a higher mortality risk at the same time⁴. Meanwhile, relative humidity influences the water regulation of ticks, affecting their questing behaviour and survival⁶⁻⁸. Furthermore, climate change could also pose indirect effects on tick populations through biotic factors, such as habitats⁴ and host populations (both density and behaviours)²,⁹. As a result, climate change has been expected to alter the distribution of tick populations, the activity of questing ticks, and human-tick exposure both geographically (expansions towards higher latitudes and altitudes)¹⁰,¹¹ and seasonally (altered phenology)¹².

*Ixodes ricinus* ticks have wide geographical and climatic ranges spanning all over Europe and some regions in northern Africa and western Asia¹³. Long-term observation of the questing activity of *I. ricinus* nymphs, for both phenology and abundance, has been an important indicator for the exposure risks to ticks and tick-borne diseases in these regions, particularly *Lyme borreliosis*¹⁴⁻¹⁷. *I. ricinus* ticks have been observed to exhibit different phenology (a seasonal variation) of questing activity across their wide range of geo-climatic distributions¹⁸, which determines the time of year when questing ticks are most or least likely to be encountered. The reported phenological patterns appeared to follow a climatic gradient across northern Africa to eastern Europe, with a unimodal pattern with a summer peak and a winter pause in freezing climates (Karelia, Russia), a bimodal pattern with spring and autumn peaks in intermediate climates (European Russia's temperate zone), and a unimodal pattern with a late winter peak in warm climates (Algeria)¹⁸. Meanwhile, the abundance (population size) of ticks in each area reflects the suitability/carrying capacity of the habitats for the populations, both land cover and climatic conditions¹⁹. Furthermore, when the suitability of the habitat is altered, long-term observations can reveal directional changes in questing activity over several years (inter-annual variation)¹⁶. However, such long-term observations that provide baseline data for evaluating how environmental changes will impact tick populations are currently lacking²⁰. Therefore, to evaluate the effects of climate change on the distribution of tick-borne diseases, it is critical to understand the impacts of environmental factors, both abiotic and biotic, on all levels of the variations of tick populations and their questing activities over a long period of time across various geographical/climatic ranges.

However, most of the previous studies on the environmental effects on *I. ricinus* nymph activities in Europe have not captured both the temporal (seasonal and inter-annual) and spatial variations. Some studies focused on short-term observations at multiple sites, allowing for the comparison of abundance estimates but not seasonal and inter-annual variations²¹⁻²⁵. While the other studies were long-term observations at a few sites, capturing seasonal and inter-annual variations and thus unable to compare with multiple sites¹⁴⁻¹⁷. Furthermore, the most widely used method for monitoring nymph activity, cloth-dragging sampling, has a low capture rate²⁶. Therefore, this method produces a high level of uncertainty in abundance estimates across different environmental conditions and field investigators²⁷, which makes comparing results from different observatories difficult.

France is a hotspot for *I. ricinus* and *Lyme borreliosis*²⁸ and has a wide range of geographical and climatic characteristics²⁹. An in-depth understanding of the characteristics of questing activity of *I. ricinus* across various geo-climatic areas is still needed to help control *Lyme borreliosis* in this area. The present study sought to investigate the activity of *I. ricinus* nymphs, both phenology and abundance, across a wide range of climatic region types in France over a long-term observation period. The sampling protocol suggested by a previous study was employed to reduce the uncertainty from the cloth-dragging method²⁷. Also, this study aimed to assess the impacts of environmental factors, such as meteorological, bioclimatic, and habitat characteristics on the variations of *I. ricinus* nymph questing activity. Long-term investigations on nymph questing activity will help us better understand the effects of climate change on the dynamics of *I. ricinus* populations and eventually *Lyme Borreliosis* risk.

## Material and methods

**Sampling sites.** Longitudinal 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 types²⁹ 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.

Each 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
<figure><img src="image_3.png" /><figcaption>**Figure 1.** 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. The map was created using QGIS version 3.8, Zanzibar (https://www.qgis.org). The climatic region types were previously classified by Joly et al.²⁹.</figcaption></figure>

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 started³⁰. 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* nymphs.** Activity of questing *I. ricinus* nymphs was observed by a cloth-dragging sampling technique³¹. 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 used²⁷. 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 identification³² 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 m². 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.
<table><caption>**Table 1.** Environmental variables (meteorological, land cover, topography, and bioclimatic variables) used to explain *I. ricinus* nymph counts per 100 m² in regression analysis. <sup>a</sup>Land cover, topography, and bioclimatic variables were transformed by the principal component analysis as coordinates of PCA dimensions before being used in the regression analysis.</caption><thead><tr><th>Type</th><th>Group</th><th>Variable</th><th>Description</th></tr></thead><tbody><tr><td rowspan="12">Meteorological</td><td rowspan="4">Temperature</td><td>T<sup>0:1</sup><sub>M</sub></td><td>One-month moving average temperature; average temperature from 30 preceding days to the sampling day (°C)</td></tr><tr><td>T<sup>3:6</sup><sub>M</sub></td><td>Previous 3-to-6-month moving average temperature; average mean temperature of preceding 3-6 months (°C)</td></tr><tr><td>T<sup>0:6</sup><sub>M</sub></td><td>Six-month moving average mean temperature; between 6 preceding months to the sampling day (°C)</td></tr><tr><td>T<sup>0:12</sup><sub>M</sub></td><td>Twelve-month moving average mean temperature; between 12 preceding months to the sampling day (°C)</td></tr><tr><td rowspan="4">Relative humidity</td><td>U<sup>0:1</sup><sub>N</sub></td><td>One-month moving average minimum relative humidity; average value from 30 preceding days to the sampling day (%)</td></tr><tr><td>U<sup>3:6</sup><sub>N</sub></td><td>Previous 3-to-6-month moving minimum relative humidity of preceding 3-6 months (%)</td></tr><tr><td>U<sup>0:6</sup><sub>N</sub></td><td>Six-month moving average minimum relative humidity; between 6 preceding months to the sampling day (%)</td></tr><tr><td>U<sup>0:12</sup><sub>N</sub></td><td>12-month moving average minimum relative humidity; between 12 preceding months to the sampling day (%)</td></tr><tr><td rowspan="10">Land cover, topography, and bioclimate<sup>a</sup></td><td rowspan="10">Topography</td><td>Daytime</td><td>Daytime duration (h)</td></tr><tr><td>Mean<sub>elv</sub></td><td>Mean elevation (m)</td></tr><tr><td>sd<sub>elv</sub></td><td>Standard deviation of elevation (m)</td></tr><tr><td>pflat</td><td>Proportion of flat area; slope ≤ 2.5%</td></tr><tr><td>p<sub>north</sub></td><td>Proportion of non-flat area facing north; aspect < 45° or ≥ 315°</td></tr><tr><td>p<sub>east</sub></td><td>Proportion of non-flat area facing east; 45° ≤ aspect < 135°</td></tr><tr><td>p<sub>west</sub></td><td>Proportion of non-flat area facing west; 135° ≤ aspect < 225°</td></tr><tr><td>p<sub>south</sub></td><td>Proportion of non-flat area facing south; 225° ≤ aspect < 315°</td></tr><tr><td>Catchment</td><td>Catchment area</td></tr><tr><td rowspan="4">Bioclimate</td><td>BIO1<sub>temp</sub></td><td>Annual mean temperature (°C)</td></tr><tr><td>BIO2<sub>diur</sub></td><td>Mean diurnal range (°C)</td></tr><tr><td>BIO5<sub>maxTemp</sub></td><td>Maximum temperature of the hottest month (°C)</td></tr><tr><td>BIO12<sub>prec</sub></td><td>Annual precipitation (mm)</td></tr><tr><td rowspan="8">Land cover</td><td>H<sub>CLC1</sub></td><td>Shannon's diversity index for level-1 CLC types</td></tr><tr><td>H<sub>CLC2</sub></td><td>Shannon's diversity index for level-2 CLC types</td></tr><tr><td>H<sub>Forest</sub></td><td>Shannon's diversity index for forest types</td></tr><tr><td>p<sub>Forest</sub></td><td>Percentage of forest-covering area (%)</td></tr><tr><td>n<sub>Forest</sub></td><td>Number of forest patches</td></tr><tr><td>ED<sub>Forest</sub></td><td>Forest edge density (m/km²)</td></tr><tr><td>pH<sub>soil</sub></td><td>Soil pH</td></tr></tbody></table>

**Environmental data.** We 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).

**Daytime duration and meteorological variables.** Daytime duration (*daytime*) from January 2013 to June 2021 at each sampling site was obtained from the corresponding latitude using *geosphere* package³³. 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*<sub>M</sub>, *T*<sub>N</sub>, and *T*<sub>X</sub>; in °C) and relative humidity (*U*<sub>M</sub>, *U*<sub>N</sub>, and *U*<sub>X</sub>; 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* package³⁴. 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 variables¹⁵. 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*<sub>M</sub> and minimum relative humidity *U*<sub>N</sub>. The interval-average variables were defined as the average values of a meteorological variable *M* ∈ {*T*<sub>M</sub>, *U*<sub>N</sub>} during
a period between t₁ to t₂ 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 rates³, 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 conditions⁶⁻⁸. 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 variables.** We 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.0³⁵. The digital elevation model (DEM) data derived from the Shuttle Radar Topography Mission (SRTM) database³⁶ was used to describe the topographic features of sampling sites, which included the mean (mean<sub>elv</sub>) and standard deviation (sd<sub>elv</sub>) of the elevation (in m above sea level), the proportion of flat area ($p_{flat}$; defined by the slope ≤2.5%³⁷), 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 database³⁸, including the annual mean temperature (BIO1<sub>Temp</sub>; in °C), the mean diurnal range (BIO2<sub>Diur</sub>; in °C), the maximum temperature of the warmest month (BIO5<sub>maxTemp</sub>; in °C), and the annual precipitation (BIO12<sub>Prec</sub>; in mm). The land cover features of each sampling site were described using the CORINE Land Cover (CLC) 2018³⁹, while the characteristics of forests were explained by the BD forêt version 2 data⁴⁰. The forest fragmentation was characterized by the percentage of forest-covering area ($p_{Forest}$), the forest edge density (ED<sub>Forest</sub>; in m/km²), 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 index⁴¹ (H) as $H = \sum_{i=1}^{S} p_i \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) database⁴².

**Statistical analysis.** All the statistical analyses were carried out using the programming language R version 3.6.0⁴³. 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_{basei}$) was defined as a summation of monthly median nymph counts $\tilde{N}_i = \{\tilde{N}_{i,t}\}$ across all 12 months $t \in \{1, 2, \dots, 12\}$ and expressed as: $N_{basei} = \sum_{t=1}^{12} \tilde{N}_{i,t}$. Subsequently, the monthly median nymph counts of each site $\tilde{N}_i$ were transformed into normalized monthly median nymph counts $\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}}{\max(\tilde{N}_i)} \qquad (1) $$

The term max($\tilde{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 < 0.05 was considered significant.

The effects of environmental variables (meteorological, land cover, topographical and bioclimatic variables) on the abundance of *I. ricinus* nymphs, were evaluated by a multivariate mixed-effects negative binomial regression, using the glmmTMB package⁴⁴. Recalled that the negative binomial distribution is a generalized form of the Poisson distribution with an additional shape parameter $\theta$ that allows the variance to be independent to the mean $\mu$. The number of *I. ricinus* nymphs per 100 m² was considered the response variable, and the sampling site served as the random effects for intercepts as displayed in Eq. (2):
$$ 
\begin{gathered}
N_i \sim \text{NegBinom}(\mu_i, \theta) \\
\ln \mu_i = \alpha_0 + \alpha_i + X_i \beta
\end{gathered}
\qquad (2) $$

where $N_i$ denoted nymph counts per 100 m² of site i, $\mu_i$ and $\theta$ represented mean of site i and the shape parameter of the negative binomial distribution, respectively. $\alpha_0$ was an intercept of the regression model, while $\alpha_i$ was the random intercept of site i. Finally, $X_i$ and $\beta$ denoted a covariate vector of site i and a vector of regression coefficients for fixed effects, respectively.

The multicollinearity among all environmental variables was evaluated using Pearson's correlation coefficient $r$. Highly correlated variables, with $|r| \ge 0.6$, were not included in the same model. The principal component analysis (PCA) was used to reduce the dimensions of multiple highly correlated land cover, topography, and bioclimatic variables, using FactoMineR⁴⁵ and factoextra⁴⁶ packages. To reduce the number of variables considered in the regression analysis, the coordinates of individual sampling sites on the PCA dimensions were subsequently explored as explanatory variables. Non-linearity and interaction effects of the fixed effects were also investigated.

Multivariable models were explored by both forward and backward selection of significant variables decreasing the Akaike information criterion (AIC). Variations of the data explained by fixed effects were monitored by a pseudo-R² statistic using the MuMln package⁴⁷. Residuals of the models were evaluated using the DHARMA package⁴⁸,⁴⁹. The performance of the best-fitted model was assessed by a resampling method. Briefly, 50% of the original dataset was randomly sampled to fit with the best model for 500 iterations, referred to as resampling models. The distributions of regression coefficients across 500 resampling models were compared against those of the model with the original dataset, called the complete model. Additionally, the coverage probability (CP) of each explanatory variable, defined as the percentage of resampling models that produced regression coefficients that were significantly different from zero ($p < 0.05$), were also calculated. The CP index is a measure of how robust each variable is to explain observed data. Finally, the effects of each environmental variable were predicted using the effects package⁵⁰.

## Results

**Characteristics of sampling sites.** During 2014 to 2021, the activity of questing *I. ricinus* nymphs was observed from 11 sample sites at 7 tick observatories in France. The characteristics of each sampling site, i.e., meteorological, land cover, topography, and bioclimatic features, appeared to be highly correlated (Fig. S1). According to the meteorological conditions (Fig. 2, Figs. S2, S3A,B), the characteristics of tick observatories can be classified as: (1) cold average temperature with the average winter temperature <2.5 °C (Saint-Genès-Champanelle and Velaine-en-Haye); (2) warm average temperature with average winter temperature ~7 °C and high summer humidity with average minimum relative humidity in summer >60% (Gardouch and Carquefou); (3) intermediate average temperature with average winter temperature between 4 and 6 °C and low summer humidity with average minimum relative humidity in summer <55% (La Tour de Salvagny, and Etiolles).

The land cover, topography, and bioclimatic characteristics of each sampling site were described in Tables S3–S5, respectively. The first two dimensions of the principal component analysis (PCA) explained 70.3% of the variation in the land cover, topography, and bioclimatic variables (Fig. 3A and Fig. S4). The first dimension differentiated each sampling site predominantly based on forest and land cover characteristics, in which the first 4 contributions were the number of forest patches $n_{Forest}$, percentage of forest-covering area $p_{Forest}$, Shannon's diversity index for forest types $H_{Forest}$, and Shannon's diversity index for level-2 CORINE Land Cover (CLC) types $H_{CLC2}$, respectively (Fig. 3A and Fig. S4). The Shannon's diversity index for forest types $H_{Forest}$ was found to be negatively correlated with the Shannon's diversity index for the CLC types, $H_{CLC1}$ and $H_{CLC2}$, (Fig. S1). The sampling sites in La Tour de Salvagny A, La Tour de Salvagny B, Gardouch Inside and External were relatively high fragmented forests, less percentage of forest-covering areas, more diverse in the CORINE Land Cover types, and less diverse in the forest types. While the sampling sites in Saint-Genès-Champanelle, Les Bordes and Etiolles were covered by a high percentage of forest areas with less forest fragmentation and high diversity in forest types.

The second PCA dimension of was driven by topographic and bioclimatic features, such as mean mean$_{elv}$ and standard deviation sd$_{elv}$ of elevation, proportion of non-flat area facing north $p_{north}$, annual precipitation BIO12$_{prec}$, and annual mean temperature BIO1$_{Temp}$ (Fig. 3A and Fig. S4). The sampling site in Saint-Genès-Champanelle was located at a high uneven elevation above the sea level (mountainous) with higher annual precipitation and lower annual temperature compared to all other sampling sites. While sampling sites of Les Bordes and Etiolles were located at flat areas with the lowest annual precipitation. The coordinates for each sampling site on the first and second dimensions on the PCA plot (Fig. 3B) were used to describe land cover, topography, and bioclimatic characteristics in subsequent regression analysis, denoted as Dim1 and Dim2, respectively. In addition, soil in the forest of Les Bordes appeared to be acidic, with an estimated pH of 4 (Table S3).

**Abundance indicators of questing *Ixodes ricinus* nymphs.** During the observation period, 631 campaigns were carried out, yielding a total of 22,912 *I. ricinus* nymphs (Table S1). Since the study design did not allow us to estimate the true abundance of questing nymphs, this study reported confirmed *I. ricinus* nymph counts, monthly median nymph counts, and baseline annual nymph count as indicators of the true abundance (Figs. 4, 5, and Fig. S3C). The baseline annual nymph count was highest in Carquefou (802 per 100 m²), Gardouch Inside (717.5 per 100 m²), and La Tour de Salvagny B (513.5 per 100 m²), while it was lowest in Les Bordes A to C (30–85.5 per 100 m²), La Tour de Salvagny A (122.5 per 100 m²), and Saint-Genès-Champanelle (138 per 100 m²) (Fig. 5).

**Phenological patterns of questing *Ixodes ricinus* nymphs.** The normalized monthly median nymph counts revealed three phenological patterns of questing nymph activity observed among the 11 sampling sites
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Monthly average meteorological conditions of 7 tick observatories: (A) average temperature *T*<sub>M</sub>; (B) minimum relative humidity *U*<sub>N</sub>. The average values were calculated using the imputed data from January 2013 to June 2021, with the exception of Les Bordes, which was calculated from January 2017 to June 2021. The points represented the median of the monthly average values over the observed period, while the error bars represent the maximum and minimum monthly average values. Background colours indicated meteorological seasons of the temperate areas of northern hemisphere: Spring, 1st March to 31st May (green); Summer, 1st June to 31st August (yellow); Autumn, 1st September to 30th November (brown); Winter, 1st December to 28th or 29th February (blue).</figcaption></figure>

(Fig. 6 and Fig. S3D): (1) Pattern 1: a unimodal pattern with a relatively long peak activity beginning in late spring, reaching its maximum in summer and fading out in autumn, with no or very low activity in winter (Saint-Genès-Champanelle, Velaine-en-Haye, and La Tour de Salvagny B); (2) Pattern 2: a bimodal pattern with an earlier nymph activity reaching its peak in spring (the first peak), fading out in summer, reappearing in autumn (the second peak), and fading out again in winter (Les Bordes A to C, La Tour de Salvagny A, and Etiolles); (3) Pattern 3: a unimodal pattern with moderate nymph activity throughout autumn and winter, reaching its peak in spring, and fading out quickly in summer (Carquefou, Gardouch Inside, and Gardouch Outside). Interestingly, questing nymphs have also been observed during winter in most sites. In addition, a decreasing trend in annual maximum nymph counts was observed more evidently in Gardouch Outside (Spearman's rank correlation coefficient = -0.8982; *p*-value < 0.05), while no apparent directional change was observed in other sites (Fig. 4 and Fig. S5).

**Environmental impacts on nymphs activity.** Several mixed-effects negative binomial regression models were evaluated to explain impacts of environmental factors on observed *I. ricinus* nymph counts. The best-fitted model was selected by the AIC value, while other model characteristics were compared by dispersion parameter, variance of the random intercepts, and pseudo-R² statistics (Table 2). Fixed effects in Model 5, predominantly meteorological variables, explained 65.1% of the total variance of nymph counts.

Model diagnostics showed no indication of lack of fit (Kolmogorov-Smirnov test: *p*-value = 0.2067), under- or over-dispersion in the residuals (*p*-value = 0.08), or residual outliers (*p*-value = 1) (Fig. S6). The strength of the best-fitted model was assessed by resampling half of the data to fit with this model. The regression coefficients yielded from 500 resampling models were similar to those of the complete model with mostly high coverage probability (Fig. S7 and Table S6).

The best-fitted model (Model 5) suggested that phenology of *I. ricinus* nymphs non-linearly responded to meteorological factors, such as 1-month moving average temperature $T_M^{0.1}$, previous 3-to-6-month moving average temperature $T_M^{3.6}$, and 6-month moving average minimum relative humidity $U_N^{0.6}$. The interaction effect
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Principal component analysis (PCA) for land cover, topography, and bioclimate characteristics of all sampling sites: (A) PCA plot for land cover, topographical, and bioclimatic variables. Colours indicate whether the variable contributes more to Dimension 1 (red), Dimension 2 (blue), or other higher dimensions (black), using the cumulative contribution of ~60%; (B) PCA plot for individual sampling sites. Coordinates of the individual PCA plot were subsequently used in the regression analysis.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Number of *Ixodes ricinus* nymphs per 100 m² obtained from 11 sampling sites. Background colours indicate meteorological seasons of temperate areas of the northern hemisphere: Spring, 1st March to 31st May (green); Summer, 1st June to 31st August (yellow); Autumn, 1st September to 30th November (brown); Winter, 1st December to 28th or 29th February (blue). Vertical dashed lines indicate 1-year interval.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Monthly median nymph counts (left), and baseline annual nymph counts (right) obtained from 11 sampling sites. Meteorological seasons of temperate areas of the northern hemisphere: Spring, 1st March to 31st May; Summer, 1st June to 31st August; Autumn, 1st September to 30th November; Winter, 1st December to 28th or 29th February.</figcaption></figure>

between the 1-month moving average temperature $T_M^{0:1}$ (weather) and the annual average temperature $BIOTemp$ (climate) explained that the nymph counts observed at different climates arrive at peaks at different 1-month moving average temperatures, where nymph counts in warmer climates peaked at lower 1-month moving average temperatures (Fig. 7A). Also, the best-fitted model explained that nymph counts were highest at low previous 3-to-6-month moving average temperatures $T_M^{3:6}$ (Fig. 7B) and high 6-month moving average minimum relative humidity $U_N^{0:6}$. (Fig. 7C). In addition, it should be noted that because the daytime duration was highly correlated with several meteorological variables, such as 1-month moving average temperature $T_M^{0:1}$, previous 3-to-6-month moving average temperature $T_M^{3:6}$, and 1-month moving average minimum relative humidity $U_N^{0:1}$ (Fig. S1); therefore, it was excluded from the best-fitted model.

The Dimension-1 PCA coordinates *Dim1*, representing the forest and land cover characteristics of different sampling sites, markedly reduced the variance of random intercepts from 1.0004 in Model 4 to 0.4224 in Model 5 (Table 2). Thus, the forest and land cover characteristics in Model 5 significantly explained model intercepts which indicates their impacts on the baseline nymph count levels. The landscape that supports more baseline nymph counts was suggested to have a moderate percentage of forest-covering areas $p_{Forest}$, number of forest patches $n_{Forest}$, and Shannon's diversity index for level-1 ($H_{CLC1}$) and level-2 CORINE Land Cover types ($H_{CLC2}$) (Fig. 7D). It should be noted that the Dimension-2 PCA coordinates *Dim2*, driven mainly by the mean and standard deviation of the elevation above sea level, was not significant and therefore not included in the final model.

Finally, the prediction from the best-fitted model replicated well the phenological patterns of observed *I. ricinus* nymphs (Fig. 6 and Fig. S8), particularly sampling sites with more observations, i.e., Saint-Genès-Champanelle, Velaine-en-Haye, Etiolles, Carquefou, Gardouch Inside, and Gardouch Outside. However, the autumn peaks of the phenological Pattern 2 (Fig. 6) and the interannual variations of the nymph counts (Fig. S8) remain imperfectly predicted by the model using meteorological variables.

## Discussion

The variations in tick activity across space and time can be classified into four levels: (1) baseline abundance variation; (2) sub-daily variation; (3) seasonal variation; (4) inter-annual variation. The baseline abundance reflects the suitability/carrying capacity of the habitats for the populations, both land cover and climatic conditions¹⁹. The sub-daily variation arises from the immediate responses of questing ticks to rapid microenvironment changes, such as local temperature, relative humidity⁵¹, solar altitude, or presence of dew drops, to maintain their water balance⁵². Such responses reflect the instantaneous probability of questing. While seasonal and inter-annual variation reflects the long-term cumulative results of the questing behaviour as well as reproduction, development, and mortality rates on the tick population³. However, previous studies on the effects of the environment on *I. ricinus* nymph activities in Europe did not capture all the aforementioned levels of variation. They largely focused on either a short-term period for multiple sites²¹⁻²⁴ or a long-term period for a few areas¹⁴⁻¹⁷. The present study attempted to explore the baseline abundance, seasonal, and inter-annual variations of *I. ricinus* nymph activity by expanding the observation across a wide range of geographical/climatic areas over a long period, involving a total
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Phenological patterns of *I. ricinus* nymph activity observed at 11 sampling sites: (A) observed and predicted normalized monthly median nymph counts of each sampling site. The value of 1 indicates the maximum nymph activity (peak), while the value of 0 indicates the absence of nymph activity. Line types distinguish between the observation (solid) and prediction from the best-fitted model (dashed). Line colours indicate different phenological patterns of the observed data: Pattern 1, a unimodal pattern with a summer peak and a winter pause (blue); Pattern 2, a bimodal pattern with a greater activity peak in spring and a smaller peak in autumn (green); Pattern 3, a unimodal pattern with a spring peak without winter pause (red). Text annotations demonstrate mean elevation above the sea level *mean*<sub>elv</sub> (m) and average annual temperature *BIO*1<sub>temp</sub> (°C). Numbers at the bottom represent the number of observations of each month throughout the study period; (B) normalized monthly median nymph counts observed at each site (circle) and the overall trend of each phenological pattern (diamond and line). The overall trend of each pattern was determined using the medians of the values from all sites that showed that pattern, as indicated by text annotations. Background colours indicate meteorological seasons of the temperate areas of northern hemisphere: Spring, 1st March to 31st May (green); Summer, 1st June to 31st August (yellow); Autumn, 1st September to 30th November (brown); Winter, 1st December to 28th or 29th February (blue).</figcaption></figure>
<table><caption>**Table 2.** The Akaike information criterion (AIC), dispersion parameter of the negative binomial regression, and the variance of the intercepts when sampling sites are considered a random effect. Asterisks indicated the interaction terms between 2 variables.</caption><thead><tr><th>Model</th><th>Description</th><th>AIC</th><th>Dispersion parameter</th><th>Variance of intercepts</th><th>Pseudo-R<sup>2</sup></th></tr></thead><tbody><tr><td>0</td><td>Null</td><td>5447.6</td><td>0.70449</td><td>1.0057</td><td>0.0000</td></tr><tr><td>1</td><td>(T<sub>M</sub><sup>0:1</sup>)<sup>3</sup> + (T<sub>M</sub><sup>0:1</sup>)<sup>2</sup> + T<sub>M</sub><sup>0:1</sup></td><td>5272.6</td><td>0.92771</td><td>1.0004</td><td>0.2495</td></tr><tr><td>2</td><td>(T<sub>M</sub><sup>0:1</sup>)<sup>3</sup> + (T<sub>M</sub><sup>0:1</sup>)<sup>2</sup> + T<sub>M</sub><sup>0:1</sup> * BIO1<sub>temp</sub></td><td>5256.5</td><td>0.95624</td><td>0.9707</td><td>0.2729</td></tr><tr><td>3</td><td>(T<sub>M</sub><sup>0:1</sup>)<sup>3</sup> + (T<sub>M</sub><sup>0:1</sup>)<sup>2</sup> + T<sub>M</sub><sup>0:1</sup> * BIO1<sub>temp</sub> + (T<sub>M</sub><sup>3:6</sup>)<sup>3</sup> + (T<sub>M</sub><sup>3:6</sup>)<sup>2</sup> + T<sub>M</sub><sup>3:6</sup></td><td>4858.0</td><td>1.98588</td><td>0.9814</td><td>0.617</td></tr><tr><td>4</td><td>(T<sub>M</sub><sup>0:1</sup>)<sup>3</sup> + (T<sub>M</sub><sup>0:1</sup>)<sup>2</sup> + T<sub>M</sub><sup>0:1</sup> * BIO1<sub>temp</sub> + (T<sub>M</sub><sup>3:6</sup>)<sup>3</sup> + (T<sub>M</sub><sup>3:6</sup>)<sup>2</sup> + T<sub>M</sub><sup>3:6</sup> + (U<sub>N</sub><sup>0:6</sup>)<sup>3</sup> + (U<sub>N</sub><sup>0:6</sup>)<sup>2</sup> + U<sub>N</sub><sup>0:6</sup></td><td>4823.9</td><td>2.13044</td><td>1.0004</td><td>0.6406</td></tr><tr><td>5</td><td>(T<sub>M</sub><sup>0:1</sup>)<sup>3</sup> + (T<sub>M</sub><sup>0:1</sup>)<sup>2</sup> + T<sub>M</sub><sup>0:1</sup> * BIO1<sub>temp</sub> + (T<sub>M</sub><sup>3:6</sup>)<sup>3</sup> + (T<sub>M</sub><sup>3:6</sup>)<sup>2</sup> + T<sub>M</sub><sup>3:6</sup> + (U<sub>N</sub><sup>0:6</sup>)<sup>3</sup> + (U<sub>N</sub><sup>0:6</sup>)<sup>2</sup> + U<sub>N</sub><sup>0:6</sup> + (Dim1)<sup>3</sup> + (Dim1)<sup>2</sup> + Dim1</td><td>4811.9</td><td>2.13050</td><td>0.4224</td><td>0.6507</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Predicted effects of environmental factors on the number of *I. ricinus* nymphs per 100 m² given the average values of all other variables: (A) 1-month moving average temperature, an averaged mean temperature between previous 30 days and sampling day (T<sub>M</sub><sup>0:1</sup>) with an interaction effect with annual temperature (BIO1<sub>Temp</sub>); (B) previous 3-to-6-month moving average temperature, an averaged mean temperature between 3 and 6 months prior to the sampling day (T<sub>M</sub><sup>3:6</sup>), (C) 6-month moving average minimum relative humidity, an averaged minimum relative humidity between previous 6 months and the sampling day (U<sub>N</sub><sup>0:6</sup>); (D) coordinates on Dimension 1 of the principal component analysis (Dim1), representing forest and land cover characteristics. Gray shaded areas indicate a 95% confidence interval. Rugs display the distribution of observed values.</figcaption></figure>

of 11 sampling sites from 7 tick observatories in France over 8 years. With the repeated removal sampling design, we controlled the uncertainty on the proxy measure of tick abundance (nymph counts per 100 m²) that arises from the variability of sampling rate across different sites and times<sup>27,53</sup>. Furthermore, all transects were covered with low vegetation, which have been shown to maximize sampling rates compared to other types of vegetation floors<sup>27</sup>, resulting in higher accuracy in estimating the proxy of tick abundance. To avoid the uncertainty caused by the variations of the environment within a sampling site, we repeated tick sampling from the same transects throughout the observation period. Furthermore, we strived to minimize the interference of tick populations and their environment during the long-term observation period by using a 1-month sampling interval. As all of the transects were open areas, ticks were free to relocate in and out of the transects, either directly or through hosts, and compensate for the monthly removals of tick populations. Nevertheless, the sub-daily variations in questing behaviour or immediate response to meteorological conditions among nymphs from different regions were not captured by this study design.

Furthermore, the activity of *I. ricinus* in other life stages (larvae and adults) was suggested to respond differently to the environment than the nymph stage, depending on their population cohorts<sup>54</sup>. However, the protocol and performance of sampling methods for estimating tick abundance were well described only for the nymph stage<sup>27,31</sup>, and to our knowledge, no studies have proposed a standard protocol and assessed the performance of the sampling method for larva and adult stages of *I. ricinus*. As a result, the sampling method used in our study
did not allow us to observe the activity of larvae and adults. To investigate potential environmental effects on the activity of these stages, more research is needed on the performance of sampling protocols for larvae and adults that give us accurate abundance estimates and allow us to compare several environmental conditions.

Phenological patterns of *I. ricinus* nymph activity observed in our study tended to be driven by climatic characteristics. In general, nymph activity was markedly reduced during two types of unfavourable temperature conditions: (1) cold conditions with monthly average temperature <3 °C; (2) hot conditions with monthly average temperature > 16.5 °C. Temperature conditions in mountainous climate (Saint-Genès-Champanelle and Velaine-en-Haye) were intensely unfavourable in winter, resulting in a near-complete halt. In contrast, conditions in summer were not as severe, resulting in a reduction in nymph activity arriving in mid-to-late summer. As a result, a unimodal pattern with a summer peak and complete winter pause was observed in this climate (Pattern 1). However, winter temperatures were relatively benign in the oceanic (Carquefou) and south-west basin (Gardouch) climates, allowing nymph activity to increase. Meanwhile, unfavourable conditions were severe during the summer, resulting in a sharp decrease in nymph activity in early summer, and therefore, a unimodal pattern without a winter pause (Pattern 3). In intermediate climates (Les Bordes and Etiolles), winter brought unfavourable conditions though not as severe as in mountainous climates, resulting in a reduction in activity but not complete cessation. While summer also brought unfavourable conditions, resulting in a significant decrease in activity from early to mid-summer, and thus a spring-autumn bimodal pattern (Pattern 2). However, despite the phenological patterns 1 and 3 being clearly distinct, we found that the existence of the autumn peaks in pattern 2 was ambiguous for the following reasons: (1) most sites with phenological pattern 2 had a low baseline abundance level. One additional nymph count due to variability of the sampling could lead to a significant false signal for the autumn activity; (2) compared to other patterns, tick activity at most sites with phenological pattern 2 was observed for a shorter period (a smaller sample size). A longer observation period is required to confirm the presence of bimodal patterns at these sites. Furthermore, we observed different phenological patterns at La Tour de Salvagny A and B, despite these two sites being approximately 2-km apart. The heterogeneity of the phenological patterns could have arisen from the following reasons: (1) THE phenological pattern 2 observed from a short-term period at La Tour de Salvagny A was ambiguous, as previously discussed; (2) different landscapes and host densities in these sites influenced differently tick population dynamics and apparent tick activities.

The overall phenological patterns observed in our study followed a climatic gradient previously reported: (1) a unimodal pattern with a summer peak and a winter pause (Pattern 1) in freezing climates (Karelia, Russia); (2) a bimodal pattern with spring-autumn peaks (Pattern 2) in intermediate climates (European Russia's temperate zone); (3) a unimodal pattern with a late winter peak (Pattern 3) in warmer climates (Algeria)¹⁸. Also, a similar phenology trend observed in one mountain region in Switzerland was reported to follow an elevation gradient where the spring-autumn bimodal phenology (Pattern 2) shifted to a late-spring unimodal phenology (Pattern 1) at a higher elevation¹⁷, presumably due to a temperature gradient. A similar climatic trend was observed in *Ixodes scapularis* nymphs in a study conducted across a large geographic area in the United States, with peak nymph activity occurring in the summer in colder regions, and late spring in warmer regions⁵⁵. However, the climatic gradient of phenological patterns observed in our study might not apply to all *I. ricinus* populations. For example, an *I. ricinus* population in a freezing climatic area (Finland) has been reported to have a spring-autumn bimodal phenology (Pattern 2)⁵⁶, rather than a summer unimodal phenology (Pattern 1) as suggested by our results. Further, *I. ricinus* populations from different climates were reported to begin their host-seeking behaviour at different temperature thresholds, either through local adaptation or phenotypic plasticity⁵⁷,⁵⁸, which may result in different phenological patterns.

However, according to the climatic gradient observed in our study, it is conceivable to hypothesize that a gradual shift in climate due to climate change towards warmer temperature and lower relative humidity will result in the following phenological changes to *I. ricinus* tick populations in France: (1) in cold climates, such as in Saint-Genès-Champanelle and Velaine-en-Haye, climate change towards warmer temperatures could shift activity peaks towards early spring, and there would be no longer a winter pause; (2) in warm climates, such as in Carquefou and Gardouch, climate change towards even warmer and drier conditions could limit the distributions of suitable habitats that could support the physiological needs of *I. ricinus* ticks, i.e., excessive thermal stress and difficulties in maintaining their water balance⁷,⁸.

Furthermore, we found that describing phenological patterns solely based on their shape, such as unimodal or bimodal, was insufficient. Rather, the timing of the phenological peak and the winter traits must be studied. To better understand the relationship between climatic gradient and phenological patterns, we recommend longitudinal studies conducted over large geographic areas to report the phenological pattern of each sampling location separately. Moreover, determining the age of nymphs through fat content analysis has been demonstrated to guide formulating hypotheses about the mechanisms for the population dynamics of each phenological pattern¹⁷.

As the computation for cumulated nymph counts using areas under the density curve¹⁶ requires strict, regular sampling intervals throughout the study period, we used annual maximum nymph counts as a proxy to explore inter-annual variations in questing nymph activity. Directional changes in the annual maximum nymph counts were observed only in Gardouch Outside but not Gardouch Inside, despite being located in the same area. We hypothesized that it could have occurred as a result of (1) a possible decrease in roe deer population outside of the experimental station, while the number of roe deer was always maintained at a higher inside the experimental station; or (2) the periods of heatwaves that impacted this area in summer 2019 and summer 2020 affected the tick survival, while the high host density in Gardouch Inside compensated for the tick mortality. Furthermore, a recent study suggested that the inter-annual fluctuation of nymph abundance could be explained by the level of tree seed production, and thus the abundance level of small mammal hosts, in previous years⁵⁹. However, retrospective data on seed production were not available for all of our sampling sites during the study period, which precluded us from exploring this hypothesis in this study.
According to the regression analysis, the relationship between observed nymph counts (response variable) and interval-average meteorological variables was non-linear, indicating that the impacts on nymph counts were not the same at different levels of these meteorological variables. For example, when the 1-month average temperature $T_M^{0:1}$ was too high or too low, observed nymph counts decreased significantly, with a peak at an optimal value of $T_M^{0:1}$, which was dependent on the bioclimatic condition (Fig. 7A). Also, nymph counts were steadily low when the previous 3-to-6-month moving average temperature $T_M^{3:6}$ was greater than 15 °C, increased when $T_M^{3:6}$ was between 2 and 15 °C, and peaked at around 0–2 °C (Fig. 7B). Finally, observed nymph counts dramatically increased when the 6-month moving average minimum relative humidity $U_N^{0:6}$ was greater than 90%, while nymph counts were relatively stable when $U_N^{0:6}$ was less than 90% (Fig. 7C).

As the number of nymphs collected in each sampling was a product of total population size of nymphs at the time of sampling (true abundance), the proportion of nymphs questing for a host, and the sampling rates, we hypothesized that the meteorological conditions influenced each of these unobserved parameters at different time lags. The population size is determined by the cumulative effects of past meteorological conditions, either directly (mortality rate or larva-to-nymph developmental rate) or indirectly (emigration/immigration from host movement). While the proportion of questing nymphs tended to respond to short-term fluctuation of the weather conditions⁴. A longitudinal observation on one site in Germany found similar effects of interval-average variables that explain the phenological pattern and inter-annual variation of *I. ricinus* activity¹⁵, which were current temperature $T_M^{0:0}$, previous 4-to-6-month average temperature $T_M^{4:6}$, and 1-month moving average relative humidity $U_M^{0:1}$. Furthermore, the best-fitted model also showed a significant interaction between the 1-month moving average temperature and a bioclimatic variable (average annual temperature $BIO_{Temp}$). This result indicated that although *I. ricinus* from colder climates start questing at lower temperature thresholds⁵⁷, their activity arrives at its peak at higher temperatures (in summer; temperature ~ 14 to 16 °C). While tick activity in warmer climates peaks at lower temperatures (in spring; temperature ~ 11 to 14 °C). Presumably, the population size of questing nymphs in cold climates is unlikely to reach its peak during a colder spring due to developmental diapause (suppression of interstadial development of engorged ticks) or behavioural diapause (suppression of host-seeking behaviour activity of active ticks)⁵⁴,⁶⁰. In addition, photoperiod has also been proposed to regulate both behavioural and developmental diapause of ixodid ticks⁵⁴; therefore, the exclusion of daytime duration from the best-fitted model, which was highly correlated with several meteorological variables, does not rule out the effects of photoperiod on tick activity dynamics.

The regression model also suggested that the baseline abundance level was associated with land cover characteristics. Sampling sites with a transition between moderately fragmented forest and non-forest areas, such as grassland or agricultural areas, support higher baseline nymph abundance (La Tour de Salvagny B, Gardouch Inside, Gardouch Outside, and Carquefou; $p_{Forest}$ ~ 17 to 52%). While landscapes with a high proportion of forest covers ($p_{Forest}$ > 90%) tended to have intermediate (Velaine-en-Haye and Etiolles) or low (Les Bordes A to C, and Saint-Genès-Chamanelle) baseline abundance level. These findings are consistent with previous studies⁶¹⁻⁶³, indicating that forest fragmentation created transitional forest areas known as ecotones, which attract several mammal hosts²³ and lead to a higher tick abundance. In addition, it has been suggested that forest fragmentation increases not only tick abundance but also the prevalence of pathogen-infected ticks⁶⁴ and the risks of human-tick exposure⁶⁵, which may lead to a higher incidence of tick-borne diseases. We also observed markedly low baseline *I. ricinus* abundance in the highly fragmented forest with a low proportion of forest-covering area (La Tour de Salvagny A; $p_{Forest}$ = 6.5%), possibly due to the lack of large foraging hosts in such a landscape⁶⁶. Furthermore, we found a remarkably low baseline abundance in all three sampling sites in Les Bordes, which could not be fully explained biologically by our available environmental factors. This could be explained by the unsuitable conditions of Les Bordes forest, where superimposed layers of clay and sand lead to a temporary perched water table in winter, and the low soil water storage capacity leads to arid soil conditions in the upper parts of the soil profile during summer months with low rainfall. Other additional factors could also contribute to explain the low tick abundance at Les Bordes site such as host abundance, competition, or predation of a forest with such high biodiversity.

Desirable model diagnostic results and repeated resampling models using half of the original data produced similar estimate values over 500 iterations, indicating that all explanatory variables in the best-fitted model are robust. Model predictions also successfully predicted baseline abundance, phenological patterns, and inter-annual variations of *I. ricinus* nymph activity. However, the prediction from the best-fitted model still imperfectly captured the autumn peak in the phenological Pattern 2, where the sampling sites with this pattern had both fewer observations and lower abundance. This has brought us to the following inconclusive explanations: (1) the model failed to replicate the apparent autumn activity peak more clearly; (2) the apparent autumn peak was observed at random, and one additional nymph collected could generate a false signal when the sample size and baseline abundance were both small. As a result, it is envisioned that a more extended observation period and more sampling sites encompassing all phenological patterns should be implemented. Finally, although the regression analysis has been primarily used to find robust associations of environmental factors and help formulate ecological hypotheses on tick activity, it is considered a phenomenological model that does not fully explain the mechanisms of the observed phenomena³. In light of growing concerns about the impact of climate change on tick-borne diseases, particularly Lyme borreliosis, the mechanistic modelling framework for *I. ricinus* population dynamics is still required to investigate these ecological hypotheses.

## Conclusions

The study design in this study allowed us to longitudinally monitor the questing *I. ricinus* nymph activity, which is comparable across multiple tick observatories. Two distinct phenological patterns were identified in colder climates (a unimodal pattern with a summer peak and complete winter pause) and warmer climates (a unimodal
pattern with a spring peak without winter pause). While a longer observation period is still required to confirm the presence of a spring-autumn bimodal phenological pattern in intermediate climates. Furthermore, the regression analysis revealed that phenological patterns and inter-annual variations in questing activity were associated with meteorological variables at different lag times. Our results also supported the notion that moderate forest fragmentation promotes higher nymph abundance. Our findings can guide formulating ecological hypotheses about the effects of climate change on the distribution, abundance, and phenology of this tick species, which can then be assessed using a mechanistic population dynamics model. Finally, the phenological patterns of *I. ricinus* activity predicted by our model could potentially be employed to develop tick-borne disease prevention tools, such as developing a seasonally dynamic human-tick exposure risk map, which could help to reduce Lyme borreliosis cases in France.

## Data availability

The datasets generated during and/or analysed during the current study are available in the “TEMPO—Réseau National d'Observatoires de la Phénologie” repository, https://doi.org/10.15454/ZSYGUM.

Received: 22 December 2021; Accepted: 31 March 2022

Published online: 12 May 2022

## References

1. IPCC. *Climate Change 2014: Impacts, Adaptation, and Vulnerability. Part A: Global and Sectoral Aspects*. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. (2014).

2. Gilbert, L. The impacts of climate change on ticks and tick-borne disease risk. *Annu. Rev. Entomol.* **66**, 273–288 (2021).

3. Ostfeld, R. S. & Brunner, J. L. Climate change and *Ixodes* tick-borne diseases of humans. *Philos. Trans. R. Soc. B Biol. Sci.* **370**, 1–11 (2015).

4. Ogden, N. H. & Lindsay, L. R. Effects of climate and climate change on vectors and vector-borne diseases: Ticks are different. *Trends Parasitol.* **32**, 646–656 (2016).

5. Randolph, S. E., Green, R. M., Hoodless, A. N. & Peacey, M. F. An empirical quantitative framework for the seasonal population dynamics of the tick *Ixodes ricinus*. *Int. J. Parasitol.* **32**, 979–989 (2002).

6. Needham, G. R. & Teel, P. D. Off-host physiological ecology of ixodid ticks. *Annu. Rev. Entomol.* **36**, 659–681 (1991).

7. Milne, A. The ecology of the sheep tick, *Ixodes ricinus* L.: Microhabitat economy of the adult tick. *Parasitology* **40**, 14–34 (1950).

8. Wongnak, P. *et al*. A hierarchical Bayesian approach for incorporating expert opinions into parametric survival models: A case study of female *Ixodes ricinus* ticks exposed to various temperature and relative humidity conditions. *Ecol. Modell.* **464**, 1–12 (2022).

9. Morellet, N. *et al*. Seasonality, weather and climate affect home range size in roe deer across a wide latitudinal gradient within Europe. *J. Anim. Ecol.* **82**, 1326–1339 (2013).

10. Sagurova, I. *et al*. Predicted northward expansion of the geographic range of the tick vector *Amblyomma americanum* in North America under future climate conditions. *Environ. Health Perspect.* **127**, 1–14 (2019).

11. Ogden, N. H., Beard, C. B., Ginsberg, H. S. & Tsao, J. I. Possible effects of climate change on ixodid ticks and the pathogens they transmit: Predictions and observations. *J. Med. Entomol.* **58**, 1536–1545 (2021).

12. Macdonald, A. J., Mccomb, S., O'Neill, C., Padgett, K. A. & Larsen, A. E. Projected climate and land use change alter western blacklegged tick phenology, seasonal host-seeking suitability and human encounter risk in California. *Glob. Chang. Biol.* **26**, 5459–5474 (2020).

13. European Centre for Disease Prevention and Control. *Ixodes ricinus*—current known distribution: March 2021. (2021). https://www.ecdc.europa.eu/en/publications-data/ixodes-ricinus-current-known-distribution-march-2021. (Accessed 20th July 2021)

14. Alonso-Carné, J., García-Martín, A. & Estrada-Peña, A. Modelling the phenological relationships of questing immature *Ixodes ricinus* (Ixodidae) using temperature and NDVI data. *Zoonoses Public Health* **63**, 40–52 (2016).

15. Brugger, K., Walter, M., Chitimia-Dobler, L., Dobler, G. & Rubel, F. Seasonal cycles of the TBE and Lyme borreliosis vector *Ixodes ricinus* modelled with time-lagged and interval-averaged predictors. *Exp. Appl. Acarol.* **73**, 439–450 (2017).

16. Hauser, G. *et al*. Influence of climatic factors on *Ixodes ricinus* nymph abundance and phenology over a long-term monthly observation in Switzerland (2000–2014). *Parasit. Vectors* **11**, 1–12 (2018).

17. Bregnard, C. *et al*. Beech tree masting explains the inter-annual variation in the fall and spring peaks of *Ixodes ricinus* ticks with different time lags. *Parasit. Vectors* **14**, 1–26 (2021).

18. Korenberg, E. I. Seasonal population dynamics of *Ixodes* ticks and tick-borne encephalitis virus. *Exp. Appl. Acarol.* **24**, 665–681 (2000).

19. Weber, M. M., Stevens, R. D., Diniz-filho, J. A. F. & Grelle, C. E. V. Is there a correlation between abundance and environmental suitability derived from ecological niche modelling? A meta-analysis. *Ecography* (Cop.) **40**, 817–828 (2017).

20. Nuttall, P. A. Climate change impacts on ticks and tick-borne infections. *Biologia Bratisl.* https://doi.org/10.1007/s11756-021-00927-2 (2021).

21. Agoulon, A. *et al*. A vegetation index qualifying pasture edges is related to *Ixodes ricinus* density and to *Babesia divergens* seroprevalence in dairy cattle herds. *Vet. Parasitol.* **185**, 101–109 (2012).

22. Boehnke, D. *et al*. Estimating *Ixodes ricinus* densities on the landscape scale. *Int. J. Health Geogr.* **14**, 1–12 (2015).

23. Perez, G. *et al*. Effect of landscape features on the relationship between *Ixodes ricinus* ticks and their small mammal hosts. *Parasit. Vectors* **9**, 1–18 (2016).

24. Hauck, D., Springer, A., Chitimia-Dobler, L. & Strube, C. Two-year monitoring of tick abundance and influencing factors in an urban area (city of Hanover, Germany). *Ticks Tick. Borne. Dis.* **11**, 101464 (2020).

25. Cat, J. *et al*. Influence of the Spatial Heterogeneity in Tick Abundance in the Modeling of the Seasonal Activity of *Ixodes Ricinus* Nymphs in Western Europe. *Exp. Appl. Acarol.* **71**(2), 115–130. https://doi.org/10.1007/s10493-016-0099-1 (2017).

26. Nyrhilä, S., Sormunen, J. J., Mäkelä, S. & Sippola, E. One out of ten: Low sampling efficiency of cloth dragging challenges abundance estimates of questing ticks. *Exp. Appl. Acarol.* **82**, 571–585 (2020).

27. Bord, S., Druilhet, P., Gasqui, P., Abrial, D. & Vourc'h, G. Bayesian estimation of abundance based on removal sampling under weak assumption of closed population with catchability depending on environmental conditions. Application to tick abundance. *Ecol. Modell.* **274**, 72–79 (2014).

28. Santé Publique France. Borréliose de lyme. (2021). https://www.santepubliquefrance.fr/maladies-et-traumatismes/maladies-a-trans-mission-vectorielle/borreliose-de-lyme/donnees/#tabs. (Accessed 20th July 2021)

29. Joly, D. *et al*. Les types de climats en France, une construction spatiale. *CyberGeo* **2010**, 2–22 (2010).

30. Lebert, I. *et al*. Distribution of ticks, tick-borne pathogens and the associated local environmental factors including small mammals and livestock, in two French agricultural sites: The OSCAR database. *Biodivers. Data J.* **8**, 1–33 (2020).
31. Vassallo, M., Pichon, B., Cabaret, J., Figureau, C. & Pérez-Eid, C. Methodology for sampling questing nymphs of *Ixodes ricinus* (Acari: Ixodidae), the principal vector of Lyme disease in Europe. *J. Med. Entomol.* **37**, 335–339 (2000).

32. Pérez-Eid, C. *Les Tiques: Identification, Biologie, Importance Médicale et Vétérinaire* (Tec & Doc Lavoisier, 2007).

33. Hijmans, R. J. geosphere: Spherical Trigonometry. R package version 1.5-10. (2019).

34. van Buuren, S. & Groothuis-Oudshoorn, K. mice: Multivariate imputation by chained equations in R. *J. Stat. Softw.* **45**, 1–67 (2011).

35. QGIS Development Team. QGIS Geographic Information System. (2021).

36. Jpl, N. A. S. A. NASA shuttle radar topography mission global 1 arc second. NASA EOSDIS Land Process. DAAC. https://doi.org/10.5067/MEaSURES/SRTM/SRTMGL1.003 (2013).

37. FAO. Elements of Topography. in *Irrigation Water Management: Training Manual No. 1—Introduction to Irrigation* (Food and Agriculture Organization of the United Nations, 1985).

38. Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. *Int. J. Climatol.* **37**, 4302–4315 (2017).

39. European Environmental Agency. Corine Land Cover (CLC) 2018, Version 2020_20u1. https://land.copernicus.eu/pan-european/corne-land-cover/clc2018 (2019). (Accessed 1 Aug 2021)

40. L'Institut national de l'information géographique et forestière. BD Forêt version 2. https://inventaire-forestier.ign.fr/spip.php?article646 (2019). (Accessed 1 Aug 2021)

41. Shannon, C. E. A mathematical theory of communication. *Bell Syst. Tech. J.* **27**, 379–423 (1948).

42. Panagos, P., Van Liedekerke, M., Jones, A. & Montanarella, L. European Soil Data Centre: Response to European policy support and public data requirements. *Land Use Policy* **29**, 329–338 (2012).

43. R Core Team. A Language and Environment for Statistical Computing. *R Foundation for Statistical Computing 2*. https://www.R-project.org (2019). (Accessed 1 Aug 2021)

44. Brooks, M. E. et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. *R J.* **9**, 378–400 (2017).

45. Lé, S., Josse, J. & Francois, H. FactoMineR: A package for multivariate analysis. *J. Stat. Softw.* **25**, 1–18 (2008).

46. Kassambara, A. & Mundt, F. factoextra: Extract and Visualize the Results of Multivariate Data Analyses. (2020).

47. Barton, K. MuMln: Multi-Model Inference. (2020).

48. Hartig, F. DHARMA: Residual Diagnostics for Hierarchical (Multi-Level/Mixed) Regression Models. (2020).

49. Feng, C., Li, L. & Sadeghpour, A. A comparison of residual diagnosis tools for diagnosing regression models for count data. *BMC Med. Res. Methodol.* **20**, 1–21 (2020).

50. Fox, J. & Weisberg, S. *An R Companion to Applied Regression*. (Thousand Oaks, CA, 2019).

51. Mejlon, H. A. Diel activity of *Ixodes ricinus* Acari: Ixodidae at two locations near Stockholm, Sweden. *Exp. Appl. Acarol.* **21**, 247–256 (1997).

52. Mejlon, H. A. & Jaenson, T. G. T. Questing behaviour of *Ixodes ricinus* ticks (Acari: Ixodidae). *Exp. Appl. Acarol.* **21**, 747–754 (1997).

53. Kjellander, P. L., Aronsson, M., Bergvall, U. A. & Carrasco, J. L. Validating a common tick survey method: CLoth-dragging and line transects. *Exp. Appl. Acarol.* **83**, 131–146 (2021).

54. Gray, J. S., Kahl, O., Lane, R. S., Levin, M. L. & Tsao, J. I. Diapause in ticks of the medically important *Ixodes ricinus* species complex. *Ticks Tick. Borne. Dis.* **7**, 992–1003 (2016).

55. Diuk-Wasser, M. et al. Spatiotemporal patterns of host-seeking *Ixodes scapularis* nymphs (Acari: Ixodidae) in the United States. *J. Med. Entomol.* **43**, 166–176 (2006).

56. Cayol, C., Koskela, E., Mappes, T., Siukkola, A. & Kallio, E. R. Temporal dynamics of the tick *Ixodes ricinus* in northern Europe: Epidemiological implications. *Parasit. Vectors* **10**, 1–11 (2017).

57. Tomkins, J. L., Aungier, J., Hazel, W. & Gilbert, L. Towards an evolutionary understanding of questing behaviour in the tick *Ixodes ricinus*. *PLoS ONE* **9**, e110028 (2014).

58. Gilbert, L., Aungier, J. & Tomkins, J. L. Climate of origin affects tick (*Ixodes ricinus*) host-seeking behavior in response to temperature: Implications for resilience to climate change?. *Ecol. Evol.* **4**, 1186–1198 (2014).

59. 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. *Parasit. Vectors* **13**, 1–12 (2020).

60. Randolph, S. E. Tick ecology: Processes and patterns behind the epidemiological risk posed by ixodid ticks as vectors. *Parasitology* **129**, 37–65 (2004).

61. Li, S., Heyman, P., Cochez, C., Simons, L. & Vanwambeke, S. O. A multi-level analysis of the relationship between environmental factors and questing, *Ixodes ricinus* dynamics in Belgium. *Parasit. Vectors* **5**, 1–11 (2012).

62. Tack, W. et al. Local habitat and landscape affect *Ixodes ricinus* tick abundances in forests on poor, sandy soils. *For. Ecol. Manage.* **265**, 30–36 (2012).

63. Perez, G. et al. Relationships between landscape structure and the prevalence of two tick-borne infectious agents, *Anaplasma phagocytophilum* and *Borrelia burgdorferi sensu lato*, in small mammal communities. *Landsc. Ecol.* **35**, 435–451 (2020).

64. Halos, L. et al. Ecological factors characterizing the prevalence of bacterial tick-borne pathogens in *Ixodes ricinus* ticks in pastures and woodlands. *Appl. Environ. Microbiol.* **76**, 4413–4420 (2010).

65. De Keukeleire, M. et al. Scouts, forests, and ticks: Impact of landscapes on human-tick contacts. *Ticks Tick. Borne. Dis.* **6**, 636–644 (2015).

66. Ogrzewalska, M., Uezu, A., Jenkins, C. N. & Labruna, M. B. Effect of forest fragmentation on tick infestations of birds and tick infection rates by Rickettsia in the Atlantic Forest of Brazil. *EcoHealth* **8**, 320–331 (2011).

## Acknowledgements

Most tick observatories were supported by the CC-EID and Climatick projects: "Adaptation of Agriculture and Forests to Climate Change" metaprogramme of the French National Research Institute for Agriculture, Food and Environment (INRAE). The observatory of Les Bordes was installed and equipped by the Centre Val-de-Loire region, the Loiret and the French National Forest Office. We also would like to thank the TEMPO French Network (https://tempo.pheno.fr/) for providing resources to support this study. The Les Bordes site belongs to the French national research infrastructure, ANAEE-F (http://www.anaee-france.fr/fr/), and is also included in the SOERE TEMPO (https://tempo.pheno.fr/). The Les Bordes site is also part of the framework of the ZAL (LTSER Zone Atelier Loire) and the GIS Coop network (https://www6.inra.fr/giscoop/), which is supported by the French Ministry for Agriculture and Food. This study was partially funded by EU grant 874850 MOOD and is catalogued as MOOD 040. The contents of this publication are the sole responsibility of the authors and do not necessarily reflect the views of the European Commission. Finally, PW is grateful for a grant from the Association Nationale de la Recherche et de la Technologie, bourse CIFRE (n° 2019/0393).
**Author contributions**

Conceptualization and design of field observations: S.B., K.C.M.; acquisition of data: P.W., M.J., S.B., A.A., L.B., A.C., J.F.C., N.D., T.H., F.H., N.K., I.L., A.Md., A.Mr., S.M., O.P., T.P., V.P., M.R.M., M.V.T., H.V., N.C., G.V., K.C.M.; statistical analysis and interpretation of data: P.W., M.J., S.B., K.C.M.; drafting the manuscript: P.W., M.J., S.B., G.V., K.C.M.; all authors revised and approved the final manuscript.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-022-11479-z.

**Correspondence** and requests for materials should be addressed to K.C.-M.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022