Modeling habitat suitability and utilization of the last surviving populations of fallow deer (Dama dama Linnaeus, 1758)

Yasin ünal ·?zdemir ?entürk ·Ali Kavgaci ·Halil Süel ·Serkan Gülsoy ·Idris O?urlu

Abstract Düzler?am? Wildlife Reserve Area (WRA) is the last natural habitat of fallow deer (Dama dama) in the world. Fallow deer is native to Turkey, however, its geographical range is currently confined to Düzler?am? WRA, Antalya. To date, a detailed habitat investigation of fallow deer distribution has not been conducted. This study is vital for the last surviving populations of fallow deer in Turkey. Therefore, we studied the habitat suitability and utilization of fallow deer in the Düzler?am? WRA. Vegetation and wildlife inventory was surveyed across a total of 304 sample areas between 2015 and 2017. Plant species were recorded according to the Braun-Blanquet method and wildlife surveys were based on footprints, feces, and other signs of fallow deer. Classification and regression tree techniques, as well as MAXENT, were used to model vegetation and fallow deer habitat. Topographic position index, terrain ruggedness index, roughness index, elevation, and bedrock formation were also calculated and included in the models. Based on our results, we drafted a habitat protection map for fallow deer. To ensure sustainability of habitats where populations of fallow deer are found in Turkey, we developed recommendations such as closuring human access of the 1st-degree Protection Area and reintroduction of the species to other potential habitats.

Keywords Ungulates·Fallow deer·MAXENT·Wildlife reserve area·Düzler?am?·Antalya·Turkey

Introduction

Many important mammal species that have existed in Anatolia in the past such as Asian lion (Panthera leo persica) (Kasperek 1986), Asian elephant (Elephas maximus) (Albayak and Lister 2012), Caspian tiger (P. tigris virgata) (Goodrich et al. 2015) and Cheetah (Acinonyx jubatus) (Chynoweth et al. 2015) (Asl?m et al. 2012) became extinct during or before the nineteenth century. Currently, there are 134 wild animal species in the critical (CR), endangered (EN), and vulnerable (VU) categories of the International Union for Conservation of Nature (IUCN) in Turkey. Among these species, fallow deer (D. dama), hyena (Hyaena hyaena), and leopard (P. pardus) are reported to be on the verge of extinction (Baskaya and Bilgili 2004; Akay et al. 2011; Avgan et al. 2016; ünal and ?ulhac? 2018; Toyran 2018).

Fallow deer has spread over a wide geographic area from European countries such as England, Germany, and the Netherlands to countries including USA, Canada, Argentina, and New Zealand (Chapman and Chapman 1980; Geist 1998; Nowak and Walker 1999; Shackleton 1999; Collins 2003; Masseti 2012). However, the size of the native population in Turkey does not exceed 500 individuals (ünal and ?ulhac? 2018). Therefore, fallow deer is a species considered to be at serious risk in Turkey despite being assessed as "Near Threatened/NT" in the IUCN Red List (Masseti and Mertzanidou, 2008). Fallow deer has only been recorded in the Düzler?am? Wildlife Reserve Area (WRA) in Turkey and is confined to a small area of about 8000-10,000 ha (Arslangündo?du et al. 2010). Additionally, fallow deer habitat is subject to anthropogenic degradation (Marzano and Dandy 2012; ünal et al. 2020).

Antalya-Düzler?am? was protected with the status of Düzler?am? WRA in 1966 (?ahin 2014; Sar?ba?ak et al. 2006) and fallow deer population size was estimated (Masseti 2007; Masseti and Mertzanidou 2008; De Marinis and Masseti 2009). These estimates showed a decline in fallow deer numbers around Antalya, Termessos National Park, within its natural distribution range, and researchers concluded the species was in danger of extinction (Sar?ba?ak et al. 2006; Masseti and Mertzanidou 2008). Düzler?am? WRA, one of Turkey’s most important protected areas, is surrounded by a variety of ecotourism facilities (De Marinis and Masseti 2009; Arslangündo?du et al. 2010; Albayrak et al. 2012). Recreational activities such as off-roading and motor crossing, trekking, hunting, and foraging for wild plants and mushrooms are common in the area (Weber and Rabinowitz 1996; Linnell et al. 2001; Madden 2008; Ripple et al., 2015). The area is home to fallow deer, as well as other species of mammals such as caracal and wild goat. Together, these factors make Termessos National Park and Düzler?am? WRA ideal locations to study human-wildlife interactions as well as fallow deer-habitat relationships.

