The study was conducted in the tree diversity experiment ORPHEE, located in southwestern France, in the Landes de Gascogne forest. The region is characterized by a mild oceanic climate, with mean annual temperature of 12.5 °C and mean total precipitations of ca. 950 mm. The ORPHEE experiment, which belongs to the global network of tree diversity experiments TreeDivNet (Verheyen et al., 2016; Paquette et al., 2018), was established in 2008 on a 12 ha clearcut of maritime pine stands, on a sandy podzolic soil. It is composed of eight blocks of 32 randomly distributed plots corresponding to all the possible combinations of one to five tree species. Five native tree species were planted in monocultures (European birch: Betula pendula, ROTH; Pedunculate oak: Quercus robur, L.; Pyrenean oak: Quercus pyrenaica, WILLD.; Holm oak: Quercus ilex, L.; and Maritime pine: Pinus pinaster, AIT.) as well as all combinations with two, three, four and five tree species. Each plot is composed of 10 × 10 trees, planted 2 meters apart, resulting in 100 trees per plot, with a plot area of 400 m². Tree species are planted according to a systematic alternate pattern (see details in the study by Castagneyrol et al., 2013; Damien et al., 2016). Four out of eight blocks are irrigated since 2015. Blocks are sprinkled at night with 42 m³ of water (i.e., 3 mm precipitation, block area of 1.28 ha), daily from May to October, which significantly increases water availability as compared to non-irrigated blocks (Maxwell et al., 2020). The site has been fenced, which may have reduced the presence or passage of large mammal in the plots.

In our study, we focused on 14 different plots with the following four tree species: Q. robur, Q. ilex, P. pinaster, and B. pendula, planted in monocultures (n = 4), two-species mixtures (n = 6), and three-species tree mixtures (n = 4). All plot compositions were replicated eight times, four in irrigated blocks and four in non-irrigated blocks, for a total of 112 plots.

Free-living ticks, being larvae, nymphs or adults, were collected in June and July 2020 on days without rain during the peak of tick activity season, using the dragflag technique (Macleod, 1932; Vassallo et al., 2000). The total number of ticks collected along four parallel 20 meter long transects in each plot (80 m² sampled area) were denoted N_{i,j,k} where i ∈ [1;14] represents the plot composition, j ∈ [1;8] represents the replication of plot composition and k ∈ [1;4] represents the k^th parallel long transect. To take into account the effect of different understorey diversity or composition on the sampling rate (Bord et al., 2014), i.e., the proportion of ticks collected by the drag in the questing tick population, each transect was passed four times, corresponding to two return trips, with a 1 m × 1 m white flannel cloth of dragged on the ground at a constant speed (walking pace). Ticks were counted after each successive sampling, stage determined and kept in a 70% ethanol tube at the end of each pass by transect. Abundance of adults, nymphs and larvae was confirmed after binocular examination in the laboratory.

Three plots could not be sampled due to the presence of dense gorse (Ulex europaeus) that would have torn the flannel. A total of 109 plots were therefore sampled, corresponding to 436 transects and 1744 sampling passes. Tick species were determined after binocular examination using the systematic key from Estrada-Peña (2017). Too few adults were collected (n = 35 including 29 I. ricinus). We therefore decided to focus only on I. ricinus nymphs and larvae.

All I. ricinus collected nymphs (n = 173) were isolated for pathogen identification. We used a microfluidic PCR technique that provides information on the presence of 35 different pathogens of the following genius: Borrelia, Anaplasma, Rickettsia, Francisella, Bartonella, Babesia and Coxiella potentially present in tick nymphs (developed by Michelet et al., 2014).

Tick nymphs were rinsed twice with distilled water. Prior to extraction, ticks were crushed and processed in a homogeniser at 5,500 rpm for 20 s with 100 μl of elution buffer BE, 40 μl of lysis buffer MG and 10 μl of liquid proteinase K and six steel balls. DNA was extracted using a Nucleospin DNA insect kit (Macherey-Nagel, Germany), and DNA from each sample was stored at −20°C in 50 μl of rehydration solution until next step.

