Materials and methods
Study site
The Veluwe area (1250 km2) is designated as a protected area (i.e., Natura 2000) in the Netherlands (52° 5′ N, 5° 48′ E) (Fig. 2). 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 ind.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 of recreational activities and red deer bodyweight and fecundity between 1985 and 2015 (Fig. 3). We quantified fecundity as a pregnancy rate.

Recreational activities

Between 1985 and 2015, intensity of recreational activities was recorded by a single observer (R. Bijlsma) with the aim to assess their effect on ground-breeding bird ecology (see Bijlsma, 2006). The recreational activities were visually observed in Mosselse zand and Otterlose zand, for a mean time (±SD) of 19.7 (±13.4) hours (range: 2-61 hours) spread across an average (±SD) of 8.9 (±5.4) days per year (range: 2-23 days) (Fig. 2). Recreationists in the area were counted whilst walking without a fixed route, but with complete visual coverage of the whole area. They 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, taking into account the time spent observing.

Red deer characteristics

During the regular hunting season between 1985 until 2015, hunters culled female red deer (n = 488) between October and June in assigned locations, within the area of the Game Management Unit (GMU)(Fig. 2). 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 bodyweight are known to affect fecundity (Borowik et al., 2016). During that period, data on the age and bodyweight were collected by Han ten Seldam and Dirk Lieftink for no specific purpose but to serve as 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. 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.

Characteristics of the study area