Among the factors affecting the extinction of a species, the loss of habitat or degradation of habitat conditions has a significant impact. Therefore, to protect fallow deer in Antalya Düzler?am? WRA, habitat factors like vegetation, climate, topography, and the relationship of the species within the habitat should be identified and investigated in detail.

Given that the fallow deer population size is estimated to be no more than 500 and that Düzler?am? RWA is the last natural habitat for the species in Turkey, the future survival of this species. This study will be the first research that focused on modeling and mapping of fallow deer habitat in is a natural environment.

Materials and methods

Study area

This research was conducted in Antalya-Düzler?am? WRA during 2015-2017. The total surface area of the WRA is 29,033 ha with an elevation range of 38 to 1645 m (Fig. 1). The climate is Mediterranean. Typically, winters bring heavy rainfall and summers are hot and dry (Gül?en and S?nmez 2016). Most of the study area was covered by red pine (Pinus brutia) and red pine-juniper mixed forests in different coverage degrees. Shrub-shaped scrub vegetation dominated areas where the natural forest floor vegetation had been destroyed or on sites where productivity was low. At some places, high surficial limestone bordered the distributions of shrublands and forests, and plants could only grow along the cracks. Two hundred and eighty-eight genera and 430 plant species of 76 families have been reported in the region including Termessos National Park and Düzler?am? (Sar?ba?ak et al. 2010). Thirty-three of these species are endemic and the rate of endemism in the region is about 7.7%.

Fig. 1 Map of the location of the study area (Düzler?am? WRA) and plots

In the northern part of the study area, streams flow only in winter and very few springs can be found. In the southern part, although water is more available, a significant part of this water is used to irrigate private lands for farming (Sar?ba?ak et al. 2010).

Data collection

A digital elevation model (DEM) of the whole study area was created by digitizing the contour curves of cell sizes of 1/25000 scale topographic maps using ArcMap 10.2. Obtained grid designs from the DEM were used to site sampling plots for quantifying the habitat use by fallow deer (Fig. 1). Inventory was carried out in two stages from January 2015 to October 2017. Firstly, vegetation was surveyed within 304 sample areas to determine plant species and frequencies. According to the stages recommended for temperate zone forest areas, plant species were identified using the Braun-Blanquet (1932) method and recorded in 400 m2(Chytry and Otypkova 2003) quadrats (Westhoff and Maarel 1973). The Braun-Blanquet method is a hierarchical system for vegetation inventory. In the second stage, fallow deer tracks were recorded using the presence-absence data recording technique in 200 plots overlapping the sample areas where vegetation inventory was made. For this purpose, we sampled plots with a radius of r = 200 cm, placed at 20 m intervals along 400 m transects (O?urlu 2003). Footprints, feces, and all other signs of fallow deer were recorded on these plots.

Elevation, slope, and aspect maps were created using DEM in ArcMap 10.2 software. A topographic position index was calculated to describe landforms and was created via an extension reported by Jennes (2006). A topographic position index approaching zero indicates a flat landscape. A negative value indicates valley or canyon areas. Positive values correspond to hilly and mountainous areas (Jennes 2006; Vinod 2017). Maps of radiation and heat index were developed based on the slope and aspect maps of the study area. Radiation index maps drafted by use of Eq. 1 (Roberts and Cooper 1989) and heat index maps were created by use of Eq. 2 (Parker 1988).

where the value of θ is measured from the north aspect, and the index ranges in value from 0 to 1.

