The park (45° 8′ 45″, 7° 36′ 2″, Figure 1) spanning 6,571 hectares at an average elevation of 386 m above sea level, with an elevation gradient of 269 m, is characterized by a temperate lowland climate and is enclosed by fencing (27). The park predominantly features deciduous forests and grasslands, which are managed as hay meadows. The park attracts approximately 2,000 visitors daily and maintains consistent wildlife management throughout the year. While a few horse farms and cultivated plots are situated within the park, access by other domestic animals is restricted (including pets). Data on the density of wild ungulates within the park have been recently established via camera trapping by the European Wildlife Observatory (27).

The camera trap images analyzed in this study originated from a separate field study conducted by Ferroglio et al. (23). This study utilized 14 sampling points, evenly distributed between open (hay meadows) and closed (deciduous forests) habitats. Camera traps at each sampling point were operated continuously for a year, from August 2020 to August 2021; for detailed information on the deployment of these traps, refer to Ferroglio et al. (23). Alongside the camera trapping, ticks were systematically collected every 2 weeks throughout the entire study period using dragging transects. As described by Ferroglio et al. (23), we implemented a 1 m² cloth, dragged to cover a 10 m² surface in front of the camera trap and a 26 m circle around it. The cloth was repeatedly checked to collect ticks, which were stored in 70% EtOH for further identification (time of storage: 1 to 2 months).

Alongside fieldwork, ticks were identified using dichotomous keys (28–30), washed to remove any EtOH residual, which would inhibit the polymerase chain reactions, and stored at −20°C for further analysis (performed after the end of fieldwork). To optimize the effort, the genomic DNA was extracted from ticks grouped into uniform pools based on specific criteria: dragging transect, sampling point, repetition, species, developmental stage, and sex for adult ticks. The extraction was performed on the entire body of the ticks in the pool, using a blackPREP Tick DNA/RNA kit (Biosense, Italy) according to the manufacturer’s instructions. Subsequent PCR tests were conducted on these samples to detect tick-borne pathogens, namely B. divergens/capreoli, B. microti-like, Theileria capreoli, A. phagocytophilum, B. burgdorferi s.l., and SFG Rickettsia. Table 1 illustrates the references and targeted genes of the primers implemented to detect pathogens’ DNA, while specific protocols are detailed in Supplementary Table S1. Our PCR protocol did not differentiate B. capreoli from B. divergens, as it did not target the specific region containing single-nucleotide polymorphisms which is typically used to differentiate the two species (31). However, given the specific features of the study area (the absence of domestic ruminants and the surrounding fence) and the documented extensive circulation of B. capreoli (6), we considered positive samples to represent B. capreoli alone and included them in the statistical analysis. Originally designed to detect B. microti, primers from Persing et al. (32) encompass the whole B. microti-like group, including vulpes-like and Munich-like clades (8), which are of interest both from a zoonotic point of view and for the wildlife species targeted in the present study, specifically the red fox. Primers for T. capreoli were designed in the current study to be species-specific. All primers for Piroplasmid species were intended to avoid the interference of co-infections with other protozoan microorganisms. Similarly, when targeting Anaplasmataceae, we selected specific primers for A. phagocytophilum, the species with the highest public health relevance, to avoid the interference of symbiotic bacteria such as Candidatus Midichloria spp. On the other hand, primers targeting B. burgdorferi s.l. and SFG Rickettsia encompass the whole group of microorganisms (e.g., B. afzelii and B. lusitaniae in the first case, and R. monacensis and R. conorii in the second), which are all relevant for public health as zoonotic pathogens and are not commonly hosted by the wild species targeted in the study (16). All PCR tests included a confirmed positive control for the target pathogen and a no-template negative control. All standard measures were taken to minimize the risk of contamination. Amplicons were analyzed by agarose gel electrophoresis (2%) and visualized by staining with GelRed Nucleic Acid Gel Stain (VWR International Milano, Italy).

For each pool and pathogen, we recorded the binary outcome of the PCR (positive-negative) and estimated the pathogen prevalence within each pool using the package PoolTestR (33) for RStudio (34). This package provides Bayesian estimates of prevalence along with 95% credibility intervals (Cr.I.), based on the number of ticks in each pool and the test outcome (33). Differences in prevalence among tick species and among developmental stages were explored through the chi-squared test or Fisher’s test, depending on the distribution matrix (35).

For each repetition and sampling point, we computed the total head count of individuals passing by the camera trap. Each instance of an animal exiting and re-entering the camera’s field view was treated as a new individual, as individual recognition was not possible. All age classes were included in the analysis. We extracted data for red deer, roe deer, fallow deer (Dama dama), wild boar, red fox, European badger (Meles meles), pine marten (Martes martes), and beech marten (Martes foina), as CT deployment was not sensitive for distinguishing animal species of smaller size, such as rodents and birds, or for detecting their presence in the whole field of view of the camera trap. We grouped data for mesocarnivores and wild ruminants.

Our database considered, as the response variable, the binary outcome of the PCR test (positive-negative). Explanatory variables involved were as follows:

Vegetation and environmental parameters were not included in the analysis for several reasons. First, the sampling points were uniformly distributed across vegetation types, ensuring homogeneity in coverage. Second, the habitat within the park exhibited overall uniformity due to its limited spatial extent and minimal elevation gradient. Finally, the influence of vegetation on the presence of tick-borne pathogens was likely indirect, as vegetation primarily affects reservoir host distribution or vector abundance rather than directly influencing the pathogen presence.

We implemented a General Linear Mixed Model with a binomial family using the lme4 package in R studio (36). We considered the sampling point as the random variable, and the variance was weighted on the season in which ticks were collected, to account for seasonality in tick and host population.

To test how well the model performed and, consequently, how useful camera trap data could be to predict the presence of tick-borne pathogens, we evaluated three parameters: (i) conditional and marginal R², to test how much the model was satisfactory in explaining the variance; (ii) accuracy; and (iii) AUC (area under curve, for which threshold interpretation was presented by studies such as Çorbacıoğlu et al. (37)), tested by splitting the database into train (70% of data) and test (30% of data) datasets.

We modeled (i) all tick pools, regardless of species and developmental stages (global model), and (ii) I. ricinus nymph pools, as particularly relevant in terms of zoonotic risk (I. ricinus nymph model) (38, 39). We ultimately created models only for the four most prevalent pathogens (B. capreoli, B. microti-like, B. burgdorferi s.l., and SFG Rickettsia).

We analyzed a total of 2,537 ticks divided into 413 pools, including 282 pools of I. ricinus, of which 112 (424 individuals) were identified in the nymphal stage (23). The two most prevalent pathogens were B. microti-like (Bayesian prevalence 13.12% with Cr.I. 11.05–15.28%) and SFG Rickettsia (10.79%, Cr.I. 9.14–12.64%), followed by B. capreoli (5.47%, Cr.I. 4.2–6.78%) and B. burgdorferi s.l. (2.57%, Cr.I. 1.83–3.45%). Finally, the prevalence was recorded below 1% for T. capreoli and A. phagocytophilum, respectively (0.23%, Cr.I. 0.02–0.71, and 0.34%, Cr.I. 0.14–0.62%).

Ixodes ricinus and H. punctata recorded at least one positive pool for each pathogen (Table 2). B. capreoli, B. microti-like, and SFG Rickettsia were found in samples from every tick species tested, including D. reticulatus, H. concinna, I. hexagonus, and the R. sanguineus complex [R. sanguineus sensu stricto, R. pusillus, and R. turanicus according to the keys implemented in this study (28, 30)]. Additionally, T. capreoli was detected in H. concinna, H. punctata, and I. hexagonus, while B. burgdorferi s.l. was identified in H. concinna and H. punctata. According to the chi-squared test and Fisher’s test, the prevalence rates of only B. burgdorferi s.l. were significantly different (p < 0.05) among developmental stages and those of B. microti-like alone did significantly vary (p < 0.05) among the tick species.

While the percentages of the positive pools of SFG Rickettsia and A. phagocytophilum were persistent throughout the study period, we observed a peak during the summer months for B. burgdorferi s.l. and the two Babesia species (Figure 2). However, some differences may be spotted while decomposing the trend according to developmental stages. For B. burgdorferi, positive larvae peaked in June, while among the other stages, positive pools were more uniformly distributed in warmer months (May–September). For B. microti-like, they peaked in all stages in June and July. Finally, for B. divergens, the peak in nymphs and adults occurred earlier in the year (from May to July) than in larvae (July). Theileria capreoli, which was almost only recorded in larvae, showed higher positivity later in the year.

Based on the best performance of the models, the number of ticks in the pool and the wild boar and mesocarnivore head counts were included in all models (Table 3), while deer species head counts were included in the models for B. capreoli, SFG Rickettsia and in the global model of B. burgdorferi s.l.

Overall, the coefficients indicated a similar effect (either dilutive or additive) of the same host species across both the global model and the I. ricinus nymph model, despite statistical significance (Figure 3). Wild ruminants and wild boar exhibited a dilutive effect (negative coefficient) on the prevalence of all assessed pathogens. In contrast, the influence of mesocarnivores varied between the global model and the I. ricinus nymph model. Overall, mesocarnivores demonstrated a statistically significant dilutive effect for B. divergens, while showing a statistically significant additive effect for B. burgdorferi s.l. and B. microti-like. Furthermore, the total number of ticks in the pool had a positive association with the presence of pathogens within the pool, except in the global model of B. burgdorferi s.l., where a detractive impact was detected.

The model performance, evaluated using R², accuracy, and AUC, was the highest for the models focused on I. ricinus nymphs, apart from those predicting B. burgdorferi s.l. Notably, SFG Rickettsia exhibited the lowest accuracy and AUC, but the highest R² among the models analyzed. The model for B. capreoli demonstrated the highest accuracy, while the B. burgdorferi model achieved the highest AUC, exceeding 0.7, marking it as the only model to surpass this threshold of good predictivity (Table 4).