DNA pre-amplification reactions was performed with 2.5 μl of PreAmp Master Mix according to the manufacturer’s instructions (Fluidigm, CA, United States). Prior to experiments, primers targeting all pathogens were pooled in equal volumes (200 nM each). The experiment was performed with 1 μL of PreAmp Master Mix, 1.25 μL of pooled primers mix, 1.5 μL distilled water, and 1.25 μL DNA for 5 μL final volume. The pre-amplifications were performed with one cycle at 95°C for two minutes, 14 cycles at 95°C for 15 s and 60°C for four minutes. Pre-amplified DNA was diluted 1:10 and stored at −20°C. Finally, a high-throughput real-time PCR was performed with the BioMark system (Fluidigm, United States) using the method explained in Michelet et al. (2014).

Pathogens detected with BioMark system were confirmed by conventional PCR or nested PCR on several samples presenting low Ct values. Primers targeting different genes or regions than those of the BioMark system were selected from the literature (listed Supplementary Table 1) and PCR were performed with Phusion High fidelity kit (Thermo scientific, MA, United States). The positive samples after gel migration were sent to sequencing at Eurofins (Germany) and assembled using BioEdit software (Ibis Biosciences). An online BLAST (NCBI; National Center for Biotechnology Information) was used to identify the sequenced organisms.

We analyzed four response variables for I. ricinus: larvae presence, nymph abundance and proportion of infected nymphs by Borrelia burgdorferi s.l. and proportion of infected nymphs by Rickettsia. Nymph abundance indicator was the cumulative number of nymphs collected after the four successive samplings ∑k=14N{i,j,k}l where i ∈ [1;14],j ∈ [1;8],k ∈ [1;4] and l ∈ [1;4] corresponding to the plot composition, the repetition number, the parallel long transect and the successive samplings, respectively. Larvae presence variable is binomial: it is set to 1 when at least one larva was found in the plot, and 0 if none was found. As we performed individual PCR’s, the proportion of nymphs infected by each pathogen was tested in the model as a vector including the number of pathogen-positive nymphs and pathogen-negative nymphs in the plot.

We selected three tree layer variables. First, canopy tree density, corresponding to the percentage of canopy individuals in the plot (here pine and/or birch, since oak trees are much smaller in ORPHEE, still being part of the understorey layer; Damien et al., 2016).

Second, tree composition corresponds to different types of tree mixtures (a total of nine type of mixtures). In order to obtain sufficient statistical power to test the effect of plot composition on ticks variables, individual plot compositions were combined into four composition types: (1) Pure pine (including monocultures of pines and mixtures of pines and oaks, which are still small enough to be considered part of the understorey); (2) Pure birch (including monocultures of birches and mixtures of birches and oaks); (3) Mixtures of pine and birch (including mixtures of pines and birches, and mixtures of pines, birches and oaks); (4) Oaks (including monocultures of oaks and mixtures of only oaks) (Table 1).

Third, tree species richness corresponds to the number of different tree species in the plot (including the oak species).

Further, we included understorey (understorey phytovolume) and microclimatic variables (presence/absence of irrigation and the saturation deficit of the plot) into the models.

To describe the understorey vegetation, we sampled one 1 m × 1 m quadrat in each transect passed for the ticks (i.e., four quadrats per plot). In each quadrat, we estimated the percentage of bare soil, litter and moss cover, the percentage cover of the main plant types found in ORPHEE (Corcket et al., 2020), i.e., grasses, heath, ferns and shrubs (mainly gorses), and the maximum height of each plant type. We considered not only plants directly rooted in the quadrat but also those whose aerial parts partially or completely covered the surface of the quadrat. The percentage cover of each type was transformed in eight categorical classes: O: 0%; A: 1–5%; B: 6–15%; C: 16–25%; D:26–50%, E:51–75%; F:75–99%; G: 100% [inspired by the Braun-Blanquet approach (Braun-Blanquet, 1964)]. The maximal height of each plant type was measured on the highest individual in the quadrat using a measuring pole. Finally, we aggregated the data at the plot level by averaging across the four quadrats to calculate the understorey phytovolume (cm³) according to equation (1):

Finally, temperature (T) and relative atmospheric humidity (RH) were measured in each plot N_i,j before each sampling (with Kestrel 3500 NV), in order to calculate the Saturation Deficit in millimeters of mercury, according to equation (2) (Vourc’h et al., 2016):

Saturation Deficit, a proxy for tick desiccation potential, is considered a relevant predictor of nymph foraging activity (Gray, 2008).

