Population development and landscape preference of reintroduced wild ungulates: successful rewilding in Southern Italy

Study area

Our study was conducted in the Cilento region (40°12′ N−15°12′E, in the Southern Apennines), defined as the Cilento, Vallo di Diano e Alburni National Park (PNCVDA) and its immediate surroundings (Fig. 1). The PNCVDA is the largest protected area in Italy, including marine, lowlands and mountainous areas. In particular, the latter are characterized by a Mediterranean landscape, structured as a mosaic of woody, natural and agricultural patches, with small rural areas.

At high elevations (with Mount Cervati reaching 1,898 m a.s.l.) the habitat is dominated by beech (Fagus sylvatica) forests, interspersed with secondary grasslands (sensus Dengler et al., 2020) originated in historical time (Blondel, 2010). These are high nature value grasslands (Veen et al., 2009; Oppermann, Beaufoy & Jones, 2012; Török & Dengler, 2018), where the natural conditions are more or less unaltered, sometimes after cessation of restation, fire, herbivory and human arable fields (Dengler et al., 2020). At lower elevations, forests of Alnus cordata and Castanea sativa dominate, followed by deciduous oaks (Quercus cerris, Quercus pubescens) and maple trees (Acer sp.). Nowadays, these forests are populated mostly by Italian hare (endemism Lepus corsicanus; Buglione et al., 2018), wild boar (Fulgione et al., 2016; Maselli et al., 2016; Maselli et al., 2014), roe deer, red deer and by the wolf (Canis lupus), the main predator in the region, expanding at a fast pace with a population currently estimated at many dozens of individuals (Buglione et al., 2020a; Fulgione & Buglione, 2022).

The study area, since the 1950s, is strongly affected by farmland abandonment (Bakudila et al., 2015) and by changes in traditional practices (de Filippo, de Luca & Nicoletti, 2002), so named non-mechanized, non-intensive or extensive agriculture (Hamadani et al., 2021) with pastures in fragmented patches. This phenomenon has determine a loss of biodiversity (Cuttelod et al., 2009; Underwood et al., 2009), primarily as a consequence of the loss endemism and those species strongly linked to the human activity (e.g., Emberiza cirlus, Alectoris graeca, Antus trivialis; (Rippa et al., 2011; Brambilla et al., 2010; Brambilla, 2015; Regos et al., 2016)). In addiction, the absence of large grazing animals determined a competition between plant species for light, inducing dominances resulting in decline of diversity of plant species in grasslands (Allan, 2022). Finally, the abandonment of agricultural areas also involves a loss of the typical Mediterranean patches, resulting in a uniformity of the landscape and a reduction of edge population (Buglione et al., 2020b; Mantero et al., 2020).

The abandonment can be considered as the main driver of landscape change given that in the PNCVDA the human activity in the forest (such as logging, construction of roads, fences or houses) are limited and under rules.

In this area, roe deer and red deer were extinct, dating back to the early 1950s, as a consequence of several factors including human persecution and habitat changes. However, roe deer and red deer were reintroduced between 2003–2006 (37 roe deer released in the Campolongo locality and 35 red deer released in the Grava della Vesola locality) (Fig. 1).

Analysis of population

In order to simulate the demography of expansion of the wild ungulates, we performed a population viability analysis (PVA). The PVA was largely adopted for evaluation of management action in endangered species, using even the supplementation scenario to examine the relative benefits of alternative management actions (Ralls, Beissinger & Cochrane, 2002; Lacy, 1993; Burgman, 2005). In reintroduced population, PVA is used to estimate, even in a short period, the mean generation and future demographic trajectory (Vonholdt et al., 2008; South, Rushton & Macdonald, 2000; Armstrong & Ewen, 2002). It is a tool to modulate reintroduction design, future reintroductions and population management.

The PVA was performed in Vortex 10.1.6.0 software (Lacy, Miller & Traylor-Holzer, 2015), putting biological parameters and ecological data for roe deer (Flajšman et al., 2018; Hewison et al., 1999; Hewison & Gaillard, 2001; Mauget, Mauget & Sempéré, 2003; Vanpé et al., 2008), as well as for red deer (Albon et al., 2008; Borowik, Wawrzyniak & Jędrzejewska, 2016; Clutton-Brock, Major & Guinness, 1985; Krzywiński, 1981; Stewart et al., 2005). For both demographic simulations, with 100 iterations each, we compared the number of individuals living currently in the study area with the number expected to be found in 16–20 years since their reintroduction turning on the option “the extinction threshold for the disappearance of each individual of one of the two sexes”.