where 202.5° represents the mean warmest aspect (SSW) and the index ranges in value from - 1 to 1.

Three additional digital index maps (roughness, terrain ruggedness, and landform position index) were created using ArcMap 10.2 software. Utilizing the geological map obtained from the General Directorate of Mineral Research and Explorations in Turkey, the bedrock types and attributes table were created. Bedrock types based on polygons were converted to raster format. Finally, maps of climate variables such as annual average temperature (°C) and annual total rainfall (mm) were obtained from the www. world clim. org database, prepared by Hijmans et al. (2005). These descriptive environmental variables were given codes with annupr and meantm before statistical evaluation (Table 1). The main criteria for including these variables in the study were high impact factors in biological diversity, vegetation structure, and distribution of the fallow deer in the study area.

Table 1 Environmental variables and codes evaluated for statistical analysis

Table 2 Multiple permutation test (MRPP analysis) results by vegetation group number

Statistical analysis and creation of model maps

Two different modeling techniques, the Classification Tree (Moisen 2008) and the MAXENT method (Phillips et al. 2006) were applied to determine potentially suitable habitats for fallow deer. Before modeling, Pearson correlation (Sokal and Rohlf 1995) was applied to reduce model variables and eliminate multicollinearity problems due to the relationships between descriptive variables. The validity of the Classification Tree Technique model was checked with training and testing ROC values. Classification Tree Technique yields explainable results by associating environmental variables with presence-absence data. While sightings or sign of deer confirm deer presence in an area, the absence of confirming records is not sufficient to conclude that deer do not use a given area. Therefore, in this study, the maximum entropy (MAXENT) approach, which uses only presence data, was used to delimit the potential distribution areas of fallow deer. The validity of the MAXENT model was checked with Jackknife statistics.

Cluster analysis (Pritchard and Anderson 1971) using the Sorensen-Wards similarity index was applied to define vegetation groups. A multi-response permutation procedure test (Clarke 1993) was used to determine the number of the most explanatory groups. Harmony between model and cluster groups was analyzed by a chi-square test (McHugh 2013). We applied Indicator Species Analysis (Dufrêne and Legendre 1997) to determine indicator plant species for the vegetation groups. Nominal categorical vegetation group data in five groups were modeled with the classification tree technique (Moisen 2008).

Plant species richness (S) at the alpha level was calculated using the following formula in each parcel (Peet 1974).

In the formula, Si represents each different species in plots. The species richness of the study area was modeled using the Regression Tree Technique (Moisen 2008). Using all models obtained in the previous stages, potential vegetation distribution, plant species richness, and fallow deer habitat suitability maps were created for the study area. Finally, a habitat protection map was created in relation to the study area based on the intersection points of these three maps. In the last stage, the boundaries obtained were evaluated in accordance with expert opinion and the necessary corrections were made. SPSS 17 (Brayman and Cramer 2011) and PC-ORD 6 (McCune and Mefford, 2011) software were used in the statistical evaluation and modeling processes. All maps were created with ArcMap 10.2 software.

Results

Vegetation surveys yielded 412 different plant species of 72 different families. Of the 412 species identified, 39 (9.5%) were endemic, of whichAristolochia lycicaandSilene fenzliiEN (Endangered),Echinops onopordum,Gladiolus anatolicusandTrigonella isthmocarpaVU (Vulnerable),Arum dioscoridis,Ballota glandulosissima,Cephalaria elmaliensis,Digitalis davisianaandStachys aleuritesLR(cd),Bupleurum subuniflorumandIberis caricaLR(nt),Scrophularia pinardiiLR(lc) were categorized in the protection group. Fabaceae was the most common family with 57 species, followed by Asteraceae with 45 species, Lamiaceae with 38 species, and Poaceae with 36 species. Those with a frequency of 5% or more made up 26% of the total plant species. The species with the highest frequency wereQ. coccifera(78.3%),Pistacia terebinthusL. (66.1%),P. brutiaTen. (61.8%),Phillyrea latifoliaL. (60.8%),Olea europaea var. oleasterL. (52.3%),Cistus creticusL. (42.1%), andArbutus andrachneL. (31.6%).

