
Modeling Impacts of Hunting on Control of an Insular Feral Cat Population^{1}
When introduced to exotic ecosystems, feral cats can inflict irrevers ible harm on native fauna. This is especially true in insular ecosystems because endemic vertebrate species often lack predator defenses. Feral cat control pro grams have been implemented on islands throughout the world with varied success. Effective and responsible management of pest populations requires knowledge of the impact of control actions. Here, we examine a feral cat control program created for the protection of the critically endangered Mariana crow on Rota Island in the Commonwealth of the Northern Mariana Islands. We apply a discrete form of the Schaefer model to a 29month time series of removal data. We use a negative log likelihood framework to determine maximum likelihood parameter estimates and Akaike’s Information Criterion for small sample size (AIC_{c}) analysis to determine the bestfitting model. The model indicated that the removal program on Rota initially reduced cat abundance from an estimated 1,218 to 952 individuals within the first 18 months and then maintained the population near 1,000 individuals for the following 11 months. Given the cur rent level of available funds, we suggest that application of uniform islandwide hunting effort may not be the optimal strategy to maximize crow protection; rather, we suggest a multifaceted, targeted approach focused on areas of high crow activity.
Feral cats have been introduced by humans to most parts of the world including at least 65 major island groups and many remote oceanic islands, both inhabited and uninhabited (Courchamp and Sugihara 1999, Nogales et al. 2004, Campbell et al. 2011). In closed insular ecosystems, introduced cats are known to be the direct cause of severe reduction or extinction of numerous populations of local vertebrate species (Iverson 1978, Taylor 1979, Moors and Atkinson 1984, King 1985, Courchamp and Sugihara 1999). The Mariana crow (Corvus kubaryi) is a critically affected avian species and is endemic to Guam, a USA territory, and Rota, a neighboring island in the Commonwealth of the Northern Mariana Islands, USA. The Mariana crow (hereafter, crow) was extirpated from Guam by the introduced brown tree snake (Boiga irregularis) in the late 1990s (Savidge 1987, Plentovich et al. 2005, U.S. Fish and Wildlife Service 2005) and is listed as critically endangered by IUCN (2016). The brown tree snake is not established on Rota (Rodda and Savidge 2007), however point count surveys showed a 93% decline in the population from 1982 to 2004 (Amar et al. 2008). Recent field evidence from recovered remains of crows using radio [End Page 57] telemetry suggests that feral cats are contributing to the decline of the Rota population (Faegre 2015). Feral cats are known to be the cause of mortality for both adults and juveniles, and it is strongly recommended that feral cat control be integrated into the crow conservation program on Rota (Zarones et al. 2015).
Feral cat control strategies often include techniques such as trapping, hunting, poisoning, exclusion fencing, disease introduction, and fertility control (Algar and Burrows 2004, Nogales et al. 2004, Robertson 2008, Campbell et al. 2011). Factors such as island size, climate, sensitivity of nontarget species, or human habitation can constrain the application of such strategies. Feral cat removal methods on Rota are limited in part because of the presence and biology of the critically endangered crow. As an omnivore and scavenger, the crow would be attracted to poison baits presented to cats and at risk of secondary poisoning via scavenged rat carcasses. Combined, these factors preclude the use of poison as a cat control agent. Fencing is also unfeasible. It would require hundreds to thousands of hectares to accommodate the crow’s mobility and large territory size (Faegre 2015). In addition, the introduction of a feline viral disease is controversial in part because some island residents own pets, but other factors also prevent its application. In the United States, no cat control disease agent is registered, a process that would require extensive regulatory review to evaluate the potential environmental consequences of its release. Live trapping has proven to be difficult on Rota because of interference from multiple nontarget species including rats (Rattus spp.), mice (Mus musculus), coconut crabs (Birgus latro), hermit crabs (Calcinus haigae), chickens (Gallus gallus), and monitor lizards (Varanus indicus). Models have shown that TrapNeuterRelease (TNR) removal strategies are successful at decreasing feral cat populations if >57% of nonneutered cats are captured and neutered annually (McCarthy et al. 2013). The same study showed that TrapVasectomyHysterectomyRelease (TVHR) caused population decline with an annual capture rate of ≥35%. Capture rates such as these could be achieved on Rota only with an extensive, islandwide trapping effort, which would require funds beyond any level of available support. Furthermore, cats would continue to prey on native species, making this method unacceptable. Therefore, two main control strategies have been employed: spotlight hunting and lethal removal via live box trapping.
Because managers on Rota are limited to only two removal techniques, it is important to fully understand how each affects the cat population and how much protective benefit each has for the crow. Here, we determine the effect of spotlight hunting on the feral cat population with the use of a Schaefer model during a 29month time period (Schaefer 1954, Clark 1990). A model parameter of particular importance is initial population size (size of population before control), because it has strong influence on all subsequent abundance estimates. It is not uncommon practice to assume that a population is at or near carrying capacity at the outset of control measures given prior absence of control effort, no recent catastrophic events, and low natural population stochasticity (Hilborn and Mangel 1997, MilnerGulland and Rowcliffe 2007). We provide data that suggest that these three characteristics apply to the population on Rota. We use model selection to determine if the data support this assumption.
The primary objectives of this study are to: (1) estimate population model parameters and test the hypothesis that the population was at or near carrying capacity at the outset of the control program, (2) provide the control program on Rota with information about the effect of hunting on cat abundance, and (3) demonstrate the utility of the Schaefer model for assessment of control projects, especially those without target species data before implementation.
materials and methods
Study Area
Rota (14° 10′ N, 145° 12′ E) is a part of the Commonwealth of the Northern Mariana Islands, USA, and is the second southernmost island of the Marianas Archipelago, a [End Page 58] chain of 15 small islands located approximately 1,500 km east of the Philippine Islands. Rota is 18 km long and 4.8 km wide. The climate on Rota is tropical, with an average annual temperature of 26°C (Amar et al. 2008). Precipitation is seasonal, with a rainy season from July to October, and averages 2,000 to 2,500 mm (MuellerDombois and Fosberg 1998). An uplifted limestone mesa known as the Sabana is located on the western half of the island. Relatively undisturbed forest occurs around the base of the Sabana cliffs and on some coastal shelves with precipitous terrain (Fancy et al. 1999). The island has two villages, Sinaplo and Songsong, located in the center and southwestern regions of the island, respectively. Around 60% of the island is covered in native limestone forest (Falanruw et al. 1989). Feral cats have been observed throughout most of the island (B.T.L., unpubl. data).
Hunting
The data used in this study were collected from late February 2012 to June 2014; however, no hunting was conducted during July 2012 and April 2014 due to permitting and logistic issues (Figure 1). Hunting occurred at night from a vehicle on all public roads that were at least 200 m away from domiciles. Spotlights were used to illuminate cat habitat; once cats were identified they were pursued and shot with either a 10/22 rifle (Ruger) or 410 shotgun (Mossburg). Removal locations were recorded using handheld GPS Navigator units (Garmin GPSMAP 60CSx). Other data collected included time of removal, [End Page 59] pelage, sex, number of fetuses if pregnant, stomach contents, locations, time of locations, and pelages of cats that were observed and not shot. There were multiple hunters employed throughout the project. The project manager started the program in February 2012 and was the exclusive hunter until July 2013. From August 2013 to June 2014, three removal technicians were employed.
Population Model
Our model is designed using methodology outlined in Hilborn and Mangel (1997). We characterized the feral cat population on Rota using a discrete form of the Schaefer model:
where N is the estimated population at time t, r is the intrinsic monthly rate of increase, K is the carrying capacity, and R is the number of cats removed at time t. We use monthly increments of time.
We use a maximum likelihood framework to accomplish several objectives: (1) estimate parameters that provide the best fit to the data, (2) compare hypotheses of initial population size with the use of Akaike’s Information Criterion corrected for small sample size (AIC_{c}), and (3) calculate Confidence Intervals (CIs) of estimated parameters using likelihood profiles. The use of maximum likelihood techniques also allows for the incorporation of different data types into a single framework. Here, we include three data types to inform the model: monthly feral cat removal rates, an independent calculation of initial population, and an independent calculation of intrinsic rate of increase.
Values of abundance N were calculated based on observational data, in our case monthly removal rates from spotlight hunting:
where q is catchability coefficient and I _{t} is number of cats removed per kilometer driven for each month (t), or Catch Per Unit Effort (CPUE). The catchability coefficient is a parameter that relates an abundance index to population size. Here, it is calculated through a maximum likelihood optimization process and can be defined as the fraction of cats removed from the population given one unit of effort. We assume that N _{t} are related to predicted population size with lognormal distributed observation uncertainty (v _{t}); negative log likelihoods (NLLs) were calculated for each time period using the same distribution. Estimates for intrinsic rate of increase r and initial population size N _{0} were calculated independently from the population model with standard deviations σ_{r} and σ _{N}_{0}. NLLs were calculated for both parameters assuming normal distributions. The completed framework is the sum of NLLs for observation error (v _{t}) across all time periods, r, and N _{0}:
where r is an independent estimate of intrinsic rate of increase, r̂ is intrinsic rate of increase to be estimated by NLL minimization, N̂_{0} is an independent estimate of initial population size, and Nˆ_{0} is initial population size to be estimated by NLL minimization. Given the data and particular values of q, r, K, N _{0}, and σ _{v}, the likelihood of that set of parameters was evaluated. Here we select the parameters that make the NLL as small as possible and call these bestfit parameters. This was done using the optim function in R version 3.1.2 (R Core Team 2014).
To reduce bias in catchability estimates because of variability among hunters, maximum likelihood estimates of parameters were initially optimized using the first 18 months of data, because there was only one hunter during that time period. For the remainder of the removal data, a second catchability coefficient and observation error standard deviation were estimated, keeping the r and K estimates from [End Page 60] the initial optimization constant, because a difference in hunting personnel should have no influence on those parameters.
Input Parameters and Assumptions
initial population size, N_{0}. An initial island population estimate was based on extrapolation of cats observed (removed and sighted) from hunting trips in March 2012; we did not use data from February 2012 because hunting began late that month, and hunting methods were still being developed during that time. We first calculated cat density along spotlighting roads. Maximum cat detection distance (the maximum distance a cat could be seen from a road) perpendicular to all spotlighting roads was measured with a range finder (Bushnell sport 450). The areas were summed to obtain a total sampling area. Total number of cats observed in the sampling area was calculated by adding the number of cats removed to the number of cats that were observed and not shot, so long as there was spatial evidence to support that they were unique. A cat was not considered unique if it was seen within 2.4 km of another cat with the same pelage that was not previously shot. A distance of 2.4 km was used because that was the maximum distance a cat was observed to travel in one 24 hr period (Leo et al. 2016). Cat per unit area was calculated for the sampling area and then extrapolated to the rest of the island. Areas from the two villages were excluded because spotlighting was always conducted 200 m away from domiciles. A coefficient of variation (CV) was calculated from the mean and standard deviation of logtransformed unique cats seen per day (n = 20 days) and was used to calculate σ _{N}_{0} in equation (3). The extrapolation was computed using ArcMap 10.1 software. We recognize that this methodology assumes that cat density is homogeneous on Rota. It is widely accepted that predator densities can be higher along roads because roads facilitate movement and increase foraging opportunity; however much of this evidence exists for large animals in areas with dense vegetation, and for species with generalist requirements such as foxes (Vulpes vulpes) (Andrews 1990, Bennett 1990, May and Norton 1996). It is not known whether feral cat densities are lower or higher along roadways than in neighboring habitats (May and Norton 1996). To address this we incorporate the initial population estimate into the likelihood framework with a CV, as a parameter to be optimized.
intrinsic rate of increase, r. We used fetus count data from cat removal necropsies (n = 44), a sex ratio calculated from removed cats on Rota, average number of litters per month per female (Hurni 1981), and average age (in months) at first reproduction (Ernest 2003) to calculate the intrinsic rate of increase (r) defined (Hone et al. 2010):
where b is monthly fecundity and α is age (in months) at first reproduction. A CV was calculated from fetuses per pregnant female and was used to calculate σ _{r} in equation (3).
catchability, q. For each time period, February 2012 to July 2013 and August 2013 to June 2014, a unique catchability was estimated. The parameter q _{1} denotes catchability for the first 18month time period, and q _{2} denotes catchability for the remainder of time.
Model Selection
We tested two models that represented competing hypotheses: that the initial cat population on Rota was at or near carrying capacity at the outset of control measures, and that the initial population was not equal to carrying capacity. We used AIC_{c} because the models contain different likelihoods. We included all constants within likelihood calculations and concluded that the model with the lower AIC_{c} is the better model.
Confidence Intervals (CIs) and Likelihood Profiles
We report CIs and likelihood profiles only for parameters of interest and those with least bias: K, r, and q _{1}. Briefly, the maximum likelihood estimates were calculated using the likelihood framework described above. Each profile was created by systematically changing the parameter of interest, then computing the values of the other parameters that minimized [End Page 61] the NLL. The parameter values used for K ranged from 500 to 5,000 in increments of 1. Parameter values for r ranged from 0.001 to 0.15 in increments of 1.0e04. Parameter values used for q _{1} ranged from 4.50e06 to 6.2e05 in increments of 6.2e08. Using the likelihood ratio test, the 95% CI is the range of parameter values for which the NLL is within 1.92 of the minimum value of the NLL.
Sensitivity Analysis
The sensitivity of model output to ranges of parameters r, K, and q _{1} was determined using the sensitivity, Hmisc, ks, and pse packages in R (Chalom et al. 2013). Values were randomly selected from a uniform distribution constrained by 95% CIs calculated from likelihood profiles. For each parameter combination, the model was run for 18 time steps, and output was recorded. The number of time steps was chosen based on the length of time the project operated with only one hunter. Model output was then plotted as the dependent variable for each parameter value. Sensitivity was determined by calculating the Partial (Rank) Correlation Coefficient (PRCC), which measures the strength of linear associations between the result and each input parameter after removing the linear effect of the other parameters.
To ensure sufficient sample size, we calculated the Symmetric Blest Measure of Agreement (SBMA) between the partial correlation coefficients of two runs with different sample sizes. We used the sample size that resulted in a SBMA value of 1, which indicates an adequate sample size (Chalom et al. 2013).
results
From February 2012 to June 2014, 530 cats were removed via spotlight hunting. From February 2012 to July 2013 (the first optimization time period) only 23 cats were removed via live box trapping. Box trap removals are not included in the model because of the relative paucity (one per month) and because trapping effort data were measured in trap nights, not kilometers driven. The amount of information gained from the trapping data did not warrant the addition of another parameter into the modeling framework.
In March 2012, 49 unique cats were seen within a sampling area of 1.83 km^{2}. The extrapolation of these values returned an initial population estimate of 2,286. The CV calculated from the mean and standard deviation of logtransformed unique cats seen per day (n = 20 days) was 0.63. Mean fetus count was 3.32 (n = 44, SD = 0.98, CV = 0.30). Using all removal data, we failed to reject the hypothesis that the sex ratio of cats on Rota was 0.5 (Exact Binomial Test, n = 417, P value = .33, 95% CI: 0.47 to 0.57). Our independent estimate for monthly r was 0.04.
The AIC_{c} analysis supported the model that equated initial population size to carrying capacity and the hypothesis that the feral cat population was at or near carrying capacity at the onset of control measures (Table 1). We therefore used that model for all subsequent analysis. Bestfit parameter estimates and 95% CIs (Figure 2) are listed in Table 2. The optimization process converged to an initial population estimate of 1,218 in February 2012. The model suggested that this number dropped to 952 by July 2013 and then increased from 945 in August 2013 to 1,027 in June 2014 (Figure 3). The SBMA calculations indicated that 300 was an adequate model run sample size to assess model sensitivity. The model was most sensitive to K and q _{1} and marginally sensitive to r (Figure 4). PRCCs were higher for K and q _{1} than for r (Table 3).
[End Page 62]
discussion
Based on the AIC_{c} analysis, we found support for the model that equated initial population size to carrying capacity, and less support for the model that estimated the two parameters separately (Burnham and Anderson 2003). The assumption that the population was at carrying capacity at the onset of control measures is a large one; however, the ecosystem on Rota possesses attributes that would permit such an assumption. Feral cat reproductive cycles are influenced by length of daylight (Hurni 1981), of which there is very low variability on Rota. Indeed, feral cats were documented (via necropsies) to be pregnant yearround. A constant rate of production should facilitate resiliency to stochastic events, supporting the hypothesis that the population existed near carrying capacity at the beginning of the project. In addition, the productive Rota ecosystem provides abundant food resources that could sustain a stable feral cat population at carrying capacity. Rota supports a high density of rodents, which are the primary diet component of the local cat population (Wiewel et al. 2009, Leo et al. 2016). Of 174 stomachs examined from cat necropsies in 2012, 150 (86%) contained prey items, indicating that it is not difficult for cats on Rota to find food resources. In summary, given the constant reproductive cycle of feral cats on Rota, access to abundant prey resources, and the almost complete absence of predators, [End Page 63]
[End Page 64] we suggest that it is reasonable to assume that the population was at or near carrying capacity upon the commencement of the removal project.
The model was marginally sensitive to r. This implies that the accuracy of the r parameter is not critical when using the model for shortterm inferences, as long as the r estimate is reasonably well estimated and supported with field data. The insensitivity is also due to the number and magnitude of timestep intervals. If the time increments were annual, then the model would be more sensitive to r. It would be prudent to further investigate the reproductive characteristics of the Rota population if the model were to be used for annual abundance projections.
The extreme sensitivity to K can be explained because we equated N _{0} to K for all model runs. We therefore would expect to see a linear relationship between model outcome and K. The model was also sensitive to q _{1}. This parameter should be taken into careful consideration if the model is to be used for future management. Catchability can vary considerably among personnel and should be estimated accordingly, especially because of its strong influence on model output.
The independent initial population estimate calculated by extrapolating the density of unique cats seen within the hunting area was coarse. The continuous nature of the jungle adjacent to spotlighting roads made it exceedingly difficult to determine an accurate measure of the total sampling area. Cats were unpredictable and difficult to detect, which produced a large CV of 0.63 for unique cats seen in March 2012. Indeed, the estimate of 2,286 initial cats was added into the likelihood framework with a standard deviation of 1,440. This nonetheless informed the model, and the initial population estimate that was produced from the optimization (1,218) appeared to be more realistic because it was more consistent with another estimate for population size, 1,512, calculated using a removal method (Hayne 1949); however, this was likely an overestimate because the removal method assumed that data collection occurred over a small period of time (no reproduction).
Even if the true cat population size on Rota were closer to 2,000 instead of 1,000, our inferences about the effectiveness of the removal program would be similar because in either case the remaining cat abundance is high. If the initial population was 1,218, then the program managed only to marginally lower and then maintain the population at approximately 1,000 individuals, a reduction of only 15%. These findings indicate that until more effective islandwide removal of cats can be implemented, effort would be best spent protecting targeted areas of high crow activity or density, perhaps with a combination of targeted trapping and hunting, as opposed to a steady islandwide application of hunting effort.
The current form of the model provides information about population parameters but needs refinement to be useful for abundance projections. Given the precision of the estimated parameters, uncertainty in model predictions will be high and envelop differences in projected abundance levels under varying levels of effort. The model would benefit from further research into the variability of cat density with habitat type, which would refine the initial population estimate. Estimates for natural mortality rates of kittens and juveniles would give a more accurate representation of recruitment. These improvements would allow for further lines of inquiry and perhaps inference into what levels of removal would be needed for sustained absence of catrelated crow mortality.
conclusions
Historically, the removal program on Rota had few criteria with which to apportion spotlight effort levels. New understanding of the effect of spotlight removal on cat abundance [End Page 65] gives reason to consider alternative management strategies different from the uniform effort approach that has been implemented thus far in the program. Spotlight removal of feral cats on Rota has clear advantages over other removal techniques; however it should be used in a strategic fashion and in concert with trapping, which can target cats known to occupy highvalue nest areas. The model showed that applying a constant, islandwide spotlighting effort is not necessarily the optimal strategy. The analysis represents a straightforward application of the Schaefer model. Although morecomplex model population forms are available, we suggest that the Schaefer model is well suited for assessing the impacts of limited predatorcontrol programs, such as the one conducted on Rota Island. We suggest that useful information on the effectiveness of the program could be obtained with survey and modeling of the response of the prey species of concern, the Mariana crow in our example.
Corresponding author (email: Brian.Leo@ colostate.edu).
acknowledgments
We thank the U.S. Fish and Wildlife Service, Pacific Island Fish and Wildlife Office for funding the feral cat control program and this research. We thank the Commonwealth of the Northern Mariana Islands Department of Land and Natural Resources Division of Fish and Wildlife on Rota, and the Guam Department of Agriculture’s Division of Aquatic and Wildlife Resources. A special thanks to the Anderson and Gallucci Laboratories in the SAFS and QERM Departments at UW, as well as the Ha Animal Behavior Psychology laboratory, specifically Aneesh Hariharn, Chloe Bracis, Jeff Rutter, Kelsey Vitense, and Anton Mates. We also thank the field techs: Alex Reyer, Daron Hartman, Brian Holt, Otton Mendiola, and Dane Horowski; and our Mariana crow collaborators: Sarah Faegre, Phil Hannon, Jen Carpenter, Jose Antonio Diaz, Andria Kroner, Kelsey McCune, Kelle Urban, and Lydia Goy.
Literature Cited
Footnotes
1. Funding for this research provided by the U.S. Fish and Wildlife Service, Pacific Island Fish and Wildlife Office. The findings and conclusions in this article are those of the author(s) and do not necessarily represent the views of the U.S. Fish and Wildlife Service. Manuscript accepted 1 August 2017.