DOI: 10.5061/dryad.0rxwdbs7t
title: Data from: Insights from a 31-year study demonstrate an inverse correlation between recreational activities and red deer fecundity, with body weight as a mediator

description: Human activity is omnipresent in our landscapes. Animals can perceive risk
 from humans similar to predation risk, which could affect their fitness.
 We assessed the influence of the relative intensity of recreational
 activities on body weight and pregnancy rates of red deer (Cervus elaphus)
 between 1985 and 2015. We hypothesized that stress, as a result of
 recreational activities, affects pregnancy rates of red deer directly and
 indirectly via a reduction in body weight. Furthermore, we expected
 non-motorized recreational activities to have a larger negative effect on
 both body weight and fecundity, compared to motorized recreational
 activities. The intensity of recreational activities was recorded through
 visual observations. We obtained pregnancy data from female red deer that
 were shot during the regular hunting season. Additionally, age and body
 weight were determined through post-mortem examination. We used two
 generalized linear mixed models (GLMM) to test the effect of different
 types of recreation on 1) pregnancy rates and 2) body weight of red deer.
 Recreation had a direct negative correlation with the fecundity of red
 deer, with body weight as a mediator as expected. Besides, we found a
 negative effect of non-motorized recreation on fecundity and body weight
 and no significant effect of motorized recreation. Our results support the
 concept of humans as an important stressor affecting wild animal
 populations at a population level and plead to regulate recreational
 activities in protected areas that are sensitive. The fear humans induce
 in large-bodied herbivores and its consequences for fitness may have
 strong implications for animal populations.
