Modelling harvest of Asian elephants Elephas maximus on the basis of faulty assumptions promotes inappropriate management solutions

Abstract

A ratio-based logistic model developed to assess elephant harvest rates, based on a study at Nagarhole Tiger Reserve in India, was recommended as a management tool to control human–elephant conflict through culling. Considering this reserve among others violates an assumption of the logistic model: isolation. Nevertheless, assuming this violation was irrelevant, we re-evaluated the model, with minor modifications, for the neighbouring Mudumalai Tiger Reserve, where we used data from 13 elephant Elephas maximus population surveys to derive bootstrapped sets of population ratios, and mortality records. We generated arrays of harvest regimes and examined which ratio outputs were closest to the bootstrapped ratios. Our results indicated that (1) model outputs corresponded best with the Mudumalai population structure when harvest regimes were extreme and unlikely, (2) there were significant differences in population structure and harvest regimes between Nagarhole and Mudumalai, and (3) only 49% of adult male deaths predicted by model outputs were recorded in official governmental records. The model provides significantly different results among reserves, which invalidates it as a tool to predict change across the entire elephant population. Variability in survey data and inaccuracies in transition probabilities are sufficiently large to warrant caution when using them as a basis for deterministic modelling. Official mortality databases provide a weak means of validation because poaching incidents are poorly recorded. We conclude that the model should be based on validated transition probabilities and encompass the entire regional population.

Introduction

Large terrestrial mammals are threatened worldwide (Yackulic et al., 2011; Ripple et al., 2014) as a result of habitat fragmentation, lack of dispersal opportunities, small population sizes, poaching and hunting (Leimgruber et al., 2003; Ricketts et al., 2005; Henschel et al., 2014; Milliken, 2014). Under these challenging conditions, modelling population variations and estimating harvest rates (Ginsberg & Milner-Gulland, 1994; Allendorf et al., 2008), when applicable, can inform the development of conservation strategies.

In 1999 CITES established the Monitoring the Illegal Killing of Elephants (MIKE) programme, based upon a multitude of assessments from individual countries. To contribute to this effort, conservationists undertake projects to better assess illegal harvest; for example, to validate records of illegal killing of elephants Chelliah et al. (2013) applied a deterministic Leslie matrix model to populations of Asian Elephas maximus and African Loxodonta africana elephants. They simulated population ratios (such as adult sex-ratio) and compared them to those elicited from survey data. Population ratios have the advantage of using broad age categories instead of narrow age intervals, which are more error prone. The model adequately replicated the population structure in the two study areas, and detected a higher poaching intensity than recorded in official databases. On the basis of these coherent results it was concluded that the model could be applied to a number of vertebrate species and could inform the resolution of human–elephant conflict through male-biased culling.

The basic assumptions of the logistic model are determinism, constant carrying capacity, constant birth and death rates, and no immigration or emigration. However, the Nagarhole Tiger Reserve in India, for which the model was developed (Chelliah et al., 2013), is connected to the Bandipur Tiger Reserve, the Wayanad Wildlife Sanctuary, and the Anechaukur and Maukal State Forests, which violates the assumption of isolation. Arguably models built with the objective of predicting population fluctuations can produce acceptable results without always strictly obeying assumptions. To investigate whether non-isolation of populations mattered in this particular case we applied the same model, with minor modifications, to the Mudumalai Tiger Reserve (Fig. 1).

Fig. 1 Location of the Mudumalai Tiger Reserve in Tamil Nadu, India.

We used population data from 13 surveys carried out by various authors (Table 1) in Mudumalai (Fig. 1), nine of which yielded sufficient details to calculate population ratios. Multiple surveys have the advantage of providing statistical variability compared with arbitrary intervals defined around a single observation, between which the ratio sets are accepted as valid. Firstly, we used arrays of harvest scenarios as entry data in the model. We then evaluated calculated ratios against bootstrapped ratio sets. From this comparison we extracted the number of elephants poached, according to the model. We compared these outputs with two databases, one belonging to the Tamil Nadu Forest Department and the other to the Wildlife Protection Society of India. Lastly, we assessed whether the model was valid for any non-isolated reserve.

Table 1 Data from elephant Elephas maximus population surveys in Mudumalai Tiger Reserve, India (Fig. 1), with source, year of survey, method, population and sex ratio.

*These data were derived from original survey data on herd composition and sex ratio, and therefore the numbers in parentheses do not always sum exactly to the totals, following rounding.

Study area

Mudumalai (321 km²; Fig. 1) is part of the Western Ghats–Sri Lanka biodiversity hotspot (Myers et al., 2000). Established in 1940, it obtained the status of Tiger Reserve in 2007, and together with Nagarhole, it belongs to the Nilgiris Biosphere Reserve, a set of protected areas surrounding the Nilgiri Mountains (Puyravaud & Davidar, 2013).