Vegetation was separated into three, five, and seven groups from the first separation lines of the dendrogram. To determine the most appropriate group separation number in these three different group applications, a Multiple Permutation Test (MRPP) was performed using the Euclidean Distance formula. T, A, andp(significance level) values obtained from the analysis are presented in Table 2.

Multiple permutation test results indicated that all vegetation discrimination groups were found to be statistically significant (P＜ 0.001). Community ecology studies have suggested that A value (within-group homogeneity) is generally less than 0.1, while values greater than 0.3 are ideal (Mielke and Iyer 1982; McCune and Grace 2002). The values closest to this criterion were identified in the group 5 and the group 7. Between the two groups, the group 5, (T: - 106.650) with a lower T value (between groups distance) represented the ideal vegetation group distinction in this study.

In groups 1-5, there were 71, 61, 19, 27, and 126 plots, respectively. The distinctive indicator species of all these vegetation groups with their indicator values are listed in Table 3. Group 1 was represented byP. brutiaforests. In some places,Quercus cocciferawas a dominant shrub together withP. brutia. At the scrub layer,P. latifoliaandP. terebintuswere dominant withQ. coccifera.Salvia tomentosa,Stipa bromoides,andTeucrium chamaedryswere dominant at the herb layer in this group.

Vegetation group 2 was represented byO. europeamaquis or sclerophyllous forests (Table 3).Q. coccifera,P. latifolia,andP. terebinthusdensely accompaniedO. europeaat the upper vegetation layer in this group. The lower vegetation layer was dominated by semi-shrub speciesP. lycia, and

Euphorbia characias,and herb species such asMicromeria myritifolia,G. setaceum,C. latifolia, andT. distachya.

Q. cocciferadominated maquis and scrublands represented vegetation group 3. Indicator species of this group are listed in Table 3.P. terebintus,J. excelsa,andJ. oxycedruswere also dominant, which differed from group 2. This community is characterized by the high frequency of ruderal plants likePoa bulbosa,B. tectorum,H. geniculatum,A. minus,Arenaria leptoclados,andCerastium brachypetalumat the herb layer. This structure may contribute to the reason for gaps and degradation in the group.

Table 3 The most powerful plant species by vegetation groups

Vegetation group 4 was also represented byP. brutiaforests. Similar to group 1,Q. cocciferawas also dominant in the scrub layer of group 4. However, typical maquis or sclerophyllous species, such asP. latifolia,C. siliqua,andA. andrachne, characteristic of group 1, disappeared in this group or were less dense, e.g.,P. terebintus. This observation was related to the elevational range of this group since group 4 represented relatively higher zones in the study area compared to group 1.

AnotherP. brutiadominated forest in the study area was group 5. This group was similar to group 1 in terms of the dominant species (P. latifolia,P. terbintus,Q. coccifera) at the scrub layer.A. andrachneandC. siliquawere also common species in group 5. However, this group was distinguished from the other groups based on specific indicator species such asC. siliqua,P. spina-christiand,S. spinosum.

C. creticus,a low scrub species of Mediterranean forests and scrublands, was especially dominant at the scrub layer of this vegetation group.

Environmental variables and vegetation groups (from 1 to 5) were modeled by the classification tree technique by assigning each group binary data (0-1). When correlation analysis was applied for descriptive variables, a high negative correlation was found between elevation and annual average temperature (r:- 0.963) and annual rainfall (r:- 0.924). Therefore, only the elevation representing the climate variables was analyzed in the modeling stages to avoid multicollinearity problems.

The model obtained as a result of the classification tree technique included 11 different open nodes. The variables that constructed the model were elevation, ruggedness index, topographic position index, and radiation index. The concordance between cluster groups and model groups was analyzed by chi-square test and it was found that all distinctions were significant atP＜ 0.001 level. As a result of this process, it was determined that the model explained in 11 different nodes yielded statistically significant results. Thus, the distribution map of potential vegetation groups was obtained on the scale of Düzler?am? WRA (Fig. 2a). The distinctive indicator species of all these vegetation groups are listed in Table 3.

