Material and Methods​​​​​​​

1. Study sites and data collection

Our study was performed in 12 nest box plots within a 10×15 km² area in southern Munich (47° 57'N, 11° 21'E) established in autumn 2009. Each plot was fitted with 50 nest boxes in a regular grid covering approximately 9 hectares. Breeding parameters of great tits (detailed below) were monitored for three years (2017-2019) within all plots except one plot that was monitored only in 2017. Each nest box was inspected biweekly during the breeding season (April-July) for nesting activity. During plot inspections, the number of recreationists was counted as a measure of human disturbance (detailed by Hutfluss & Dingemanse, 2019).

When hatchlings were 10-12 days old, both parents were captured in the nest box using a spring trap. Birds were given a unique, numbered band if not previously banded. They were then tested for their exploration behaviour (detailed below), weighed, measured morphologically (Moiron et al. 2019), aged (based on plumage characteristic as first year breeder or older; Dingemanse et al. 2020), and screened for ticks (protocol detailed below). Most birds were captured in May and a few in June. These captures were used to calculate breeding density, defined as the number of breeding pairs producing first clutches per hectare.

1. 2. Exploration behaviour

Exploration behaviour was measured using a novel environment test adapted to the field using established protocols (see Stuber et al. 2013). Birds were transferred to a holding box attached to a cage (61L × 39W × 40 H cm) fitted with a mesh front and three perches, representing the novel environment. The holding box was covered with a cloth bag and the subject allowed to acclimatize for one minute. After acclimatization, the focal bird was released into the novel environment without handling and recorded for two minutes with the field observer located out of sight. Videos were subsequently scored by dividing the cage into six equal sections and three floor sections as described in Stuber et al. (2013). The exploration score was calculated as the sum of movements between all sections within the first two minutes of entering the novel environment (Abbey-Lee and Dingemanse 2019; Dingemanse et al. 2020).

1. 3. Tick burden and collection

Ticks tend to concentrate around the eyes and beak of great tits (Heylen and Matthysen 2008; Fracasso et al. 2019). A regional patch examination protocol (reviewed in Lydecker et al. 2019) was developed to standardize the screening: We screened 1) around the eyes/ears and along the margin of the beak on both sides of the bird, 2) underneath the beak, 3) along the top margin of the beak, and 4) on the top of the head. We subsequently calculated the total number of ticks carried by each captured bird. In 2018-2019, all captured birds were screened; in 2017 only a sample. In 2018 and 2019, a sample of ticks were collected using fine tweezers and stored in 99% ethanol as part of another study. Those were subsequently morphologically identified to life-stage and species according to published taxonomic keys (Filippova 1977; Hillyard 1996; Estrada-Peña et al. 2014).

1. 4. Recreation pressure index

During each plot inspection, all observers recorded each recreationist seen and their specific location (detailed by Hutfluss and Dingemanse 2019). Each recreationist was connected to the first observation to avoid double counts. The probability to observe a recreationist during a plot inspection is biased by various factors (Hutfluss and Dingemanse 2019). To obtain an unbiased index of recreation pressure, the binary probability to observe recreationists was calculated using all inspections conducted from 2010-2019 (n = 3724 inspections). This probability was calculated using a binomial generalized linear mixed effects model (GLMM) fitting fixed effects for plot inspection duration, the number of observers, and starting time (in hours from sunrise) (Hutfluss and Dingemanse 2019). We fitted random intercepts for each unique combination of plot and year (termed plot-year, see Abbey-Lee et al. 2016; Araya-Ajoy et al. 2016; Araya-Ajoy and Dingemanse 2017) to acquire an average value of recreation pressure for each plot in each year, and for date of observation to control for date-specific environmental effects (Supplemental Table 1). We extracted best linear unbiased predictors (BLUPs) for each plot-year between 2017 and 2019 (n = 34) and used them in subsequent models as a recreation pressure index. The usage of BLUPs has been criticised when associated uncertainty is not taken forward; some have proposed to estimate a posterior distribution of possible BLUP values and thereby take forward uncertainty in subsequent analyses (Hadfield et al. 2010; Houslay and Wilson 2017). However, taking forward uncertainty in BLUP-values can result in biased estimates, whereas utilizing average BLUPs as fixed effects result in less precise but unbiased estimates (Dingemanse et al. 2020). Therefore, we here present the estimated effect of the average BLUP values.

1. 5. Data Preparation

The scaled mass index (Peig and Green 2009) was used as our measure of body condition. Following Peig and Green (2009), this index was calculated separately for each sex according as follows:

1     SM_i=M_i×L₀L_i^b_SMA

where, M_i and L_i are the mass and tarsus measurements of individual i respectively; L₀ is the arithmetic mean of all tarsus lengths of either all females or all males (dependent on the sex of individual i); and b_SMA is the scaling exponent calculated as:

2    b_SMA=b_OLSR

where b_OLS is the slope of linear regression of log-transformed body mass as a function of log-transformed tarsus length of either all females or all males (dependent on the sex of individual i) and R is the Pearson’s correlation coefficient of the regression. We used all birds recorded from 2010-2019 (n = 4273 records) to calculate the scaling coefficient (including any repeated measures).

We calculated individual-mean exploration scores using a linear mixed effects model (LMM) fitting exploration behaviour as the response variable (Supplemental Table 2). We fitted test sequence as a fixed effect covariate to control for sequence effects ( Dingemanse et al. 2012, 2020), and included random intercepts for individual and plot-year. The model included data from all years (2010-2019) (n = 4251) and assumed a Gaussian error distribution. We extracted an average BLUP-value per individual, which we used in subsequent analyses as a measure of an individual’s average (experience-corrected) exploration score.

1. 6. Statistical Analysis

Tick infestation of great tits was analysed in two parts using GLMMs. First, infestation probability was estimated using all recorded captures from 2017 through 2019 (n = 784 records). Second, tick burden (i.e. the absolute number of ticks observed) was modelled using only infested individuals (n = 520 records). For both response variables, the model set up was the same. Fixed effects were fitted for recreation pressure, breeding density, body condition, sex, age, and individual-mean exploration score. Sex and age were fitted because females (versus males) and older (versus younger) birds often have higher tick burdens (Heylen and Matthysen 2008; Heylen et al. 2013). Random intercepts were fitted for year, plot, plot-year, field observer, nest box, and individual. Prior to analysis, all fixed effects were mean-centred, and variance standardized, such that the statistical intercepts of our models reflected the value for the average individual in the average environmental condition.

All statistical analyses were performed in R (version 3.5.3) (R Core Team 2019). Models were run using the functions lmer and glmer from the package lme4 (Bates et al. 2015) assuming either binomial (infestation probability) or Poisson (tick burden) error distributions. Mean estimates and their 95% credible intervals (CI) were estimated based on 5,000 simulations using the sim function from the arm package (Gelman and Su 2016). Residual errors were calculated according to Nakagawa & Schielzeth, (2010). Adjusted repeatability values were calculated as the proportion of variance unexplained by the fixed effects that was explained by the focal random effect. Estimates with 95% CIs not overlapping zero were treated as statistically significant in the frequentist’s sense. Models were checked graphically for fit (Supplemental Fig. 1 & 2).

R-script for processing the data has been included to prepare all data (as described in Data Preparation) and run models (as described in Statistical Analysis).