Habitat connectivity and resource selection in an expanding bobcat (Lynx rufus) population

Overview

We implemented population-level habitat selection using a dataset of verified bobcat sightings (camera trap, personal observation, incidental trapping) collected by the Ohio Department of Natural Resources, Division of Wildlife (ODOW) between 2010 and 2019 (out of a dataset spanning 1960–2019) and logistic regression to evaluate macro-scale land cover predictors of occurrence and predict current Ohio-wide habitat suitability. More than 95% of the observations in the original dataset were collected after 2000, and 91% were recorded after 2010. We implemented individual-level (within home range) habitat selection analyses using a GPS telemetry dataset collected between 2012 and 2014 in two subpopulations (Prange & Rose, 2020), and a weighted-distribution resource selection approach (Lele & Keim, 2006; Lele, 2009) using fine-scale predictors, such as distance to roads and streams, road density and habitat types. Lastly, we used the habitat suitability map developed based on the population-level selection analysis as the resistance (or permeability) layer to parameterize a state-wide habitat connectivity model.

Long-term bobcat sighting data

We used bobcat sightings recorded by the ODOW between the years 2010 and 2019 (Fig. 1). The dataset included multiple types of verified records, such as camera trap images, observations of live animals by wildlife experts, roadkill animals and incidentally trapped animals. The majority of the sightings were reported by the general public, who were encouraged to report all bobcat sightings through the yearly regulations booklet and other ODOW announcements via phone, email, or an online form (starting in 2017). The original dataset contained >1,500 reported sightings, and we eliminated all roadkill incidents, which would introduce a bias towards roadways, which is the largest source of bobcat mortality in Ohio (Bencin et al., 2019). After eliminating entries where geospatial data were missing, 975 unique observations remained.

Telemetry data

Telemetry study areas were located in the southern and eastern regions of the state, and were selected based on the presence of two recognized subpopulations ((Anderson, Prange & Gibbs, 2015; Prange & Rose, 2020); Fig. 1). For the individual-level resource selection analysis, we used a telemetry dataset for 28 bobcats collected between 2012 and 2014. Telemetry data collected was previously described in (Bencin et al., 2019) and (Prange & Rose, 2020). Bobcats were fitted with Telemetry Solutions Quantum 4000 (Telemetry Solutions, San Francisco, CA, USA) or Tellus GPS System (Followit, Lindesberg, SWE) collars that were programmed to locate individuals twice daily at 12-h intervals on a system that rotated through a 24-h period. GPS collars frequently located individuals only once daily, and collar performance varied. We removed several animals from analysis because of (1) a small number of fixes (<30) and (2) dispersal behavior. As such, we focused the analysis on 20 resident animals with >30 GPS fixes (maximum number of fixes was 630); 13 individuals (seven females; six males) occurred in the eastern subpopulation and seven individuals (four females; three males) in the southern subpopulation.

Habitat selection analyses

Permits/ethics statement

The program administrator for Ohio’s Division of Wildlife, Wildlife Management and Research Group approved the telemetry study along with the agency’s executive administrator for wildlife research, and the wildlife federal aid coordinator (state approval codes: WFSR12 and WFPR18). Bobcat capture and handling techniques were carried out in accordance with the American Society of Mammalogists guidelines. Personnel were trained and supported by a professional USDA APHIS (Animal and Plant Health Inspection Service) trapper, and were Safe-Capture International (Snohomish, WA, USA) certified.

Connectivity modeling

We investigated the habitat connectivity for bobcats using circuit theory-based software Circuitscape 4.0.5 (McRae et al., 2008; McRae, Shah & Mohapatra, 2013). The strength of a circuit theory approach to modeling connectivity is the ability to simultaneously evaluate contributions of multiple dispersal pathways, as opposed to simpler least-cost path metrics, which yield single connectivity pathways (McRae, 2006; McRae et al., 2008). In Circuitscape, landscapes are represented as conductive surfaces; high conductance/low resistance values are assigned to landscape features that are most permeable to movement or best promote gene flow, and high resistance/low conductance is assigned to barriers to movements.