According to the Köppen–Geiger climate classification system the regional climate is at the confluence of four climatic types: tropical monsoon, tropical savannah, temperate dry winter warm summer and temperate dry winter hot summer (Kottek et al., 2006; Peel et al., 2007). Rainfall occurs mostly during the south-west (May–August) and north-east (September–December) monsoons. Mudumalai receives 600–3,000 mm of rainfall annually and the forest types vary, from moist-deciduous and semi-evergreen towards the west to deciduous and dry deciduous in the east (Prabhakar & Pascal, 1996).

Methods

Mortality records

We obtained records of elephant deaths from the offices of the Tamil Nadu Forest Department in the Nilgiris District, and the Wildlife Protection Society of India in New Delhi. The records included the date of discovery of the carcass, the sex of the elephant, estimated age, location, possible cause of death, and wildlife offence, if any. If age and sex were not recorded, an elephant was assumed to be an adult male if the reported cause of death was poaching, as males with tusks are the main target of ivory poachers (Daniel et al., 1987; Baskaran et al., 2010). The same deaths were sometimes recorded in both databases. If n₁ and n₂ are the number of carcasses recorded by the Tamil Nadu Forest Department and the Wildlife Protection Society of India, respectively, and n₃ is the number of carcasses recorded in both databases, assuming independence three probabilities can be defined: P₁ = n₁/N, the probability a carcass will be tallied by the Tamil Nadu Forest Department; P₂ = n₂/N, the probability a carcass will be tallied by the Wildlife Protection Society of India; and P₃ = n₃/N = P₁ * P₂, the probability a carcass will be tallied by both. The unknown total number of carcasses can be calculated as N = (n₁ * n₂)/n₃. Both databases contained records for the period 1986–2010, and therefore, for the purposes of comparison, we ran the model for this time span.

Elephant population data

We compiled data from 13 population surveys conducted during 1985–2009 (Table 1) by various researchers using a variety of methods: line transects, synchronized block count, known/unknown (i.e. the ratio of individually identified elephants to those that could not be identified), visual and photographic methods. In some cases the total number of elephants sighted and the age structure of the population were reported (Arivazhagan & Sukumar, 2008; Ramesh et al., 2012), whereas in others only the density was reported (Varman & Sukumar, 1995; Baskaran et al., 2010; Kumara et al., 2012). For the latter cases we extrapolated the population estimate to the area of the sampled reserve. Nine studies reported the population age structure following Sukumar & Santiapillai (1993), with calves (< 1 year), juveniles (1–5 years), subadults (5–15 years) and adults (> 15 years). We merged the calves and juveniles categories into a single category: juveniles. Normality of distributions was appraised using the Shapiro–Wilk test. We calculated three ratios (Table 2) from the population surveys: the adult sex ratio (asr), the male adult to subadult ratio (mas) and the proportion of adult males in the population (pam). The mas ratio was a minor modification of the male adult to young adult ratio of Chelliah et al. (2013), as no survey data were available on young adults versus older adult males. The ratios are correlated, and therefore for a given year the three ratios were bootstrapped together to identify confidence intervals from their distributions. Each set of three ratios was taken at random to form samples with 22 measurements, representing the 22 years of surveying. A total of 5,000 samples were produced. The mas distribution was normal and the asr and pam were close to normal. We used non-parametric statistics (2.5% quantile, median and 97.5% quantile) of the bootstrapped distributions to identify ratio sets for use as a reference against which the calculated ratios were compared.

Table 2 Adult sex ratio (asr), male adult to subadult ratio (mas) and proportion of adult males ratio (pam) recorded during surveys of the elephant population in Mudumalai Tiger Reserve (Fig. 1).

The population dynamics model

The population model was a discrete time-step, logistic growth model based on a Leslie matrix age-structured model (Jensen, 2000), combined with Williamson's two-sex model (Williamson, 1959). Details on the model and the rationale for using ratios are provided in Chelliah et al. (2013). We obtained survivorship data from Chelliah et al. (2013; Table A.1), Sukumar (1989) and Sukumar et al. (1998). We used a mean of 0.2 calves per adult female per year for the entire reproductive span (Arivazhagan & Sukumar, 2008), with an equal number of male and female births. The carrying capacity K was fixed at 634, the mean population at Mudumalai over the time span of the field study. We used popbio v. 2.4 in R (Stubben & Milligan, 2007) to compute the stable age distribution of the Leslie projection matrix. The model was simulated in R v. 3.1.2 (R Development Core Team, 2014).

Harvest rates

