Does predation risk affect spatial use in an introduced ungulate species? The case of a Mediterranean mouflon alpine colony

Predation risk is known to affect the spatial use of prey species, imposing a trade-off between feeding requirements and predation avoidance. As a result, prey species can leave high-quality forage areas to use sub-optimal, but safer, habitat patches, defined as “refuge areas.” In this study, we describe changes in the spatial use of an introduced ungulate species, the Mediterranean mouflon Ovis aries musimon, following the recolonization (in 1996) of wolves Canis lupus into the Albergian Hunting Estate (Italian Western Alps). Since 1988, we monitored the mouflon population by spring counts from vantage points. We georeferenced all observations and recorded the size and structure of the spotted groups. Finally, we identified available refuges by selecting patches characterized by (i) the presence of rocks and (ii) high values of steepness and ruggedness. We found that mouflons significantly reduced the average distance from refuge areas over the years, with the yearly average distance from refuges being 56% lower after wolves recolonized the area (i.e., 93.8 ± 32.1 vs. 213.1 ± 40.9 m). The analysis of orographic parameters showed that mouflons used patches with higher values in elevation, slope, ruggedness, and a significant difference in all three parameters when comparing years pre and post wolf return. Both sexes were significantly affected, but ewes were particularly sensitive and selected patches closer to refuge areas (75.8 ± 30.3 m) than males (131.0 ± 53.6 m). Our results suggest that the presence of new predators can alter the distribution of an introduced species such as the Mediterranean mouflon, triggering the resurgence of anti-predation behavior.

Predation risk may also have a critical influence on prey spatial selection, affecting their distribution (Festa-Bianchet 1988; Ripple and Beschta 2004;Bongi et al. 2008), gregariousness (Lima and Dill 1990), and vigilance levels (Hochman and Kotler 2007). Under pressure by predators, ungulates are likely to modify their feeding behavior and eventually leave high-quality forage areas to settle in predator-safer patches (Festa-Bianchet 1988;Walker et al. 2007): an adaptive response, function of what is commonly defined as a species "landscape of fear" (Ripple and Beschta 2003). Refuge areas are assumed to be of critical importance for the long-term survival of prey populations (Festa-Bianchet 1988; Ripple and Beschta 2004;Grignolio et al. 2007). The characteristics of these areas depend on the ecology of different ungulate species. Close habitats such as forests are selected because they reduce the probability to be spotted by a predator (e.g., Red deer Cervus elaphus -Creel et al. 2005; Roe deer Capreolus capreolus - Bongi et al. 2008), while open environments are often preferred by gregarious species, as they increase the chances of seeing a possible threat from a distance (e.g., O. canadensis -Risenhoover and Bailey 1985;Wakelyn 1987). Some species select rough, hardly reachable areas where attacks by predators are less likely to occur (e.g., rocky slopes, Nubian ibex Capra nubiana - Kotler et al. 1994;C. ibex -Grignolio et al. 2007;steep slopes and cliffs, O. canadensis -Bleich 1999;O. dalli dalli -Corti and Shackleton 2002) or less likely to be successful (Baruzzi et al. 2017). Other species use habitats close to human activities which are more likely avoided by predators (White-tailed deer Odocoileus virginianus -Hebblewhite and Merrill 2009). Moreover, sex-related differences in habitat use are also described, with sexual segregation reported to increase reproductive success, according to two leading hypotheses, reviewed in Ruckstuhl and Neuhaus (2002). The "reproductive strategy-predation risk hypothesis" predicts that males select areas with more abundant forage but higher predation risk while females maximize the survival of their offspring by using safer habitats to the detriment of nutritional quality (Geist 1971;Mooring et al. 2003;Grignolio et al. 2007). In contrast, the "forage selection hypothesis" (Main et al. 1996) postulates that female ungulates tend to select higher-quality sites to forage, segregating from males in order to meet higher nutritional demands of lactation (Fattorini et al. 2019).The Mediterranean mouflon O.aries musimon was initially introduced to the islands of Corsica (France) and Sardinia (Italy) at the beginning of Neolithic (Vigne 1992). Since the late eighteenth century, mouflons were also introduced to continental Europe for hunting purposes (Uloth 1972). These wild sheep prefer hills or low mountain ranges where they favor open habitats interspersed with forest and rocky areas, if available (Pfeffer 1967;Apollonio and Meneguz 2003). In the Alps, where they were introduced in the first half of the twentieth century (Pfeffer 1967), mouflons show a preference for meadows and the shrub/meadow interface (Pfeffer 1967;Apollonio and Meneguz 2003). Stray dogs (Perco 1977;Rossi et al. 1988;Cugnasse 1992;Nasiadka et al. 2021) and wolves (Meriggi and Lovari 1996;Poulle et al. 1998;Ståhlberg et al. 2017;Nasiadka et al. 2021) are the main predators of the Mediterranean mouflon in Europe. Since the last decade of the past century, wolves have been successfully recolonizing the Alps, where they had been eradicated at the beginning of the twentieth century (Fabbri et al. 2007), well before the first European mouflons were introduced. Accordingly, the wolf-mouflon interaction was a totally new event in the area.
The aim of this contribution is to describe the changes in spatial use of mouflon in the Western Alps (Italy), using a 25-year-extended series of count data of mouflons, collected with the same method, under the supervision of one of the authors (PGM). Counts were conducted during the birth season in early spring, when predation risk is likely to be stronger as adult mouflons are weakened by winter nutritional deprivation and lambs are more vulnerable to predation (Main et al. 1996). We expected that, according to the "landscape of fear" hypothesis (Ripple and Beschta 2003), wolves ranging in the study area would (i) negatively influence the mean yearly distance of mouflons from available refuge sites, and (ii) force mouflons to select sites with different characteristics, compared to years preceding wolf return. Finally, we investigated whether, according with the "reproductive strategy-predation risk hypothesis" (Ruckstuhl and Neuhaus 2000), spatial behavior would be more affected by predation risk in females than in males.