For the purpose of this study, we decided to adopt a simple approach and used the habitat suitability map resulting from the population-level habitat selection analysis as the input resistance layer. Previous research has shown that habitat suitability can be a good surrogate for explaining functional connectivity and gene flow across landscapes (Wang et al., 2008; Chan, Brown & Yoder, 2011; Gherghel & Martin, 2020). Connectivity modeling using Circuitscape requires nodes, either points or regions between which connectivity is to be modeled. For this analysis, we selected five regions as nodes representing areas that bobcats could disperse to and from in Ohio. These areas were selected at the county level based on whether they represented (1) core areas for genetically distinct eastern and southern subpopulations (Anderson, Prange & Gibbs, 2015): one county in southern Ohio and four counties in eastern Ohio; and (2) counties with recent verified sightings in southwestern (one county), northwestern (one county) and northeastern (one county) Ohio. For this analysis, we were primarily interested in exploring two questions. First, the two genetically distinct subpopulations in the telemetry study were separated by ~200 km, and the area between the subpopulations is predominantly forested with no natural or significant man-made barriers to dispersal. As such, we wanted to evaluate the potential for admixture between the two subpopulations (a separate genetics study of 120 roadkill animals collected throughout southeastern Ohio in 2019–2020 is underway). Second, outside southeastern and southern Ohio, land cover is dominated by agriculture and urban and suburban development (Fig. 1), yet bobcats have been sighted in many areas of the state, indicating dispersal through presumably suboptimal habitat. As such, we wanted to evaluate corridors for dispersal through suboptimal habitat, which could further inform management of the population (e.g. viability of lethal management (i.e. trapping or hunting) in areas that serve as conduits for dispersal).

We modeled the relationship between habitat suitability and permeability to movement using three functions: a linear function (suitability values were inversely proportional to the resistance values (i.e. cells with highest habitat suitability values had the lowest resistance to movement, and vice versa)) and two non-linear functions. The non-linear functions took an exponential decay form using the transformations suggested by Keeley, Beier & Gagnon (2016), with an exponent c = 2 and c = 8 determining the shape of the curve (c = 8 assumes that resistance to movement decreases sharply at low suitability values, while c = 2 provides an intermediate relationship between the other two).

$$R=100-99\ast {\displaystyle \frac{1-e^{-c\ast H}}{1-e^{-c}}}$$where R = resistance, H = suitability, and factor c determines the shape of the curve (c = 2 and c = 8 in our case).

We used the resistance maps under the three scenarios to model connectivity as described above, and visualized the connectivity map using a quantile distribution of map values; this visualization approach performs best at showcasing areas of high and limited connectivity, and highlights corridors with high flow in otherwise low permeability landscapes (McRae et al., 2008). To test which relationship between resistance and habitat suitability performed best and produce a final connectivity map informed by biological data, we identified 46 bobcat locations in our dataset that were outside the high suitability regions of the state that could potentially be locations of dispersing individuals. For each of these locations (1’s), we generated three random points within a 5-km buffer (0’s), which represent available locations in the landscape. For both bobcat and random locations, we extracted the quantile (1–10) associated with the connectivity predictions under the three relationships between resistance and suitability. We used a Generalized Linear Mixed Effect Model (GLMM) framework and tested three models; each model had the quantile value for that perspective resistance-suitability relationship as a fixed effect, and location ID as a random effect (to account for variation within each location). The goal of this procedure was to evaluate which connectivity map was a better predictor for occurrences of potential bobcat dispersers. We ranked the three models using AICc, and extracted the AICc weight for each model in the model set. Lastly, we used the AICc weights to calculate a weighted average of the three maps, and produce a model-averaged connectivity map. GLMM’s were implemented via package ‘lme4’ (Bates et al., 2015), and model selection procedure via package ‘MuMIn’ (Barton, 2020) in program R 4.0 (R Core Team, 2021).

Population-level habitat selection