We constructed an array of male (hm) and female (hf) harvest rates. High harvest rates result in population oscillations (Higgins et al., 1997). Values of hm were 0.01–0.40 in increments of 0.001. Values of hf were 0.01–0.1 in increments of 0.0023. In total, 1,600 harvest regimes combined 40 hm and 40 hf. We began the simulation with the population in its stable state and then applied the various harvest regimes. Starting from natural conditions, the application of a particular harvest regime (a single combination of hm and hf) leads to a change in ratios. Significant levels of poaching of male elephants in southern India began during the late 1970s and early 1980s (Sukumar, 1989, 2003), and therefore we modelled a harvest regime for 1980–2010, in six time steps. We used three sets of ratios provided by bootstrapping (Table 3) to compare the outputs to moderate (2.5% quantile), intermediate (median) and severe (97.5% quantile) scenarios of harvesting. The calculated ratios were compared to bootstrapped ratios by means of the composite index, D:

$$D = \sum \left[ {\left\vert {\displaystyle{{asr_o - asr_c} \over {asr_o}}} \right\vert + \left\vert {{\displaystyle{{mas_o - mas_c} \over {mas_o}}} } \right\vert + \left\vert { {\displaystyle{{\,pam_o - pam_c} \over {\,pam_o}}} } \right\vert} \right],$$

where o indicates observations, and c calculations. The metric D assigned the same relative importance to each ratio. The smallest distance (D _min ) between the sets of observed and calculated ratios provided the harvest regime.

Table 3 Bootstrapped ratios for the adult sex ratio (asr), male adult to subadult ratio (mas) and proportion of adult males ratio (pam) recorded during surveys of the elephant population in Mudumalai Tiger Reserve (Fig. 1).

Results

There were 126 elephant deaths recorded in Mudumalai by the Forest Department (1984–2011) and 34 by the Wildlife Protection Society of India (1983–2010), with 12 records overlapping. A total of 98 adult deaths were recorded (38 males, 53 females and 7 individuals whose sex was not recorded). Poaching was the major cause of death and peaked during 1996–2000, with the killing of 11 male elephants.

Nine surveys carried out in Mudumalai during the study period provided sufficient data to calculate population ratios (Table 2). The population (Shapiro–Wilk test, W = 0.96, P = 0.77) and number of individuals in various age categories (calves, juveniles, subadults and adults; data not provided) were normally distributed. The mean bootstrapped asr was as high as 13.1, and the median 8.8, indicating a bias towards females (Table 2).

The simulated stable age distribution ratios were asr = 2.3, mas = 0.51 and pam = 0.148. The median harvest regime of adults with signature ratios closest to the observed median ratios was hm = 0.23 and hf = 0.02, with a minimum of hm = 0.20 and hf = 0.01 and a maximum of hm = 0.27 and hf = 0.02 with 95% confidence (Table 4). These harvest regimes corresponded to the removal of 61 male elephants in 30 years, with a minimum of 56 and a maximum of 66 in a 95% confidence interval. The calculated ratios were relatively far from the median bootstrapped ratios (D _min  = 1.096), with a low sex ratio (calculated asr = 6.867 versus bootstrapped asr = 13.098) and a low median proportion of male adults to subadults (calculated mas = 0.309 versus bootstrapped mas = 0.815). As the model did not capture the bias in sex ratio sufficiently well, and underestimated the proportion of adult males to subadults, we ran the model with a harvest of subadults and then a harvest of all males with the same intensities as previously (Table 4). The model that best predicted the ratios (D _min  = 0.343) assumed all males were harvested, with a harvest regime of hm = 0.30 and hf = 0.01. In this scenario 69 male elephants were removed, with a minimum of 59 and a maximum of 75 in a 95% confidence interval.

Table 4 Comparison of bootstrapped ratios, quantile (%), with calculated ratios (adult sex ratio, asr; male adult to subadult ratio, mas; proportion of adult males ratio, pam) for the elephant population in Mudumalai Tiger Reserve (Fig. 1). D _min is the shortest distance between calculated ratios and bootstrapped ratios. The harvest regime of males (hm) and females (hf) indicates the number of males and females poached, under the given calculated set of ratios.

During 1986–2010 the Tamil Nadu Forest Department and the Wildlife Protection Society of India recorded 13 and 19 adult male deaths, respectively, with an overlap of six records. Correcting with probabilities indicates that 41 male deaths should have been recorded. By comparison, the model adjusted for the reduced time period estimated 53 male elephants poached for a median harvest regime.

Discussion

Population models have been used to estimate levels of harvest for plants, fishes and terrestrial animals (Roughgarden & Smith, 1996; Jensen, 2000; Milner-Gulland & Akçakaya, 2001; Berry et al., 2008). However, deterministic models are based on assumptions that are probably never entirely met in reality. For this reason it is important to assess whether assumption violations can be ignored, particularly when the culling of a threatened species is recommended as a management tool (Sukumar, 1991). We therefore attempted to calculate the harvest regime of Asian elephants, using nearly the same methods as Chelliah et al. (2013), in Mudumalai Tiger Reserve, which is connected to Nagarhole Tiger Reserve by the Bandipur Tiger Reserve.