Table 4 Relationship between fallow deer distribution and other wild mammals detected in the area

Plant species richness was calculated for each plot and was modeled by using the regression tree technique. The model was described in six different nodes. The variables explaining the tree model were topographic position index, elevation, terrain ruggedness, and roughness index. After obtaining the estimation values for each node of the obtained model, a spreading process was used to develop the distribution map of potential plant species richness for the WRA (Fig. 2b).

Fig. 2 Model maps of Düzler?am? WRA boundaries a Map of potential vegetation groups, b Map of potential plant species richness, c potential distribution map of fallow deer (classification tree technique), d potential distribution map of fallow deer (MAXENT technique)

According to the population survey, 63 (31%) of the 200 sample areas had evidence of the presence of fallow deer. Other mammalian wildlife species found in the sample areas were jackal (Canis aureus), fox (Vulpes vulpes), European hare (Lepus europaeus), rock marten (Martes foina), wild boar (Sus scrofa), wild goat (Capra aegagrus), badger (Meles meles), caracal (Caracal caracal) and weasel (Mustela nivalis) (Table 4).

The model based on the classification tree technique showing the potentially suitable habitats of fallow deer included 6 different nodes. The variables that shaped the model were elevation, roughness, and topographic position index. The training set ROC value of this model obtained for the target species (fallow deer) was 0.806, and the ROC value for test data was 0.710. In terms of the ROC values, the description level of the classification tree model obtained for fallow deer proved accurate (Araújo et al. 2005). Following confirmation tests of the obtained tree model, the potential habitats of fallow deer were mapped according to the classification tree model technique (Fig. 2c).

According to the MAXENT modeling approach, elevation, bedrock formation, topographic position index, and terrain ruggedness were found to be important (Fig. 3). The area under the ROC curve (AUC) was used for evaluating the obtained model and this value was 0.887. The AUC value approaching 1 (0.887) indicates that the model was successful (Baldwin 2009). Therefore, by using this model, the most suitable potential habitats for fallow deer were mapped (Fig. 2d). Maps drawn by the classification tree technique (Fig. 2c) and MAXENT modeling (Fig. 2d) shared about 70% overlap in area.

Fig. 3 MAXENT model’s Jackknife statistics for fallow deer

Fig. 4 Habitat protection map of fallow deer in Düzler?am? WRA borders

A positive relationship between 5 plants species and fallow deer distribution was confirmed (Table 5). Vegetation group 5 had a statistically significant (P＜ 0.05) positive correlation with fallow deer distribution. The distribution of fallow deer was not related to that of other mammal species on the study area (Table 4).

Table 5 Relationship between the distribution of fallow deer and plant species-vegetation groups

A potential habitat protection map was drafted for fallow deer conservation by evaluating all the model maps (habitat suitability maps of fallow deer, potential vegetation groups map, and potential plant species richness map) together (Fig. 4). The area shaded red on the map (5.196 ha) was determined to be the highest priority protection area for fallow deer. It contained the most suitable habitats for fallow deer and also aligned with the 5th vegetation group, which had a positive correlation with the distribution of fallow deer. The brown-shaded regions were partially suitable in the potential habitat suitability models, with high plant species richness, and thus classified as 2nd-priority protection areas (7.007 ha). Green regions (16.441 ha) were determined to be almost incompatible with the habitat requirements of fallow deer and were therefore identified as 3rd-degree protection areas.

Discussion

MAXENT modeling showed a negative relationship between terrain ruggedness and deer habitat suitability. Fallow deer prefer forest openings in flat landscapes in the Mediterranean region (Masseti et al. 2008). According to the topographic position index values, we concluded that fallow deer preferred valleys, flatlands, and areas characterized by gentle slopes. In addition, when we considered elevation, fallow deer preferred landscapes up to approximately 300 m a.s.l., after which habitat conditions gradually became unsuitable.

