Study area
We assessed habitat selection and activity of free-ranging caracals on the Cape Peninsula and within a small nature reserve (False Bay Nature Reserve, Fig.1a), both of which are isolated from other populations by dense residential and commercial development within the municipal limits of the City of Cape Town (Fig. 1a). The City of Cape Town has an average population density of 1,530 people/km2 (Statistics South Africa Census 2011, statssa.gov.za/). Within the Cape Peninsula, we primarily sampled individuals in Table Mountain National Park (TMNP), which encompasses approximately 320 km2 of fragmented, natural habitat, which we split into two regions (Fig. 1).
The ‘urban-dominated region’ of the study area comprises five habitat patches that encompass 187 km2 (northern half, Fig. 1), 78.4% of which is bordered by urban area. We sampled two additional individuals in an urban-dominated region within the small fragmented False Bay Nature Reserve (comprising <3 km2) surrounded by an extremely densely populated region of Cape Town (Cape Flats human population densities: 9,000–17,000 people/km2). The ‘wildland-dominated region’ (southern half, Fig. 1) encompasses 133 km2, 46.2% of which is bordered by urban area. This region consists largely of contiguous protected national park areas and small patches of vineyards (3.9%). In both regions, land uses include low to high density urban development, light industrial and commercial areas, golf courses, vineyards, eucalyptus stands, and pine plantations (urban-dominated region only). The fynbos biome dominates the natural areas of both the urban- and wildland-dominated regions and is comprised of extremely diverse, low-growing, dense shrubland with a fire-dependent ecology [49]. Across the study period, 35.2% of the available wildlife habitat in the study area intensively burned. Extensive coastal sand dune fields are also a common natural feature throughout the study area.
Capture and GPS-collaring
We sampled 29 caracals spanning the period 2014-2016. Individuals were sampled throughout the urban-dominated region, while in the wildland-dominated region, our permits constrained sampling to the southern-most location (Cape of Good Hope Nature Reserve, Fig. 1a) furthest (15.8 km) from urban development in our study area. Captured animals weighed a minimum of 6 kg so that the collar weighed no more than 3% of the animal’s body weight. Individuals were chemically immobilized, and we recorded age class, sex, weight, and morphological measurements. Individuals were classified as subadults (<2 years) or adults (>2 years) based on body size, weight, tooth wear and eruption, and reproductive status [50]. Individuals were fit with Tellus 1C collars (FollowitTM, Lindesberg, Sweden) that collected GPS locations at three-hour intervals throughout the 24 hr cycle. To assess fine-scale movement patterns and the use of potential movement corridors through urban areas, we increased the GPS-sample rate to 20-minute intervals every 10th day for 24 – 36 consecutive hours (resulting in a target of 72-96 consecutive locations). Collar data were downloaded via that FollowIt™ GEO web interface or obtained directly from collars after units dropped-off. Collars were equipped with both a timed drop-off mechanism and a rot-off cotton spacer to ensure eventual drop-off.
Landscape covariates
Caracals primarily used undeveloped habitat patches with intact shrubland vegetation (hereafter ‘natural areas’), and to a lesser extent, areas within the fully transformed urban matrix (consisting of commercial, industrial, residential, and altered open areas). Therefore, to capture the influence of the urban matrix on caracal habitat selection, we used distance measure (URBAN.DIST) that fell along a continuum of negative to positive values, where 0 indicated the line (in our GIS layer, Supplemental Table S1) demarcating the urban matrix-natural area boundary (hereafter ‘urban-natural boundary’). Negative values represented the distance inside the urban matrix measured from the urban-natural boundary whereas positive values represented the distance into natural areas measured from the urban-natural boundary. We also tested whether the availability of COVER influenced URBAN.DIST selection (COVER x URBAN.DIST).
To control for the relative influences of both natural and anthropogenic landscape features on caracal habitat selection (modelling approach described below), we tested the influence of a suite of natural features including topography (SLOPE, ELEVATION), distance from freshwater (WATER.DIST) and coastline (COAST.DIST), and whether individuals were located within sand dune fields (DUNES). We used three different measures of vegetation: i) Cover (COVER, binary), ii) normalized difference vegetation index (NDVI, continuous), and iii) habitat burn intensity (BURN, normalized burn ratio 2 [NBR2], continuous). Vegetative COVER was an indicator variable (high cover: intact fynbos, vineyard vines, pine plantations, etc., low cover: bare ground, beaches, manicured grass lawns, etc.).
Because of frequent wildfires, we attempted to sample NDVI and NBR2 covariates within a two-week time point of each GPS-point (matching the frequency of Landsat Spectral 8 updates, usgs.gov) although in some cases, daily cloud cover precluded our ability to do so. Thus, for all locations, we extracted the nearest time point for which information was available (Supplemental Methods). The median time difference between location date and vegetation index date was seven days for 3-hour fix intervals (mean = 23 days, SD = 27) and six days for 20-min fix intervals (mean = 19 days SD = 26).
Other anthropogenic landscape features included pine plantations (PINE), eucalyptus stands (EUCALYPTUS), vineyards (VINEYARDS), and distance from arterial roads (ROAD.DIST, high traffic primary and secondary roads [51]). In the wildland-dominated region, ROAD.DIST and URBAN.DIST were highly colinear (VIF > 20). We thus used only URBAN.DIST in model-fitting.
Each categorical variable was converted to indicator covariates (VINEYARD, EUCALYPTUS, PINE, DUNES, COVER). All covariates were standardized (mean = 0, standard deviation = 2, [52] with the exception of URBAN.DIST in the urban-dominated region, which was standardized by two standard deviations because the URBAN.DIST distribution for available locations (see modelling approach below) was disproportionately skewed towards negative values. Importantly, this method of standardizing URBAN.DIST did not change the selection coefficients that resulted from models; rather this standardization approach to aid in our ability to create representative plots of the relative probability of selection for URBAN.DIST (e.g., Fig 2). We calculated variance inflation factors (VIF) to ensure that no two covariates were strongly colinear (VIF ≤ 3.0 in all cases, [53]). Further details for each covariate are reported in the Supplemental Methods and Supplemental Table S1.
Movement-explicit step selection functions
We analyzed movement-explicit habitat selection using step-selection functions (SSFs [54]). To ensure we used only movement relocations, we filtered consecutive GPS-locations to have ≥50-meter step lengths (for both 3-hour and 20-min fix intervals) because, for both fix intervals, a 50-meter step length was the validated criteria for classifying GPS-clusters that formed when caracals in our study area were explicitly engaged in feeding or resting [44].
Caracals in the urban and wildland-dominated regions were subset for analyses because: i) there were stark contrasts in the degree of fragmentation by urbanization in the two-regions and sampling approach (described above), and ii) the distribution of the distance of available points from the urban matrix-natural boundary also varied substantially (Supplemental Table S2). Caracals were thus split into three subgroups: i) five adults in the wildland-dominated region, ii) seven subadult males in the urban-dominated region, and iii) 14 urban adults in the urban-dominated region. Two subadult females in the urban-dominated region were included in the adult group because, unlike subadult males, they moved and used the landscape similar to adult females in their respective areas. Specifically, their median step length and their average distance from urban-natural boundary were similar, and for the one subadult female for which we were able to assess home range, her home range extensively overlapped the adult female in the same area (Serieys, unpubl.data).
We performed SSFs on both 3-hour and 20-min interval data separately using a 1:20 match-case control empirical design. Strata were created where each individual’s end ‘used’ location (t) was paired with 20 ‘available’ locations [55,56]. ‘Available’ locations were created using random vectors originating from the location immediately preceding ‘used’ location t (i.e., location t-1). Random vectors were drawn based on the empirical distribution of turn angles and step lengths between consecutive locations derived from data on all individuals that were the same sex and age group as the focal individual but excluding the focal individual to avoid circularity [56]. For both fix intervals, three resource-independent movement parameters (STEP, LOG.STEP, and DIR.PERSIST) were calculated for all steps within each stratum to control for inherent biases in animal movement that influence habitat use [6,55,57,58]. For the fine scale 20-min interval data subset, we additionally calculated terms to control for potential tendency for movement along the same topographic gradient (SLOPE.PERSIST). See Supplemental Methods for more details about these behavioral parameters.
To correct for autocorrelation inherent in movement datasets, we calculated robust standard errors for selection coefficients using generalized estimating equations (GEE [59–61] by specifying intra-group (‘cluster’) correlation [59,61]. To specify clusters, we used a destructive sampling approach described by [61] which required that we first estimate the time lag past which temporal autocorrelation between successive caracal locations is no longer significant (α = 0.05). We fit an intercept-only mixed-effects model [57] and plotted the autocorrelation function of the model (Supplemental Fig. S1). Autocorrelation decayed after a lag of 54 hours.
We collected data during 31 collaring events from 25 unique individuals. Prima et al. [61] recommend a minimum of 20 clusters per analysis group. Thus, we destructively sampled our data (by removing 54 hours of data between successive clusters) resulting in 1-6 clusters per collaring event (31 adult and 27 subadult male clusters in the urban-dominated region, and 24 adult clusters in the wildland-dominated, range: 13–92 days/cluster). Because 20-min data collection occurred every 9th-10th day, we assigned clusters according to the individual ID and consecutive period of collection, resulting in a total of 277 clusters. To assess whether habitat selection varied with differing degrees of exposure to human activities, we divided our data into diel period based on local sunset and sunrise times. We defined night as the period when the sun is at least 9o below the horizon.
We estimated selection (ß) coefficients using the coxph function in the survival package [62] for R statistics software [63]. A primary objective was to evaluate the relationship between proximity to the urban matrix and habitat selection. Thus, we built our modelling framework around the URBAN.DIST covariate. We had two a priori expectations regarding the relationship between URBAN.DIST and caracal habitat selection: i) the relationship may be non-linear, and ii) the relationship may differ between caracals in the wildland-dominated vs. urban-dominated regions. Thus, for all subgroups, we first evaluated whether the relationship was best-captured using a simple linear, quadratic-transformed, log-transformed URBAN.DIST, or a segmented linear spline regression approach.
The segmented linear spline regression approach split the distance from the urban-natural boundary into two (or three) covariates with different slopes on either side of a breakpoint [64]. We selected the optimal breakpoint (hereafter ‘knot’) using a grid search approach by comparing quasi-likelihood under independence criterion (QIC) scores of candidate spline models (Supplemental Fig. S2) and selecting that with the lowest QIC. QIC is a criterion suggested by [59] as best for comparing model performance in GEE-based analyses.
We tested candidate knot values at one meter intervals ranging between the minimum and maximum standardized distances from the urban-natural boundary. At the scale of one meter intervals, multiple well-supported knots were identified (∆QIC < 2) but the difference in distances among the well-supported knots was trivial (< 10m). Thus, we averaged well-supported candidate knot values to obtain a single value. To validate the use of segmented regressions, we performed the segmented regression approach on each caracal subgroup and compared the segmented model QIC with models with a simple linear, log-transformed, or quadratic-transformed distance from the urban-natural boundary. The segmented approach performed best for caracal subgroups in the urban-dominated region, whereas the quadratic transformed distance performed best for caracals in the wildland-dominated region.
Next, for each caracal subgroup in the urban-dominated region, we split their data according to diel period, and performed the same grid search approach to determine whether the optimal knot shifted for each diel period. Subadults and adults thus had three knot values each (day, night, composite [day and night] dataset), ranging between 44-153 m from the urban-natural boundary. These values indicate that there exists a narrow (natural area) buffer zone around the urban matrix that caracals perceive as equal to the urban matrix itself (Fig. 2a-c). We hereafter refer to the urban matrix and the surrounding narrow, natural buffer zone as the ‘urban interface.’
For caracals in the urban-dominated region, we next considered that the relative probability of selecting certain landscape features was dependent on where individuals encountered this feature (urban interface vs. natural area). We considered this possibility for candidate covariates PINE, BURN.INDEX, COVER, ELEVATION, WATER.DIST, and COAST.DIST (hereafter ‘split’ covariates). We predicted the potential for split PINE, BURN.INDEX, and COVER because these features are indicative of available vegetative concealment that may be more strongly selected (i.e., COVER) or avoided (i.e., PINE, BURN.INDEX) within the urban interface where concealment would mitigate detection risk. Within the urban interface, we predicted that the relative importance of ELEVATION would diminish substantially. WATER.DIST within the urban interface may be linked with increased opportunity to find urban-associated prey (e.g., Egyptian geese, Alopochen aegyptiaca [44]). Finally, we considered a split potential for COAST.DIST because we predicted that caracals would avoid coastline closely bordered by residential or commercial development. We did not consider a split potential for NDVI, DUNES, VINEYARDS, EUCALYPTUS, or ROAD.DIST (hereafter ‘simple’ covariates), because with rare exception, roads and eucalyptus were primarily located within, or in close proximity to, the urban interface, while vineyards were always directly adjacent to the urban interface. We could not think of a biological reason why selection for DUNES or NDVI would differ with proximity to the urban-natural boundary.
To adapt our modelling approach to incorporate split covariates in urban models, we created an indicator variable I, where I =1 when an individual was within natural areas, and I = 0 when individuals were located within the urban interface. This allowed us to model parameters αn representing the selection estimate when individuals were within the urban interface, and αn + ßn representing the selection estimate when individuals were within natural areas. For caracals in the urban-dominated region, we thus adapted our model as follows:
w(x) = exp(ß1x1 + (α2 + ß2I)x2… ßnxn+ (αn + ßnI)xn)
In contrast, for caracals in the wildland-dominated region, we calculated the relative probability of selection w for landscape covariate x using coefficients ß via the standard exponential model [65]:
w(x) = exp(ß1x1 + … ßnxn)
Finally, for each of the three caracal subgroups, diel period-, and fix interval-specific data subset, we modeled the relative probability of a caracal selecting a particular location as a function of NDVI, BURN, ELEVATION (or SLOPE, wildland-dominated fine-scale data subset), COVER, EUCALYPTUS, VINEYARD, PINE, COAST.DIST, WATER.DIST, URBAN.DIST, ROAD.DIST and interaction URBAN.DIST x COVER). For caracals in the urban-dominated region, we tested an additional interaction URBAN.DIST x ELEVATION to assess whether caracals select for proximity to the urban-natural boundary simply because the urban-natural boundary tends to be at low elevations in this region of the study area. For each subset, we fit models consisting of every combination of these terms and compared the QIC of all resulting models. In cases where >1 candidate best-fit model was identified with a ∆QIC < 2, we selected the most parsimonious model (the model with the fewest coefficients). In cases where there were two models with the same number of terms and ∆QIC < 2, we selected that with the lowest QIC to report in figures and results, and report both competing best-fit models in the Supplemental Tables. For all top models, we report the selection coefficient ß, and 95% confidence intervals calculated using the robust standard errors (Supplemental Tables S3-S13).
Cross validation
To evaluate the robustness of all top (best-fit) models we used 5-fold cross validation following [32,66]. We built case-control SSFs using covariates retained in the top model fit to 80% of randomly selected strata. We used this model to predict w(x) scores for used and available locations from the 20% of withheld strata. For each stratum comprising used (n=1) and available (n=20) movement locations, we calculated the w(x) score which was then used to assign each location a rank of 1 to 21 (lowest to highest). We next used spearman rank correlations to compare the frequency with which the withheld used locations received each of the 21 possible rankings. We iterated this process 100 times for each top model. We report the mean and standard deviation values of Spearman correlation coefficients (rs) across all iterations (Table 1). Models with strong predictive power exhibit mean rsvalues close to 1. We also report mean rsvalues expected under random habitat selection (i.e., if models have no predictive power [66]).
Activity in high risk areas
We assessed relative activity patterns to test whether caracals in the urban-dominated region changed activity patterns mitigate risk of detection in areas with high (lethal and nonlethal) human activity [27]. We additionally assessed whether subadults in the urban-dominated region alter activity patterns to avoid interactions with adults. We thus defined risk as a function of whether: i) individuals utilized the urban-dominated region (as opposed to the low-risk wildland-dominated region), and ii) individuals were of the subadult male demographic For risk scenario (i), we compared the activity data of four adults from the wildland-dominated region, and 12 adults from the urban-dominated region (for which we had 20-min data). For risk scenario (ii), we compared the activity data of 12 adults to the data of six subadult males in the urban-dominated region. We assessed effect size using log risk ratios (“RR”,[27]) calculated from 20-min fix interval data that we subsampled to one-hour intervals to reduce autocorrelation (as in [67,68]). RR is a ratio of the mean nocturnality at sites of high human disturbance to the mean nocturnality at sites of low human disturbance. We defined each caracal observation as “active” if individuals moved ≥50 meters in the 20-minutes preceding each one-hour observation (otherwise, the time point was classified as “inactive”). This selected movement threshold is consistent with the threshold described above to isolate movement locations from the larger GPS dataset, and follows the approach previously described for classifying active periods based on GPS-movement data [67].
Our dataset was skewed towards more daytime than nighttime locations because 20-min data were collected for an interval of 36-hours beginning at 02:00. We thus used bootstrap resampling to control for individual effects and maintain even sampling among individuals. We sampled without replacement 24 daytime observations and 24 nighttime observations for each caracal for 1000 replications. For each replicate, we calculated the log proportion of active observations that occurred during the night for populations in both the urban-dominated and wildland-dominated regions. We report a weighted mean log RR with weights assigned according to the inverse variance of each replicate, the proportion of replicates for which nocturnality was greater in the population in the urban-dominated region, and the 95% bootstrap confidence interval based on the 25th and 95th greatest RRs. To interpret our results of these analyses, we back-transformed mean effect sizes and converted to unlogged risk ratios to assess percent shift towards nocturnality in areas of high risk. We assessed significance using confidence intervals– those that overlapped zero were not significant.