Habitat suitability analysis

To collected data for spatial elaborations, we gathered direct and indirect (i.e., footprints, excrement, barking of trees, shrub browsing) presence information ascribable to the roe deer and red deer, by using diurnal and nocturnal surveys, with walking transects and spotlight method (Parkes, 2001).

To define the potential distribution for roe deer and red deer, presence data were employed in spatial elaboration by using Maximum Entropy Distribution Model (MaxEnt, Phillips, Anderson & Schapire, 2006) in MaxEnt 3.4.1 software (http://www.cs.princeton.edu/).

A total of 13 spatially environmental variables were selected: elevation, aspect, slope, roads, waterways, broadleaf forests, mixed woods, conifer forests, scrublands, grasslands, tree plantations, agricultural fields, urban area.

Elevation was obtained as a Digital Terrain Model (DTM) from the National Geoportal (http://www.pcn.minambiente.it/mattm/), while both aspect and slope were calculated from the DTM, using “Aspect” and “Slope” GDAL functions in QGIS 3.10 (https://www.qgis.org). The remaining variables were obtained as categorical vector files from Corine Land Cover 2018 classes (EEA, 2018). Then, in the model we considered the Euclidean distance from each point of presence to each categorical variable. To do this, the GDAL function Proximity (raster distance) was implemented in QGIS 3.10, by setting distance units to “geographical units”.

Pearson’s correlation coefficient was calculated to avoid collinearity between the predictive variables and, after the test, distance from urban areas variable was excluded from analysis.

MaxEnt model settings that best fit our data were selected using the R package “ENMeval” (Muscarella et al., 2014). The software implements features of five classes (linear, quadratic, product, threshold, and hinge) and the following alternative combinations of features were tested: linear, linear + quadratic, hinge, linear + quadratic + hinge, linear + quadratic + hinge + product and linear + quadratic + hinge + product + threshold (Muscarella et al., 2014), as well as the regularization values between 0.5 and 4, with 0.5 increments. The model with the lowest Akaike Information Criterion (AICc) was chosen out of the 48 combinations, corrected for a small sample size (Warren & Seifert, 2011). Predictive performance was evaluated by calculating the area under the receiver operating characteristic curve (AUC) (Swets, 1988).

The habitat suitability maps was obtained converting habitat suitable area in discrete maps using “10th percentile training presence” (Pearson et al., 2006) and “maximum training sensitivity plus specificity” as thresholds (Liu, Newell & White, 2016). The total area included within these thresholds was used as suitable habitat area and constituted the clipping mask in future analyses.

The evaluation of suitable area for the two species was also useful to select where to install the camera traps.

Animal distribution and density

To estimate distribution and density of roe deer, red deer and livestock, we used data from camera trapping, that has been proven to be a reliable method to obtain density estimation of wild and livestock populations, compared to other more expensive and time-consuming techniques (Palencia et al., 2021; Romani et al., 2018).

Locations for cameras were selected by dividing the study area in 10 × 10 km UTM quadrants and selecting 12 quadrants which included suitable habitat for roe deer and red deer, gradually distant from the known point of reintroduction of the species (Fig. 2). Within each quadrant, one smaller cell of 2 × 2 km was selected randomly, and four camera traps were placed as close as possible to the four corners of each cell, based on feasibly reachable locations (Tables S1 and S2).

All cameras had a detection radius (r) = 20 m, except for Uovision UV575, (r = 12 m). Detection zone and radius were obtained from manufacturer manuals for each camera trap model used (according to Rowcliffe et al., 2008). Detections were considered independent if they were at least ≥ 10 min apart.

Camera traps were affixed to trees, around 1.50 m from the ground, pointed in the direction most likely to record animal activity and preferentially oriented in order to minimize sun glares. The devices were set to video mode, with a recording length of 30 s and with a refractory period varying from one to ten minutes. PIR sensitivity was set to normal in locations with open vegetation, and to low in locations with surrounding vegetation, to reduce the risk of accidental camera activations due to movement of branches and/or leaves. The camera traps were left active for an average of 15 days, before being moved to a different cell.

In analyzing the video trapping data, we evaluated every sign that allowed to distinguish individuals, such as the horns shape, marks of predation or fighting, peculiarities of the coat color.

Camera trapping took place from spring to fall 2021 and was carried out using: 8 Scout Guard SG-2060-X cameras, with a detection zone θ = 0.87 rad, 4 LTL Acorn 5310 with a detection zone θ = 2.09 rad, 2 Browning Recon Force 4K with a detection zone θ = 0.77 rad, 2 Stealth Cam DS4K, with a detection zone θ = 0.79 rad, 2 Victure HC300, with a detection zone θ = 1.57 rad, 2 Victure HC520, with a detection zone θ = 1.57 rad, and 4 Uovision UV575, with a detection zone θ = 1.75 rad. Density of roe deer, red deer and livestock was estimated according to the Random Encounter Model (REM) (Rowcliffe et al., 2008) in R (R Core Team, 2020), using “camtools” functions (Palencia et al., 2021), taking into account the camera trap models, with the following equation: ${\displaystyle D=\frac{y}{t}\frac{\pi }{vr\left(2+\theta \right)}}$ where D represents the density (individuals/km²), y/t is the trapping rate (expressed as the ratio of independent capture events - y - by the number of trap days - t -, times 100), r (radius of the detection zone) and θ (the angle of the detection zone) are the camera trap detection parameters, v represents the animal movement speed (Rowcliffe et al., 2008).

When different models of cameras have been used for the same quadrant, radius and angle of detection zone were average. Averaged daily movements for ungulate species and livestock were obtained from the literature: 1.99 km/day for roe deer (Pépin et al., 2004; Romani et al., 2018), and 2.88 km/day for red deer (Pépin et al., 2004; Pépin, Morellet & Goulard, 2009). Average daily movements for livestock were 4 km/day for cattle (Pauler et al., 2020), 6 km/day for sheep (Broom, 2021), and 7.2 km/day for horses (Hampson et al., 2010).

For each species, we estimated the average density for the study area as well as local densities for each surveyed cell, along gradually distant areas from the release locations. The density of each cell was interpolated using the IWD Interpolation feature in QGIS and clipped to the extent of suitable habitat for each species. To obtain the distribution of density, we fitted GLM models (using Pseudo-R²) selecting an “Inverse gaussian” family for roe deer and “Gamma” family for red deer, using the R package “lme4”. We considered density as a response variable, and distance from release points as a predictor variable. The latter was calculated using the algorithm Proximity (raster distance) implemented in QGIS 3.10, by setting distance units to “geographical units”. This allowed also to understand whether our deer populations are represented exclusively by the reintroduced individuals or have been enriched by other. Finally, to obtain the area affected by the cumulative presence of grazers, the livestock density was merged to roe deer and red deer density rasters.

Landscape changes

To infer the potential effect of reintroduced wild herbivores and livestock on the landscape, we compared the land use and land cover from 2000 and 2018, considering Level III Corine Land Cover (Bossard, Feranec & Otahel, 2000; EEA, 2018), a third-level hierarchical classification system with 44 classes (i.e., artificial surfaces, agricultural areas, forest and semi natural areas), in a scale of 1:100,000, with a minimum cartographic unit (MCU) of 25 ha, a geometric accuracy better than 100 m.

To obtain areas where land use changes occurred, we compared vector layers for CLC classification from 2000 and 2018 using the overlay function in QGIS Layers, cleaning polygons smaller than 0.25 ha, as these are generally misalignment artifacts or represent areas too small to have a significant impact on the landscape (Büttner & Kosztra, 2011). The resulting polygons were categorized in different land use change classes according to the type of change (de Filippo, de Luca & Nicoletti, 2002). We considered Spontaneous and Inverse successions, as they are opposing processes and are strongly linked to herbivore activity.

Spontaneous Succession (SS), indicator of loss of semi-open habitats (open area surrounded by the forest), was defined as the transformation of natural areas (class 3) into more advanced succession stages according to natural dynamics (grasslands in bushes, bushes in woods, etc.).

Inverse Succession (IS), indicator of landscape opening areas, was defined as transition from natural advanced forms to less advanced states of succession, for example from woods to bushes (e.g., 311 in 322 or 321 or 312 in 322 or 323).

The polygons were subsequently rasterized to obtain files comparable with the density distributions (100 × 100 m resolution).

To infer on the combined effect of wild and domestic grazers’ pressure, we compared the mean grazing density in areas where Inverse and Spontaneous succession occurred, as well as where they not occurred. We tested the differences using Welch’s t-test, to account for different sample sizes as the number of pixels where the changes occurred different.