2.1. Study area
Bruny Island is located off the south-east coast of lutruwita/Tasmania, Australia (-43.27S, 147.34E). The island is 362 km2 and is divided into two distinct continental land masses, North and South Bruny, which are connected by a narrow sand isthmus called “The Neck”. Cat control took place on the isthmus, where there is a seabird breeding colony, and adjacent areas to the north and south (which we refer to hereafter as North and South). Whalebone Point (WBP), located in far South Bruny, also supports a seabird colony but was not subject to cat control. We used it as a control site (Fig. 1). The vegetation at both sites is a mosaic of coastal heathland, dry coastal forest, saltmarsh and wetland (sourced TASVEG, the state government GIS vegetation layer (Department of Primary Industries, 2013)). The climate is classified as mild temperate, with four distinct seasons.
The principal terrestrial predators on Bruny Island are the invasive alien feral cat (1.7-4.8 kg, unpublished data for Bruny cats) and the native eastern quoll (0.7-2.0 kg). Other native predators that occur on the adjacent Tasmanian mainland, the Tasmanian devil (Sarcophilus harrisii) and spotted-tailed quoll (D. maculatus), are absent from the island. Cats have been present on the island since the mid-1800s and occur across the entire island. There are two native and two invasive rodents present, which are prey of cats and may prey on seabirds. The rakali/water rat (Hydromys chrysogaster), at 0.8 kg (up to 1.3 kg) Australia’s largest native rodent, forages along the coast and may take some terrestrial prey, particularly seabirds in the breeding colonies. Swamp rats are an omnivorous native rodent, present mainly in denser vegetation and unlikely to prey on seabirds. Black rats and house mice are numerous, particularly in the seabird colonies and may prey on seabirds. At The Neck and WBP seabird colonies, adult shearwaters are present from September to late March. Chicks fledge in early May and leave the island soon after. Penguins are present at The Neck all year around.
2.2. Field study and study design
The field study was repeated in 2017 and 2019, with paired time periods in each year: during the short-tailed shearwater breeding season “BS” and after the shearwater breeding season “PBS”. We deployed a total of 113 cameras stations in 2017 and 111 in 2019, for 24 to 87 consecutive days, respectively, in a grid array with cameras spaced on average 150 meters apart. We placed 93 and 91 camera stations across The Neck shearwater colony and the adjacent areas North and South in 2017 and 2019, respectively (Appendix S1). We placed 20 camera stations on Whalebone Point shearwater colony in 2017 and 2019 (Appendix S1).
Each camera station was set up with a fresh lure secured in a vented PVC canister, suspended 60 cm above the ground. We used lures to maximize the time that cats spent in front of the cameras to aid individual identification. The canister was filled with food attractants for predators (sardines, liver treats, and tuna oil), omnivores (truffle oil and walnut oil) and herbivores (peanut butter and rolled oats).
We monitored the lure stations using passive infrared-triggered cameras (Reconyx Hyperfire Professional PC600 and PC800) programmed to record three consecutive images each time the sensor was triggered, with one second interval between pictures, and no break between sets. Cameras were strapped to a tree or post ~60 cm above ground and ~2.5 m from the lure. The cameras operated 24 hours a day, with most images obtained at night as all species were nocturnally active.
Following the first paired survey in 2017, BICMP undertook an intensive control program of cats at The Neck from July 2017 until July 2019 with the intention to reduce cat numbers around the shearwater colonies with a trapping effort of 4,185 cage trap-nights spread over two years. Cats were controlled on the isthmus and adjacent areas (North and South) as well as 7 km further north along the east coast at Cape Queen Elizabeth, which also has a shearwater colony. The combined area of cat control was much larger than the study area and allowed more effective reduction of cat numbers on and around The Neck. This program followed many years of opportunistic cage trapping and removal of feral cats at The Neck seabird colony during the shearwater breeding season but was much longer in duration, larger in geographic extent and more intensive in terms of the numbers and frequency of trapping. The intensive control program was undertaken using large wire cage traps (710 x 305 x 305 mm) baited with either tuna or chicken. Twenty two cats were trapped and removed (euthanized) during the study period on The Neck seabird colony, only one was removed from the South of The Neck, and none were removed from the site North of The Neck. An extra 20 cats were removed from Cape Queen Elizabeth during the same period.
2.3. Individual cat identification
Individual feral cats were uniquely identified from camera images through a combination of coat pattern and morphological characteristics (McGregor et al., 2015). All discrimination of individuals was done by one experienced observer (Author). We first grouped the images of feral cats as marked or unmarked (black) individuals. Although some black cats had small white patches on the neck or chest, these were not always visible depending on orientation of the image, and so all black cats were considered unmarked to avoid double-counting. The marked portion were tabby cats with naturally unique coat markings. We identified individual cats from this group based on matches in unique markings. Due to poor image quality, some images were considered as unidentifiable.
2.4. Statistical analysis
We estimated the effect of lethal control on the population density of feral cats over the two paired surveys using two different methods. First, we applied spatial mark-resight (SMR) models, a form of spatially explicit capture-recapture (SECR) that is used when only a part of the population is uniquely identifiable (McClintock et al., 2009; Efford and Hunter, 2018). Second, we used dynamic N-mixture models (Royle, 2004) for open populations (Dail and Madsen, 2011) to estimate relative abundance, detection probability and population dynamics of four species that may have been affected by the reduction in the cat population.
2.4.1. Population density estimates of feral cats – spatial mark-resight models
We estimated the density of cats for each of the eight deployments (two sites, two years, paired breeding/post-breeding season - BS/PBS) using a spatial mark-resight model, performed using the “secr” package (Borchers and Efford, 2008) in R version 4.0.2 (R Core Team, 2020). Such models have been used elsewhere in Australia to estimate abundance and density of feral cats (McGregor et al., 2015; Cunningham et al., 2020). To estimate the contribution of unmarked individuals to the overall population, the model assumes that marked and unmarked individuals have identical sighting probabilities (McClintock et al., 2009).
We constructed a set of capture histories for each of the eight deployments by dividing each 3-12 week survey period into a series of five-day intervals for each camera, in which each the number of separate individual cats detected was counted. This approach provided a sufficient number of days within each session to achieve high detection probabilities, while allowing a sufficient number of survey periods to construct encounter histories.
To aid estimation of model parameters, a ‘mask’ has to be chosen as the likely distance from an individual’s home-range centre at which its probability of detection is essentially zero (Balme et al., 2009). In our case, the mask serves several purposes: i) to distinguish habitat sites from non-habitat sites (i.e., the sea) within the outer limit, and ii) to define a region for which a post-hoc estimate of population size is required (i.e., North and South of The Neck). We chose a mask width of 2000 m around the cameras with a spacing of 100 m between the cells, as this is the estimated maximum width of home ranges of male cats in similar habitat type on the island (personal observation from GPS data), clipped to exclude the ocean. This resulted in an area of 20.14 and 7.92 km2 in The Neck and Whalebone Point, respectively.
To estimate the detection function, we first assessed which functional form, half-normal HN or negative exponential EX, best fitted the observed data. The candidate forms differ primarily in the length of the distribution tail, that is, the probability they assign to very distant detections. Models were compared using information-theoretic multi-model inference (Burnham & Anderson 2002), and the detection function with the lowest AIC was used as the basis for further modelling. Next, we constructed biologically plausible models, influenced or not by the seabird breeding season, on the shoulder of the detection function (‘sigma’) and compared these with the null model. Finally, we constructed models to compare estimates of cat densities across: i) The Neck (treatment site where cat removal occurred) and Whalebone Point (control site), and ii) across the area to the North of The Neck seabird colony, The Neck colony, and the area to the South of The Neck colony. We divided The Neck area into three zones to nuance the cat densities of the seabird colony itself on The Neck isthmus compared to the adjacent land to the north and south of the isthmus (Appendix S2).
2.4.2. Estimates of relative abundance of mesopredators – open N-mixture models
We employed the extended Dail-Madsen open population model, also called an open N-mixture model, (Dail and Madsen 2011; Hostetler and Chandler 2015) with a Poisson distribution to estimate variation in demographic parameters of four native and invasive mesopredator and prey species: eastern quoll, swamp rat, black rat, and house mouse. The European rabbit and the water rat were not considered in our analysis because there were too few detections of these species. The Dail-Madsen model relies on temporally and spatially replicated detection histories, which are counts in the N-mixture model (counts of detections need to often exceed one). We defined a detection as independent if it was separated by at least 30 minutes from the next detection of that species at that camera site, as is commonly done in similar studies (e.g., Brook et al. 2012, Cunningham et al. 2019). We created species-specific detection histories for each camera, allowing us to assess factors that may affect species-specific detection. Missing data during a sampling occasion resulted from cameras malfunctioning and was accounted for in the detection histories. The Dail-Madsen model can be described as follow:
where Nit is the number of individuals at site i on survey occasion t, Git is the number of gains (recruits) between seasons, Sit is the number of survivors that do not permanently emigrate, and yit is the observed count at site i on survey occasion t. The four structural model parameters are initial abundance (lambda; λ), recruitment rate (gamma; γ), apparent survival (omega; ω) and detection probability (given presence; p), i.e., the probability of detecting Nit on a single survey. Apparent survival represents animals in site i at time t who survived in that site since time t – 1, and recruitment denotes gains at site i since time t – 1. Note that because individual animal identities are not recorded, no distinction is possible between either deaths and emigrants for the animals departing (apparent survival), or between births and immigrants for new arrivals (recruitment) (Dail and Madsen, 2011).
The four parameters (λ, γ, ω and p) can be modelled as functions of covariates. Due to the large number of covariate combinations possible among the four demographic parameters in the models, our approach was to do the modelling in four stages. In stage 1, we built a set of models to estimate p (see below for model selection). Using the model with the lowest AICc resulting from stage 1, we then built a set of stage 2 models to estimate λ. Using the same procedure, we then estimated γ and ω during stage 3 and 4 of the modelling stages, respectively.
In each of the four modelling stages, model candidates were compared using multi-model inference in an information theoretic framework (Anderson and Burnham, 2004), and all models with ∆AICc > 6 were discarded as non-competitive fits to the data (Richards, 2015). Of the remaining competitive models, the more complex versions of models that had simpler, nested versions with smaller AICc values were also excluded to yield the final model set, and the top-ranked model was used to provide estimates and visualise fitted values (Richards, 2008, 2015). Although we used the top model with the lowest AICc in the successive stages of the best model construction, we used the best model set to describe the relative contribution of each of the covariates for each parameter.
To fit the open N-mixture model, we used the unmarked package (Fiske and Chandler, 2011) in R (R Core Team, 2018), which provides a unified modelling framework for hierarchical models. Data were modelled using maximum likelihood methods with the function pcountOpen, specifically written to use the Dail and Madsen (2011) model in unmarked (Chandler and King, 2011). We estimated seasonal abundance for each species at each site using empirical Bayes methods, using the function ranef from unmarked and demographic parameters from the best-supported model.