description: Study site The Veluwe area (1250
 km<sup>2</sup>) is designated as a protected area (i.e.,
 Natura 2000) in the Netherlands (52° 5′ N, 5° 48′ E). It has a temperate
 maritime climate and its geology mainly consists of ice-pushed ridges,
 fluvioglacial deposits, and wind-blown sands. The area is elevated 50-100m
 above sea level, holding a large freshwater aquifer. Dry coniferous forest
 (i.e., mainly Scots pine (<em>Pinus silvestris</em>)),
 deciduous forest (i.e., mainly oak (<em>Quercus sp.</em>) and
 beech (<em>Fagus sp.</em>)) make up the majority of the
 habitat, interspersed with dry to moist heathland, containing common
 heather (<em>Calluna vulgaris</em>) and grasses (Ten Houte de
 Lange, 1977). Apart from red deer, three other ungulate species occur on
 the Veluwe, namely: wild boar (<em>Sus scrofa</em>), fallow
 deer (<em>Dama dama</em>), and roe deer (<em>Capreolus
 capreolus</em>) (Broekhuizen et al., 2016). Hunting levels in the
 Veluwe area can be up to 9.48 ungulates.km-2.year-1 (Ramirez et al.,
 2021). Data collection and preparation
 We performed a quantitative study using a correlational research
 design to test our hypotheses. In order to do this, we collated data on
 recreational activities and red deer body weight and fecundity between
 1985 and 2015. We quantified fecundity as a pregnancy rate.  
 Recreational activities Yearly, between 1985
 and 2015, the relative intensity of recreational activities was recorded
 by a single observer (R. Bijlsma) to assess their effect on
 ground-breeding bird ecology (see Bijlsma, 2006). The recreational
 activities were visually observed in one open landscape area (30 ha) split
 into two zones: Otterlose zand (open for public; 15 ha) and Mosselse zand
 (closed for public; 15 ha). Observations were done for a mean time (±SD)
 of 17.2 (±9.3) hours (range: 6-45 hours) spread across an average (±SD) of
 8.0 (±3.9) days per year (range: 3-21 days). For the majority of the time
 (86%), observations were done between March and August (i.e., the peak of
 recreational activity in the Veluwe) (Bijlsma, 2006). We assumed that the
 recreational activities observed in the two zones were representative of
 the Veluwe area, especially to assess whether the relative recreation
 intensity (hereafter recreation intensity) significantly differed between
 various years. Recreationists were counted (with complete visual coverage)
 whilst walking without a fixed route but covering the whole area.
 Observations were divided into six different categories, namely: hikers,
 horse riders, vehicles (cars and quads), dogs (off-leash), and cyclists
 (including All-Terrain Bikes, ATBs). Using this data, we calculated the
 annual number of recreationists for each type observed per hour.
 Red deer characteristics and counts During
 the regular hunting season between 1985 and 2015, hunters culled an
 average of 82.3 (SD = 31.3) female red deer per year (total n = 488)
 between October and June in assigned locations, within the area of the
 Game Management Unit (GMU). The GMUs predetermine the number of animals to
 cull in each age class before the hunt, to reduce population numbers and
 promote tree regeneration. Age and body weight are known to affect
 fecundity (Borowik et al., 2016). During that period, data on the age and
 body weight were collected by Han ten Seldam and Dirk Lieftink for no
 specific purpose but to serve as a baseline of red deer demographic data.
 Age was defined through post-mortem examination. For red deer younger than
 2.5 years, age was determined through incisor and molar changes. The age
 of older individuals was estimated from the wear of the teeth in the lower
 jaw (Lowe, 1967). After evisceration, female red deer were weighed (n =
 261) and the presence or absence of a foetus was recorded.
 Red deer counts were performed twice (i.e., during dawn or dusk)
 in three days, with one day in between. A single count took 3 hours (from
 30 minutes before sunrise until 2.5 hours after, or from 2.5 hours before
 sunset until 30 minutes after). Every year, count sites within fixed
 locations were determined by the coordinators of the GMU, in consultation
 with the hunting keepers, based on the occurrence of red deer in the area.
 These sites were divided over the area, with an average of one count site
 per 400 ha. Each count started at a central location, where all
 participants got instructions. Two people were counting per site, a game
 warden or field expert of the GMU and an independent counter from a
 different GMU. Every team started and ended their counting simultaneously
 and used the same form to register the counts. Animals were counted
 through visual observations, using a binocular. Each count team noted the
 number of animals in different categories (i.e., stag, yearling-buck,
 hind, yearling-doe, and fawn) and the time and location of the animals
 when counted. Afterward, all teams discussed their counts with adjacent
 teams to avoid double counts. Characteristics of the
 study area Because of our correlative study design, we
 collated data of population variables, to control for their effects on
 body weight and fecundity (i.e., confounding effects). Red deer body
 weight and pregnancy rates can be influenced by the density of red deer
 (Putman, et al., 2019; Carpio et al., 2021), and the density of competitor
 ungulates (Barrios-Garcia and Ballari, 2012; Borowik et al., 2016). We did
 not control for roe deer density, as roe deer seem to be more affected by
 the presence of red deer than vice versa, possibly because roe deer are
 more selective in their foraging compared to red deer (Borkowski et al.,
 2021). Similarly, we did not control for fallow deer density, as there
 seems to be little competition between the sympatric fallow deer and red
 deer (Bartos et al., 2002). Therefore, we collated data on red deer
 density (t<sub>-3</sub>, t<sub>-2</sub>,
 t<sub>-1,</sub> and t<sub>0</sub>) and wild boar
 density estimated by hunters and game managers (Hušek et al., 2021), using
 a yearly census according to a fixed protocol (see Supplementary Materials
 I) (Game Management Unit Gelderland, 2019). Census of red deer took place
 in March and April at count sites (1 per 400 ha) within fixed subareas,
 when animals got active and vegetation cover did not limit visibility.
 Wild boar were counted between May and June at fixed locations (1 per 200
 ha) within the area of the GMU, using bait. Each count took three hours
 (from 2.5 hours before sunset until 30 minutes after). Bait was used to
 attract the wild boar to the counting sites and binoculars and infrared
 cameras were used to count the animals. The team composition and working
 method used for counting wild boar was identical to the aforementioned
 method for red deer. In the form used for counting wild boar, each team
 noted the date, start- and finishing time, and the number of animals of
 different age and sex classes (i.e., boar, sow, pig of the sounder,
 piglet). If the group of animals left the count site, the direction and
 time they left were noted to avoid a double count by adjacent teams.
 Later, all data of the counts was registered by the coordinators and the
 total number of individuals per year was recorded in a database.
 Furthermore, we collated data on the densities of semi-domesticated herds
 of Sayaguesa and Scottish Highland cattle that were present and monitored
 in the area since 2002. Moreover, we collated culling data of red deer
 from the Game Management Unit (GMU), to account for their effect on
 chronic stress (Vilela et al., 2020) and red deer density (Stewart et al.,
 2005; Carpio et al., 2021), and therefore bodyweight and fecundity (Putman
 et al., 2019; Bötsch et al., 2020). Additionally, we
 collated environmental variables that are known to affect body weight and
 fecundity via food availability (Borowik et al., 2016), such as habitat
 availability, the presence of supplementary feeding sites, and mean annual
 temperature and precipitation (Rodriguez-Hidalgo et al., 2010). We
 assessed changes in habitat availability (Stankowich, 2008) using ArcMap
 (10.7.1) and LGN1 – LGN7 Landsat images (WUR-Alterra, 1980-2012). We
 quantified the available habitat for red deer by including habitat types
 related to foraging (i.e., grasslands, deciduous forest, coniferous
 forest, mixed forest, and heath; Gebert and Verheyden-Tixier, 2001). In
 addition, we recorded the presence and absence of supplementary feeding
 sites in the area. Nevertheless, even though chronic stress can reduce the
 intake of food (Zanette et al., 2013) and therefore affect the body weight
 and fecundity of animals, our study design did not allow us to take this
 into account. Moreover, we collected mean annual temperatures and
 precipitation from the nearest weather station (i.e., Deelen at 8.4 km
 from the study site), as it can affect energy expenditure, body weight,
 and the production of biomass (i.e., food) (Rodriguez-Hidalgo et al.,
 2010). Data analysis Data has been
 explored in SPSS (IBM SPSS statistics 28) following the protocol by Zuur
 et al. (2010). Because the different types of recreation were strongly
 correlated, we used a Principal Component Analysis with a varimax rotation
 to extract two components (96.9% of the total variance) with an eigenvalue
 larger than 1. We characterized the first component as a ‘non-motorized’
 axis, due to its strong correlation with the number of dogs, hikers,
 cyclists (including ATBs), and horse riders. The second component was
 characterized as a ‘motorized’ axis based on its strong correlation with
 the number of vehicles (including quads). Since the vehicle data contained
 a lot of zeros (60.9%), we transformed it into a presence-absence
 variable. We took into account the following control
 variables: age, available habitat, presence of feeding sites, red deer
 density (t<sub>-3</sub>, t<sub>-2</sub>,
 t<sub>-1,</sub> and t<sub>0</sub>), wild boar
 density, cattle density, mean annual temperature and precipitation and the
 annual number of red deer culled (i.e., shot). The density of red deer and
 densities of wild boar and cattle were square-root transformed to reduce
 skewness. The continuous independent variables were standardized and
 assessed for multicollinearity (r &gt; 0.7; Dormann et al., 2013).
 Correlation between categorical and continuous variables was visually
 assessed through overlap in boxplots. Multicollinearity was assumed in the
 absence of overlap and therefore one of the variables examined was not
 included in the initial model. The initial model did not include wild boar
 density and cattle density because of their strong correlation with the
 non-motorized variable. The annual number of culled red deer and red deer
 density (t<sub>-1</sub>) were excluded because of their strong
 correlation with red deer density (t<sub>0</sub>). Red deer
 density (t<sub>-2</sub>) was omitted due to its strong
 correlation with red deer density (t<sub>-3</sub>). Lastly,
 the mean annual temperature in years with feeding sites present was
 significantly lower than when absent, suggesting a spurious correlation,
 leading to the exclusion of the feeding sites variable from the initial
 model. After selection, the following control variables remained: age,
 density of red deer (t<sub>-3</sub> and
 t<sub>0</sub>), available habitat, mean annual precipitation,
 and mean annual temperature. We constructed two
 Generalized Linear Mixed Models (GLMMs) in R (v.4.2.1.; R Core Team 2021)
 to assess the direct and indirect effect via body weight of two
 recreational components, motorized and non-motorized, on red deer
 pregnancy rates. For the direct effect of recreation on pregnancy rates,
 we used a Binomial GLMM with a log-link function with the glmmTMB package
 (v.1.1.4; Brooks et al., 2017). In addition to the control variables
 mentioned above, body weight was included as an independent variable. Age
 was square-root transformed to reduce skewness. To obtain odds ratios and
 95% Confidence Intervals the model parameters function (parameters
 package, version 0.18.2; Lüdecke et al., 2020) was used. The inclusion of
 the random factor year in the pregnancy model, aimed at addressing the
 collection of data from multiple deer within a single year, proved to have
 little added value as the variance of this random factor was close to
 zero. To assess the effect of recreation on body weight we used a Gaussian
 GLMM, which included the linear- and quadratic term of age (Putman et al.,
 2019) in addition to the control variables mentioned above and the random
 factor year. We used the “drop1” protocol of Zuur et
 al. (2009) to select the models with the lowest value for the Akaike
 Information Criterion with a correction for small sample sizes. We then
 assessed the fit of both final models using residual diagnostics and
 tested the pregnancy model for overdispersion and zero inflation using the
 DHARMa package (v.0.4.6; Hartig and Lohse, 2022). Finally, we used
 piecewise Structural Equation Modelling in R (Lefcheck, 2016) to perform a
 pathway analysis. To integrate the final models for pregnancy rates and
 body weight into a unified Structural Equation Model (SEM), we utilized
 the psem function from the "piecewiseSEM" package in R
 (Lefcheck, 2016). This function enabled the calculation of standardized
 path coefficients, providing a measure of the strength and direction of
 relationships between variables within the SEM. By standardizing the path
 coefficients, we obtained a consistent measure of the magnitude of
 effects, enabling meaningful comparisons across different variables. Our
 analysis focused on assessing the direct effects of recreation intensity
 on body weight and pregnancy rates. Additionally, we explored the indirect
 and total effects of recreation intensity on pregnancy rates, considering
 body weight as a mediator.
