Monitoring, modeling and harvest management of non-native invasive green iguanas on Grand Cayman, Cayman Islands

The green iguana (Iguana iguana) was most likely introduced as a pet and became overabundant during the last 20 years on Grand Cayman. Because negative impacts were unmanageable (e.g., damage to buildings and other infrastructure), a harvest management strategy was developed and implemented, and over 874,252 green iguanas were removed between October 2018 and August 2019. Distance sampling surveys were conducted to estimate abundance in February 2019 and annually in August 2014–2019. Abundance estimates were used to develop a Bayesian state-space logistic model, generate the posterior distributions of population and harvest management parameters, and make future predictions of abundance for August 2020–2030. Abundance increased over fivefold between August 2014 and 2018, from an average of 254,162 to 1,319,939 green iguanas. However, after harvesting for 5 months, abundance declined to an average of 600,113 green iguanas in February 2019; and after 11 months, abundance declined to an average of 103,020 green iguanas in August 2019. Maximum population growth rate averaged 1.552, carrying capacity averaged 1,611,013, equilibrium abundance averaged 805,506, maximum sustainable total harvest averaged 628,491, and maximum sustainable harvest rate averaged 0.776. With harvest rates between 0.600 and 0.800, predicted abundance averaged 28,751 green iguanas for August 2020–2030. However, harvest mortality may have unforeseen outcomes due to the release from density dependence and overcompensation through high survival and fecundity rates. Because natural resource managers have partial control over harvesting and incomplete understanding of green iguana population dynamics, monitoring and modeling are essential to assess population response and guide harvest management decisions.