According to our hypotheses we were interested in testing the effects of tree layer variables, understorey vegetation variables and microclimatic variables on tick response variables. A first sorting of the explanatory variables was done with univariate tests. Variables with a p value < 0.1 were selected as candidates for the full models.

Tick nymphs, larvae and associated pathogens were analyzed separately. For nymph abundance, we built one negative-binomial generalized linear model (glmer.nb) suitable for overdispersed data (Breslow, 1990). We modeled the presence of larvae and proportion of infected nymphs with binomial generalized linear models (glm, binomial family).

Nymph abundance, presence of larvae, and proportion of infected nymphs were modeled as a function of tree composition or canopy tree density, irrigation, tree species richness, understorey phytovolume and their interactions, and saturation deficit as fixed effects. Spatial autocorrelation of plots in the block was accounted for by including Block as random factor.

All continuous predictor variables were scaled and centered prior to modeling to make their coefficients comparable (Schielzeth, 2010). We first built full models fitting the most complex model consistent with the experimental design including all fixed effect factors, interactions and random factors. Then we applied model simplification by sequentially removing one-by-one, the less significant terms to allow a non-singular fit (Barr, 2013). Once optimally simplified, we selected the best fitted model with the lowest AIC. Statistical analyses were performed using the R software version 4.0.5 (R Core Team, 2021).

Larvae were present in 56 out of the 109 sampled plots and their abundance was on average (±sd, n = 1458) 13.38 ± 0.37 per plot and 26.03 ± 0.46 in the plots where they were present. However, because larval distribution was highly aggregated due to clustering in larval nests, we analyzed only presence/absence data.

The probability of presence of the larvae significantly differed between tree compositions (Figure 3). Presence of larvae was highest in pine plots (pure and mixed with birch), intermediate in birch plots and lowest in oak plots. We also found a significantly positive effect of irrigation, with higher probability of presence in irrigated plots, irrespective of plot composition (i.e., lack of significant interaction) (Table 2 and Figure 3). The probability of larvae presence was significantly higher in irrigated pure and mixed pine plots than in non-irrigated oak plots (Figure 3).

We collect I. ricinus nymphs in 65 out of 109 plots. Nymph abundance was on average (±sd, n = 173) 1.59 ± 0.02, with a median of 1 nymph per plot and 2.66 ± 0.02 in the plots with at least one nymph was present. According to the AIC of compared models, three explanatory variables were retained as significantly influencing the abundance of tick nymphs (Table 2). Nymph abundance increased with the density of canopy trees (Figure 4) and decreased with the phytovolume of understorey vegetation (Figure 5), being significantly higher in irrigated plots, irrespective of their composition (no interaction effect). We did not find any effect of the saturation deficit or tree species richness.

The control probe was positive for 172 of 173 collected nymphs, so one DNA extraction failed. Infected nymphs were found in 40 of the 65 plots with nymphs, 72 nymphs were infected with at least one pathogen and 6 nymphs presented a co-infection with several pathogens. Among the 72 infected nymphs, the majority of associated pathogens were from the B. burgdorferi s.l. complex (n = 35), and Rickettsia spp. (n = 32, with 69% of Rickettsia helvetica and 31% of Rickettsia monacensis), the other pathogens were not found frequently enough to test for a significant relationship with covariates.

The Borrelia complex included five different species: 80% of B. lusitaniae, 11% B. afzelii, 3% B. valaisiana, 3% B. garinii and 3% B. burgdorferi s.s. No B. spielmanii, B. bissetti, and no B. miyamotoi (relapsing fever Borrelia) were found. B. lusitaniae, which was the more frequently detected species, is known to be vectored by ticks with reptiles as hosts. Two other reptile-borne pathogens were found: 10 nymphs carried Rickettsia monacensis and one carried Candidatus cryptoplasma. The proportion of nymphs infected by B. burgdorferi s.l. was significantly influenced by the richness of tree species and also by irrigation, with no interaction between these factors (Table 2). The occurrence of nymphs infected by B. burgdorferi s.l. pathogen was lower in irrigated plots and decreased with increasing tree species richness, being highest in non-irrigated pure plots and lowest in irrigated three-species mixtures (Figure 6).

None of our variables were significantly related to the proportion of ticks infected with a Rickettsia pathogen (data not shown).