We found that habitat selection by bobcats in Ohio at the population level was predicted by broad scale habitat variables and distance to main roads (Table 1). The predictive ability of our models was high: average predictive ability based on a k-fold cross validation analysis = 0.94 (Supplemental Material) and AUC ROC values for candidate models used for model averaging was >0.95 (Table 1). Several competing models were within two AICc units from the top model and suggested that the proportion of forest within 50 km² had a strong, positive relationship with habitat suitability (Fig. 2). The proportion of pasture within 50 km² and the distance to high traffic roads also had a positive association with habitat suitability, but not as strong as forest (Fig. 2). The proportion of natural herbaceous vegetation at finer scales (15–50 km²) was included in different models within the two AICc model set, in this order, and had a strong, positive association with habitat suitability (Fig. 2; only relationship to natural herbaceous at 30 km² shown). The highest habitat suitability for bobcats was in eastern and southeastern Ohio, which coincides with the areas of forest cover, and bobcat population expansion (Fig. 3). There was high suitability along the southern border with Kentucky and the southwestern Ohio border with Indiana. Central and northwestern Ohio are dominated by cropland, and had low suitability, except for patches of habitat along river corridors; these patches may act as stepping stones for dispersing individuals but are unlikely to be occupied permanently. Northeastern Ohio, which is a mosaic of forest, agriculture and human development, had intermediate suitability (Fig. 3).

Individual-level habitat selection

The number of bobcat GPS telemetry locations ranged between 31 and 630 locations per individual (Supplemental Material). Females had smaller mean 95% UD home ranges (27.27 km², range = 7.70–77.82 km²) compared to males (69.79 km², range = 27.85–165.66 km²) (Kruskal–Wallis chi-squared = 6.477, df = 1, p-value = 0.011). Home ranges were larger in the southern study area (mean = 79.23 km², range = 23.35–165.55 km²) compared to the eastern study area (mean = 28.73 km², range = 7.70–8.88 km²) (Kruskal–Wallis chi-squared = 5.841, df = 1, p-value = 0.015).

The RSF fitted using a weighted-distribution approach with an equal number of locations and random availability points with individual home ranges showed that at the home range scale, bobcats exhibited the strongest selection for areas with low road density (Fig. 4A) and those farther away from main, high-traffic roads (Fig. 4C). Although bobcats also preferred areas with higher canopy cover and a low cover of impervious surfaces, these variables had lower selection strength (for or against; Table 2). Distance to streams and distance to all roads were also weak predictors of selection (Figs. 4D, 4F), relative to road density and high-traffic roads. Land cover played an important role in habitat selection within home ranges, with the strongest selection for forest habitat (Fig. 4B). Bobcats showed weaker selection for open habitat (pasture, grassland, Fig. 4B); although the percent of impervious habitat in a given grid cell was negatively associated with bobcat use, bobcats exhibited selection towards low-intensity developed habitat (Fig. 4B).

Habitat connectivity

We found that the three relationships between habitat suitability and connectivity performed equally well (similar AICc weights) when predicting high connectivity area for dispersing individuals (Table 3). As such, the model-averaged maps were equally weighted in the final connectivity map (Fig. 5). When examining the statewide habitat connectivity for bobcats, several interesting patterns emerged. First, habitat was highly permeable to movements between the core areas of two genetically distinct subpopulations—eastern (region 1) and southern (region 2) (Fig. 5). Connectivity between southern and southwestern Ohio (region 3) was lower and concentrated along the Ohio River Valley; this is likely an outcome of the high proportion of agricultural lands and urban development outside of the river valley (Cincinnati, Dayton and suburbs). Recent sightings in southwestern Ohio could be either animals from neighboring Indiana or Kentucky or dispersers from southeastern Ohio (e.g. one GPS-collared subadult male captured in region 2 traveled to the Cincinnati area during the telemetry study). A pattern of high connectivity was evident between eastern and northeastern portions of the state (region 5; Fig. 5). This region is characterized by a mosaic of forest (including the Cuyahoga Valley National Park), agricultural lands, and low-intensity urban development, which allows for many conduits of dispersal. One of the most interesting outcomes of the connectivity analysis was the potential dispersal corridors in central Ohio, to and from the northwestern region (Fig. 5). Within the intensive row-crop agriculture landscape, which is not conducive to bobcat dispersal (Reding et al., 2013), riparian forest cover along water courses could act as conduits for dispersal in this landscape, although overall connectivity was lower compared to the eastern part of Ohio. These areas coincide with several recent verified sightings in central and western Ohio (Fig. 1), and are potential movement corridors across the landscape. Western and northwestern Ohio had relatively low connectivity, spread across a larger area; this area is relatively homogeneous, with agricultural lands interspersed with small forest patches. Lastly, there were several areas with no connectivity in central Ohio; these areas lack forest cover and have intensive agricultural use (Figs. 1, 5).

