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.