Age structure changes indicate direct and indirect population impacts in illegally harvested black rhino

Overharvesting affects the size and growth of wildlife populations and can impact population trajectories. Overharvesting can also severely alter population structure and may result in changes in spatial organisation, social dynamics and recruitment. Understanding the relationship between overharvesting and population growth is therefore crucial for the recovery of exploited species. The black rhinoceros (Diceros bicornis; black rhino) is a long-lived megaherbivore native to sub-Saharan Africa, listed as Critically Endangered on the IUCN Red List of Threatened Species. Since 2009, the targeted illegal killing of rhino for their horns has escalated dramatically in South Africa. Given their slow life trajectories, spatial structure and social dynamics, black rhino may be susceptible to both direct and indirect impacts of overharvesting. Our study compared black rhino demography before and during extensive poaching to understand the impact of illegal killing. The population exhibited significant changes in age structure after four years of heavy poaching; these changes were primarily explained by a decrease in the proportion of calves over time. Population projections incorporating both direct poaching removals and decreased fecundity/recruitment were most similar to the observed demographic profile in 2018, suggesting that indirect impacts are also contributing to the observed population trajectory. These indirect impacts are likely a result of decreased density, through processes such as reduced mate-finding, population disturbance and/or increased calf predation. This study illustrates the combined effect of direct and indirect impacts on an endangered species, providing a more comprehensive approach by which to evaluate exploited populations.