Our results show that forest diversity and composition can influence tick populations across different stages and their level of infection by pathogens. The probability of presence of larvae was lowest in oak plots, intermediate in birch plots and greatest in pure and mixed pine plots, suggesting a strong tree species identity effect. The abundance of nymphs was significantly affected by canopy tree density, being greater in denser forest plots (e.g., 100% pines, 100% birches or 50% pine–50% birch mixtures), suggesting a stronger effect of forest structure than tree identity. For a tick to be present in an environment, it must (i) have been deposited by a host for nymphal and adult stages or have been laid by a female for the larval stages. (ii) survived in a favorable environment. Overstorey trees can thus have an impact on young tick population by controlling local microclimate (Renaud et al., 2011) and host presence. For example, dense, shaded forests as a result of canopy packing (Jucker et al., 2015; Williams et al., 2017) can be optimal habitats because they provide host animals with food resources while maintaining sufficient moisture for tick survival. A 19-year study in a mature oak-maple mixt forest showed a bottom-up, cascading effect of acorns yield on increasing mouse population levels, resulting in increased abundance of Ixodes scapularis (Say, 1821) nymph and associated Lyme disease risks (Ostfeld et al., 2018), and similar results were recently found with I. ricinus (Bregnard et al., 2020, 2021). High moisture is important for survival of tick larvae and nymphs (Knülle and Devine, 1972; Tukahirwa, 1976; Davey et al., 1991) and Gray (1998) explains that I. ricinus are usually found in environments with 80% atmospheric humidity. The questing activity of tick nymphs occurs higher in the understorey vegetation than that of larvae, so that they can attached to a passing animal, but this presents a greater risk of desiccation during this questing period (Schulze and Jordan, 2005). Denser forests with more trees and with a higher tree canopy packing could therefore limit this risk of desiccation and thus favoring the survival of tick nymphs. This explanation for the positive correlation between nymph abundance and proportion of canopy trees in our experimental plots is also consistent with our result of a positive effect of irrigation.

The presence of larvae and their questing activity are related to the microclimate closer to the ground (Yoder et al., 2006; Ehrmann et al., 2017). Leal et al. (2020) suggest that forest types where the most larvae are found differ between tick species. Indeed, several studies have found greater abundance of I. ricinus larvae in deciduous forests (Walker et al., 2001; Tack et al., 2012; Ruyts et al., 2018). Due to the much slower growth of these species, the plots with only oaks in ORPHEE are currently characterized by small trees of about 2 m height and thus an open canopy. Our results in deciduous plots are therefore to be distinguished from the scientific literature on mature deciduous forests. These conditions on ORPHEE were probably less favorable for the maintenance of suitable humidity in the understorey and for the visit of host animals. In contrast, pine plots characterized by dense and tall trees, may have favored larval survival by providing thick needle litter as a protective barrier against desiccation. Birch plots are in an intermediate position in terms of likelihood of larvae presence, probably because they are composed of high trees but with lower canopy closure and perhaps lower quality litter for ticks.

The proportion of nymphs vectoring Borrelia pathogens significantly decreased with increasing tree species richness, with the highest frequency in pure pine plots and the lowest in mixtures associating pine, birch and oaks. Our results thus support the finding by Ehrmann et al. (2018) that a reduction in the prevalence of tick-borne diseases is associated with greater forest diversity, confirming the beneficial role of tree species mixing for natural risk reduction (Jactel et al., 2017). Ruyts et al. (2016) did not find a tree species diversity effect on the prevalence of tick-borne pathogens but they found that the risk of Lyme borreliosis differed between different forest compositions. They showed that B. afzelii was more often found on ticks in pine plots while B. garinii and B. burgdorferi s.s. were more prevalent in oak plots. A negative effect of tree species richness on infected nymphs could be explained by an increase in ecological niches for host species, resulting in a change in the ratio of competent to incompetent reservoir hosts, with diluting effects on Borrelia (Levi et al., 2016; Ehrmann et al., 2018). Empirical multitaxonomic studies (e.g., Ampoorter et al., 2020) and meta-analyses (e.g., Stemmelen et al., 2021) have indeed shown greater faunal diversity in mixed forests than in pure forests. Additional research is now needed to inventory potential tick host animals in pure and mixed forests and better assess their level of competence for tick-borne pathogens.