Materials and methods
The studied mouflon herd dwells in the Albergian Hunting Estate (7,170 ha, 45°2′0″ N; 7°3′0″ E), on the Italian side of the Western Alps. The area has mean elevation of 1,852 ± 457.2 m a.s.l. (minimum = 980; maximum = 3,049) and mean slope of 28 ± 10.9° (minimum = 0.2; maximum = 70.7). Mixed coniferous forest (Larix decidua and Pinus sylvestris) is present at lower altitudes, whereas shrubs (Rhodondendron ferrugineum, Juniperus spp., and Vaccinium spp.) and extensive meadows are dominating above the timber line, at approximatively 2,000 m a.s.l. Mouflons were introduced in 1962, with 12 founder individuals originating from current Croatia (Rossi et al. 1988). The study area is also home of the Northern chamois R. rupicapra, Alpine ibex C. ibex, roe deer C. capreolus, red deer C. elaphus, and wild boar Sus scrofa (Rossi et al. 1988). Potential predators are present, i.e., the golden eagle Aquila chrysaetos, the red fox Vulpes vulpes, and the wolf Canis lupus, the last one returned in 1996 after being locally extinct for approximately 70 years (Fabbri et al. 2007). Thenceforth, one to three wolf packs have been ranging in the area, with two more packs colonizing the area from 2001 (Marucco and Avanzinelli 2010). The mean number of wolves in each pack, measured at the end of winter, was reported to be 4.2 ± 1.8 (Marucco and Avanzinelli 2010).
Throughout the study, the herd was subject to regulated recreational hunting under a system of annual quotas (Supporting Table 1) defined following count-based annual assessments of the mouflon population. From 1988 to 2012, we recorded the size and structure (sex and age composition) of the mouflon population by spring counts carried out from 9 vantage points, scanning 85 subzones (Cruveille and Tufféry 1981). The herd was at high density from 1988 to 1995 (phase H -mean density value: 12.0 ± 1.0 animals/100 ha; mean absolute abundance: 618.5 ± 49.2 individuals) while a sharp decline occurred between 1996 and 1999, from 641 to 146 individuals (phase D). Then, the herd persisted at relatively lower density from 2000 to 2012 (phase L -mean density value: 1.9 ± 0.4 animals/100 ha; mean absolute abundance: 132 ± 43.5 individuals), without any signs of recovery. During counts, observations were reported by subzone on paper maps (scale 1:10,000), then georeferenced as the centroid of the subzone where each individual or group was spotted (accuracy ± 50 m). A point shapefile was built with an open-source GIS software (Qgis 2.0.1 -Quantum GIS Development Team, 2013).
In order to study the changes in spatial use of mouflons, three orographic rasters, i.e., elevation, slope, and terrain ruggedness index -TRI (Riley et al. 1999), were derived from the digital elevation model of the study area (dem -50-m resolution -PrjCRS ED50 32 N) supplied by the Regional Cartographic Catalogue (http:// www. geopo rtale. piemo nte. it/ cms/). TRI provides an objective quantitative measure of topographic heterogeneity. Information regarding land cover were derived from the Local Forestry Plans with a native resolution of 1:10,000 (PrjCRS ED50 32 N) (Piedmont Regional Cartographic Catalogue -http:// www. siste mapie monte. it/ monta gna/ sifor/).
We defined available "refuge areas" as patches characterized by at least one of these characteristics: (i) presence of rocky slopes or rocky meadows, (ii) a slope value above 36° (third quartile of the slope values distribution in the study area), and (iii) a TRI value above 29 (third quartile of the TRI values distribution in the study area). In this way, we derived a raster map of the "refuge areas" which we used to calculate a raster of the distance of mouflons localizations from the nearest refuge site (QGIS 2.0.1 -"Proximity tool").
We considered the yearly average distance (YAD) from a potential "refuge area" as an index of the herd response to predation pressure (QGIS 2.0.1 -"Point sampling tool"). We obtained YAD values for each spotted group by using the following formula: where d i is the distance of a localization from the closest refuge area, n i is the number of mouflons spotted in the single localization, and n TOT represents the total number of mouflons counted each year. The use of n TOT enabled to normalize the index by the mouflon population size along the years.
We used this index to assess (i) the YAD values from "refuge areas" in H and L phases, (ii) the trend of YAD values from 1988 to 2012, and (iii) the "sexual segregation hypotheses" (Main et al. 1996;Ruckstuhl and Neuhaus 2000), by comparing the YAD trend over the years in (1) ewes ≥ 1 years old and lambs vs. (2) males > 1 years old.
To evaluate the population response to predation risk, we also considered the orographic characteristics (i.e., elevation, slope, and TRI) associated to all localizations to assess: (i) the mean values of each characteristic in phases H and L, and (ii) the trend of these values from 1988 to 2012.
Between 1996 and 2000 (phase D), the study area was affected by two exceptionally severe winter seasons, particularly 1995/1996, which led to mass mortality due to starvation (PGM personal communication). To avoid a confounding effect of unusual weather on the herd spatial use (in particular the snow cover), we excluded this period from our comparative analysis. We checked all data series for normality, performing a Shapiro-Wilk test. According to the test results, we used either a Mann-Whitney nonparametric test or a Student's t-test to compare these parameters before and after the return of wolves in the study area (H and L phases) and to test for the "sexual segregation hypothesis." Accordingly, we used either Pearson's or Spearman's correlation tests to evaluate these parameters trend from 1988 to 2012. All statistical analyses were run, and figures were created using R version 3.0.2 (R Core Team 2013).

Discussion
To our knowledge, this is the first study in Europe investigating the effects of wolf presence on the Mediterranean mouflon spatial use. So far, other studies have dealt with habitat and spatial use of mouflon, but all of them referring to areas free from wolves (Dubois et al. 1994;Bon et al. 1995;Cransac and Hewison 1997;Ciuti et al. 2009;Darmon et al. 2011;Marchand et al. 2015;Bourgoin et al. 2018).
"Refuge areas" are critical for the persistence of mouflons, even in areas where wolves are absent, but predators such as eagles and humans are present (Bon et al. 1995;Ciuti et al. 2009). In wild sheep, "refuge areas" have been characterized as steep slopes of broken, rocky terrain (O. aries -Pfeffer 1967;O. canadensis -Smith et al. 1991, Bleich 1999O. dalli -Corti and Shackleton 2002). We thus have considered slope as a critical feature of a "refuge area" as well as the topographic ruggedness index (TRI), another important factor in identifying "refuge areas" (Frair et al. 2005;Laporte et al. 2010;Ciuti et al. 2012). Our data support our working hypotheses. First, we showed that mouflons significantly reduced their average distance from refuge areas over the years, according to Pfeffer (1967) and Bon et al. (1995) who suggested that this species, under predation pressure, moves to broken, steep terrain leaving patches with higher forage quality. The same behavior has been observed for other species (Festa-Bianchet 1988;Berger et al. 2001;Walker et al. 2007;Grignolio et al. 2019) and it is consistent with the "landscape of fear" (Ripple and Beschta 2003). Our hypothesis is also confirmed by the fact that YAD from refuge areas decreased by 56% (from 213 to 94 m) in the low-density phase L, following wolves recolonization in 1996, supporting the idea that the first reaction of the herd was a change in spatial use, through a shift to safer areas (Creel et al. Fig. 2 Variation in distance to refuge areas and orographic parameters between phases H (high density) and L (low density). Clockwise, betweenphases variation in: YAD (yearly average distance) from refuge areas; slope; TRI (terrain ruggedness index); elevation. All comparisons are statistically significant; * = p < 0.05  . As a result, mouflons appeared to have modified their distribution in response to the predator settlement, with individuals exploiting patches closer to refuge areas being more likely to avoid predation, during phase L.
Our expectations were supported also by the orographic parameters analysis. The herd changed its spatial use over the years, selecting patches with higher values in elevation, slope, and ruggedness during phase L. Changes in elevation, as a reaction to predators because of lamb protection, have been reported by Festa-Bianchet (1988) referring to O. canadensis ewes during the pre-weaning period.
Remarkably, our results highlight that mouflons of both sexes (and not only ewes with lambs) dwelt closer to refuge areas over the years, although ewes were more affected than males. These findings are consistent with our assumptions and agree with the key predictions of the "reproductive strategypredation risk hypothesis" (Ruckstuhl and Neuhaus 2000): (i) males select sites of higher quality (forage abundance), even with greater predation risk; (ii) females select sites of reduced predation risk, even with lower forage abundance. However, we were unable to directly test these predictions as we did not measure forage availability across the study period.
Our results are also consistent with studies reporting a development of anti-predatory behavior with increasing predation risk (Berger et al. 2001;Ale and Brown 2009). As mentioned above, the studied mouflons underwent a severe population crash between 1995 and 1998, due to a mass mortality event caused by winter starvation (PGM, personal observation). A recovery was expected in the following years, as observed in other wild sheep populations (Boussès et al. 1994;Moorcroft et al. 1996). However, the herd settled at low numbers. The return of a major predator in 1996 (Marucco et al. 2005) is likely to have affected the viability of the mouflon herd in our study area, preventing a recovery to previous population numbers. Mouflon herds on the French side of the Western Alps decreased in numbers soon after the establishment of wolf packs in their home-range (Poulle et al. 1997;Espuno 2004). In Europe, mouflons are among the favorite prey of the wolf, and they have been reported to be positively selected even when in sympatry with other ungulates (Meriggi and Lovari 1996;Poulle et al. 1998;Ståhlberg et al. 2017;Nasiadka et al. 2021). Possible reasons for this sensitivity to predation, compared with sympatric native ungulates, may be: (i) the limited skills when moving on snow covered terrain (Cruveille and Tufféry 1981); (ii) the apparent weakness of efficient anti-predatory behavior by ewes with lambs (Ferrier M. personal communication); and (iii) the fact that the mouflon is the only ungulate, out of the six species present in the study area, that makes unusually loud and persistent vocalizations from birth to 7 months of age (Apollonio and Meneguz 2003).
However, although our results support our hypothesis that the presence of a major predator had a strong impact on the spatial use of prey, they do not explain if predation has affected the proper recovery of this mouflon herd.
Other variables could also have had a role in driving the described pattern, namely (i) intraspecific competition, (ii) hunting pressure, (iii) density-dependent factors, and (iv) climate changes.
The population size of red deer and chamois significantly increased among the years in the Albergian Hunting Estate (PGM, unpublished data). These species showed a significant dietary overlap with the mouflon herd in this study (Bertolino et al. 2009). However, there is no clear indication in the literature about competition between mouflon and red deer, with available studies suggesting that mouflons appear to outcompete red deer in Hungary (Nahlik and Dremmel 2017) and Ukraine (Smagol et al. 2019). Similarly, competition with the chamois is reported to be favorable to the mouflon in the Alps (Pfeffer and Settimo 1973;Chirichella et al. 2013).
Hunting has always been conducted in the study area, both before and after the wolf returned in the area. Remarkably, the only change in hunting activity was that the number of mouflons harvested decreased together with the population size, with the number of culled mouflons reduced to less than 10 animals/year in the period after the predator return (mean annual harvested = 6.8 individuals; PGM unpublished data). Furthermore, it has often been shown that over-hunted ungulates populations may recover relatively rapidly if hunting is forbidden or substantially reduced (Coulson et al. 2004). Accordingly, an opposite pattern would have been expected in response to hunting reduction, with animals less threatened and thus more likely to expand their range.
Population density could also have influenced the herd spatial use. The population crashed by 77% in only 4 years (decreasing phase D), after two particularly severe winters. Lower densities might have prevented aggregation of individuals and, consequently, reduced the dilution effect provided by large groups (Hamilton 1971). As a result, an increase in perceived predation risk might have pushed the herd to use safer patches (Kotler et al. 1994;Bleich 1999;Corti and Shackleton 2002;Grignolio et al. 2007;. However, other authors have shown the opposite, with ungulate gregariousness increasing with population density (Pépin and Gerard 2008;Vander Wal et al. 2013). In addition, according to the theory on density-dependent habitat selection and the ideal free distribution (Fretwell and Lucas 1970), it is expected that at low density, mouflon would select patches with higher forage abundance, because of lower intraspecific competition, even in the presence of wolves (van Beest et al. 2016). Finally, according to what was observed in an Alpine ibex population, where modifications in ibex behavioral patterns (including changes in spatial use) following wolf recolonization of the area could not be explained by density changes only (Grignolio et al. 2019), density alone does not explain the lack of recover after phase D.
Finally, there is increasing evidence that climate changedriven warming temperatures are resulting in wild mountain ungulates to move at higher elevation, where more nutritious pastures are present (Büntgen et al. 2017;Lovari et al. 2020). According to the forage selection hypothesis (Main et al. 1996), ungulate females are expected to select patches with higher nutritive forage to meet the energetic demands of lactation (Fattorini et al. 2019). If this was true, the stronger decrease in YAD for females might have been the result of the joint action of predation risk and global warming, resulting in an increase in elevation in the ewe's chase for more nutritious grounds. To disentangle the effects of these two factors (Main 2008), future studies should assess pasture productivity changes in respect to climatic factors (e.g., mean temperature) over the period of study.
In conclusion, our study provides insights on the role played by wolf predation on the Albergian mouflon population. Rather than a mere direct predation impact, we suggest a mixed direct and indirect effect, the latter implying the change in habitat use through a shift to safer patches. This result has been often reported for other ungulates (Kotler et al. 1994;Bleich 1999;Corti and Shackleton 2002;Grignolio et al. 2007;, but to the best of our knowledge, not clearly showed so far for the Mediterranean mouflon (but see : Pfeffer 1967;Bon et al. 1995;Ciuti et al. 2009).
To obtain a broader view of the determinants impacting on the dynamics and the viability of free-ranging ungulate populations, we recommend that indirect effects, such as those described in this study, be considered for the effective management and conservation of both alien and native ungulate populations. As wild ungulate populations increase across Europe, they also expand their distribution range generating large-scale human-wildlife conflicts in terms of e.g., crop damage and collisions (Apollonio et al. 2010;Carpio et al. 2020). The results of our study suggest that predator recovery at European scale may help limiting and regulating wild ungulate population growth. In this way, it could be considered a sustainable management tool particularly on natural or protected areas (Pascual-Rico et al. 2020).