Introduction
Large terrestrial mammals are threatened worldwide (Yackulic et al., Reference Yackulic, Sanderson and Uriarte2011; Ripple et al., Reference Ripple, Estes, Beschta, Wilmers, Ritchie and Hebblewhite2014) as a result of habitat fragmentation, lack of dispersal opportunities, small population sizes, poaching and hunting (Leimgruber et al., Reference Leimgruber, Gagnon, Wemmer, Kelly, Songer and Selig2003; Ricketts et al., Reference Ricketts, Dinerstein, Boucher, Brooks, Butchart and Hoffmann2005; Henschel et al., Reference Henschel, Coad, Burton, Chataigner, Dunn and MacDonald2014; Milliken, Reference Milliken2014). Under these challenging conditions, modelling population variations and estimating harvest rates (Ginsberg & Milner-Gulland, Reference Ginsberg and Milner-Gulland1994; Allendorf et al., Reference Allendorf, England, Luikart, Ritchie and Ryman2008), 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. (Reference Chelliah, Bukka and Sukumar2013) 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., Reference Chelliah, Bukka and Sukumar2013), 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 km2; Fig. 1) is part of the Western Ghats–Sri Lanka biodiversity hotspot (Myers et al., Reference Myers, Mittermeier, Mittermeier, da Fonseca and Kent2000). 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, Reference Puyravaud and Davidar2013).
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., Reference Kottek, Grieser, Beck, Rudolf and Rubel2006; Peel et al., Reference Peel, Finlayson and McMahon2007). 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, Reference Prabhakar and Pascal1996).
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., Reference Daniel, Desai, Sivaganesan and Ramesh Kumar1987; Baskaran et al., Reference Baskaran, Udhayan and Desai2010). The same deaths were sometimes recorded in both databases. If n1 and n2 are the number of carcasses recorded by the Tamil Nadu Forest Department and the Wildlife Protection Society of India, respectively, and n3 is the number of carcasses recorded in both databases, assuming independence three probabilities can be defined: P1 = n1/N, the probability a carcass will be tallied by the Tamil Nadu Forest Department; P2 = n2/N, the probability a carcass will be tallied by the Wildlife Protection Society of India; and P3 = n3/N = P1 * P2, the probability a carcass will be tallied by both. The unknown total number of carcasses can be calculated as N = (n1 * n2)/n3. 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, Reference Arivazhagan and Sukumar2008; Ramesh et al., Reference Ramesh, Sankar, Qureshi and Kalle2012), whereas in others only the density was reported (Varman & Sukumar, Reference Varman and Sukumar1995; Baskaran et al., Reference Baskaran, Udhayan and Desai2010; Kumara et al., Reference Kumara, Rathnakumar, Ananda Kumar and Singh2012). 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 (Reference Sukumar and Santiapillai1993), 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. (Reference Chelliah, Bukka and Sukumar2013), 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, Reference Jensen2000), combined with Williamson's two-sex model (Williamson, Reference Williamson1959). Details on the model and the rationale for using ratios are provided in Chelliah et al. (Reference Chelliah, Bukka and Sukumar2013). We obtained survivorship data from Chelliah et al. (Reference Chelliah, Bukka and Sukumar2013; Table A.1), Sukumar (Reference Sukumar1989) and Sukumar et al. (Reference Sukumar, Ramakrishnan and Santosh1998). We used a mean of 0.2 calves per adult female per year for the entire reproductive span (Arivazhagan & Sukumar, Reference Arivazhagan and Sukumar2008), 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, Reference Stubben and Milligan2007) 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., Reference Higgins, Hastings and Botsford1997). 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, Reference Sukumar1989, Reference Sukumar2003), 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, Reference Roughgarden and Smith1996; Jensen, Reference Jensen2000; Milner-Gulland & Akçakaya, Reference Milner-Gulland and Akçakaya2001; Berry et al., Reference Berry, Gorchov, Endress and Stevens2008). 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, Reference Sukumar1991). We therefore attempted to calculate the harvest regime of Asian elephants, using nearly the same methods as Chelliah et al. (Reference Chelliah, Bukka and Sukumar2013), 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. (Reference Chelliah, Bukka and Sukumar2013) 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., Reference Karanth, Nichols, Hedges and Hedges2012), the variability among surveys (Jathanna et al., Reference Jathanna, Karanth, Kumar, Goswami, Vasudev and Karanth2015) and a short time span of observations since Daniel et al. (Reference Daniel, Desai, Sivaganesan and Ramesh Kumar1987) preclude trend analysis. The corollary is that a single survey (as used by Chelliah et al., Reference Chelliah, Bukka and Sukumar2013) 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, Reference Sukumar1989), 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, Reference Martin and Vigne1989).
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., Reference Chelliah, Bukka and Sukumar2013), 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., Reference Baskaran, Balasubramanian, Swaminathan, Desai, Daniel and Datye1993; Daniels, Reference Daniels1993; Desai & Baskaran, Reference Desai and Baskaran1996; Baskaran, Reference Baskaran2013; Lakshminarayanan et al., Reference Lakshminarayanan, Karanth, Goswami, Vaidyanathan and Karanth2015). 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., Reference Traill, Bradshaw and Brook2007) but further pressure could lead to its collapse (Maisels et al., Reference Maisels, Strindberg, Blake, Wittemyer, Hart and Williamson2013). 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, Reference Sukumar1991; Chelliah et al., Reference Chelliah, Bukka and Sukumar2013) should be examined more broadly than on the basis of a simple population model, incorrectly applied.
Acknowledgements
This study was supported by a grant from the United States Fish and Wildlife Service (ASE-0485 for the project Elephant Habitats of the Nilgiri Biosphere Reserve: Location, Threats and Management. We thank the Tamil Nadu Forest Department for granting us permission to carry out the study. Pratheesh C. Mammen, Wilfred Lamuel and J. Duraimurugan helped with data collection from Mudumalai. Tito Joseph collated data from the wildlife crime database of the Wildlife Protection Society of India. We thank Henri de Parseval, Anuradha Reddy and two anonymous reviewers for their incisive comments.
Biographical sketches
Jean Philippe Puyravaud works on landscape-level planning for the conservation of threatened species. Priya Davidar is an ecologist interested in the large-scale patterns of species distributions. Rajeev Srivastava has worked as a forest officer in tiger reserves in southern India, and is dedicated to forest and wildlife conservation. Belinda Wright works to combat poaching and the illegal wildlife trade.