Ruyts et al. (2016) suggested that the dilution effect would depend on the species of B. burgdorferi s.l. In our study, 80% of B. burgdorferi s.l. were in fact B. lusitaniae, which appears to be carried and multiplied by lizards in several countries around the Mediterranean rim (Michelis et al., 2000; Younsi et al., 2001; Sarih et al., 2003; Tomassone et al., 2017). As ectothermic animals, lizards often prefer environments that allow them significant sun exposure for active thermoregulation (Smith and Ballinger, 2001). This may explain the negative effect of irrigation on the proportion of nymphs infected by Borrelia burgdorferi s.l. pathogens, either directly by increasing soil humidity or indirectly by increasing the density and thus shading of understorey vegetation (Corcket et al., 2020). We do not have yet direct data from the ORPHEE experiment but surveys in the Aquitaine region reported the presence of three lizard species in forests: the common wall lizard [Podarcis muralis (Laurenti, 1768)], the viviparous lizard [Zootoca vivipara (Lichtentein, 1823)] and the western green lizard (Lacerta bilineata, Daudin, 1802) (Berroneau, 2015). All are known as I. ricinus potential hosts (Scali et al., 2001; Richter and Matuschka, 2006; Wu et al., 2019; Gwiazdowicz et al., 2020), and two of them (common wall and western green lizards) have already been involved in the transmission of B. lusitaniae in I. ricinus populations (Richter and Matuschka, 2006; Mendoza-Roldan et al., 2021). To date, no studies on the dilution effect of lizard hosts have been tested on B. lusitaniae; it would therefore be interesting to inventory them along a gradient of tree diversity and also to classify them according to their level of competence.

Finally, no variable explained the proportion of I. ricinus nymphs infected with Rickettsia spp. pathogens. The ecological drivers of distribution and prevalence of Rickettsia spp. remain poorly understood in the literature and even less the effect of forest composition or diversity.

Our study showed a distinct negative effect of the phytovolume of understorey vegetation on nymph abundance, independently from canopy diversity or composition. Several studies found opposite results with the abundance of nymphs positively correlated with high density of understorey vegetation, a condition considered favorable for visit by potential tick host animals (Tack et al., 2013; Ostfeld et al., 2018). However, other studies found a negative relationship between high density of understorey vegetation and nymph abundance, suggesting on the contrary that these environments would be less visited by hosts or would not be associated with good tick survival rates (Stein et al., 2008). Similarly, a denser understorey vegetation can also dilute tick infection prevalence (Ehrmann et al., 2018). A dense understorey can provide a suitable habitat for ticks but also for small mammals; this increases the likelihood of ticks finding a host, but at the same time decreases the likelihood that many ticks can make contact with the same host species. Thus, a lower average tick load and a lower co-feeding probability on hosts could result in dilution of Borrelia prevalence.

Bord et al. (2014) also highlighted that the dragflag technique, although very effective for sampling free-living ticks in forest, is less effective when the understorey is too dense, thus resulting in underestimation of tick abundance. In our study, the negative effect of understorey phytovolume on nymph abundance is, however, still present at low levels of understorey phytovolume. In addition, the phytovolume does not take into account plant species identity, an important factor for host attraction. In ORPHEE, the main understorey plant species were Molinia caerulae, Ulex europaeus, Pteridium aquilinum and some Ericaceae spp. (Corcket et al., 2020), which are not or only little palatable for many herbivorous hosts and whose density or height may limit the passage of large ones.

All Rickettsia and Borrelia pathogens detected in our study are transmissible to humans and of medical importance (Steere et al., 2004; Fournier and Raoult, 2009). R. helvetica and R. monacensis are belonging to the spotted fever group (Merhej and Raoult, 2011) and B. afzelii, B. valaisiana, B. garinii and B. burgdorferi s.s. are bacteria involved in the Lyme disease in humans (Cutler et al., 2017). Although the transmission of B. lusitaniae to humans has been demonstrated, the link with Lyme disease has yet to be investigated (Stanek and Reiter, 2011). Previous studies found B. lusitaniae associated with medical symptoms in humans after a tick bite (Collares-Pereira et al., 2004; Franca et al., 2005; Lopes de Carvalho et al., 2008). In our study, the proportion of infected nymphs did not correlate with nymph abundance, which confirms that it is important to distinguish tick abundance vs. tick-borne pathogen dynamics and to study pathogens at the species level, due to their host and home range specificity.