We collated data of population variables, to control for their effects on bodyweight and fecundity. Red deer bodyweight 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). Therefore, we collated data on red deer density (t-3, t-2, t-1 and t0) 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 (Fig. 2), 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. 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 from the Game Management Unit (GMU), to account for their effect on chronic stress (Vilela et al., 2020) and red deer density (Carpio et al., 2021), and therefore bodyweight and fecundity (Putman et al., 2019; Bötsch et al., 2020).
Additionally, we collated environmental variables that were known to affect bodyweight 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) (Fig. 4). In addition, we recorded the presence and absence of supplementary feeding sites in the area. 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, bodyweight 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 (Table 1). 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 to 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 t0), wild boar density, cattle density, mean annual temperature and precipitation and the annual number of red deer culled. 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 absence of overlap and therefore one of the variables was not included in the initial model. After selection, the following control variables remained: age, density red deer (t-3 and t0), available habitat, mean annual precipitation and mean annual temperature.
We constructed two Generalized Linear Mixed Models (GLMMs) in R (v.4.2.1.) to assess the direct and indirect effect via bodyweight 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, bodyweight 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. To assess the effect of recreation on bodyweight 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 Akaike Information Criterion with a correction for small sample sizes (Table 3). 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. We assessed the direct effects of recreation intensity on bodyweight and pregnancy rates, and the indirect and total effect of recreation intensity on pregnancy rates, with bodyweight as a mediator.
Results
In total 228 red deer out of 261 were pregnant (mean = 79.3 (± 26.5%) pregnant per year) with an average bodyweight of 61 kg (SD = 8; range: 41-83 kg) (Table 2). The relative density of red deer in the area was on average 2.4 deer.km-2 (SD = 0.4 ; range: 1.8-3.5). Between 1985 and 2015, red deer pregnancy rate decreased by 10.2%, while average female bodyweight decreased with 8 kg.
In the final models, non-motorized recreation had a strong direct negative effect on red deer pregnancy rates (Fig. 5a, Table 4) (β = -0.59), and an even stronger negative effect on bodyweight (Fig. 5b) (β = -1.69). Motorized recreation did not affect pregnancy rates and bodyweight. Bodyweight had the strongest relative (positive) effect on pregnancy rates (Fig. 5c) (β = 1.24). Except age, none of the other control variables, density red deer (t-3 and t0), available habitat, mean annual precipitation and mean annual temperature, affected red deer pregnancy rates and body weight.
Based on a pathway analysis with standardized beta’s (Fig. 6), the indirect effect of non-motorized recreation on pregnancy rates, with bodyweight as a mediator was 29% of the total effect (-0.27). Therefore, the direct effect of non-motorized recreation on red deer pregnancy rates (-0.19) is much larger than the indirect effect via bodyweight (-0.08).
Discussion
Our study investigated the effect of fear of humans on the fitness of a large-bodied ungulate over a 30-year period. It is important to understand the fitness consequences of the fear of humans as a ‘super predator’ on other species for the management and conservation of these species (Darimont et al., 2015; Larson et al., 2016; Naidoo and Burton, 2020). In particular, it enhances our understanding of the mechanisms that can result in trophic cascades (Chitwood et al., 2022), because large-bodied ungulates can structure landscapes such as forests (Ramirez et al., 2018; 2023), affecting forest-related plant and animal species (Svenning et al., 2015). Even though it is known that predators can affect the fitness of small-bodied prey species (Weterings et al., 2022a), large-bodied ungulates may be expected to fear humans more strongly than other (apex) predators, such as wolves, thereby exhibiting strong responses to humans (see Ciuti et al., 2012; Zbyryt et al., 2018; Crawford et al., 2022), especially hunters (Cromsigt et al., 2013). However, compared to hunting, species responses to recreation are generally considered less strong (Gundersen et al., 2021), in particular if recreation is restricted to trails (Reimers and Colman, 2006). Besides, repeated exposure to humans can also lead to habituation to fear (Blumstein, 2016). Nevertheless, in line with our first hypothesis, we found a direct negative correlation between the intensity of non-motorized recreation and the pregnancy rates of red deer. Besides, recreation correlated negatively with bodyweight, whilst bodyweight was strongly positively related with the pregnancy rates of red deer (e.g. Creel et al., 2014; Borowik et al., 2016; Creel, 2018; Putman et al., 2019). This substantiates our second hypothesis that recreational activities can have an indirect effect on pregnancy rates with bodyweight as a mediator. Even though we assumed that repeated recreation could lead to chronic stress, ultimately affecting the pregnancy rate of red deer, we did not measure chronic stress. Boonstra (2013) suggested that chronic activation of the stress axis is an evolved adaptive response to variation in risk. Red deer possibly adapted to the recent increase (i.e., last couple of decades) in disturbance by recreational activities, as adaptive evolution can operate on a generation-to-generation time scale (Bonnet et al., 2022). Nevertheless, both Boonstra (2013) and Creel (2018) argue that acute (unpredictable) stressors (e.g., a risky encounter) can have stress-mediated costs, in particular if they would impact an animal’s stress physiology for a longer duration. It is therefore probable that frequent encounters with recreational activities could result in chronic stress in red deer in the Veluwe (see also Zbyryt et al., 2018).
Even though recreational activities can affect the fecundity of species (Larson et al., 2016; Bötsch et al., 2020), the knowledge about the mechanisms that affect the responses of ungulates to human activities like recreation is far from complete. Recent studies suggest that the responses of ungulates to non-risky cues of recreationists could have consequences on the population level if ungulates incorrectly perceive these cues as risky (i.e., over-response) (see assessment mismatch hypothesis; Smith et al., 2021). Habitat characteristics can modify risk or affect risk assessment of ungulates (Weterings et al., 2022b). Especially in open habitats like the heathland areas in the Veluwe animals have ample visibility and show increased perceptions of risk when disturbed compared to closed habitats (Stankowich, 2008). Next to this, population dynamics are suggested to be driven by a response mismatch (Smith et al., 2021) if ungulates display a generalized antipredator response to all non-risky cues initiated by humans. This could be the case when ungulates associate humans with the potential presence of a dog or hunting activities (Stankowich, 2008), which can both be associated with increased levels of risk (Cromsigt et al., 2013; Darimont et al., 2015).
Our results showed standardized effect sizes ranging from small to medium (Cohen, 1992), which corroborates Stankowich and Blumstein (2005). Possibly, these effect sizes are relatively high because of the presence of hunting in the Veluwe, which could exacerbate species’ responses to recreation (Thiel et al., 2007) especially in areas with a high level of disturbance (Ciuti et al., 2012). Moreover, hunting in the Veluwe is also done with high-powered rifles that is more difficult to associate with a predictable cue (Cromsigt et al., 2013; Ramirez et al., 2021). Besides, recreation in protected areas can be more substantial than recreation in areas that are not protected, because areas that are protected attract visitors (Reinius and Fredman, 2007). This is also seen happening in protected areas around the world where recreation has increased to great numbers (Vallecillo et al., 2019). As a result, recreation in protected areas with a high intensity of use can have a strong impact on species ecology (see e.g. Sundersen et al., 2021). While zonation of recreationists in less-sensitive areas, could be used to manage ungulate populations (Widén et al., 2022), our results primarily emphasise the importance of regulating recreational activities in sensitive locations such as protected areas. In addition to animals, this could also benefit vegetation and soils, as recreation is an important source of pressure on the environment (Ballantyne and Pickering, 2013). Notably, visitors should be more actively informed by protected area managers about the effects of recreation on wild animal populations and their surroundings, as most visitors are unaware about the consequences of their activities (Gruas et al., 2020). This requires protected-area management to be knowledge driven, evidence based, site specific and effective (Pressey et al., 2015), especially when related to management of recreational pressure and its corresponding effects on nature values.
Corroborating our third hypothesis, motorized recreation showed less effect on red deer fecundity than non-motorized recreation. On the one hand, this could apply to other species as well (Stankowich, 2008; Larson et al., 2016). For example, coyote (Canis latrans ), bobcat (Lynx rufus ) (George and Crooks, 2006) and wolverine (Gulo gulo ) (Krebs et al., 2007) showed no responses to motorized recreation as opposed to non-motorized recreation. On the other hand, non-motorized recreation can affect species ecology in many ways (Papouchis et al., 2001; Newsome et al., 2004; Stankowich, 2008; Larson et al., 2016).
Besides recreation, we took various variables into account to control for their effect on fecundity and bodyweight (see Tarlow and Blumstein, 2007), in particular, mean annual temperature and precipitation, food supplementation, available habitat and population density. However, we did not find an effect of the mean annual temperature and precipitation. Possibly, this is because in different periods throughout the year temperature and precipitation have a different effect on the quality and availability of food, red deer bodyweight and fecundity (Schmidt et al., 2001) that could be negated at a large temporal scale. Additionally, we did not find density-dependent effects on fecundity and bodyweight. Possibly density dependence in red deer reproduction and body condition operates only in resource restricted areas (Putman et al., 2019), in contrast to our Veluwe area (Ramirez et al., 2021). Moreover, density dependence is not shown in populations at equilibrium (Clutton-Brock et al., 1985). The relative densities of our population did not vary that much between 1985 and 2015 and fluctuated around an ‘equilibrium’ of 2.4 deer.km-2, making it more difficult to detect density dependence. Similarly, the low variation in available habitat made effects of this variable on fecundity and body weight difficult to detect.
In conclusion, recreational intensity was negatively correlated with fecundity and body weight of red deer. Our results support the regulation of the number of recreationists in sensitive nature areas, especially non-motorized recreationists such as hikers, horse riders and dog-walkers.
Acknowledgements
We are grateful to Rob Bijlsma for providing data on numbers of recreationists in the Veluwe area and critically reviewing a draft version, Han ten Seldam from Natuurmonumenten for providing data on red deer fecundity, age and body weight, Erik Koffeman and Arend Heinen from the Game Management Unit Gelderland for providing data on red deer and wild boar population numbers and the number of red deer shot, and Peter van Geneijgen for providing data on the Sayaguesa cattle. This study was funded by Nationaal Regieorgaan Praktijkgericht Onderzoek SIA (SVB/RAAK.PRO 02.048 to MW) and Van Hall Larenstein University of Applied Sciences.
Competing interests Statement
The authors declare that they have no conflict of interest.
Ethics approval
Red deer were culled and made available by hunters. Not being part of an animal experiment as referred to in the Dutch Act on Animal Experiments, an ethical assessment was not needed. Besides, for this type of retrospective study formal consent is not required.
Data Accessibility Statement
During submission the data base and R scripts will be provided as supplementary data files. After acceptance, the data base and R scripts will be deposited in the Dryad Digital Repository under the reference number [identifier number].
Author Contributions
MW originally formulated the idea. MW, EE, BS, HK designed the study. GJS collated and contributed culling data, data on competitor density and the number of road kills. EE, BS and HK performed statistical analyses. MW, EE, BS and HK wrote the manuscript, while GJS provided editorial advice. All authors reviewed the manuscript.