Introduction
Overharvesting of wildlife populations is a central concern of conservation biology. Overharvesting can occur as a result of poorly-regulated legal harvest (e.g. sport hunting; [1]), subsistence removals (e.g. bush-meat snaring; [2]), population control strategies (e.g. culling; [3]) and illegal killing for profit (e.g. ivory poaching; [4] such practices, but the extent and consequences of such activities can be difficult to assess, particularly for illegal activities where information must be gathered indirectly. Overharvesting affects both population size and growth as a result of the direct mortalities suffered, and can have profound consequences on population trajectories [5]. Continued removals over time, as is often seen with illegal killing, can also profoundly alter the structure of a population [6]. Large mammal populations exhibit defined age and sex structures, often with class-specific survival rates [7]; as a result, changes in population structure can severely alter temporal dynamics and the population-level response to stochastic environmental variation [8]. Processes that also alter the structure of populations can therefore have severe consequences on the growth and trajectories of endangered species. In addition to the direct mortalities of overharvesting, changes in population structure can also indirectly affect fecundity, recruitment and reproduction [8,9]. Overharvesting can induce behavioural changes that may carry fitness costs as a result of decreased energy intake or increased energy expenditure [10]. Human-induced selection-intentional or unintentionalhas the potential to cause rapid phenotypic and evolutionary change in exploited populations that may reduce viability [11,12]. Overharvest may also cause changes in spatial organisation [13,14], social dynamics/interactions [15] and recruitment [16,17], which can significantly affect population growth [9]. Understanding the relationship between overharvesting and the indirect effects on population growth and viability is therefore crucial for the recovery of exploited species. The black rhinoceros (Diceros bicornis; black rhino) is a long-lived megaherbivore native to sub-Saharan Africa. Once widespread across the continent, the population declined by 97% during the 1900's as a result of overharvesting, and although subsequent conservation initiatives increased numbers to the current population of 5,250 animals, the black rhino remains Critically Endangered on the IUCN Red List of Threatened Species [18]. The Kruger National Park (Kruger) hosts one of the largest black rhino populations in South Africa, which grew from 81 re-introduced founders, beginning in 1971, to an estimated 627 (95% CI: 588-666) animals in 2009 [19]. Analysis of the population vital rates in 2009 showed 6.75% growth per annum, high adult and calf survival and an average inter-calving interval (ICI) of 2.45 years [19]. In 2009, however, the targeted illegal killing of rhino (poaching) for their horns began in earnest within Kruger. Rhino horn is one of the most sought-after commodities in the illegal wildlife trade and can fetch prices in excess of US $35,000 per kilogram [20]. The demand originates primarily in Asia, with rhino horn used as a traditional medicine to treat ailments such as fevers, headaches and cancer [21]. More recently, demand has extended to include rhino horn as a status symbol among the wealthy and a financial investment relying on species extinction for value increase [21].
Black rhinos are elusive mammals and are considered an asocial species [22], but increasing evidence suggests that their social dynamics and interactions are complex [23] and potentially mediate population vital rates [24]. Black rhino exhibit a polygynous mating system [25,26], structured around well-defined spatial organisation. Males are territorial and compete for social and breeding dominance from 8-10 years of age [27], and their territories typically overlap the home ranges of several adult females [28]. Female black rhino usually exhibit overlapping home ranges [28], and may engage in female-based philopatry by sharing part of their home ranges with adult female offspring. Black rhino have a gestation time of 15 months, with age at first calving around seven years [22,25]. Rhino calves are dependent on their mothers for milk for the first year after birth and protection for the first two years, as they are vulnerable to predation primarily by spotted hyenas (Crocuta crocuta) and possibly lions (Panthera leo; [29]). Reproductive senescence has been suggested to occur between 30-35 years of age, with a similar lifespan in the wild [22,30].
In recent years, the proportion of black rhino killed in Kruger relative to their population size has increased, raising concerns regarding the impact of poaching on the stability and viability of the population. Given their slow life-trajectories, spatial structure and potentially complex social dynamics, black rhino may be susceptible to both direct and indirect effects of overharvesting. Our study compared the population demography of black rhino in Kruger before and during extensive poaching to infer the direct and indirect impacts of illegal killing on the population stability and trajectory. More specifically, we: (i) compared pre-poaching age structure with the annual population structure seen over five years after the poaching upsurge, (ii) determined operational sex ratios before and during extensive poaching, and (iii) investigated calf proportion changes between the pre-poaching and poaching period. Finally, we simulated population growth under different poaching and fecundity scenarios and compared the outputs to the current population size and structure to determine the possible contribution of indirect poaching effects to population demography.

Study site
The Kruger National Park (24˚0 0 41@ S, 31˚29 0 7@ E) is South Africa's largest protected area, encompassing more than 19,500 km 2 ( Fig 1A). While mean temperature and rainfall varies considerably along a north-south gradient, summers are typically hot and wet and winters are warm and dry. Black rhino historically occurred throughout this semi-arid savanna ecosystem prior to their local extinction in 1936 following an extended period of targeted hunting [31]. Black rhino re-introductions occurred within the southern region (south of the Olifants River; Fig 1A) of the park [32], with 81 animals re-introduced from Zimbabwe and KwaZulu Natal, South Africa, in the 1970's and 1980's. Black rhino have remained primarily within this area to date. Landscape types within the park are comprised of mixed species woodland on granite and gneiss deposits in the west, and wooded savanna on nutrient-rich basalts in the east [33].

Data
We used observation data of black rhinos collected during aerial surveys in southern Kruger from 2009 to 2018 (SANParks; Fig 1B). Ethics approval was not required as no animals were handled for this study; the data utilised was collected non-invasively by South African National Parks (SANParks) as per the organisation's Wildlife Management Policy. Surveys occurred at the end of the dry season each year, eliminating any effects of month or season on age class frequencies. These observations included age class classification for all individuals and sex classification for adult animals. Black rhinos were aged by standard age classes A-F for rhinos (A: 0-3 months, B: 3-12 months, C: 1-2 years, D: 2-3.5 years, E: 3.5-7 years, and F: 7 years and older; [34]). As we were interested in any changes in functional composition of the black rhino population for this study, we did not require the fine-scale resolution of every age class. In order to increase observation numbers per category, remove possible uncertainty in A and B class assignments and make the data more functionally relevant, we collapsed the age classes into the following categories: 'calf' (0-2 years; A, B and C classes), 'sub-adult' (2-7 years; D and E classes) and 'adult' (7+ years; F class). Poaching of black rhino in Kruger began in 2010, with 6-43 animals killed annually in subsequent years (SANParks; Fig 1C). Annual mortality data included cause and estimated date of death of each rhino carcass found, and sex and age group were recorded when known. Only rhino deaths recorded as poached were included in the analyses as natural deaths were assumed to be reflected in annual survival rates.

Temporal demographic comparison
We compared age category frequencies in each year from 2013-2018 (poaching) to 2009 (pre-poaching). For each comparison, we used chi-squared tests to compare the frequency of individuals in each age category in the poaching year to the frequency obtained pre-poaching. Following the method of Jones et al. (2018), we then calculated the standard residuals (SR) between the observed (O) and expected (E) frequencies for each age category between each poaching year and the pre-poaching frequencies, using the formula SR = (O-E)/ p E. Negative and positive standard residuals indicate observed frequencies that were lower and higher than expected, respectively. To test for evidence of sex-bias in the poaching of black rhino, we calculated the ratio of adult females to adult males (operational sex ratio) for all years and compared the frequency obtained in each poaching year to the pre-poaching ratio using chisquared tests. Finally, we tested whether the ratio of adult females per dependent calf had changed over the study period. We calculated the number of adult females per dependent calf ('calf' category; 0-2 years) in the population for each year. In order to estimate the total number of adult females per year, given variable proportions of unknown sex adults per year, we assigned sex to unknown sex adults based on a 1:1 ratio. We compared this ratio in each poaching year with pre-poaching ratios using chi-squared tests. All analyses were performed in R v3.5.3 [35].

Population projection
In order to visualise the possible black rhino population trajectories over this time period, we used age-based population projection models to simulate population growth from 2009 to 2018 under three different poaching scenarios, using different age and sex ratios of animals killed and, in some scenarios, adjusting for the likely death of dependent calves that were never found (Table 1). Rhino calves left without the protection of their mothers are quickly killed by predators and are consumed in their entirety; rhino calf carcasses are never found in Kruger but clearly occur (pers comm.). For all poaching scenarios, pre-poaching population size and survival rates from Ferreira et al. (2011) were used in the population projections.
In scenario 1, annual adjustments were made to include recorded poaching incidences according to the ratio of known-sex carcasses in the SANParks rhino carcass database ('recorded'). We calculated the sex ratio of known-sex poached carcasses across all years from 2009 to 2018. We then assigned sex to the unknown carcasses each year based on that ratio and added the expected males and females to the known males and females, respectively, to obtain predicted annual totals for each sex (Fig 1C). In scenario 2, annual adjustments were made to include recorded poaching incidences, assigning an equal adult sex ratio to unknown sex individuals and removing the predicted dependent calves killed but never found ('no sex bias + calves'). In scenario 3, annual adjustments were made to include recorded poaching incidences, assigning equal sex and adult:sub-adult ratios, and removing the predicted dependent calves killed but never found ('no sex/age bias + calves').
Age group of illegally killed black rhinos is difficult to assign for the sub-adult and adult categories, in part due to the absence of horns on carcasses and the speed at which tissue is removed from bone; more than 92% of carcasses with recorded ages were categorised as adults (SANParks). Thus for the recorded poaching scenario, we disregarded age group and ran the population projection with all carcasses as adults. For the latter two scenarios which included dependent calves that were likely killed but never found, we used the number of illegally killed adult females per year in each scenario multiplied by the annual proportion of adult females expected to have a dependent calf at foot. For the annual poaching deductions, individuals were removed from the projection vectors by random sampling across the applicable age years. An annual birth rate of 0.33 was used, as a conservative estimate based on the predicted ICI of 2.45 years calculated in 2009 [19]. We compared the predicted population size, age structure, operational sex ratio and dependent calf to adult female ratio under the three different poaching scenarios with the corresponding metrics observed in 2018 to explore the extent to which direct poaching mortalities accounted for the population changes seen.
Finally, to explore the extent of indirect effects, we repeated the population projections under the same poaching scenarios and initial demographic parameters, but decreased the reproductive rate of females. We used ICIs of 4, 5 and 6 years as our proxies for indirect effects; this would encompass decreased calf survival (pre-natal or post-natal) as well as increased time between successive pregnancies. Inter-calving intervals are influenced by many social, spatial and demographic factors likely impacted as a result of overharvesting and therefore are a suitable proxy for combined indirect effects. We compared the resulting trends in predicted population size, age structure, dependent calf to adult female ratio, and operational sex ratio to the corresponding metrics observed in 2018. Population projections were performed in R v3.5.3 [35].

Temporal demographic comparison
The number of black rhino observations in southern Kruger ranged from 122-233 per year. The age category structure of the black rhino population in the pre-poaching period was 22.1% dependent calf, 13.6% sub-adult and 64.3% adult ( Table 2). Comparison of the age group frequencies between the pre-poaching and poaching years showed no significant difference until 2017. Both 2017 and 2018 age structures were significantly different to those observed in 2009 ( Table 2). The population showed an increase in the proportion of sub-adults and a decrease in the proportion of dependent calves with time, and a marginal increase in the proportion of adults (Fig 2). The operational sex ratio of the black rhino population showed no significant differences between pre-poaching and poaching years (range = 0.78-1.25; Chi-squared: 0.007-2.741, all p > 0.05; S1 Table). There was, however, a decreasing trend in the proportion of adult females from 2013-2015, followed by an increase in females in 2016 that remained fairly constant until 2018 (Fig 3A). The number of adult females per calf was significantly higher in both 2017 and 2018 compared to 2009 (Table 3). Annual ratios showed a fairly consistent increasing trend over time (Fig 3B).

Population projection
The 2018 population estimate from the annual rhino survey was 291 (95% CI: 151-441) black rhinos (SANParks internal report). Simulating the direct poaching mortalities under the three scenarios using the pre-poaching demographic parameters resulted in population size projections from 461 (no sex bias + calves) to 589 animals (recorded poaching; Fig 4A). Age structure projections for all poaching scenarios were significantly different to that observed in 2018   Table) were significantly different to those observed in 2018 (dependent calf: adult female = 2.59; operational sex ratio = 1.17) for all poaching scenarios. When including the proxy for indirect poaching effects, increased ICIs of five and six years produced the most similar population size trends to those observed in 2018 (observed: 151-441; ICI 5: 334-408 animals; ICI 6: 303-367 animals; Fig 4). Similar age structures were simulated under the 'no sex/age bias + calves' scenario for ICI 5 and 6 (ICI 5: chi-squared = 1.23, p = 0.5397; ICI 6: chi-squared = 0.71, p = 0.7026; S1A Fig) to those observed in 2018; all other simulations produced significantly different age group proportions (S2 Table). Simulated operational sex ratios converged on those observed in 2018 under the 'no sex/age bias + calves' scenario with ICI's from 4-6 (ICI 4: chi-squared = 3.264, p = 0.07082; ICI 5: chi-squared = 3.3021, p = 0.06919; ICI 6: chi-squared = 2.9854, p = 0.08402; S1B Fig); all other simulations were significantly different (S3 Table). Dependent calf to adult female ratios in all poaching scenarios were similar to those observed in 2018 for ICI 4, ICI 5 & ICI 6 (S4 Table). changes were primarily explained by a decrease in the proportion of dependent calves and an increase in the proportion of sub-adults over time. The ratio of adult females per dependent calf in 2018 was almost double that seen in the pre-poaching period, despite a similar proportion of adults in the population. Population projections incorporating both direct removals and indirect effects (decreased fecundity/recruitment) were most similar to the observed demographic profile in 2018 when poaching was simulated with no sex or age bias, the likely dependent calves of poached females were removed and inter-calving interval was increased to 5 or 6 years.
In our study, both population age structure and the number of adult females per dependent calf changed significantly in 2017 and 2018 compared to the equivalent pre-poaching values in 2009. This suggests that the effect of the direct removals took some years (or some proportion of population removals) to manifest in a detectable change in these metrics. Operational sex ratios did not change significantly, indicating a lack of intentional or unintentional sex bias in poached black rhinos. As there was no survey data prior to the onset of poaching in 2009, neither stochastic population variation nor sampling variation could be definitively estimated in this system. However, annual stochastic variation in population composition in a long-lived mammal with a long gestation period such as black rhino is expected to be low. Furthermore, by combining traditional rhino age classes into functional classes that spanned multiple years, we reduced potential sampling variation across surveys. While overharvesting is generally expected to affect both ends of a population age distribution, in our study it appeared to disproportionately impact the dependent calf age category. This is unlikely purely the result of direct removals, even if the dependent calves that die after their mothers have been killed (and whose carcasses are never found) are considered. For the calf age category to be more affected than the adult category to which their mothers belong, additional factors (indirect effects) must be impacting this black rhino population. A decrease in the proportion of dependent calves following the removal of older individuals has also been seen in elephant populations subjected to poaching pressure [6], and is considered to reflect suppressed recruitment as a function of indirect impacts on breeding and/or calf survival. Demographic Allee effects [36] encompass a range of mechanisms under which a population experiences decreased growth as a result of decreased population size [37]. Cooperative interaction mechanisms include any process where fitness is reduced by reducing conspecific interactions, such as fewer mating opportunities [38,39] or increased mortality due to insufficient defence strategies by smaller groups [40]. Mate-finding, in particular, has been the focus of many studies; if a proportion of individuals are not mated to their full potential, fewer offspring per individual over time can lead to substantial declines in population growth [39]. Given the large expanse of southern Kruger and the uneven distribution of black rhinos across this landscape, it is possible that interactions between individuals are spatially or temporally reduced. Both male and female black rhino maintain well-defined home ranges and are slow to disperse into or recolonize vacant areas [23,41]. The regular removal of individuals by poaching is therefore likely to increase the distances between animals, and/or decrease home range overlap, potentially for many months or years. Linklater and Hutcheson [23] showed that female black rhino did not change their spatial organisation when a neighbouring/overlapping male was harvested from the population, while males decreased their range when a neighbouring female was removed. While the behavioural mechanisms that govern receptivity and mating are poorly understood, the spatial organisation of black rhino is likely an important component of their breeding system and reduced interactions between the sexes at low population density may delay mating and thus increase inter-calving intervals. Hrabar and Toit [42] found that maternal success increased with population density of black rhino in Pilanesberg, South Africa, supporting the importance of cooperative mechanisms within black rhino populations.
Overharvesting of individuals can disrupt the established social dynamics of a population; these perturbation effects are particularly likely to impact populations that exhibit highly structured socio-spatial systems [43]. Increasing the rate of social re-organisation in a population can increase the frequency and intensity of encounters between new individuals [16]. For example, increased male turnover has been demonstrated to increase infanticide thereby decreasing recruitment in large carnivores such as leopards (Panthera pardus; [44]), lions (Panthera leo; [1]) and brown bears (Ursus arctos; [45]). Increased male turnover may also increase natural deaths from territorial fights [44] and decrease adult survival. Furthermore, females may be reluctant to mate when dominant male tenure has been disrupted and a new dominant male has not yet established himself. Frequent turn-over of males as a result of continued poaching in the same area may therefore inadvertently slow breeding rates until social stability is regained. In black rhinos, intra-specific fighting is frequent and often results in injury and/or death [29]; avoiding conflict until the local male hierarchy has stabilised would likely increase survival of both the female and her future offspring.
The change in dependent calf recruitment may, however, also be influenced by increased predation. While adult black rhino are seldom at risk, calves and sub-adults are vulnerable to predation by large predators [29,46]. Observations of black rhino calves with mutilated pinnae and/or missing tails in the Hluhluwe-iMfolozi Complex, South Africa, were attributed to hyena attacks, and subsequent calf disappearance suggested that predation significantly contributed to calf mortality [29]. These findings supported early observations of black rhino calf vulnerability to hyena predation [46]. The hyena population in southern Kruger increased from 3,667 (95% CI: 3443-3891) in 2008 to 7,339 (95% CI: 6998-7680) animals in 2015 (SANParks internal report), and while the drivers of this population growth are unknown, the increase may be related to food availability. The number of rhino carcasses (both black and white) available as a result of poaching increased from 50 to 826 per year from 2009 to 2015 (SANParks, unpubl. data), providing a dramatic increase in easily accessible meat for scavengers annually. Numerous observations of black rhino calves and adults without one or both ears in southern Kruger further supports the potential role of hyena predation on black rhino calf mortality. Lion attacks on black rhino have been recorded on occasion [47,48], but are generally considered infrequent and are therefore unlikely to be a major threat to dependent calf recruitment. Increased predation may also be linked with reduced black rhino density, as single adult black rhino are likely to be less of a deterrent to hyenas and thus a more appealing target for calf predation. Increased calf mortality by hyena predation can, to some degree, also be considered an indirect overharvesting impact, if the increase in rhino carcasses drove the population explosion.
In conclusion, while there has been much focus on the prevention of rhino poaching, little research has been directed towards understanding or mitigating the potentially devastating effects of such frequent removals on the reproductive dynamics of black rhino. The impact of reduced density is considered to be the proximate cause of drastically reduced breeding rate, reproductive malfunction and ultimately, the wild population demise of the Sumatran rhino (Dicerorhinus sumatrensis; [49]), further supporting the role of Allee effects as a primary driver of rhino population dynamics. Future work should attempt to untangle the drivers of recruitment declines in black rhino so that both direct and indirect overharvesting impacts may be addressed, and ultimately ensure that black rhino populations under such extreme pressures do not follow a similar trajectory.