Introduction
The green iguana (Iguana iguana) was most likely introduced as a pet on Grand Cayman, Cayman Islands (Echternacht et al. 2011). First sightings in the wild occurred in the early 1990s (Burton, personal observation), with green iguanas becoming overabundant during the last 20 years (Rivera-Milán et al. 2014a, 2016avan Veen et al. 2015; Department of Environment [DoE] 2018[DoE] , 2019. Although quantitative data about negative impacts are lacking for Grand Cayman, green iguana overabundance caused economic losses (e.g., damage to buildings and other infrastructure in urban areas), posed health and safety and risks, and triggered negative interactions with native plants and animals (Engeman et al. 2005;Krysko et al. 2007;Meshaka et al. 2007;Falcón et al. 2012;van Veen et al. 2015;Vuillaume et al. 2015;Moss et al. 2017). When eradication is not economically feasible due to low detection probability and high population growth rate, natural resource managers may decide to harvest an overabundant population above the maximum sustainable yield to reduce and maintain abundance at acceptable levels during management cycles (Caughley and Sinclair 1994;Mack et al. 2000;Baxter et al. 2008;Runge et al. 2009;Alisauskas et al. 2011;Koons et al. 2014;Pasko and Goldberg 2014;Anderson et al. 2016). For the harvest management of green iguanas on Grand Cayman, natural resource managers prescribed harvest rates above 0.600, with 2.3 million harvested by August 2020, and abundance levels below 50,000 by August 2027 (DoE 2018). However, natural resource managers were uncertain about the effectiveness of harvest management at low abundance due to diminishing returns. In addition, natural resource managers had an incomplete understanding of green iguana population dynamics, and recognized that harvest mortality could have unforeseen outcomes due to the release from density dependence and overcompensation through high survival and fecundity rates (Courchamp et al. 2003;Krysko et al. 2007;Meshaka et al. 2007;Zipkin et al. 2009;Alisauskas et al. 2011;Koons et al. 2014;Pasko and Goldberg 2014;Cortez and Abrams 2016;McIntire and Juliano 2018). Under these circumstances, monitoring and modeling were deemed essential to assess population response and guide harvest management decisions (Nichols and Williams 2006;Lyons et al. 2008;Runge et al. 2009;Lindenmayer and Likens 2010;Lindenmayer et al. 2011;Pasko and Goldberg 2014;Rivera-Milán et al. 2014b, 2016b.
Not accounting for individuals wounded but not retrieved for processing (i.e., crippling loss), registered hunters shot with air guns or trapped with snares and other devices over 874,252 green iguanas between October 2018 and August 2019 (for more information, see https://doe.ky/green-iguana-cull-updates/). Assuming the peak of reproduction occurred between April and July, when courtship and other reproductive behaviors were observed (e.g., dominant males defending territories and nesting females), distance sampling surveys were conducted to estimate abundance in February 2019 and annually in August 2014-2019, using methods that accounted for changes in detection probability (Buckland et al. 2001(Buckland et al. , 2015Royle 2004;Thomas et al. 2010;Thomas and Marques 2012;Burton and Rivera-Milán 2014;Rivera-Milán et al. 2014b, 2016b. To assess population response and guide harvest management decisions, distance sampling abundance estimates in August 2014-2019 were used to develop a Bayesian state-space logistic model, generate the posterior distributions of population and harvest management parameters (i.e., maximum population growth rate, carrying capacity, equilibrium abundance, maximum sustainable total harvest and maximum sustainable harvest rate), and make future predictions of abundance with harvest rates from 0.100 to 0.800 for August 2020-2030, using methods that accounted for the variation associated with observation and process errors (Millar and Meyer 2000;Kéry and Schaub 2012;Hobbs and Hooten 2015;Rivera-Milán et al 2014b, 2016b.

Monitoring
Distance sampling was used to estimate density (i.e., D = number of individuals/ha) and make abundance inferences for the whole island (i.e.,N = number of individuals in 19,600 ha; Fig. 1). Six-minute counts were conducted by teams of three observers at a minimum of 157 and a maximum of 212 points established along and away from roads using systematic sampling with random starts (Rivera-Milán et al. 2014a, 2016a. Because on-road points may not be representative of green iguana abundance at interior habitats, between 50 and 75 points (32-35%) from the total of 157-212 points were located away from roads. Off-road points were separated by a minimum distance of 400 m and on-road points were separated by a minimum of 800 m. From the 157-212 points that were visited once, random subsamples of 50 points were visited twice within 2-week survey periods in February 2019 and annually in August 2014-2019.
Survey effort accounted for the number of points and the number of visits per point (Buckland et al. 2001). The surveys included paved and unpaved roads in urban and rural areas, forest and farm trails, as well as remote areas with little or no disturbance from human activity. Accessibility to mangrove forests, wetlands and farms was limited, particularly in the central and eastern portions of the island (Fig. 1). The surveys also included areas of special interest, such as the Queen Elizabeth II Botanic Park and Colliers Wilderness Reserve in the east, where green iguanas may have negative interactions with the endangered endemic blue iguana (Cyclura lewisi: Burton and Rivera-Milán 2014, Conley et al., in press); and Safe Haven Golf Course and urban developments in the west, where green iguanas caused economic losses and were harvested before and during 2014-2019 (Rivera-Milán et al. 2014a, 2016a. Detection probability of the ith single green iguana or cluster of green iguanas was modeled as a function of distances r i and vector z i representing the following categorical and continuous covariates: sampling period (February 2019 and August 2014-2019; categorical 0-6), age class (hatchling [small] or juvenileadult [medium-large]; categorical 0 or 1), time of day (08:30-16:30; continuous in hours and minutes), point location (on road and off road; categorical 0 or 1), disturbance and visibility (none-low or moderatehigh; categorical 0 or 1), cluster size (C 2 individuals within 10 m of each other, displaying similar behavior; continuous 2, 3, 4, …), and teams with three observers (2 of them mainly measuring detection distances with binocular rangefinders and the third mainly recording the data; categorical 0, 1 or 2). Points were surveyed at varying times during the morning and afternoon. A quadratic effect was also considered for time of day, because detection probability was expected to increase as the morning warms up and decline as the afternoon cools down. Density was estimated aŝ where n is the total number of detections; s is the average cluster size, which was used when cluster detection was not size biased; k is the total number of points; and detection of available green iguanas within right-truncation distance w iŝ rĝðr; z i Þdr. Expected cluster sizeÊðsÞ was used instead of s when cluster detection was size biased; or cluster size was included as a covariate in the models (Buckland et al. 2001(Buckland et al. , 2015. Detection probability had two components; that is,P da ¼P dPa , whereP d was the probability of an observer detecting an available single green iguana or a cluster of green iguanas (''detectability''), andP a was the probability of a single or a cluster being available for detection (''availability''). Repeated point counts were used to estimate availability and its standard error, which were included as a multiplier in the density estimator (Royle 2004;Thomas and Marques 2012;Burton and Rivera-Milán 2014). The fit of uniform, half-normal and hazard-rate detection models was evaluated with Cramer-von Mises family tests (Buckland et al. 2015). Model selection was based on minimization of Akaike Information Criterion (AIC). The half-normal and hazard-rate detection models included the above mentioned categorical and continuous covariates. After the 6-min counts, the observers moved around point centers to verify detection distances and cluster sizes, and reach a consensus about disturbance from traffic and other human activity, as well as visibility from vegetation cover and other visual obstructions. The basic assumptions of distance sampling were: (1) green iguana distribution was independent of point location; (2) certain detection of singles and clusters at point centers; (3) accurate measurement of detection distances and cluster sizes; and (4) no responsive movement, away or toward the observers. To relax the third and fourth assumptions, detection distances were grouped as: 0-15, 16-30, 31-45, 46-60, 61-90, 91-120, 121-180, 181-240, 241-340, 341-440, [ 440 m (Rivera-Milán et al. 2014b, 2016b. DISTANCE version 7.2 was used for data analysis (Thomas et al. 2010). Nonparametric bootstrapping was used to estimate standard errors (SEs). Results are presented as means with bootstrapped SEs and lognormal 95% confidence intervals (CIs),

Modeling
Details of the Bayesian state-space logistic model are presented in previous publications (Rivera-Milán et al. 2014b, 2016b. Briefly, annual changes in population state N t were modeled with where r max is the maximum intrinsic rate of population growth, K is the population carrying capacity, N t is the true unknown state of the population and H t is the total number of green iguanas harvested in time period t. where h t is the harvest rate, which was generated as part of the Markov chain Monte Carlo (MCMC) algorithm using the uniform distribution (e.g., h * uniform [0.100, 0.600] or [0.601, 0.800]). The unknown population state was reparameterized as a proportion of population carrying capacity N t /K to reduce autocorrelation of MCMC samples. Based on this reparameterization, the state dynamics were projected forward in time according to Model prediction error e was assumed to be lognormally distributed with mean 0 and standard deviation r process . Assuming linear density dependence, the following parameters were derived as where h msy is the maximum sustainable harvest rate, N eq is the population equilibrium abundance, and H msy is the maximum sustainable total harvest. Because natural resource managers opted for a simple strategy, with harvest rates above 0.600 and abundance levels below 50,000 for age and sex classes combined, the model was simplified by assuming that harvest mortality occurred after the reproductive peak (presumably April-July) and that age and sex classes had equal probability of being shot or trapped in time period t. In addition, additive natural and harvest mortality was assumed, although the model formulation allowed for compensation through density dependent population growth. Uniform prior distributions were used for maximum population growth rate (r max * uniform [0.001, 2.000]), population carrying capacity (K * uniform [200,000, 2,000,000]), the mean of the initial population proportion on the logarithmic scale (P 0 * uniform [-10, 0]), process and initial population proportion SD (r process and r P 0 * uniform [0, 10]).
MCMC methods were used by running JAGS version 4.3.0 (Plummer 2003) with R2JAGS version 0.5.7 (Su and Yajima 2015). Trace plots and node summary statistics were used to check for convergence of the MCMC algorithm (Gelman et al. 2004). After a burn-in period and thinning of three Markov chains, results from 8000 samples are presented as means and MCMC standard deviations (SDs) with 2.5%, 50.0% and 97.5% percentiles. Model-based predictions were used to estimate outcome probabilities. That is, the probabilities of having abundances below 50,000 green iguanas with harvest rates between 0.100-0.600 and 0.601-0.800 in 2020-2030.
Abundance estimates for the age classes separated and combined are presented in Table 2  As expected from the logistic model, there was a linear relationship between estimated population growth rate and abundance in August 2014-2019 (Fig. 2b). Population and harvest management parameters generated with the Bayesian state-space logistic model are presented in Table 3. With harvest rates between 0.100 and 0.600, future abundance was predicted to be 1,063,962 green iguanas (SD = 605,525,with 2.5%,50.0% and 97.5% percentiles = 28,427,1,096,332 and 2,186,599,respectively); and the probability of having abundance below 50,000 green iguanas was 0.036 for August 2020-2030. However, with harvest rates between 0.601 and 0.800, future abundance was predicted to be  28,751 green iguanas (SD = 65,407, with 2.5%, 50.0% and 97.5% percentiles = 532, 11,227 and 166,706, respectively); and the probability of having abundance below 50,000 green iguanas was 0.626 for August 2020-2030 (Fig. 2a).

Discussion
Age class and cluster size were the most important detection covariates. Not surprisingly, juvenilesadults were more detectable than hatchlings. For both age classes, large clusters were more detectable at longer distances than small clusters or singles. However, the reliability of distance sampling abundance estimates depended on the ability of observers to meet basic method assumptions. Teams with three observers did not likely miss green iguanas available for detection at point centers during 6-min counts. Green iguanas did not show much responsive movement to the approaching observers. Moving green iguanas were not included in the analysis, unless detection distances were measured to their initial locations. Distance categories were used to relax the assumptions pertaining to distance measurements. Although cluster detection was size biased, expected cluster size was used instead of average cluster size in conventional distance sampling, and cluster size was included as a covariate in multiple-covariate distance sampling. Points were located along and away from roads to cover low, medium and high abundance areas across the island, independently of green iguana spatial distribution. Therefore, the basic assumptions of distance sampling were likely met, and reliable abundance estimates were obtained for monitoring and modeling green iguana population dynamics.
Although a number of simplifying assumptions were made (e.g., linear density dependence and equal harvest mortality for age and sex classes after a reproductive peak), the Bayesian state-space logistic model was pragmatic, problem-oriented, used monitoring data efficiently, and sufficed to represent green iguana population dynamics under varying harvest rates (Starfield 1997;Nichols and Williams 2006;Lindenmayer and Likens 2010;Lindenmayer et al. 2012;Rivera-Milán et al. 2014b, 2016b. The model generated useful baselines for population and harvest management parameters, and facilitated future comparisons between estimated and predicted abundances, despite the lack of data on sex-specific and agespecific demographic rates. For example, with harvest rates above 0.600, the population declined sharply between August 2018 and 2019; and simulations showed that abundance can be maintained below 50,000 green iguanas during the modeled time horizon in August 2020-2030. However, with harvest rates below 0.600, the population would likely rebound and become overabundant again in a short time period (e.g., doubling time t = ln[2]/1.552, or about 5 months, based on mean r max ).
Maximum population growth rate denoted the potential rate of increase under ideal environmental conditions (e.g., abundant food and space to reproduce). The model assumed additive mortality from natural and anthropogenic disturbances, but allowed for compensation through density dependent population growth (Runge et al. 2009;Rivera-Milán et al. 2014b, 2016b. That is, when unharvested, the green iguana population would be expected to increase rapidly from low abundance, slow down and fluctuate around carrying capacity, with density dependent factors limiting and regulating further growth (e.g., see Runge et al. 2009: Fig. 1). At low abundance, the age and sex classes can survive and reproduce more or less equally through resource partitioning; but at high abundance, dominant territorial males and nesting females would likely partition resources unequally through contest competition. Because of green iguana life-history characteristics (e.g., rapid growth rate of hatchlings, and early maturity and high fecundity rate of females; see Meshaka et al. 2007), harvest mortality may result in overcompensation due to a release from density dependence, increasing rates of survival and fecundity with more foraging and nesting resources available in less crowded space (e.g., see Zipkin et al. 2009: Fig. 1).
Harvesting to reduce and maintain abundance below 50,000 green iguanas is a management action but not necessarily a management objective in itself, particularly when the ultimate goal is to minimize economic losses, health and safety risks, or to minimize negative interactions with native flora and fauna (for a review of harvest theory and management, see Caughley and Sinclair 1994;Runge et al. 2009). To that end, the success of harvesting green iguanas, or any other non-native invasive species, should be measured in terms of the cost and benefit accrued from the action taken (e.g., see Engeman et al. 2002). However, for r-selected species, such as green iguanas, harvesting can be ineffective when directed towards decreasing survival instead of decreasing fecundity rate, assuming population closure to immigration and additional introductions (e.g., see Meshaka et al. 2007). Moreover, apart from airports, golf courses and other areas of special interest (e.g., see Engeman et al. 2005), the cost may exceed the benefit of harvesting green iguanas on islands larger than Grand Cayman, unless there are well-established recreational and commercial incentives in addition to bounty programs and contracts to service providers (for a review of harvest incentives, see Pasko and Goldberg 2014).
In conclusion, because natural resource managers have partial control over harvesting and incomplete understanding of green iguana population dynamics in a changing island environment, abundance predictions should be taken only as suggestive of population response to harvest management. Reduction and maintenance of abundance below 50,000 green iguanas by August 2027 seemed feasible with adequate funding for sustained harvesting. The initial population ''knockdown'' accomplished between August 2018 and 2019 should be followed by sustained harvest rates above 0.600 until August 2020. Subsequently, natural resource managers may decide to harvest the population during shorter management cycles (e.g., September-January), allocating harvest effort (i.e., the number of hunters and shooting hours/day and days/cycle) to specific areas, based on estimated and predicted abundances in February and August 2020-2030. The harvest management strategy implemented can be adapted and allow learning from monitoring and modeling changes in green iguana abundance before and after the peak of reproduction on Grand Cayman.