Population-level selection

Unsurprisingly, the highest probability of bobcat occurrence (i.e. suitability) was in the forested southeastern part of Ohio. Eastern and southern Ohio regions are home to two genetically distinct (as of 2012) sub-populations (Anderson, Prange & Gibbs, 2015), and the current population is thought to have expanded from these two regions (roughly, areas 1 and 2 in Fig. 5); preliminary data suggest that as of 2020, the two subpopulations are no longer distinct (V Popescu, 2021, unpublished data). The earliest modern-day records of bobcats in Ohio occurred in these regions, which are in close proximity to well-established populations in Kentucky and West Virginia, and they have continued to have the most consistent sightings over the past several decades (Fig. 1). Ohio bobcats selected for areas with highest forest cover (Fig. 2A), which corroborates other bobcat populations in Midwestern US states with similar history of recovery (e.g. Illinois and Iowa; (Nielsen & Woolf, 2002b; Woolf et al., 2002; Reding et al., 2013)). These findings are very different than heavily forested portions of their range (e.g. New England), where bobcat habitat suitability was largely driven by snowpack depth (Reed et al., 2017b). However, we also found that bobcats selected for higher proportions of natural herbaceous vegetation, as well as pasture (Figs. 2B, 2C). As a species with high flexibility in habitat selection at the population scale (bobcats have a large distribution in North America overlapping a broad range of habitat types; (Roberts & Crimmins, 2010; Dyck, Wyza & Popescu, 2021)), bobcats can take advantage of existing habitat resources and be successful in a variety of natural, rural and suburban settings (Kamler & Gipson, 2000; Linde et al., 2012; Beattie, 2020).

Our models predicted that northern and northeastern Ohio had low-medium suitability. While our dataset did not include many bobcat occurrences in these regions, bobcats are known to thrive in heterogeneous, human-dominated landscapes, and recent observations of animals in these areas are increasing. The increasing number of observations in peripheral areas is concomitant with an increase in observations in SE Ohio, and could also be an artifact of greater availability and affordability of monitoring technology (i.e. trail cameras), which also likely led to increased reporting of animals. Bobcats in other US Midwest regions (Tucker, Clark & Gosselink, 2008; Reding et al., 2013) and the Eastern US (Beattie, 2020) were found to be resident of heterogeneous mixed-use and suburban landscapes, which may include lower suitability areas. Given that the population has likely not yet reached equilibrium and continues to expand, it is likely that bobcats can take advantage of resources in this region composed of a mosaic of forest, agriculture and suburban development. Moreover, bobcat densities in heterogeneous landscapes (including suburban and forest interspersed with open natural habitat) may reach higher levels compared to forested landscapes; early seral habitats likely have higher prey availability compared to mature forest, which results in smaller home ranges (Tucker, Clark & Gosselink, 2008; Prange & Rose, 2020) and potentially higher density. Because bobcat occupancy and density are linearly and positively related (Clare, Anderson & MacFarland, 2015), it is important to continue to collect bobcat occurrence data and update the population-level selection models periodically to reflect future expansion in areas currently predicted as low-medium suitability.

Habitat suitability in the central, western and northwestern part of Ohio was low due to the landscape characterized by intensive agriculture with little natural habitat remaining. Bobcat observations in these areas were rare, only occurred in the last 3–4 years of the dataset and likely reflect dispersal events. These observations are important from two perspectives; first, they showcase the ongoing expansion process, which makes the case for monitoring these areas closely to evaluate residency (e.g. recording females with kittens) and for updating the population-level selection in the future; second, there is a need to understand the selection of dispersal habitat and if, or how, it differs from selection of residential habitat (Reding et al., 2013). The Ohio bobcat population is part of the broader US Midwestern and Eastern population; while genetic differences exist between populations at a regional level (Anderson, Prange & Gibbs, 2015), long-distance dispersal movements do occur, even in unsuitable landscapes that pose high resistance to movement. Bobcats may make use of forested riparian corridors or forest patches embedded in agricultural landscapes, thus maintaining the integrity of such features could be critical for the regional population.

Individual-level selection