The classification tree model indicated that fallow deer found the most suitable habitats at elevations between 171.5 and 296.5 m a.s.l. and in areas where the topographic position index value was less than- 126. These findings confirm the observations that fallow deer is drawn to valleys and canyons for resting or hiding after feeding in the early morning and evening hours. From the classification tree model we inferred that fallow deer almost lost all potentially suitable habitat, especially in places where elevation was higher than 368.5 m.

In the MAXENT model, a rock formation proved to be a suitable habitat for fallow deer. Also, areas composed of sandstone and volcanic sediment accounted for the most positive contribution to the obtained model. Even though it is not known whether these two rock formations indirectly affect the habitat preferences of the species; these formations are abundant in the region, especially at low elevations, where there the landform was flat or slightly inclined. We observed that areas covered with limestone were not preferred by fallow deer as limestone was found in mostly steep areas at the upper elevations.

While a habitat protection map for fallow deer was our primary goal, we also produced a map of plant species richness. The plant species richness model and map indicate that species diversity was the greatest at elevations higher than 380 m a.s.l. with topographic position index less than 138 and terrain ruggedness greater than 0.022. Therefore, the 2nd zone seen in the habitat protection map should not be ignored. Protective management of the 2nd zone will also serve the conservation of plant species diversity and especially endemic plants.

We inferred a positive correlation between the 5th vegetation group and the distribution of fallow deer, but no significant relationship was found between the other vegetation groups. The most powerful indicator species identified for the 5th vegetation group wereA. andrachneL.,C. siliquaL.,E. manipulifloraSalisb.,M. communisL., andP. brutiaTen., respectively. We consider these indicator species to be important food sources and to act as thermal and camouflage cover for fallow deer. These five plants are also common in the floristic structure of group 1. Thus, fallow deer preference towards group 5 instead of group 1 may be due to a combination of floristic structures.

In addition to natural areas, we observed fallow deer frequently using agricultural areas such as orange, pomegranate, and olive fields. They often fed on trees, fruits, and leaves of other plants. In this study, we found no other relationship with other mammalian wildlife species choosing the same or similar habitats as fallow deer (ünal et al. 2020). These results confirm that fallow deer, confined to a small area of about 12,000 ha in the 1st-degree protection zone in Duzlercami WRA in Turkey, are under heavy anthropogenic pressure. The reason why the fallow deer cannot expand beyond this small area in Düzler?am? WRA despite the human pressure is that they are surrounded by steep rocky habitat. For this reason, fallow deer sought cover in this habitat during the day in the closed red pine habitat and frequently approach agricultural areas in search of food (Honda 2009).

In conclusion, we offer the following recommendations to protect and increase fallow deer population size in its natural habitat. First, we recommend closure to human access of the 1st-degree Protection Area (Fig. 4) to prevent human activities such as off-roading and motor crossing, trekking, hunting, and mushroom collecting. Second, new fallow deer reintroduction projects should be prepared to introduce it to other potential habitats outside the Düzler?am? WRA to ensure survival of the species. We also recommend that the General Directorate of Nature Conservation and National Parks strictly enforce the ban on fallow deer hunting and levy heavy fines to deter deer hunting (Alkan and Ersin 2018). Additionally, we think that some projects can be prepared to raise native people’s awareness about the importance of fallow deer and protecting its natural habitat.

AcknowledgementsWe thank fallow deer team (Ahmet Koca, Mehmet ?irin Yelsiz, Mevlüt Zenbilci, Kür?at Bal, Abdülkadir Ery?lmaz, Hasan Uysal, Esra Bayram, Tayfun ?sa Aygül, Eda D?nmez, Dilek Arslan) and we gratefully thank Ibrahim Ozdemir, Ahmet Mert, Serkan Ozdemir, Aslan Merdin, Gül?ah Y?lmaz for analysis, comments, and all contribution.

FundingThis work was supported by the Scientific and Technological Research Council of Turkey (TüB?TAK) with a Grant number of 214O248.