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.

Study site

The Veluwe area (1250 km²) 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 (Pinus silvestris)), deciduous forest (i.e., mainly oak (Quercus sp.) and beech (Fagus sp.)) make up the majority of the habitat, interspersed with dry to moist heathland, containing common heather (Calluna vulgaris) and grasses (Ten Houte de Lange, 1977). Apart from red deer, three other ungulate species occur on the Veluwe, namely: wild boar (Sus scrofa), fallow deer (Dama dama), and roe deer (Capreolus capreolus) (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_-3, t_-2, t_-1, and t₀) 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_-3, t_-2, t_-1, and t₀), 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 > 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_-1) were excluded because of their strong correlation with red deer density (t₀). Red deer density (t_-2) was omitted due to its strong correlation with red deer density (t_-3). 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_-3 and t₀), 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.