We considered population data from 13 surveys carried out in Mudumalai (Fig. 1; Table 1), of which nine yielded sufficient detail to calculate population ratios. In 22 years of surveys no clear trend could be detected for the Mudumalai elephant population or for individual age categories. Poaching pressure varies, and the lack of a trend may be related to sampling variation. We therefore considered that the population (634 elephants) and the various age categories were normally distributed and were appropriately described by their arithmetic means. Against the simulated reference of a stable age distribution (asr = 2.3, mas = 0.51 and pam = 0.148), the population structure (Table 2) was biased towards females, with a deficit of subadult males.

Under the scenario of only adult males being poached we found that the closest combination of simulated ratios to the bootstrapped median of observations indicated a harvest regime of hm = 0.23 and hf = 0.02 over 30 years (1980–2010), where a total of 61 adult elephants (56–66) would have been poached. By comparison, Chelliah et al. (2013) calculated harvest regimes of hm = 0.13 and hf = 0.01 over 30 years, which was below our 95% confidence interval. The observed mas ratio in Mudumalai was as high as 0.815 but the calculated mas ratio was low (0.309). We ran the model with harvest of subadults and then harvest of the entire male population, including calves. The shortest distance between model outputs and observations (D _min  = 0.343) was obtained for the median harvest regime of the entire male population, without distinction of age. Under this scenario a total of 221 male elephants would have been killed (185–245) during a 30-year period, with a harvest regime of hm = 0.30 and hf = 0.01.

The databases of the Forest Department and the Wildlife Protection Society of India recorded 26 cases of poaching of adult male elephants in 25 years (1985–2010). Adjusting the model for this time period, the estimated harvest was 53 males, under the assumption that only adult males were harvested. This suggests that only 49% of adult male deaths were reported. A probabilistic correction of the collated databases increased the observed number of deaths to 41 adult males, amounting to 77% of the model output. Poaching incidents were significantly underrecorded in the Forest Department database, with only 13 records of adult male deaths.

We consider the model is not appropriate for a reserve that is connected to other reserves, for several reasons. Firstly, the model did not reproduce the Mudumalai population structure well, unless it assumed a poaching scenario with significant juvenile mortality. However, we have no evidence of such a scenario. Secondly, despite attempts at standardization (Karanth et al., 2012), the variability among surveys (Jathanna et al., 2015) and a short time span of observations since Daniel et al. (1987) preclude trend analysis. The corollary is that a single survey (as used by Chelliah et al., 2013) is probably not sufficient as a reference to validate model outputs. Lastly, the fact that modelled harvest numbers were higher than those recorded is hardly a validation of the model: records of elephant deaths are poor not only because not all carcasses are detected but also because poaching incidents, if detected, are often deliberately unreported. More generally, other factors could also explain the lack of correspondence between the model and observations: allocation of individuals to age categories on the basis of imprecise height measurements, transition probabilities calculated based on the false assumption of a stable age distribution (Sukumar, 1989), the postulation that harvesting of male elephants began only 30 years ago, when in fact wild elephants have been captured in India, and males killed for ivory, from time immemorial (Martin & Vigne, 1989).

The adult sex ratio was higher in Mudumalai (13.1 females to one male) than in Nagarhole (4.8 females to one male; Chelliah et al., 2013), with less than 5% probability that they belonged to the same population. These significant differences are not unexpected, as Nagarhole and Mudumalai have different habitats and different connectivity patterns to other reserves. However, many authors consider these two reserves to be within the range of a single Asian elephant population (e.g. Baskaran et al., 1993; Daniels, 1993; Desai & Baskaran, 1996; Baskaran, 2013; Lakshminarayanan et al., 2015). As the model outputs for Nagarhole and Mudumalai were statistically different we must avoid applying the logistic model in this form to single connected reserves. When the model is utilized for a single connected reserve it produces a ratio structure that may not match that of the regional population. Proposed culling on the basis of a balanced sex ratio in a single reserve could have deleterious effects on the entire population.

The regional Asian elephant population surrounding the Nilgiris Biosphere Reserve is the largest population of the species. It may exceed the minimum viable population size of c. 4,000 adult individuals (Traill et al., 2007) but further pressure could lead to its collapse (Maisels et al., 2013). Population data are scattered, error-ridden and incomplete. To avoid biases in models of population dynamics, assumptions must be respected as closely as possible. We caution that the proposal to cull the Endangered Asian elephant (Sukumar, 1991; Chelliah et al., 2013) should be examined more broadly than on the basis of a simple population model, incorrectly applied.