Axis deer (Axis axis) are invasive species that threaten native ecosystems and agriculture on Maui Island. To mitigate negative effects, it is necessary to understand current abundance, population trajectory, and how to most effectively reduce the population. Our objectives were to examine the population history of Maui axis deer, estimate observed population growth, and use species-specific demographic parameters in a VORTEX population viability analysis to examine removal scenarios that would most effectively reduce the population. Only nine deer were introduced in 1959, but recent estimates of >10,000 deer suggest population growth rates (r) ranging between 0.147 and 0.160 even though >11,200 have been removed by hunters and resource managers. In VORTEX simulations, we evaluated an initial population size of 6,000 females and 4,000 males, reflecting the probable 3F:2M sex ratio, with annual removal rates of 10%, 20%, and 30% over a 10-year period. A removal rate of 10% resulted in a positive growth rate of 0.103 ± 0.001. A 20% removal rate resulted in only a slightly negative growth, while a 30% removal rate resulted in –0.130 ± 0.004. By increasing the ratio of females removed to 4F:1M in the 30% harvest scenario, the decline nearly doubled, resulting in –0.223 ± 0.004. Effectively reducing axis deer will most likely require an annual removal of approximately 20–30% of the population and with a greater proportion of females to increase the population decline. Selective removal of males may not only be inefficient, but also counterproductive to population reduction goals.
Herbivorous mammals have been introduced on oceanic islands throughout the world, often with devastating consequences for native biota (Coblentz 1978, Courchamp et al. 2003). The adverse ecological effects of non-native ungulates in the Hawaiian Islands have been well-documented in more than 58 studies (Leopold and Hess 2017). Management actions have eradicated some feral ungulate species from enclosed portions of mostly federal lands in Hawai‘i, resulting in measurable ecosystem recovery (Hess 2016). However, management has been difficult to achieve throughout larger landscapes, especially for wild ungulate species that have never been domesticated. Ungulate populations have persisted or have increased in some cases despite large numbers of removals over long periods of time, and because of conflicting mandates to protect native species while providing hunting opportunities (Banko et al. 2014, Leopold and Hess 2017). Strategies for more effective large-scale control or eradication programs of wild ungulates based on species specific mating systems and population dynamics have not been developed. One such wild ungulate species introduced to Hawai‘i is axis deer (Axis axis) which have a polygynous mating system and have not been domesticated (Hess 2008).
Axis deer were first introduced to the island of Moloka‘i as a gift to King Kamehameha V in 1868. Populations were then established for hunting on Lāna‘i in 1920, and on Maui in 1959 and have reached high levels of abundance, becoming problematic to ranching, agriculture, and conservation of natural areas in Hawai‘i (Tomich 1986). Axis deer populations are female-biased, comprising as much as 70% of individuals, and 95% or more of dams give birth to fawns after reaching maturity as early as 6–10 months of age resulting in rapid annual population growth rates ranging from 20 to 30% (Anderson 2003).
Deer abundance reached an estimated 7,500–11,000 on Maui by 2013 since their introduction in 1959 (Tom Gieder, Wildlife biologist, unpublished data). Deer on Maui have caused more than one million USD of damage each year to vegetable crops, sugarcane fields, vineyards, ranches, golf courses, and ornamental plants at resorts (Hess 2016). Although resource managers generally remove all ungulates from federal lands throughout Hawai‘i, deer are difficult to manage on some inaccessible state-owned lands, while other lands are managed for sustained-yield hunting (Hess 2016). Despite lifting hunting seasons and bag limits for deer on public lands and commercial removals of large numbers, deer remain abundant on Maui, with no apparent lasting measurable decrease (Hess et al. 2015). It is unclear what level of removals would be necessary to make a sustained reduction in abundance. Moreover, it is likely that population reductions may be achieved more efficiently by focusing efforts on females rather than by lifting all restrictions on hunting. Hence, population viability modeling may demonstrate expected outcomes by examining the ability of different management scenarios to effectively reduce the axis deer population.
We hypothesized that increasing the proportion of females removed would favor reduced population growth rates and provide more efficient control strategies. Our objectives were to examine the population history of axis deer on Maui, and retrospectively estimate observed population growth, and then use species-specific demographic parameters in a population viability analysis to examine removal scenarios that would most effectively reduce the population. We ran simulations where 0%, 10%, 20%, and 30% of the estimated annual population abundance was removed each year. We then adjusted the ratios of males and females removed to explore how sex bias may affect population change.
MATERIALS AND METHODS
Observed Population Growth
We compiled dates of axis deer introduction, abundance observations over time, and records of removals through large scale management actions on the island of Maui (1,883 km2; 20° 48′ N, 156° 20′ W) to calculate plausible rates of growth (λ) in a simple population model (Eberhardt 1987) (Table 1). Beginning in 1959, three males and six females were introduced to Maui (Kramer 1971). The State of Hawai‘i Division of Forestry and Wildlife (DOFAW) reported estimates in the 1990s and 2000s, and the Maui Invasive Species Committee (MISC) continued to monitor the axis deer population. Most population estimates collected from 1960 to 2013 lacked methodological details. However, recent population estimates were determined by line-transect aerial surveys. T. Gieder (written comm.) estimated densities ranging from 4.6 to 58.1 deer/km2. The total geographic range on Maui had not been defined, but Gieder (written comm.) used a mean density of 9.1 deer/km2 and estimated a total population of more than 10,000 deer within a 1,100 km2 area. The highest reported abundance of 12,000 deer was in 2011, before nearly 10,000 deer were removed (Table 1). We plotted lower and upper abundance estimates or confidence intervals where available and calculated annual population growth rates (λ) by solving Eberhardt’s (1987) simple model for population projections:
TABLE 1
Reported Axis Deer (Axis axis) Abundance Estimates and Removals Since Their Introduction to the Hawaiian Island of Maui in 1959
We searched for minimum values of λ corresponding to 0.1% increments, that when projected from the founding population number (N1), simultaneously satisfied abundance observations, number of removals (R), and yielded terminal abundance (N2) greater than most recent abundance estimates. The annual multiplicative growth rate λ was then approximated by the exponential growth rate term (r):
This simple modeling approach relied on an assumption that the rate of population growth was fixed over time and the population had not reached density-dependent limitation, which has probably not yet happened given the short history of the species on Maui.
VORTEX Population Modeling
We modeled the axis deer population on Maui under various management scenarios of population removal using the population viability analysis program VORTEX (Lacy 2000, Lacy and Pollak 2014, Lacy et al. 2018). VORTEX provides a matrix population modeling framework that can simulate the extinction process of wildlife populations by comparing the relative effects of potential management scenarios on population growth or persistence (Reed et al. 1999, Fantle-Lepczyk et al. 2018). The program can also be used to estimate the number of removals necessary to control abundance or eradicate invasive species (Pruett-Jones et al. 2007, Licht 2014). A standard life table calculates exponential growth rate (r) by solving the Euler equation:
in which lx and mx are the age-specific mortality and fecundity rates, respectively for age class x to x + 1 and the summation is over all age classes (Lacy 1993, Lacy et al. 2018). The program then estimates the effect of environmental variability, demographic stochasticity, and genetic stochasticity on wildlife populations. The simulations move through a series of discrete and sequential annual population processes, such as reproduction, dispersal, mortality, and harvest (Lacy and Pollak 2014). Optional parameters include the random processes of catastrophes, inbreeding depression, and density dependence.
The population viability analysis provides annual mean population abundance, mean population growth rates, probability of extinction, and estimate of time to extinction. The program allows users to input annual harvest or removals, which in this case could represent either public hunting or culling by management agencies. We modeled several simulations (described below) to determine the baseline rate of population growth of axis deer, the percent of mortality that would need to occur for the population to remain stable, and the percent of additional mortality (by culling) to cause population decline. Furthermore, we investigated how the disproportionate removal of sexes affected the outcome of control efforts. We chose 1,000 iterations of each scenario for improved precision (Lacy et al. 2018). We assumed one closed population and modeled axis deer over a 10-year period. We chose this duration because most preliminary simulations yielded meaningful results (sharp decline in growth or reaching carrying capacity) within 10 years.
TABLE 2
Demographic Parameters Applied in VORTEX Population Viability Analysis Simulations of Axis Deer (Axis axis) on the Hawaiian Island of Maui
Baseline Model Parameterization
We primarily used demographic data from the Hawaiian Islands to parameterize the VORTEX model (Table 2). Where data were lacking from Hawai‘i, we cited research on the demography of axis deer in their native habitat or other locations (Schaller 1967, Ramesh et al. 2012).
Reproduction Parameters — Axis deer are polygynous and an estimated 95% of females >1 year of age breed and 27% of males breed, starting at 2 years of age (Graf and Nichols 1966, Pariwakam 2006). Gogan et al. (2001) examined necropsies during the axis deer control effort at Point Reyes National Seashore, California, and estimated the maximum lifespan and age of reproduction at 10 years of age for both sexes. In VORTEX, we set the mean number of progeny to 1.5 ± 0.5 (SD) fawns per brood in a normal distribution. Occurrences of twinning have been documented in Hawai‘i (Graf and Nichols 1966) and Texas (Fuchs 1977). However, twinning was rarely observed in zoo births (Crandall 1964, Schaller 1967). Based on observed high population growth rates substantiated by simulations, some twinning likely occurred. Furthermore, gestation is approximately 235 days (Chapple et al. 1993), and if fawns die early, dams can give birth a second time within a year (Crandall 1964).
Sex Ratio of Population — Anderson (2003) reported a 30–50% male population based on observations and removals. Gogan et al. (2001) reported the age and sex structure of axis deer removals during control effort at Point Reyes National Seashore. Males accounted for approximately half of all the removals; however, females accounted for more than twice the deer <1 year of age. In our simulations, we opted for a 40% male population, which was the midpoint of Anderson's (2003) observations. Female-biased populations are not unusual in cervid populations, because bucks are often sought for trophies (Jenks et al. 2002). There was scant data regarding sex ratio at birth of axis deer; Graf and Nichols (1966) noted a likely equal ratio of males to females on the islands of Moloka‘i and Lāna‘i, which we used in our simulations.
Inbreeding Depression and Mortality Rates —The impact of inbreeding depression has not been reported in axis deer. We presumed some effect because of the small founding population in Hawai‘i, but because of the reported high densities on Maui, we assumed no deleterious effects to population growth. Nevertheless, we chose a mean inbreeding coefficient and estimate of lethal equivalents for mammals of 3.14 (F) and 50%, respectively, as was reported for mammals by Ralls et al. (1988). Most axis deer mortality estimates come from their native range on the Indian subcontinent, where large predators are present. In their native range, Schaller (1967) reported that 48% of fawns were depredated during their first year. Mortality of yearlings and adult bucks were about 35% (Schaller 1967). Given that there are no large predators in Hawai‘i, we presumed mortality was lower. Although feral dogs (Canis lupus familiaris) are common in Hawai‘i the effect of feral dog predation is unknown. Dogs have depredated larger white-tailed deer Odocoileus virginianus (Huegel et al. 1985) and we presumed similar depredation occurs in Hawai‘i. We assumed predation rates to be lower in Hawai‘i than India, so we reduced Schaller’s (1967) mortality rate for each age class. After running several preliminary simulations, we reduced percent mortality estimates to 35 ± 0.5 (SD) for age 0–1 and 25 ± 0.5 (SD) for ages 1–2 years and >2 years for both sexes.
Carrying Capacity and Starting Population Estimate — We reviewed reported density estimates of axis deer on Maui to estimate carrying capacity. The highest density estimate on Maui was 58 deer/km2 (T. Gieder, written comm.). However, density estimates vary by area depending on land ownership, foraging resources, and barriers that may impede immigration. Given the uncertainty of aerial survey estimates, we chose a more conservative median density estimate of 20 deer/km2, which equated to a carrying capacity and best estimate standard deviation of 22,000 ± 200 SD deer in the 1,100 km2 range on Maui (T. Gieder, unpubl. data.). Based on aerial survey data, Gieder (written comm.) estimated an overall abundance of approximately 10,000 deer in a 1,100 km2 area of Maui, which is what we used in population simulations. We used a function in VORTEX that accounts for estimated birth and death rates to determine a “stable age distribution” for each of the 10 age classes, while maintaining a sex ratio of 40% males in each simulation (Table 3).
TABLE 3
Age Classes of a Starting Population of 10,000 (60% Female and 40% Male) Axis Deer (Axis axis) on the Island of Maui in VORTEX Population Viability Analysis Simulations
Sensitivity Evaluation
We performed a sensitivity analysis in VORTEX to evaluate the importance of demographic or population parameters that lacked certainty. These evaluations examined a range of input values, set by the user, to estimate how a parameter will affect population growth. The sensitivity parameter is defined as:
where ΔX is the change in the observed response of the parameter under examination. We evaluated the sensitivity of 10 parameters including the percent of females and males breeding, the mean number of progeny, mortality of males and females at different lifestages, and sex ratio at birth, and levels of inbreeding depression by adjusting lethal equivalent values. Parameters were evenly spaced across the range of values, according to a “Latin Hypercube” design (Lacy et al. 2018). The parameter space was sampled 1,000 times using a random selection of a value from each parameter. We provided values of the baseline model, the range of values evaluated, and resulting minimum and maximum values in population growth for each model parameter evaluated in the sensitivity analysis (Table 4).
TABLE 4
Demographic Parameters for Axis Deer (Axis axis) on the Island of Maui and in Their Native Range that Lack Certainty in the Published Literature
RESULTS
Observed Population Growth Rates
A growth rate (r) of 0.25 was required to achieve the estimated 85–90 deer in 1968. Abundance increased in the early 2000s to an estimated 12,000 deer in 2009. Overall, estimated minimum and maximum population and growth ranged between 0.147 and 0.160 (Figure 1).
VORTEX Management Scenarios
Starting with 6,000 females and 4,000 males and without hunting or control, the mean baseline growth rate (r) was 0.208 ± 0.001 (SE) and the population reached carrying capacity in 3 years. With annual removal of 10% of the starting population (600 female and 400 male removals), the population increased at a mean growth rate of 0.103 ± 0.001 (SE) and reached a carrying capacity at 10 years (Figure 2). By applying annual removals of 30% (1,800 females and 1,200 males the first year), the mean growth rate declined to –0.130 ± 0.004 and after 10 years the population was 2,759 ± 15 deer. After female removals of four for every male in another 30% annual removal scenario (2,200 females and 800 males the first year), the population growth rate declined sharply to –0.223 ± 0.004 and mean population of 1,086 ± 6 deer at 10 years (Figure 2). These parameter values and five scenarios of Maui's axis deer population trajectories modeled are presented in Tables 4 and 5.
Figure1.
Estimated population growth (r) of axis deer (Axis axis) since their introduction to the Hawaiian island of Maui in 1959. Simple population projections were used to approximate reported abundance (y-axis) after accounting for removals using methods from Eberhardt (1987). Abundance and removal data are from Table 1.
DISCUSSION
Modeling observed population growth indicated that Maui axis deer had life history characteristics that maintained steep population growth despite the large numbers of animals that were removed. The growth rate (r) may have initially exceeded 0.25 if historical abundance estimates were accurate. Reproduction parameters needed to achieve this rate of growth include some combination of adult female-biased sex ratios, low female mortality, and more than occasional twinning. Data on perinatal mortality was notably lacking; however, we assumed it was low because no diseases or parasites that could affect mortality have been documented on Maui (F. Duvall, Hawai‘i Division of Forestry and Wildlife, written comm.). We lacked substantive data on twinning and first-year survival, but our sensitivity analysis results corroborate that both parameters have the highest influence in our VORTEX simulations are important for actual population growth.
Figure2.
Simulated mean population growth rates of axis deer (Axis axis) from the Hawaiian island of Maui during four levels of removal intensity using program VORTEX. Simulations of 1,000 iterations were initialized with 10,000 individuals (6,000 females and 4,000 males). The five scenarios include projected abundance without harvest (solid black line), a population with an annual harvest of 10% and 20% the annual population estimate and two 30% harvest scenarios where three females were removed for every two males (gray dashed line) and another where four females were removed for every male (solid gray line).
TABLE 5
Mean Growth Rates of Axis Deer (Axis axis) on the Island of Maui Estimated in 1,000 Iterations in Program VORTEX
Figure3.
Range of stochastic growth rates of each parameter that lacked certainty modeled in a VORTEX sensitivity analysis of axis deer (Axis axis) on the Hawaiian island of Maui. The black dot represents the baseline scenario growth rate. Parameters with the longest lines have the greatest effect in VORTEX modeling scenarios.
Deer on Maui may have initially encountered favorable habitat conditions in the absence of intraspecific competition and predation that facilitated high fecundity and survival of young animals. In our VORTEX simulations we applied a modest population density of 20 deer/km2 to estimate a carrying capacity of 22,000 deer. However, Maui's axis deer population may reach abundances much higher than our simulations. The core population in central Maui occupies a 300 km2 area without barriers. Densities of deer in areas of East and West Maui are higher than the core population (T. Gieder, written comm.), but those areas are smaller and separated by barriers on private and public lands. In areas of Nepal, axis deer reached densities exceeding 200/km2 (Moe and Wegge 1994). Thus, in favorable foraging conditions without consistent hunting or culling, the Maui axis deer population could possibly reach an abundance previously undocumented and cause further economic damage and degradation of native ecosystems.
The VORTEX simulations demonstrated that reducing axis deer population growth requires a sustained and rigorous management effort of annual harvest. Reducing the population required annual removals of 30% of each annual estimate could be a challenging target for resource managers. Furthermore, our modeled time to a population reduction to approximately 1,086 in 10 years may be unrealistically rapid, because the effort to removal animals typically increases as the population declines (Gogan et al. 2001, Banko et al. 2014, Judge et al. 2017). Judge et al. (2017) monitored the eradication effort of mouflon sheep from a 65 km2 area in Hawai‘i Volcanoes National Park from 2003 to 2017. Removals averaged 6.5 sheep/staff day but the effort expended to remove each sheep increased nearly 15-fold during the last 3 years of the program (Judge et al. 2017). Eradication was not achieved until a smaller 26 km2 subunit was completely enclosed with fence in 2012. Similarly, population reduction of axis deer in areas of Maui may only be accomplished with a consistent effort in fence-enclosed units. Our management simulations also emphasized the importance of determining baseline abundance estimates to set ungulate management targets and set benchmarks for monitoring the success of control efforts.
Given that 95% of adult females breed at only 1 year of age emphasizes the importance of targeting females of all ages during control efforts. Because of their polygynous breeding system, removing males may not reduce productivity. Consistently biased removal of adult males was identified among the factors leading to the failure to reduce axis deer abundance in Argentina (Gürtler et al. 2017). Gogan et al. (2001) also found that simulated removals of only male axis deer and fallow deer Dama dama resulted in populations of both species reaching carrying capacity within a decade, whereas the removal of only females led to the eventual extirpation of both species. Some state game management agencies have implemented “earn-a-buck” programs whereby hunters were required to harvest either young or female white-tailed deer without antlers before harvesting an adult male, which has reportedly increased the harvest of female deer, improved the adult sex ratio, and proven beneficial for forest regeneration and biodiversity (Van Deelen et al. 2010, Boulanger et al. 2011).
A more common management response to overabundant ungulates is to remove all restrictions on hunting. However, unrestricted hunting may lead to a disproportionate removal of males, which are sought for trophies (Stephens et al. 2008, Hess et al. 2015). While severe reductions of either sex can ultimately cause population decreases, in practice, reducing a large proportion of either sex is exceedingly difficult as the population declines (Judge et al. 2017). In such cases, an insufficient number of male removals would likely result in a counterproductive increase in the proportion of females in a population and a consequent increase in the per capita rate of growth, whereas moderately disproportionate removals of females would be more likely to reduce abundance. Therefore, limiting male removals may be a more effective way to reduce abundance or eradicate deer populations. Regardless of whether removals are conducted individually by hunters or by culling, population reduction could be achieved more efficiently by the disproportionate removal of females. The removal of a substantial proportion of males may not only be inefficient, but also counterproductive to any population reduction goals.
Limitations of this research are due to the lack of a comprehensive understanding of total abundance and its change over time, data on the age and sex structure of the population, the proportion of pregnant females in the population, fetal sex ratio, proportion of females bearing twins, and perinatal survival. Precise estimates for the overall abundance of axis deer on Maui would be challenging and expensive to obtain given their evasive behavior and size and complexity of the island. In addition, annual removal data were variable and fragmentary, and abundance was probably underestimated in most cases by unknown amounts. These deficiencies in data most likely resulted in underestimates of overall productivity as well as the true population growth rate. Nonetheless, a repeatable index of abundance could be useful if the design were sensitive enough to reflect changes over time due to removals or mortality from periodic drought-induced starvation. Data on the age and sex structure of the population could be obtained by well-designed observational surveys that are geographically representative to capture annual changes due to recruitment and mortality events. Data on the proportion of pregnant females in the population and proportion of those bearing twins could be obtained from necropsies of animals harvested by hunters or from management removals providing that these data are representative of the entire population. Perinatal survival estimates could be obtained by comparing data from observational surveys of age and sex composition to necropsy data. A well-designed comprehensive population-monitoring scheme is important to determine if management actions result in intended effects. The sex ratio of the population within different age groups and how this relates to population change is a key parameter of any such monitoring.
ACKNOWLEDGMENTS
We thank the staff of Maui office of the Hawai‘i Division of Forestry and Wildlife, the Maui Invasive Species Committee, and the Leeward Haleakala- Watershed Restoration Partnership for providing data. We thank reviewers for helpful suggestions. Model output data are available at https://doi.org/10.5066/P9QXKE7Y (Judge and Hess 2019). Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
