Integrating spatial analysis and questionnaire survey to better understand human-onager conflict in Southern Iran
Study areaQatruiyeh National Park, established in 2008, is a core zone in the Bahram-e-Goor Protected Area (established in 1972) at the border of Fars and Kerman provinces in southern Iran (Fig. 1). It covers 310 km2 and is part of the Zagros Mountains. It is a semi-desert with temperate arid climate, vegetated mainly with Zygophyllum eurypterum and Artemisia sieberi20. There are seven villages in the vicinity of the protected area, where pastoralism is the main source of livelihood21.Figure 1Location of the study area. The software ArcGIS. Version 10.2. was used to generate figure. DEM map was downloaded from the WorldClim database (http://www.worldclim.org).Full size imageOne of the major threats for the Persian onager populations in this area is increasing construction of new roads and increasing road traffic. The Sirjan-Yazd (Hassan Abad-Meshkaan) asphalt road, which passes through the Bahram-e-Goor Protected Area, was recently converted into a highway and represents a substantial threat to Onagers (Fig. 1). This road has two lanes in each direction. The day-time speed limit on this road is 110 km/h and 90 km/h at night. Most vehicles on this road are heavy trucks, which pass at high speed (more than 90 km/h), with high traffic volumes at night. During winter, late autumn and summer of drought years, when fodder is scarce, onagers frequently cross the road to access gardens and agricultural fields, which causes high onager mortality due to vehicle collisions. In this research, we used spatial randomization of vehicle collisions and crossing locations to test the predictive ability of resistant kernel and factorial least-cost path predictions of movement18. We also conducted questionnaires with residents from local communities to determine the most important factors influencing human-onager conflicts in the Bahram-e-Goor Protected Area.Human-onager conflict assessmentQualitative data collectionWe administered a questionnaire through a personal interview to 200 randomly chosen farmers residing near onager populations in the Bahram-e-Goor Protected Area in Fars province. Data were collected through a questionnaire between May and August 2018 (Table S1). Ethical clearance was obtained from the DOE (under permit 32–239). All participants were given a printed descriptive summary of the research (if participants were illiterate, the document was read to them). Prior informed consent was obtained orally from all participants. In this research, we followed legal requirements of ethical issues.We calculated the sample size needed by using the family size in rural areas around Bahram-e-Goor Protected Area using the Daniel method22 (Table S1) as described below (Eq. 1):We randomly conducted 200 questionnaires in total.$$N=frac{ {Z}^{2 }P (1-P) }{{d}^{2}}$$
(1)
In this equation, Z is the Z statistic for a level of confidence, P is expected prevalence or proportion (if the expected prevalence is 20%, then P = 0.2), and d is precision (if the precision is 5%, then d = 0.05). In this research, we used d = 0.5 and p was selected according to family sizes in each district of rural areas22.All interviewees were adult males. We collected information on interviewees’ demographic and socioeconomic background (occupation, property, age, and income) as well as their knowledge and opinion on how to prevent onager crop-raiding.We used logistic regression to analyze the significance of sociological factors related to crop damage by onagers. Our dependent variable was “Have you had any of your crop raided by onager during the last year? (Binary response: 1 = Yes, 0 = No)”. Our independent variables included: (1) traditional solutions for reducing Persian onager damages (Response: 1 guarding dogs, 2: fencing around agricultural land, 3: use of traditional barriers (a plastic cuff with a bell on it), 4: scarecrow, 5: turn on the lights at night , 6: Bird-Scarer (Kalaghparan in Persian); (2) which of these solutions could be effective in reducing Persian onager damages (Responses included: 1: fencing around Persian onager habitat, 2: fencing around farmland, 3: give fodder and provide water for Persian onager, 4: buying fodder from local people by DoE, 4: capturing and relocating Persian onager); (3): do you agree with Persian onager hunting? (Binary response: 1 = Yes, 0= No); (4): what is the role of the Persian onager in the wild? (Response 1: distributing seed of plants, the rangelands are restored, 2: it attracts tourists in the region, 3: beauty of nature: God’s creature with a right to live (Intrinsic value), 4: none) (5): age (response: 1: 50 Years), (6): education (response: 1: Incomplete Elementary (lower than 5th grade of elementary), 2: Complete Elementary (5th grade of elementary), 3: Incomplete High school, 4: Associate Degree, 5: Bachelor of Science (BSc), 5: Master of Science (MSc) or Higher), (7) Experience of Persian onager observation in nature: Have you ever seen a Persian onager in the wild? (Response scale: 1 = Yes, frequently, 2 = Yes, several times, 3: Yes, a few times 4: No, never, 5: only seen the Asiatic wild ass carcass), (8) the presence of a Persian onager around your village damages your farms and gardens. How do you feel about this statement? (Response scale) 1: completely disagree, 2: Somewhat disagree, 3: I do not agree or disagree, 4: I agree somewhat, 5: completely agree.All statistical tests were conducted in IBM SPSS Statistics (V. 23.0). Independent variables in the logistic regression analysis were coded as showed in Table S1.Naïve Bayes classificationNaïve Bayes Classification uses a group of simple classifiers based on probabilities, which are applicable to the types of random independent variables in our study. This approach is a supervised machine learning algorithm based on the Bayes Theorem that is used to solve classification problems by following a probabilistic approach. We used the e1071 library23 in R version 3.5.324 for Naïve Bayes classification of onager crop-raiding under this scheme. We considered: Yes (local communities with experience of crop-raid damages), or No (local communities without experience of crop-raid damages during the last one year) as a dependent variable, as a function of the independent variables described in logistic regression section, except we also included farm land area (1: 5 ha) as an additional variable.We categorized data into two groups (testing and training) to determine whether the model performed correctly based on training data. Subsequently, 70% of the data were used to test and run the model along with training confirmation. The Naïve Bayes Classifier was trained to anticipate each attitude in the test data. We calculated the randomness of our results using the Mclust library25 in R version 3.5.324.Onager vehicle collisionsA 25-km section of the 99-km Hassan Abad-Meshkaan road (the area with the highest wildlife-vehicle collision reports) was monitored by motorcycling and walking daily from August to October 2017 (3 weeks). Every morning, we inspected for mammal roadkill within a 30-m buffer on each side of the road, and all carcasses of mammals were recorded using a handheld GPS (Garmin GPS Map 62S). To avoid double-counting, we removed the carcasses after recording. We also obtained collision location data during 2004–2018 from the DoE.The crossing data for onager were obtained from a variety of sources including opportunistic direct observation, environmental guard’s information, and monitoring by LED portable flashlight at night (summer and autumn seasons of 2017 and 2018).Habitat connectivity analysisHabitat suitability modelingA total of 103 presence points were obtained from DoE (2015) in the study area, including Bahram-e-Goor Protected Area, as well as nearby surroundings. To minimize spatial autocorrelation, a 1-km radius was used to eliminate points around each presence location using the SDM toolbox26. The remaining 90 presence points were used in the modeling.A habitat suitability map for onager was developed using MaxEnt software version 3.3.3k27 to create a resistance map for connectivity modeling28. We used 10,000 pseudo-absence points29. For the training data set, 75% of the presence points were randomly chosen to train and the remaining 25% were used to test the model30. We used the area under the ROC curve (AUC) to evaluate model performance. MaxEnt models were completed with 10 bootstrapped replicates.Environmental layers included in MaxEnt modeling included (1) elevation (digital elevation model [DEM]), (2) slope, (3) land cover, (4) distance from agricultural lands, (5) distance from roads and (6) distance from villages. All layers had a 30 m × 30 m resolution (Table 1).Table 1 Environmental variables used for habitat modeling of the Persian onager in the study area.Full size tableSlope was calculated from the DEM layer. Land cover for 27 vegetation classes in the study area was reclassified to 10 classes based on similarities between classes in the original landcover map and due to the importance of agricultural lands (5% of the study area) to onagers. Distance from agricultural lands, roads and villages were included as predictor variables, and were calculated with the Euclidean distance tool in the Spatial Analyst extension of ArcGIS 10.2. We checked for multi-collinearity among variables and correlation was 3 were used as a threshold to exclude variables32. VIF ranged from 1.2 to 1.8 for all variables. Therefore, all variables were retained for habitat modeling.Resistance surface for connectivity analysisTo estimate landscape resistance, we converted the habitat suitability maps to resistance maps using a negative exponential function (R = 1000(−1×HS)) where R represents the cost resistance value assigned to each pixel and HS represents the predicted habitat suitability derived from the suitability models described above33. We used 1000 as the base of our exponential decay function such that areas with > 0.3 habitat suitability would have low-cost resistance. We rescaled the resistance values to a range between 1 and 100 by linear interpolation, such that minimum resistance (Rmin) was 1 when HS was 1, and maximum resistance (Rmax) was 100 when HS was 033.Connectivity corridor network simulationWe used the universal corridor network simulator (UNICOR)34 to predict movement core areas and corridors for Onagers. UNICOR’s key features include a driver-module framework, connectivity mapping with thresholding and buffering, and graph theory metrics. UNICOR produces two kinds of connectivity predictions: (1) resistant kernels16 and (2) factorial least-cost paths15. The factorial least-cost path analysis implanted in UNICOR simulator uses Dijkstra’s algorithm34 to solve the single-source shortest path problem from every mapped species occurrence location on a landscape to every other occurrence location34. The analysis produces predicted least-cost path routes from each source point to each destination point. The resistant kernel algorithm calculates the resistance cost weighted dispersal kernel around each source point up to a user-defined dispersal threshold, and then sums these, producing an incidence function of the rate of organism movement through every pixel in the landscape as a function of the number and density of source points, the dispersal ability of the species, and the resistance of the landscape.According to observation and reports of experts in the DoE, the maximum dispersal of threshold for movement of Onagers is about 100 km. We thus specified a dispersal threshold of 100,000 cost units for the resistant kernel analysis35. We calculated the factorial least-cost path network without dispersal the threshold35 to provide a broad-scale assessment of the regional pattern of potential linkage and to map corridors. The buffered least-cost paths were then combined through summation15 to produce maps of connectivity among all pairs of presence points.Evaluating congruence between crossing points and predicted connectivityWe used a spatial randomization testing procedure to evaluate congruence between the locations where onagers were observed crossing the road and resistant kernel values of predicted connectivity18. Spatial randomization testing of this kind is recommended in cases where there is spatial dependence among observations, and produces an unbiased estimate of the probability of the observed outcome given the data18.We compared the median value of predicted connectivity (resistant kernel) for the 104 actual onager crossing locations with the distribution of median values of 1 × 107 random samples of 104 locations along the highway within the study area. For each combination of resistance surface and connectivity modeling approach, we calculated the ranking of the median of observed values within the distribution of the medians of the 1 × 107 random samples. More
