Study area
Data was collected in the Ecological Park of Funchal, on Madeira Island, Portugal (Fig. 1). Madeira is the largest island of the Madeira archipelago, encompassing Madeira, Porto Santo, and the Desertas Islands. The island is an autonomous region of Portugal, situated around 900 km from the mainland. A census in 2023 recorded the island’s population at over 250,000 people, of which ca. 100,000 were living in Funchal (INE, 2023). The island is characterised by a subtropical climate, with mean daytime temperatures between 20 and 26\(\:^\circ\:\)C, and varying levels of rainfall, accumulating from 2 mm in July to 115 mm in December (Weather Atlas, 2024). The island features an altitudinal progression of distinct vegetation belts, including xerophytic, sclerophyllous flora; the UNESCO-protected Laurisilva forest, and high mountain shrubland (Capelo et al., 2005).
The Ecological Park of Funchal (7.5 km2) is located North of Funchal capital city, between 470–1818 m asl and was established in 1994. It is an Important Bird Area, a Site of Community Importance, a Special Protection Zone, and is listed in the Natura 2000 Network (Silvestre, 2011). It is also a popular tourist destination, attracting many outdoor enthusiasts. Its vertebrate fauna includes over 25 confirmed breeding bird species, including both Madeiran (e.g., Regulus madeirensis and Fringilla maderensis) and Macaronesian (e.g., Anthus berthelotii) endemism. Manx shearwaters (Puffinus puffinus) are severely depleted, and the first two active nests were discovered in 2005 in the Santa Luzia valley, within the Ecological Park of Funchal (Nunes et al., 2010). In addition to free-ranging cats, vertebrate diversity includes three species of bats, one reptile (Teira dugesii) and non-native rodents (Rattus rattus, R. norvegicus and Mus musculus), rabbits (Oryctolagus cuniculus) and ferrets (Mustela furo) (Gonçalves et al., 2024; Soto et al., 2023 and 2024). Cat-human dependency exists on a continuous spectrum, generating ambiguity in exactly when a free-ranging cat becomes “feral”. Here, we refer to “free-ranging cat” as any owned or unowned cat that has some degree of outdoor access (Crowley et al., 2020; Farnworth et al., 2010).
Camera-trapping surveys
Camera-trapping took place on two occasions: August-November 2021 and September 2023. The 2023 survey replicated the design of Soto et al. (2023), who divided the park into eleven 1 km2 grid cells (Fig. 2). The chosen area for each grid reflected average home-range sizes for cats in comparable conditions (López-Jara et al., 2021; Zhang et al., 2022). Three camera-traps were deployed in each grid, totalling 33 camera locations. As in Soto et al. (2023), two out of three camera locations were strategically deployed in known cat pathways, whereas the third location was randomly selected using QGIS, Version 3.30.2 (Sherman, 2023). Cameras were largely deployed at the same locations in both surveys. However, three camera locations (i.e., locations 3.3, 5.2 and 5.3—see tables S1 & S2) had to be moved between surveys due to changes in the landscape (e.g., landslides and management interventions), that rendered the original sites inaccessible. Cameras were placed at least 300 m apart within each grid cell, elevated 50 cm above the ground, and attached to wooden stakes or suitable tree trunks. The cameras were locked to their attachment points with a braided steel cable lock, and a sign was attached explaining the purpose of the survey to deter theft (Clarin et al., 2014). Most of the cameras were from the model Browning Dark OPS Pro X (29 cameras), but additional models included Browning 2020 Patriot (three cameras) and Bushnell Core (one camera). As in Soto et al. (2023), cameras were set to take three photos when triggered by an infrared motion sensor (at normal sensitivity), with minimum interval between photos and a 30 second interval between photo sets. To avoid false triggers and improve detectability, we cleared obstructing vegetation within a 5 m radius around the camera. Cameras were baited with a mixture of chorizo and oil, identical to the bait used by Soto et al., (2023). In 2023, all 33 cameras were left simultaneously active in the field, drawing a contrast to Soto et al., (who had a limited number of camera traps and therefore employed a more clustered study design). However, sampling effort was equally distributed across surveys, amounting to 582 and 588 trap-nights in 2021 and 2023, respectively.
Camera detections and identification
Photos were manually analysed using Timelapse Image Analysis software, version 2.3.0.8 (Greenberg, 2023). Following Soto et al. (2023) and several other camera-trap studies (e.g., Garvey et al., 2022; Khan et al., 2023), we considered detections to be independent when animals were recorded more than 30 minutes apart, or, in the case of cats, when different individuals were captured (Wearn & Glover-Kapfer, 2019).
Wherever possible, cats were identified to the individual level from their pelage and morphological features, such as wounds or notched ears, in line with previous studies (Elizondo & Loss, 2016; Hohnen et al., 2023). Any cat record that could not be identified at the individual level (due to poor photo quality) was excluded from the cat density or abundance estimates. Although we considered kittens in both density and abundance estimates, this likely did not constitute recruitment during the survey nor violate model assumptions of a closed population. Kittens are altricial for the first three weeks of life, so would not be able to roam on their own, as all kittens in our images were observed to be doing (Spotte, 2014). Previous studies excluding kittens from density estimates had a longer duration (several months to years), with more cause to assume that a kitten could be born and reared during the sampling period (Juhasz et al., 2022).
Explanatory variables
In both 2021 and 2023, we measured several abiotic and biotic variables that were hypothesised to influence cat abundance and activity (Table S6), including the cover of (1) open/rocky area, (2) native shrubland, (3) non-native shrubland, (4) native trees, and (5) non-native trees. Each of these metrics were calculated for a circular area with a 300 m radius around each camera trap using orthophotos from the Regional Directorate for Spatial Planning in Madeira (IRIG, 2018). We also measured the (6) distance to a water source (water-dist), described as the two main water courses running through the park: the Ribeira das Cales and the Ribeira de Santa Luzia. This was done using QGIS, with GIS layers supplied by coordinating authorities of the protected area (IRIG, 2018). We recorded the (7) altitude and (8) slope at each camera trap location using a digital elevation model with a 10 m resolution (IRIG, 2018). Relating to human influences on cat ecology, we measured the (10) distance to urban area (urb-dist) and (11) distance to human resource subsidies (food-dist), defined as any area where cats could access anthropogenically sourced food (e.g., picnic benches, rubbish bins, rental accommodations, an astronomy observatory, and restaurants). These two variables were measured using the same orthophotos (IRIG, 2018). We also considered the activity of rodents given by the number of camera-trap detections of either Mus musculus, Rattus rattus or Rattus norvegicus.
Data analysis
Cat population density and turnover
Efficiency of cat capture was calculated as the number of cat detections for all cameras divided by the total number of trap-nights and multiplied by 100. We calculated turnover rate for the cat population with the same formula used by Dul’a et al. (2021): (number of individuals recorded in 2021 but not in 2023/total number of individuals recorded in 2021) \(\:\times\:\) 100 (Dul’a et al., 2021).
As per Soto et al. (2023), cat density (number of cats/km2) was estimated for the 2023 survey using Spatially Explicit Capture–Recapture (SECR) models. This approach allows to incorporate a spatial aspect to traditional statistical capture-recapture models by including trap locations, and therefore, the location of each detection (Broekhuis & Gopalaswamy, 2016). Capture probability is dependent on the distance of the animal from the detector (Borchers, 2012), allowing for a non-uniform distribution of sampling effort across a study area (Efford & Fewster, 2013). Individuals are assumed to move randomly around the sampling area, whereby the sum of their activity centres constitutes the estimate of density (Green et al., 2020). SECR models were fitted using the ‘secr’ R package (Efford, 2023a). We made a data frame of capture records for cats with the trap ID, session (survey), occasion of detection and cat ID. To generate a capture histogram, we combined this capture record with an array of cameras, detailing trap ID and their XY coordinates, including usage (reflecting recording effort). The ‘count’ detector function was used, as our data consisted of a vector of counts per detector. We assumed that cats in the Ecological Park followed a homogenous Poisson spatial distribution throughout the sampling period. When fitting SECR models to this data, we selected the Halfnormal (HN) detection function for the data, as it reached a prompt plateau from a plot of effective sampling area (ESA) and ensured comparability with the 2021 survey (Figs. S2 & S3).
To represent the area in which the target species could potentially occupy, we used two different approaches and ran SECR models with two different habitat mask areas. For the best suitability to our 2023 data and study design, we applied the suggest.buffer() function to our capture histogram, specifying a radius for the habitat mask around detectors (Efford, 2023b). This buffer value was tested for suitability from the ESA plot. With this method, we selected a 1700 m-buffer around each camera, resulting in a habitat mask area of 39.24 km2. However, to be more comparable with the 2021 study, we also created SECR models with a 1200 m-buffer and a 26.01 km2 habitat mask area, as was used by Soto et al. (2023). SECR density estimates tend to be robust against a range of study designs and detection probabilities (Efford & Fewster, 2013), but the estimated population size of cats in the study area is heavily influenced by the area of the habitat mask, so we presented estimates using both masks in acknowledgment of this fact.
Using the SECR models, we also calculated Mean Maximum Distance Moved (MMDM), which was the maximum average distance moved by an individual between detections, and the spatial scale of movement, which reflected how far an animal moved within their home range during the sampling period. Finally, we estimated free-ranging cat population size by multiplying cat density by the habitat mask area.
Direct and indirect drivers of cat activity and abundance
We investigated the direct drivers of cat activity and abundance, and indirect drivers as mediated through rodent activity using Piecewise Structural Equation Models (SEMs) from the “piecewiseSEM” package from R (Lefcheck, 2023). SEMs comprise pathways representing hypothesized causal relationships between variables, so can be used to quantify indirect effects that may be overlooked by any single model. Piecewise SEMs further allow analyses to be carried out for smaller sample sizes (Lefcheck, 2016). We generated one piecewise SEM for cat activity and another for cat abundance.
Prior to running the Piecewise SEMs, we applied a z-score standardization to the variables regarding cover of open/rocky area, native and non-native shrubland and trees, and run a Principal Component Analysis (PCA) with those variables. The scores from the PC1 were retained as an explanatory variable named “vegetation simplification” (describing a gradient from tree to shrubland cover, Fig. S5) as a predictor variable, reducing habitat dimensionality (e.g., Ferreira et al., 2017). We then examined inter-variable correlation and removed variables that were highly correlated (r > 0.70) from our analysis. In this case, both altitude and slope were highly correlated with urb-dist, as the city of Funchal is located at the base of a valley, so neither altitude nor slope were included in our models.
Each SEM included two Generalized Linear Mixed Models (GLMMs): one relating rodent activity with open/rocky area, food-dist, urb-dist, water-dist and vegetation simplification, and another relating cat activity/abundance with open/rocky area, food-dist, urb-dist, water-dist, vegetation simplification and year of sampling. A random variable regarding the grid identity was retained to account for the spatially nested placement of the camera-traps within the grids. We applied Shipley’s test of directed separation to ensure no important relationships were missing from the SEM basis set (Shipley, 2000) and evaluated models by reporting standardised estimates of relationships and R2 for response variables. Indirect effects of environmental variables on cats as mediated through rodent activity were obtained by multiplying the corresponding standardised estimates (Lefcheck, 2016).
Cat activity drivers in 2021 and 2023
Given that the year of sampling was an important predictor of cat activity, we further examined the relationship between cat activity and the above-mentioned environmental variables (i.e., rodent activity, open/rocky area, food-dist, urb-dist, water-dist and vegetation simplification) separately for each of the sampling years and compared whether cat activity was still being affected by the same factors. Cat activity was modelled using GLMMs. We generated four GLMMs: namely, cat activity and abundance in 2021 and in 2023. These models consider both fixed effects (environmental covariables), and random effects (variation among grids) that may influence cat abundance/activity. We considered 11 environmental covariables (detailed above). We also included an offset accounting for the number of camera trap-nights (log10 x) which was further excluded from the analysis as it did not provide any additional explanatory power. We applied the “dredge” function to GLMMs to evaluate all combinations of fixed effects, followed by “model.avg” to help tackle model selection uncertainty. Both functions are found within the “MuMIn” R package (Barton, 2023). The residuals of the final models were evaluated with the “DHARMA” package (Hartig, 2022; Figs. S7 and S8). All data analyses carried out in R version 4.3.2 (R Studio Team, 2023).