Integrating spatial analysis and questionnaire survey to better understand human-onager conflict in Southern Iran

Study area

Qatruiyeh 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 km² and is part of the Zagros Mountains. It is a semi-desert with temperate arid climate, vegetated mainly with Zygophyllum eurypterum and Artemisia sieberi²⁰. There are seven villages in the vicinity of the protected area, where pastoralism is the main source of livelihood²¹.

One 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 movement¹⁸. 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 assessment

Qualitative data collection

We 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 method²² (Table S1) as described below (Eq. 1):

We randomly conducted 200 questionnaires in total.

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 areas²².

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: < 30 Years, 2: 30–50 Years, 3: > 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 classification

Naï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 library²³ in R version 3.5.3²⁴ 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: < 1 ha, 2: 1–5 ha, > 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 library²⁵ in R version 3.5.3²⁴.

Onager vehicle collisions

A 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 analysis

Habitat suitability modeling

A 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 toolbox²⁶. The remaining 90 presence points were used in the modeling.

A habitat suitability map for onager was developed using MaxEnt software version 3.3.3k²⁷ to create a resistance map for connectivity modeling²⁸. We used 10,000 pseudo-absence points²⁹. For the training data set, 75% of the presence points were randomly chosen to train and the remaining 25% were used to test the model³⁰. 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).

Slope 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 <|0.7| between all pairs of variables. In addition, multi-collinearity among variables was further examined using USDM³¹ package version 3.5.3²⁴ and variables with VIF (variance inflation factor) > 3 were used as a threshold to exclude variables³². VIF ranged from 1.2 to 1.8 for all variables. Therefore, all variables were retained for habitat modeling.

Resistance surface for connectivity analysis

To 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 above³³. 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 0³³.

Connectivity corridor network simulation

We used the universal corridor network simulator (UNICOR)³⁴ 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 kernels¹⁶ and (2) factorial least-cost paths¹⁵. The factorial least-cost path analysis implanted in UNICOR simulator uses Dijkstra’s algorithm³⁴ to solve the single-source shortest path problem from every mapped species occurrence location on a landscape to every other occurrence location³⁴. 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 analysis³⁵. We calculated the factorial least-cost path network without dispersal the threshold³⁵ 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 summation¹⁵ to produce maps of connectivity among all pairs of presence points.

Evaluating congruence between crossing points and predicted connectivity

We 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 connectivity¹⁸. 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 data¹⁸.

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 × 10⁷ 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 × 10⁷ random samples.

Source: Ecology - nature.com