Study region
The study region, Swabian Alb in Baden-Württemberg, is part of the German Biodiversity Exploratories. Spanning 420 km2, it features diverse landscapes and land-use types. We selected 25 of 50 forest plots covering the full range of management intensities [24]. In these forests, land-use is defined by a silvicultural management intensity index (SMI), which combines the three main components of a forest stand: stand age, stand growth, and main tree species [25] (Additional file 1: Supplementary Material Table S1) [26]. For details on the study region and plot selection, see Weilage et al. [17].
Large mammal camera trapping
Two camera traps per plot recorded wildlife from spring 2023 to spring 2024 (9,609 camera days) following a standardized protocol [27, 28]. Images were processed with Agouti [29], and artificial intelligence-based species identification was manually verified. Human activity was excluded. A detailed description of the method can be found in Weilage et al. [17].
Small mammal camera trapping
Custom-built small mammal camera traps [30] recorded wildlife during tick collection periods in 2023–2024. Two traps per 300-m2 area operated 7 days per season (984 camera days), capturing 15 s (s) videos with a 30 s delay between triggers. Species identification followed field guides [31]. For details, see Weilage et al. [17].
Tick collection, morphometric species identification, and age measurement
Ticks were collected in spring (May), summer (August), and autumn (October) of 2023, as well as in spring (May) of 2024 using the flagging method on 300 square meters per plot. They were identified to their developmental stage and species using morphological keys [12, 32] under a stereomicroscope (Motic® SMZ–171, Motic Europe, S.L.U., Barcelona, Spain). Ticks were stored in 50-ml falcons (sterile, cat. no. AN79.1, Carl Roth GmbH + Co. KG, Karlsruhe, Germany) with a blade of grass at +7 °C in the refrigerator until morphometric age measurements, which were conducted within a maximum of 1 week after collection.
The same specimens were also used in a previous study [17], which provides detailed information on the study design and sampling protocol but focused on landscape and host-related drivers of tick density and B. burgdorferi s.l. prevalence, while the present study investigates the influence of tick age and land-use on B. burgdorferi s.l.
To determine the morphometric age of ticks, body length (BL) and width (BW), as well as scutum length (SL) and width (SW) of live I. ricinus nymphs were first measured using an already established method [19, 20] with a Keyence VHX-900F digital microscope (Itasca, IL, USA) at 200× magnification and then incorporated into a specific formula by Uspensky et al. [19] resulting in the alloscutal/scutal index. The formula can be found in Springer et al. [20]. Previous studies have classified the values of this index into eight distinct subgroups [19], which were further assigned to three overarching categories: old (IV), middle-aged (III), and young (II) (Additional file 1: Supplementary Material Table S2) [18]. Afterwards, ticks were stored at −20 °C until further examination.
DNA extraction of ticks and molecular analyses for Borrelia spp.
DNA extraction and molecular analysis via real-time polymerase chain reaction (PCR) followed established protocols and are described in detail in Weilage et al. [17, 33]. DNA was extracted using the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany), and the PCR targeted the p41 flagellin gene of Borrelia burgdorferi sensu lato. The analytical sensitivity of the assay was ten genome copies per 10 µl, as determined by a validated standard dilution series [33]. Additionally, larvae were homogenized in 200 µl of phosphate-buffered saline (PBS), while adults were homogenized in 300 µl of PBS, as not mentioned in the original text.
To identify Borrelia genospecies and STs in samples with a cycle threshold (CT) value of ≤ 41, MLST was performed, targeting eight housekeeping genes nifS, pyrG, clpX, pepX, uvrA, rplB, clpA, and recG, using the GoTaq® G2 Hot Start Green Master Mix (Promega GmbH, Walldorf, Germany) with slight modifications [34] to the original protocol [11, 34]. Additionally, the previously modified protocol was further adapted through several procedural adjustments, as detailed in Additional file 1: Supplementary Material Table S3. PCR products were visualized using the ultraviolet products (UVP) GelSolo Simplified ultraviolet (UV) Gel Documentation System (Analytik Jena, Germany). Sequencing was performed using the forward and reverse primers specific to each gene from the previous amplification (Eurofins Genomics, Ebersberg, Germany). Sequences were analyzed with Bionumerics software (version 7.6.1; Applied Maths, Austin, TX, USA) and compared with GenBank sequences through BLASTn (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Aligned sequences were assigned to allelic ST profiles in the MLST database (http://pubmlst.org/Borrelia) with newly identified STs submitted to the curators of the pubMLST platform. For samples that did not amplify PCR products for all housekeeping genes, genospecies identification was based on at least one of the following genes: clpA, pepX, recG, and rplB.
Statistical analysis
Confidence intervals (CIs; 95% CI) were calculated using the Clopper and Pearson method in GraphPad software (Graph Pad Software, SanDiego, CA). Statistical analysis was performed using R-software (version 4.1.2. for Windows; RStudio, Boston, MA) with the lme4 package [35]. A generalized linear model (GLM) with a quasibinomial error distribution and a logit link function was developed to evaluate tick age in dependence of (i) SMI (independent continuous variable), (ii) season (independent categorical variable: spring, summer, or autumn), and (iii) Borrelia infection status (independent binary variable; prevalence: positive: 1, negative: 0). This test was conducted for three broader I. ricinus nymph age groups (young, middle-aged, and old) and for eight more defined age groups according to Uspensky [19]. Furthermore, a Mann–Whitney U test was conducted to evaluate whether the mean morphometric age of Borrelia-infected ticks differs significantly from that of uninfected individuals. Borrelia sequence type diversity was log-transformed (log(x + 1)) to meet the assumptions of linear modeling.
In addition, a proportional odds logistic regression model (polr) was applied, using the MASS package in R, to examine the relationship between morphometric age groups (dependent variable) and variables such as SMI, canopy openness, mean tree diameter at breast height (dbh), tree species richness, shrub cover, and dead wood volume (dwv) [26, 36,37,38,39] (see Table 1). The ordinal regression model with included random effect of the plot (ordinal package using the clmm command) proved the random effect of the plot to be marginal (random effect variance: 2.398 × 10−9); therefore, it was not included in the final model.
For analyzing B. burgdorferi s.l. prevalence in I. ricinus nymphs in relation to tick age, season, and land-use (SMI), a generalized linear mixed model (GLMM) with binominal error distribution and link logit function was implemented using the lme4 package in R [35]. The model was designed to assess whether (i) tick age (independent categorical variable), (ii) season (independent nominal variable: spring, summer and autumn), and (iii) SMI (independent continuous variable) have an impact on Borrelia infection status (dependent binary variable; Borrelia abundance; positive: 1, negative: 0). The test was conducted twice, once for three broader nymph age groups (young, middle-aged, and old) and a second time for eight more defined age groups according to Uspensky [19].
Post hoc Tukey’s tests were performed to assess pairwise differences between seasons and tick age groups, adjusting for multiple comparisons. The significance threshold was set at P ≤ 0.05.
The Borrelia concentration (cp/µl) was determined by interpolating CT values against a standard curve based on known concentrations, akin to probit analysis. Unlike probit analysis, the standard dilution series was established once and later used for comparison with all positive samples, rather than being included in each PCR reaction. We used a bacterial culture of Borrelia afzelii, which was serially diluted in tenfold steps from an initial concentration of 2.0 × 105 cells/µl down to 2.0 × 10−1 cells/µl, resulting in the following dilution levels: 105, 104, 103, 102, 101, 10⁰, and 10−1.
Furthermore, a linear regression model (LM) was used to investigate how the Borrelia genospecies diversity is influenced by various predictors. The model was intended to evaluate whether (i) the Shannon diversity of specific host species [independent continuous variables: H_Pre (diversity of predator hosts), H_small (diversity of small mammal hosts), and H_ large (diversity of larger mammal hosts)], (ii) the relative abundance indeces (RAI) of specific host species (RAI_Pre (relative abundance index of predators), RAI_small (of small mammals), and RAI_large (of large herbivors), independent continuous variables), and (iii) total mammal host species richness (S_all, independent continuous variable) significantly affect the diversity of Borrelia genospecies at distinct locations.
To further explore factors influencing Borrelia ST diversity in ticks, another LM evaluated the impact of the same set of predictors on ST diversity.
RAIs for each species were computed as events per 100 camera days using the R script outlined in Rovero and Zimmermann [40]. The Shannon index was then derived from RAIs.
We selected and averaged the best-fitting models (ΔAICc < 2) using the MuMln package [41]. We report conditional averaged model results. All included candidate models are in the supplement.