description: # Data from: Insights from a 31-year study demonstrate an inverse
 correlation between recreational activities and red deer fecundity, with
 body weight as a mediator
 [https://doi.org/10.5061/dryad.0rxwdbs7t](https://doi.org/10.5061/dryad.0rxwdbs7t) Date of data collection: 1985-2015 Geographic location of data collection: Netherlands, Western Europe Recommended citation for this dataset: Weterings et al. (2024), Data from: Insights from a 31-year study demonstrate an inverse correlation between recreational activities and red deer fecundity, with body weight as a mediator. Dryad. ## DATA & FILE OVERVIEW #### File 1: Database\_Annual #### File 2: Database\_Animal In **database annual**, all variables are presented per year (e.g., mean temperature in degrees Celsius). In database_animal, each row represents an individual red deer. Variables are presented on an annual level, in the year the red deer was shot. Names used in the manuscript might differ slightly from the names used in the article. **Year** = the year in which the red deer was shot **Pregnant** = presence (1) or absence (0) of a foetus **Weight** = bodyweight in kg **Age** = age in years, determined through dental eruption and wear patterns *Data per year* **Feeding sites** = presence (1) or absence (0) of supplementary feeding sites **Red_deer** = total number of red deer **D_RD** =  Density of red deer, calculated from red deer counts (see ‘Red_deer’). **Red_deer1,2,3** = total number of red deer in former years (t = 1 up to -3 years) **D_RD1, D_RD2, D_RD3** = Density of red deer in former years (see ‘Red_deer1,2,3) **WB** = total number of wild boar **D_WB** = Density wild boar, calculated from wild boar counts (see ‘WB’) **Cattle** = total number of semi-domesticated Sayagueas and Scottish Highland cattle, present and monitored in the area since 2002 **D_C** = Density cattle, calculated from the number of cattle (see ‘Cattle’) **Vehicle_bin** = presence (1) or absence (0) of vehicles. Since the vehicle data contained a lot of zeros, we transformed it into a binary variable. **RD_shot** = total number of red deer shot; between 1985 and 2015, hunters harvested female red deer (n = 488) between October and June in assigned locations **Tr_col** = number of traffic collisions **An_temp** = mean temperature per day in 0.1 degrees Celsius from the nearest weather station (i.e., Deelen at 8.4 km 156 from the study site) to calculate the mean temperature per year in degrees Celsius. **An_rain** (precipitation) = mean precipitation in millimetres per day (in 0.1 millimetres) from the nearest weather station (i.e., Deelen at 8.4 km 156 from the study site) to calculate the mean precipitation per year in millimetres. **Forage** (available habitat) = available habitat in km^2 (including habitat types related to foraging (i.e., grasslands, deciduous forest, coniferous forest, mixed forest and heath)) **Nmotor_Recn** = Non-motorized recreation (i.e., dogs, hikers, cyclists (including ATBs) and horse riders)