Habitat selection at the individual animal level based on GPS telemetry of 20 adults (11 females and nine males) in two areas of Ohio that coincided with genetically distinct populations revealed complex responses to various cover types, topography, and human disturbance. An earlier study of habitat selection using a slightly larger version of the same dataset used here and Manly’s selection ratios (Prange & Rose, 2020) showed that adult bobcats selected for forest habitat in the southern population and heterogeneous landscapes resulted from reclaimed coal strip-mining in the eastern population (Fig. 1). Reclaimed strip mine lands provide a mosaic of forest and natural herbaceous vegetation that may provide more foraging opportunities and prey items compared to closed-canopy forest habitat (Rose & Prange, 2015; Prange & Rose, 2020). Our findings corroborate current studies, and also adds other evidence for home range scale selection, particularly related to artificial features. Bobcats in our study area selected for habitats with less road density, and farther away from high traffic roads (interstate, state and US routes), and showed a positive association with lower-traffic roads (including unpaved roads). These patterns of habitat use in relation to roads have been found in other telemetry studies across the bobcat range in North America (New Hampshire (Broman et al., 2014; Reed et al., 2017a); California (Poessel et al., 2014)). The information provided by the individual-level selection complements the population-level selection findings and has value when considering habitat connectivity, as discussed below. High traffic roads were a negative predictor at both population and individual-level selection analyses. While bobcats persist in human-dominated landscapes, including suburban, these findings, along with avoidance of high road density areas, suggest that areas away from high-traffic roads and with low road density are important for population persistence. High traffic roads not only affect habitat suitability statewide, but also bobcat activity at home range level, and have negative effects on population demography via road mortality (Litvaitis et al., 2015; Bencin et al., 2019).

Habitat connectivity

Our assessments of habitat connectivity corroborate those of Reding et al. (2013), who studied bobcat habitat resistance to movement using genetic data. They found that in forested areas, the landscape was permeable to movement, movements were determined by landscape composition (i.e. availability of forest), and that Euclidean distance did not explain genetic connectivity. For example, we found high connectivity between the eastern population core (region 1 in Fig. 5) and northeastern Ohio (region 5), suggesting that the paucity of observations in this area is not due to a lack of dispersal potential, but likely stems from lower overall suitability and the fact that the population is yet to reach equilibrium and achieve its full distribution in Ohio. We also arrived to similar conclusions concerning low connectivity in the intensive, row-crop agriculture and urban landscape of central and western Ohio. These areas have both low habitat suitability, and low permeability to bobcat movements (relative to the forested areas; Hughes et al., 2019), as well as urban centers, which have been shown to act as barriers to movement in some cases (Reed et al., 2017a). For example, Hughes et al. (2019) found that intensive row-crop agriculture regions of Iowa pose a barrier to movement of dispersing animals, thus limiting the expansion of the Iowa bobcat population.

Although Ohio is one of the most road-dense states in the United States, and roads can be a significant source of mortality for bobcats (Litvaitis et al., 2015; Bencin et al., 2019), we did not include roads as potential barriers in the connectivity models. Bobcat dispersal movements are typically performed by 1- and 2-year old individuals (Nielsen & Woolf, 2003; Johnson, Walker & Hudson, 2010), and are biased towards males (in Iowa, 65% of males and 26% of females dispersed; Hughes et al., 2019). These studies found that bobcats can travel hundreds of kilometers through a variety of habitats, and roads, including high traffic interstate routes, were not a barrier or pose resistance to bobcat dispersal movements (Nielsen & Woolf, 2003). However, dispersing bobcats, particularly males, are susceptible to higher road mortality; an ongoing analysis of 300 roadkill individuals in Ohio between 2010 and 2020 showed that 60% of animals were males and roadkill proportion was twice the stable stage distribution for 1-year old animals (M Dyck, 2021, unpublished data).

The habitat resistance-suitability relationships are still the subject of ongoing research. As such, imposing resistance-to-movement values in environmental layers is highly controversial, and different approaches may yield vastly different results (see Reed et al. (2017b) for an example on habitat connectivity modeling for bobcats in New Hampshire using different data sources). Optimization techniques based on genetic information can be used to determine the best environmental layers for modeling connectivity (Peterman, 2018). However, habitat models have been shown to perform well in explaining dispersal and gene flow across landscapes (Wang et al., 2008), and in the absence of information on habitat use by dispersers.