Multiscale habitat relationships of a habitat specialist over time: The case of ocelots in Texas from 1982 to 2017

Evaluating temporal trends in habitat and behavioral responses is critical for conservation, yet long-term monitoring studies are rare. We used a 35-year dataset (1982 – 2017) to assess multiscale habitat use and selection by an endangered carnivore, the ocelot ( Leopardus pardalis ), in South Texas, USA. We used a time series of remotely sensed imagery to map changes in availability of woody cover, habitat critical to ocelots that has diminished due to anthropogenic development and increased road infrastructure. Our objectives were to characterize habitat relationships, predict high-quality habitat, and assess behavior with changing environmental conditions. We fit functional response (third order) and individual-specific resource selection (second order) functions to assess multiscale habitat use of vegetation cover and roads. Within home ranges, ocelots used woody cover greater than availability. Ocelots used areas near roads in proportion to availability, with minor exceptions. We observed changes in habitat use by ocelots across time with higher proportions of woody and non-woody cover used in later time periods. Average availability of woody cover decreased in the study area between the 1980s and 2010s (0.44 in 1985 to 0.39 in


INTRODUCTION
Human-mediated habitat change has been identified as the main threat to biodiversity worldwide (Tilman et al., 2017;Vitousek et al., 1997). Understanding relationships between wildlife and their habitat is critical for conservation of biodiversity (Cramer & Bissonnette, 2005) and improved management of natural resources (Morrison, 2001). Managers need to have a clear definition of habitat for a species of interest (Morrison, 2001) and understand animal-habitat relationships, particularly for species of conservation concern (Elith & Leathwick, 2009;Holbrook et al., 2017). Identification of habitat use patterns is critical in efforts to conserve endangered species as it provides fundamental information such as attributes important for survival and reproduction (Manly et al., 2002). Changes in movement, survival, and reproductive success of individuals may result from increased human activity on the landscape (Chen & Koprowski, 2016). As landscapes are altered through time, it is crucial to understand underlying demographic trends in animal populations and resulting changes in behavioral patterns (Anderson et al., 2012;Shoemaker et al., 2018).
Evaluation of temporal changes in habitat, and subsequent behavioral responses by wildlife is particularly important for federally threatened or endangered species. The US Endangered Species Act requires the designation of "critical habitat" and suggests monitoring "critical habitat" over time is essential. Habitat and species conservation strategies hinge on knowledge of space use by animals (Aarts et al., 2008). Highly mobile species further complicate habitat monitoring due to large-scale occupancy and use patterns, yet emerging threats to habitat can often only be understood from broad-scale monitoring across space and time (Simons-Legaard et al., 2016).
Long-term, multidecadal studies and monitoring of habitat can provide crucial information for conservation planning. Remote sensing techniques have been widely used in vegetation mapping and assessments of spatial distributions of wildlife (Mata et al., 2018;Xu et al., 2014). Remotely sensed repositories are increasingly being used to map land cover, monitor temporal trends, and compare land cover between landscapes (Hansen et al., 2013;Savage et al., 2018;Vogeler et al., 2018). Satellite imagery has been widely used in building habitat models for many animal species, including assessments of functional landscape connectivity and designation of critical habitat for recovery (Morzillo et al., 2011;Roever et al., 2013). Despite these varied applications, long-term studies are often lacking in animal ecology and conservation (e.g., Bartel & Sexton, 2009;Simons-Legaard et al., 2016).
Anthropogenic land development is a major cause of species endangerment, with linear infrastructure, such as roads and utility lines, representing a widespread and significant contributor (Chen & Koprowski, 2016;Forman & Alexander, 1998). Linear infrastructure, particularly roads, can affect species in many ways (D'Amico et al., 2016;Fahrig & Rytwinski, 2009;Pattison et al., 2020) including direct impacts from roads through vehicle collisions and as barriers to movement (Forman & Alexander, 1998). Indirectly, roads can influence the loss of functional connectivity across the landscape (Latham et al., 2011), affecting gene flow within and among populations (Baigas et al., 2017;Forman & Alexander, 1998;Jahrsdoerfer & Leslie Jr., 1988). Moreover, roads may alter animal behavior through increased avoidance and reduced habitat use (Dwinnell et al., 2019), ultimately leading to a reduction in potential habitat that exceeds the footprint of the road itself. By understanding behavioral responses of animals to roads and traffic, researchers can gain insights into the effects of road networks on wildlife and contribute to effective mitigation efforts (Chen & Koprowski, 2016;Haddad, 1999).
As landscapes are altered through time from increasing disturbance, it is crucial to understand resulting changes in habitat relationships at the individual and population level (Anderson et al., 2012;Shoemaker et al., 2018;Squires et al., 2020). The responses of populations toward anthropogenic influences are often inconsistent, usually associated with local environmental conditions (Mumma et al., 2019). These context-dependent responses can be broadly described as functional responses in habitat use or selection (Holbrook et al., 2019;Mysterud & Ims, 1998). Functional responses in habitat use or selection are a means to assess how animals alter habitat relationships across spatial and temporal variation in resource availability. For example, boreal woodland caribou (Rangifer tarandus caribou) responses to anthropogenic features varying across seasons (Mumma et al., 2019). Indeed, functional responses can help advance the understanding of animal-habitat relationships as well as inform conservation efforts.
In this study, we assessed functional responses and resource selection to assist conservation planning and understand carnivore habitat relationships across broad time scales. We applied our approach to the ocelot (Leopardus pardalis), a carnivore that is federally endangered and considered a habitat specialist in the United States (Tewes, 1986). More specifically, ocelots use dense woody vegetation, specifically native thornshrub, with 95% vertical and >75% horizontal canopy cover, and demonstrates avoidance for open land cover types at coarse spatial extents and higher orders of selection (i.e., second and first order; Johnson, 1980;Harveson et al., 2004;Jackson et al., 2005;Horne et al., 2009;Veals, 2021). In the United States, ocelots are restricted to extreme South Texas in two remnant populations along the United States-Mexico border and occur on mostly private lands (Haines et al., 2006;Janečka et al., 2011Janečka et al., , 2016. Both populations occur within or adjacent to the Lower Rio Grande Valley (LRGV), a rapidly expanding urban area resulting in increased infrastructure (Tiefenbacher, 2001), including the expansion and development of transportation networks . With the growth of human population centers in the LRGV, ocelots are facing growing pressure to survive in an increasingly modified landscape. This pressure is exacerbated by the development of road networks, with vehicle collisions estimated to account for 35%-40% of ocelot mortality-the highest source of mortality for ocelots in this region (Blackburn et al., 2021;Haines et al., 2005). Expansion of road networks will likely lead to a continued transportation-related ocelot mortality and decrease in accessible habitat in Texas (Blackburn et al., 2021;Haines et al., 2005).
The goal of our study was to quantify multiscale spatiotemporal trends in ocelot habitat use and predicted habitat from the 1980s to present as well as assess how habitat use changed with respect to availability through time. By employing a time series of satellite imagery and ocelot telemetry data, we evaluated long-term changes in ocelot habitat across their current distribution in the USA. Our objectives were to (1) evaluate and compare functional responses of ocelots relative to surrounding environmental variables across multiple decades, (2) quantify habitat availability over 35 years and evaluate resource selection by ocelots across South Texas, (3) assess spatiotemporal consistency and variation of ocelot habitat across time, (4) characterize the distribution of consistently high-quality ocelot habitat across landownership to inform future conservation efforts, and (5) provide recommendations for ocelot habitat conservation and road mitigation. Land use changes are impacting wildlife globally (Tilman et al., 2017;Vitousek et al., 1997) and monitoring long-term trends in habitat and animal space use will serve an important role for developing conservation strategies for threatened and endangered species (Cramer & Bissonnette, 2005;Elith & Leathwick, 2009;Morrison, 2001).
We addressed these objectives by implementing a multiscale and multidecadal approach to assess effects of cover and road networks on ocelot habitat use and provided a time series of predicted ocelot habitat for conservation and monitoring. We quantified habitat use of ocelots at the third order (within home range; Johnson, 1980) and second order (home range placement; Johnson, 1980) with respect to roads and habitat in a highly fragmented landscape over several decades. Within home ranges, we hypothesized that ocelots would show stronger disproportionate use for woody cover at lower availabilities and, as the resource became more abundant, use would decline and become proportional to availability (Holbrook et al., 2019). This pattern of use is consistent with a habitat specialist. In contrast, we hypothesized ocelots would decrease use of areas closer to roads, such that as roads became more abundant, ocelots would increasingly avoid areas closer to roads within home ranges. Regarding home range placement, we expected ocelots to strongly and consistently select for areas with higher woody cover, and that woody cover would largely drive the predicted distribution of ocelots through time (Harveson et al., 2004;Horne et al., 2009;. Vehicle-induced mortality of ocelots continues to occur (Blackburn et al., 2020(Blackburn et al., , 2021Haines et al., 2005). Therefore, we expected the impact of roads to increase through time as indicated by ocelots avoiding areas closer to high-traffic roads in more recent years. Finally, given the high proportion of private landownership across Texas, we hypothesized most high-quality habitat to occur within private lands .

Study area
Our study focused on the three southernmost coastal counties of South Texas: Kenedy, Willacy, and Cameron. This area includes differing land use practices and vegetation characteristics. We buffered known ocelot locations from the dataset by 8 km to define the available study area for the analyses to include all areas with documented ocelot populations and the landscape between ( Figure 1; Janečka et al., 2011; Tewes, 2019).
One ocelot population resides in and around Laguna Atascosa National Wildlife Refuge (hereafter refuge population), and the other occupies private ranchlands to the north (hereafter ranch population; Tewes & Everett, 1986;Janečka et al., 2016;. The refuge population resides in the southernmost county, Cameron County, within the LRGV (Janečka et al., 2016;Tewes & Everett, 1986). This area consists of salt flats, marshes, chaparral, and thornshrub-grasslands (Lonard & Judd, 1985). The ranch population occupies private ranchlands to the north in Kenedy and Willacy counties (Janečka et al., 2016;Tewes & Everett, 1986). The ranchlands have a similar vegetation composition to the refuge as well as native woodlands dominated by mesquite (Prosopis spp.) and live oak (Quercus virginiana) with varying levels of understory cover (Leonard, 2016). The refuge is enclosed by developed areas containing a network of highways with an average of 600-11,000 vehicles/day (Texas Department of Transportation [TXDOT], 2019). In contrast, the ranches in Kenedy and Willacy counties are contiguous rangelands with few roads bordered by two highways with comparable traffic volumes to the refuge (TXDOT, 2019). These two populations are separated by $30 km characterized by agriculture fields, cattle pastures, and high-traffic roadways.

Ocelot spatial data
We used telemetry data from radio-collared adult ocelots of both populations (n = 78 ocelots) collected between 1982 and 2017 from previous studies (Harveson et al., 2004;Horne et al., 2009;Laack, 1991;Leonard, 2016;Leonard et al., 2020;Lombardi, 2020;Shindle & Tewes, 2000;Tewes, 1986). Data from the refuge population were collected with very high-frequency (VHF) collars from 1982 to 2001 (n = 50; 28 males and 22 females). Data for the ranch population were collected using VHF and GPS collars from 1982 to 2017 (n = 28; 12 males and 16 females). Locations were collected at least once a day using ground-based triangulation with multiple observers for VHF collars and GPS collars had a fix rate between four and eight locations per day. We sampled used locations from individual ocelots (n = 52 with VHF and n = 26 with GPS data). Some individuals were collared multiple times with at least 2 years between sampling events and ≤15% overlap in home range area (n = 7). We treated each sampling event as distinct observations given the change in availability across time and space. This resulted in 85 sampling events, which we used to characterize ocelot resource use. We built 95% minimum convex polygons to approximate use areas while extensively sampling availability across orders of selection (Hebblewhite & Merrill, 2008;Holbrook et al., 2017). The dataset spanned 35 years ; therefore, the areas we classified as available to ocelots changed based on the period of their spatial data. We subset individuals by sex then grouped them into one of four periods based on when spatial data were collected: 1985 groups (14 males and 15 females), 1995 groups (16 males and 9 females), 2005 groups (4 males, 6 females), and 2015 groups (10 males and 11 females). Thus, we had eight groups of ocelots (4 time periods Â 2 sexes).

Resource variables
We determined land cover types using unsupervised classification methods (Perotto-Baldivieso et al., 2009;Xie et al., 2008) in ERDAS Imagine 2018 (Hexagon Geospatial, Norcross, GA, USA) from Landsat imagery with 30-m spatial resolution for the study area. Each output image was classified into four general land cover classes: woody vegetation, non-woody vegetation, other (i.e., sand dunes, mudflats, agriculture, developed), and water (adapted from Mata et al., 2018;. To ensure land cover data matched the temporal scales of the ocelot telemetry data, we classified four images, one every 10 years (1985, 1995, 2005, and 2015) to match the four periods. We downloaded Landsat images through the US Geological Survey Global Visualization Viewer (2018). Images before 2015 were from Landsat Satellite 5 Thematic Mapper and the 2015 image was from Landsat Satellite 8 Operational Land Imager and Thermal Infrared Sensor. We clipped all images to the same area of interest and spatial extent to ensure all rasters aligned correctly. After classifying imagery, we performed desk accuracy assessments to refine classifications with ≥85% image accuracy (Jensen, 2016). We generated 200 random points per image using ArcMap version 10.6.1 (ESRI, Redlands, CA, USA), compared visual observation and image classification for each point, and calculated an overall accuracy assessment using a confusion matrix (Congalton, 1991;Foody, 2009). If overall accuracies were <85%, classified images were reviewed and reclassified until the 85% threshold was achieved (Mata et al., 2018). Finally, experts with experience in the study area (n = 3) reviewed our classifications to validate accuracy assessments.
We created focal raster layers for woody and non-woody vegetation cover for each of the four periods. We conducted moving window analyses on each raster layer to characterize proportion land cover (Hagen-Zanker, 2016). We determined mean percent land cover for each pixel with a sex-specific focal radius (150 and 250 m for females and males, respectively). We calculated sex-specific focal radii based on the average daily step lengths sampled across the GPS telemetry dataset. We generated eight raster layers for proportion of woody cover and proportion of non-woody cover, which characterized the vegetation composition around a focal cell that male and female ocelots could encounter on the landscape. We also characterized the spatial variation in woody cover by calculating the SD of proportion of cover for each cell using the same sex-specific focal radii.
We classified roads within our study area based on annual average daily traffic (AADT) records from the TXDOT. We classified roads into soil, gravel, or paved. We further classified paved roads by traffic volume since 1999, which was the earliest year reported for the study site (TXDOT, 2019). Low-traffic roads were classified as those with <1000 AADT (<1 car/min), medium as 1000-5000 AADT, and high as >5000 AADT. Traffic volume data were not collected prior to 1999, therefore through much of our study (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998) we did not have traffic volume estimates. However, when reviewing trends in traffic volume from 1999 to 2019 and extrapolated backwards to pre-1999, roads classified as high volume still would have had substantially greater traffic volume than roads classified as low or medium volume, relatively. Additionally, when analyzing satellite imagery from the period of this study, there did not appear to be a considerable amount of change in the road structure or number of roads on the landscape that ocelots would have encountered (adapted from Blackburn et al., 2021). For each of the five road types, we calculated the log-distance from each raster cell in the study area to each road type using ArcMap to account for potential non-linear responses. Preliminary analyses indicated that soil and gravel roads were not evenly distributed across the landscape; thus, we removed them from analyses to maintain reasonable assumptions of availability across ocelots.

Functional response analyses-Third order of selection
To address how ocelots modified habitat use across spatial and temporal variation in availability, we characterized habitat use within home ranges (third order of selection; Johnson, 1980) by implementing Approach 1 as described by Holbrook et al. (2019) for additive habitat use. We tested the null hypothesis that the additive difference between use and availability is constant. We assessed how the sample of use changes with the sample of availability on the additive scale. For each used location, we randomly generated three available locations within the corresponding home range. This ensured appropriate coverage of availability for all ocelots while sampling such that the ratio of use: availability for each ocelot far exceeded the ratio necessary to generate stable coefficients for a patchy, heterogeneous landscape (Northrup et al., 2013). We extracted the proportion of woody cover, proportion of non-woody cover, and distance to low-, medium-, and high-traffic roads for each used and available location. We first characterized habitat use and availability for each individual ocelot by calculating a mean value for all covariates (Holbrook et al., 2017(Holbrook et al., , 2019. We then built linear regression models to assess functional responses in habitat use for male and female ocelots by time period. We evaluated the slope and intercept of our models to indicate the occurrence of a functional response in habitat use in relation to the surrounding availabilities of a resource. A linear model with a slope coefficient (AE90% CI) that did not overlap 1 (i.e., slope is significantly different from 1, p < 0.1) represented a functional response in habitat use (Holbrook et al., 2019). Statistical deviations from proportional habitat use (slope = 1, and y-intercept = 0) could indicate a diversity of responses (positive or negative) that deviate from proportional habitat use (Holbrook et al., 2019;Squires et al., 2020). We then plotted response curves while also displaying the interquartile range (IQR) of availability across all groups (sexes and time periods) for each variable; the IQR aided in contextualizing response curves and availability for each variable.
We extracted and then plotted the slopes and intercepts from the model for each of the eight groups. If the slope estimate was statistically significant, the 90% CI would not overlap 1, and if the y-intercept estimate was statistically different, the 90% CI would not overlap 0. We then examined our model estimates and the 90% CIs to determine if there were temporal changes in spatial functional responses for male and female ocelots. If the 90% CIs did not overlap between groups of ocelots for either the slope or the intercept across time, this result was taken as evidence of temporal differences in habitat use. We conducted analyses in program R (v3.6.2; R Core Development Team, 2019).

Habitat selection and mapping-Second order of selection
To adequately sample resource availability across the study area (second order), we increasingly sampled the most variable resource layer (SD of woody cover) to determine the sampling intensity that approximated the true variation across our study area. Per time period, 5000 locations were randomly distributed across the study area with a minimum density of one point/30 m 2 . This ensured appropriate coverage of availability for all ocelots while sampling such that the ratio of use:availability for each ocelot far exceeded the ratio necessary to generate stable coefficients for a patchy, heterogeneous landscape (Northrup et al., 2013). Based on our 35 years of experience sampling ocelots, we felt that our study area adequately captured what was available to the two ocelot populations for home range establishment.
We assigned resource variables (Table 1) to used and available locations in ArcMap. We calculated means and 95% CIs at used and available locations for every ocelot using the package Rmisc (Hope, 2013) in program R (version 3.6.1; R Core Development Team, 2019). We compared use versus availability of these variables across the four time periods for each sex to understand patterns of habitat use as these variables likely influence ocelot behavior (Haines et al., 2006;Harveson et al., 2004;Horne et al., 2009).

Resource selection functions
We built resource selection functions (RSFs; Manly et al., 2002) at the second order (home range placement; Johnson, 1980) with logistic regression (e.g., logit link, with a dependent variable of use vs. available) using program R to understand multivariate resource selection by ocelots. We treated resource selection as individual-specific, such that we estimated an RSF for each ocelot. Because of the expansive time series and two populations, we evaluated each ocelot's unique pattern of selection.
For each individual dataset, we fit a global model with scaled covariates, using the scale function in R, including proportion of woody cover, proportion of non-woody cover, heterogeneity of woody cover, and log-distance to low-, medium-, and high-traffic roads. We only included variables in the global model if they were not correlated (jrj < 0.7; Dormann et al., 2012). We then estimated sex-level RSF coefficients for each period using a two-staged approach (Fieberg et al., 2010;Marzluff et al., 2004), where we weighted each coefficient based on the inverse of the variance (DeCesare et al., 2012;Murtaugh, 2007). We then averaged RSF coefficients for each sex within each period.
Habitat mapping, validation, and temporal monitoring We used the group-level RSF coefficients to map the relative probability of use (w[x]) per pixel for ocelots across the study area at 30-m resolution: w x ð Þ ¼ exp fβ 1 woody cover ð Þ þ β 2 non-woody cover ð Þ þ β 3 heterogeneity of woody cover ð Þ þ β 4 log of distance to low-traffic road ð Þ þ β 5 log of distance to medium-traffic road ð Þ þ β 6 log of distance to high-traffic road ð Þ We reclassified each map into 10% equal-area bins to balance the heavily skewed distributions in our predictions (Boyce et al., 2002). To assess if our predicted probability of use maps captured ocelot resource use, we examined the percent of ocelot telemetry locations across the 10 quantile bins for each sex and time period. We calculated Spearman rank correlation coefficients to assess the association between the quantile bins and the percent of ocelot locations. If the RSF models predicted areas of high ocelot use, we expected a strong Spearman rank correlation coefficient (ρ; Boyce et al., 2002). We then assessed spatiotemporal consistency in our predictions of ocelot habitat across our study area. We applied the approach of Holbrook et al. (2017) to identify the cut-point that captured 90% of the cumulative ocelot use (from high to low) for each sex and time period. This allowed us to develop binary characterizations of ocelot use as well as assess spatiotemporal patterns in how-high use areas changed across years within each sex. For each year and sex, we mapped the predicted area that captured 90% of ocelot telemetry locations. We then overlaid these maps and assessed if a given pixel was mapped once, twice, three times, or four times. If a pixel was mapped four times, it indicated that it was classified as high use for all four time periods (i.e., for time period 1985, 1995, 2005, and 2015). The pixels that were never mapped as high-use areas were characterized as non-habitat. We then used our spatial representations of high-use areas (i.e., 1 decade, 2 decades, 3 decades, and 4 decades) to examine how high-use areas for ocelots compared to contemporary landownership (circa 2015).
Finally, we assessed variation in the predicted probability of use across time periods using a pixel-based correlation analysis. We calculated a Spearman rank correlation coefficient (ρ) per pixel across time periods using our equal-area binned maps characterizing the relative probability of use by male and female ocelots. This produced an independent raster with values ranging from À1 to 1. A value of ρ = À1 represented strong negative correlation and therefore a decreasing trend in predicted ocelot use, while a value of ρ = 1 indicated a strong positive correlation and an increasing trend in predicted ocelot use across time. When ρ = 0, there was no correlation, which indicated variability with respect to predicted ocelot use across time periods.
T A B L E 1 Description of resource variables from raster layers used in habitat use analyses, South Texas, USA, 1982-2017.

Proportion woody cover
Raster created from moving window analysis based on sex-specific focal radii a to determine proportion of woody cover around a focal cell. Characterized vegetation composition around a focal cell.
Proportion non-woody cover Raster created from moving window analysis based on sex-specific focal radii to determine proportion of non-woody cover around a focal cell. Characterized vegetation composition around a focal cell.
Heterogeneity of woody cover SD of proportion woody cover for focal cell based on sex-specific radii. Measure of heterogeneity of woody cover around a focal cell.

Log-distance to low traffic
Log-distance from center of focal cell to low-traffic road type (<1000 AADT b ).
Log-distance to medium traffic Log-distance from center of focal cell to medium-traffic road type (1000-5000 AADT).

Log-distance to high traffic
Log-distance from center of focal cell to low-traffic road type (>5000 AADT).

Functional response analyses-Third order of selection
We found evidence that ocelots demonstrated functional responses to vegetation cover and road-related attributes at the third order. Several groups of ocelots had statistically significant slope estimates for cover and road variables (Table 2). Females exhibited a pattern of proportional use of woody cover, rather than a functional response, for 1995, 2005, and 2015. During the 1985 time period, however, females displayed a pattern consistent with relaxed selection (slope < 1); as the proportion of woody cover became less available within their home range, selection for woody cover increased by female ocelots. Males exhibited proportional use of woody cover in the 1985 group, but in the 1995 and the 2005 groups, they displayed patterns consistent with relaxed selection. Males in 2015 demonstrated increasing use and consistent selection of woody cover (slope > 1; Figure 2); regardless of availability, woody cover was always selected. Ocelots in later periods had higher proportions of woody cover available within their home ranges on average (0.49 in 1985 to 0.74 in 2015 for males, 0.55 to 0.69 for females). We found that both sexes exhibited proportional habitat use for non-woody cover across all decades, except for males in 1995 and 2015 ( Figure 3). Males in 1995 demonstrated relaxed selection, while in 2015, males exhibited consistent selection for non-woody cover. Most groups of ocelots demonstrated proportional use of areas farther from roads, regardless of traffic volume and with few exceptions. For areas farther from low-traffic roads, females in 2005 demonstrated relaxed selection, while males in 1985 exhibited a trade-off (Appendix S1: Figure S1). This trade-off indicates increasing use of areas farther from low-traffic roads, but males in 1985 selected areas closer at low availabilities and then avoided these closer areas at high availabilities. Finally, females in 2015 demonstrated consistent selection for areas farther from medium-(Appendix S1: Figure S2) and high-traffic roads (Appendix S1: Figure S3). There were differences in availability within ocelot home ranges for habitat and road variables across groups (sexes and time periods) based on the observed IQRs. Ocelot home ranges had higher availabilities of woody cover (IQR = 0.47-0.56) than non-woody cover (IQR = 0.21-0.24). Within ocelot home ranges, a broader range of woody cover was available than of non-woody cover. Ocelot home ranges had a greater range of availability of areas away from low-traffic roads (IQR = 5.37-8.21 km) than medium-(IQR = 3.36-4.65 km) or high-traffic roads (IQR = 3.27-4.66 km).
We observed temporal changes in habitat use across availability by male and female ocelots at the third order. Female use of woody cover in 1985 was significantly different compared to 1995 and 2015 (90% CIs around slope estimates did not overlap; Appendix S1: Figure S4). Females exhibited a pattern of proportional use for 1995, 2005, and 2015; however, during 1985, they displayed a pattern consistent with increasing use where selection was relaxed. Habitat use of woody cover was significantly different for males in 2015 from 1995 and 2005 (Appendix S1: Figure S4). Males in 1995 and 2005 exhibited increasing use but relaxed selection, while males in 2015 demonstrated increasing use and consistent selection. We discovered female ocelots consistently displayed a pattern of proportional habitat use for non-woody cover. Male use of non-woody cover was significantly different in 1995 from 2015 (Appendix S1: Figure S5). Males exhibited a pattern of proportional use for 1985 and 2005 however, during 1995, they demonstrated increasing use and relaxed selection, while males in 2015 increasingly used and consistently selected for non-woody cover as it became available within their home ranges. When evaluating how functional responses changed through time for roads, females in 1985 were significantly different in their use from 2005 in use of areas near low-traffic roads. Otherwise, male and female ocelots were relatively consistent in their use of areas near low-(Appendix S1: Figure S6), medium-(Appendix S1: Figure S7), and high-traffic roads (Appendix S1: Figure S8) in proportion to availability across time.

Habitat selection and mapping-Second order of selection
Ocelots demonstrated use of areas with a high percent of woody cover (≥0.48) across sexes and time periods (Figure 4) at second order. Ocelot use was greater than availability for areas with a greater percentage of woody cover on average than was available to an individual on the landscape. Both females (n = 37) and males (n = 39) used a greater percent of woody cover than was available, with a few exceptions. One female in the 1995 group and one male in the 1985 group (Figure 4) used woody cover less than available on the landscape. From 1985 to 2015, there was a continued decline in the proportion of woody cover available across our study site (0.44 in 1985 to 0.39 in 2015). However, there were no consistent changes in ocelot use patterns for males or females ( Figure 4); most ocelots (97.4%) used areas with a greater proportion of woody cover at the landscape scale. Ocelot use of non-woody cover was more variable (Appendix S1: Figure S9). A slight majority (54%) used non-woody cover T A B L E 2 Intercept and slope estimates with 90% CIs from linear models comparing the mean use of a resource as a function of availability of male and female ocelots for a given period (1985, 1995, 2005, and 2015 F I G U R E 2 The response (and 90% CI bands) to proportion of woody cover by female and male ocelots as a function of availability, as predicted using the regression model for each time period (1985, 1995, 2005, and 2015). The dashed line represents when use is equal to availability (proportional line) and the gray box represents the interquartile range of resource availability within home ranges for all ocelots (male and female) across periods.
greater than availability, while a few individuals used non-woody cover in proportion to (15%) or less than availability (31%). Ocelots tended to use areas farther from high-traffic roads relative to availability at the landscape scale across periods. Most females (n = 35) and males (n = 38) avoided areas closer to high-traffic roads. Females used areas that were 8.4 AE 2.9 km (SD) away from high-traffic roads and males used areas 7.98 AE 3.55 km away, compared to an average availability of areas 4.20 km away. Four females (three in the 1985 group and one in the 2015 group) and six males (five in the 1985 group and one in the 1995 group; Figure 8), however, used areas closer to high-traffic roads. Two females in the 1985 group demonstrated use consistent with availability (Appendix S1: Figure S10). Similarly, a majority of ocelots (81%) used areas farther from medium-traffic roads relative to availability (Appendix S1: Figure S11). However, there were temporal differences in ocelot use of areas around low-traffic roads (Appendix S1: Figure S12). For both sexes, a majority of ocelots in the 1985 and 1995 groups used areas closer to low-traffic roads, while all ocelots in the later periods of 2005 and 2015 used areas farther from low-traffic roads.
Most ocelots across all groups had individual models that showed positive selection (β > 0) for percent woody cover (95.3%) and avoidance of high-traffic roads (87.1%, Table 3). We estimated average effect sizes of β coefficients across groups of ocelots from the global models. All eight groups showed positive selection for woody cover and seven groups selected areas farther from high-traffic roads (Table 4). Average effect size of woody cover was positive for male and female ocelots across time (β woody = 0.423-13.479), while average effect size of distance from high-traffic volume roads was positive (indicating selection farther away from roads) for all groups except males in 1995 (β roads = À0.318-3.778; Table 4). Female ocelots in 1985Female ocelots in , 1995Female ocelots in , and 2005 as well as male ocelots in 1995 selected areas farther from medium-traffic roads (Table 4). All groups, except females and males in 2015, selected areas closer to low-traffic roads (Table 4). Additionally, we observed positive selection by all groups for areas with higher proportion of non-woody cover and greater heterogeneity of woody cover (Table 4); however, effect sizes were greatest for woody cover.
We used our RSF coefficients to map the relative probability of use per pixel for male and female ocelots F I G U R E 3 The response (and 90% CI bands) to proportion of non-woody cover by female and male ocelots as a function of availability, as predicted using the regression model for each time period (1985, 1995, 2005, and 2015). The dashed line represents when use is equal to availability (proportional line) and the gray box represents the interquartile range of resource availability within home ranges for all ocelots (male and female) across periods. across the study area during the four periods ( Figure 5). In general, across all groups, probability of use was positively associated with woody cover and negatively with high-traffic roads.
Across all groups, 63%-99% of telemetry locations occurred in the upper 80th and 90th quantiles of our equal-area characterization of ocelot habitat (Appendix S1: Figure S13). Consequently, there was a high correlation between the number of telemetry locations across quantile bins (ρ = 0.84-0.99). Ocelot telemetry locations were more likely to occur in the upper quantile bin of our maps, indicating satisfactory model fit.
We mapped areas of consistent high-use over time for male and female ocelots ( Figure 6). For females, approximately 47% of the study area was classified as high use for any one time period, 14% high use for any F I G U R E 4 Average use by female (a) and male (b) ocelots of woody cover (%) compared to average availability of woody cover on the landscape of the respective time frame in South Texas, USA. The y-axis depicts proportion of woody cover and the x-axis depicts individual ocelot identification codes. Triangles depict average use of woody cover for an individual ocelot. The solid line is the average woody cover available on the landscape with dashed lines representing a 95% CI. Log-distance to low traffic 3 0 0 9 4 2 0 9 Log-distance to medium traffic 11 7 5 2 8 13 1 2 Log-distance to high traffic 12 9 6 10 8 15 4 10 combination of two time periods, 2% for any three time period combination, and 0.15% for all four time periods. For males, approximately 45% of the study area was classified as high use for any one time period, 12% for any two time periods, 2% for any three time periods, and 0.0002% for all four time periods. Both sexes were predicted to consistently (in at least three time periods) use areas with high proportions of woody cover around the refuge and the private ranchlands in Willacy and Kenedy counties. Additionally, high-traffic roads were consistently predicted as non-habitat for both populations. We determined trends in predicted probability of use across time using a pixel-based correlation analysis (Figure 7). We discovered strong, positive correlations (ρ = 1) in predicted ocelot use across time for areas in the northern ranchlands of our study area as well as some areas in the refuge with high proportions of woody cover for both sexes. We observed strong negative correlations (ρ = À1) in predicted use for areas with low proportions of woody cover, indicating a decline in ocelot use of these areas over time. Some areas showed no correlation across probability of use, which indicates variation in predicted use by ocelots. Overall, areas farther from high-traffic roads with high proportions of woody cover were increasingly used by ocelots across time.
Finally, in our assessment of ocelot habitat by landownership, we discovered most high-use areas occurred on private lands compared to public lands in Texas (Figure 8). Across our characterizations of high-use areas, ≥79% of high-quality ocelot habitat occurred on private lands for both sexes.

DISCUSSION
Our study is the first to assess multiscale habitat use and quantify trends in spatial functional responses by ocelots across multiple decades. Ocelots demonstrated functional responses to vegetation cover and road-related attributes at the third order. At the second order, ocelots exhibited strong and consistent trends to use and select for areas of high woody cover and farther from high-traffic roads across time. The extensive anthropogenic disturbances in extreme southern Texas, which are projected to continue increasing within the LRGV , and high road mortality rates (Blackburn et al., 2021;Haines et al., 2005), warrant concern for ocelot persistence in the United States. Evaluation of temporal changes in habitat and behavioral responses by wildlife are particularly important for federally endangered species such as the ocelot. Collectively, this work provides a long-term lens with respect to ocelot ecology and habitat monitoring, which will serve an important role for developing conservation strategies for ocelots in the United States.
Ocelots demonstrated functional responses to vegetation cover at the third order. We hypothesized that ocelots would show stronger use for woody cover at lower availabilities, and as the resource became more abundant, use would decline and approximate proportional use. However, this prediction was not observed for most groups of ocelots. Further, habitat use varied across sexes and periods. We observed significant differences across time for male and female ocelot use of woody cover and male use of non-woody cover. Ocelots used areas within their home T A B L E 4 Effect sizes (mean AE SE) across groups of ocelots estimated from the two-stage group-averaged global models. ranges with greater proportions of woody cover than non-woody cover. Further, ocelots across time had a broader range of woody cover (0.47-0.56) available than non-woody cover (0.21-0.24) based on the IQR. We saw a shift in ocelot use of both cover types in 2005 and 2015, where ocelots used higher proportions in later periods. This shift in habitat use may coincide with a shift from VHF to GPS technology and more consistent sampling across the diel period. The temporal differences observed for non-woody cover may stem from differences in sampling effort and monitoring methods.
Ocelots are known to rest in dense woody cover during the day in Texas but show strong nocturnal activity patterns (Leonard et al., 2020). Since their diet is predominantly small mammals and lagomorphs in Texas (Booth-Binczik et al., 2013;Laack, 1991;Murray & Gardner, 1997), of non-woody cover was always lower than woody cover, increased ocelot use of non-woody cover in later time periods is likely more realistic given the improved technology of GPS-based tracking.
Ocelots are considered a forest interior species and use larger patches of dense woody vegetation communities in South Texas (Horne et al., 2009; and throughout their geographic range (Emmons, 1988;Wang et al., 2019). Average availability of woody cover decreased from 0.44 in 1985-0.39 in 2015 within our study area. While this decrease was minimal in an absolute sense, it follows broader habitat loss trends for the LRGV, which has a reduced and more patchy distribution of native thornshrub Tewes & Everett, 1986;Tremblay et al., 2005). Despite decreasing availability of woody cover over time, ocelots continued to use areas with a high proportion of woody cover, confirming that over time ocelots show consistent habitat use at the second order. Additionally, we observed strong and consistent selection for areas of high woody cover at the individual and group level with our RSFs. While we observed selection by all groups for greater proportions of F I G U R E 5 (Continued) F I G U R E 6 Predicted high-use areas for ocelots across time in South Texas, USA, from 1982 to 2017. We determined the binary cut-point for predicted high-use areas by identifying the cumulative percentile of the probability of use that captured 90% of ocelot telemetry locations (as described in Holbrook et al., 2017). We assessed the number of time periods where a 30 Â 30 m pixel was assigned high use. For instance, 1 Decade indicated a pixel was classified as high use for any one of the four time periods (1985, 1995, 2005, and 2015) for females (left) and males (right), and 2 Decades indicated that a pixel was classified as high use for any of the two time periods. Similarly, 3 Decades were pixels classified as high use for any three time periods, and 4 Decades indicated a pixel was classified as high use for all four time periods. Non-habitat represented pixels that were never classified as high use across any time frame.
F I G U R E 7 Calculated Spearman rank correlation coefficient (ρ, ranging from À1 to 1) per pixel for the probability of use for female (left) and male (right) ocelots in South Texas, USA from 1982 to 2017. This map characterizes spatial tendencies to increase, decrease, or remain variable with respect to predicted ocelot use across our four time periods. A value of ρ = À1 represented strong negative correlation and therefore a decreasing trend in probability of use (depicted in dark green), while a value of ρ = 1 indicated strong positive correlation and an increasing trend in probability of use across time (depicted in light pink and red). When ρ = 0, there was no correlation and this was interpreted as variability in probability of use across time (depicted in yellow). non-woody cover, habitat use was not as consistent. Effect sizes for RSFs were greatest for woody cover, and a majority of ocelots used woody cover greater than availability at this spatial scale. This indicates that ocelots may not use suboptimal areas and are seemingly confined into what remains of conventional habitat in the study area.
We hypothesized ocelots would decrease use of areas closer to roads, such that as roads became more abundant, ocelots would increasingly avoid areas closer to roads within their home ranges. However, most groups of ocelots demonstrated proportional use of areas farther from roads regardless of traffic volume, with minor exceptions. Male and female ocelots mostly used areas near roads of all traffic volumes in proportion to availability, which did not support our hypothesis that ocelots would avoid roads. This could be the mechanism resulting in continued ocelot-vehicle collisions. We observed all groups using areas farther from roads compared to availability across time based on the IQR for medium-and high-traffic roads. However, this could be an artifact of VHF data from earlier periods influencing the IQR for distance to roads (7.11 km in 1985 to 11.02 km in 2015 for high-traffic roads). Physical accessibility when sampling animal locations with VHF technology is a known source of potential bias (Hebblewhite & Haydon, 2010) and may be a contributing factor to the observed trends in ocelot use of areas near roads.
Roads influenced home range placement, with ocelots demonstrating a consistent temporal trend to use and select for areas farther from medium-and high-traffic roads. However, ocelots used and selected for areas closer to low-traffic roads; this is likely because of proximity of low-traffic roads to high-quality habitat Grilo et al., 2015). These patterns in habitat relationships could have been a function of where animals were trapped in combination with a saturation of territories, where territorial interactions affected behavioral patterns. This, however, is unlikely given ocelots are relatively rare in the United States (Janečka et al., 2016;Tewes, 2019) as well as the continued loss of habitat and high mortality rates from vehicle collisions since the 1980s (Blackburn et al., 2021;Haines et al., 2005).
Male and female ocelots selected for areas further from high-traffic roads when placing their home ranges on the landscape. At the third order, we found proportional use of areas farther from roads by ocelots, indicating that higher order selection processes (i.e., second order, home range placement) are likely dictating patterns of ocelot use at lower orders (i.e., third order, within home range). This pattern has important implications for ocelot conservation where individuals with greater availability of roads within their home ranges may be more susceptible to vehicle collisions. Since ocelots did not avoid roads as much as expected at third order, crossings structures and exclusion fencing should funnel movement into safer areas (i.e., higher proportions of woody cover and lower road densities; Blackburn et al., 2020). Strategically placing crossing structures within suitable habitat that connect to the broader landscape are a key factor in crossing success (Clevenger & Huijser, 2011;Forman et al., 2003).
Evaluation of temporal change in habitat and subsequent behavioral responses by wildlife is particularly important for federally threatened or endangered species. The extent of predicted consistent habitat never exceeded 47% (1515 km 2 ) of the study area for a given decade and never exceeded 14% (451 km 2 ) across 2 decades, illustrating the confined nature of ocelot habitat. Monitoring trends in woody cover, we saw a decrease in availability across the study area since the 1980s, yet ocelots strongly and consistently selected for areas with greater proportions of woody cover than average availability. We observed an increase in use of certain areas within the study area (i.e., northern ranchlands and refuge) across time. As woody cover decreased in other portions of the study area, ocelots increasingly used these strongholds of dense woody cover. These refugia of woody cover coincided with greater distances from high-traffic roads. We did observe some differences across time and breeding populations, likely stemming from different sampling methods. However, the F I G U R E 8 The percent area classified as high use by ocelots for time periods (1985, 1995, 2005, and 2015) overlapped by contemporary land ownership of public versus private land in South Texas, USA. The y-axis represents the percent area across ownership and the x-axis depicts the number of time periods the area was classified as high-use ocelot habitat. For instance, the x-axis label "4" illustrates the area that was classified as ocelot habitat across all time periods and is mapped (4 Decades) in Figure 6. spatiotemporal consistency in patterns of habitat use we observed reinforces previous work on ocelot spatial ecology (Harveson et al., 2004;Horne et al., 2009;Leonard, 2016). Importantly, we observed ocelots consistently avoided areas closer to paved roads, which were identified as non-habitat over time. This has important implications for the future of ocelot conservation in South Texas. Increased road density and expanding lanes to accommodate growing human populations in the valley could affect the availability of quality habitat for ocelots (Lombardi et al., 2021;. Of the area predicted as high-quality ocelot habitat over time, ≥79% occurred on private ranchlands. Texas is comprised of 96% privately owned land (Texas A&M Natural Resources Institute, 2014; United States Bureau of the Census, 1991). Large portions of the northern private ranchlands remained consistently high quality or were predicted to increase in relative quality for ocelots over time, indicating private ranchlands will likely act as important refugia for ocelots in South Texas (Tewes, 2019). Land use changes are impacting wildlife globally and monitoring long-term trends in habitat and animal space use will serve an important role for developing conservation strategies for threatened and endangered species.
Globally, biodiversity conservation has relied heavily on protected areas, which are usually owned by the public (Kamal et al., 2015). The Laguna Atascosa Wildlife Refuge has maintained consistently high-quality ocelot habitat over time, showing the importance of such protected areas for ocelot conservation. However, areas surrounding the refuge decreased in use over time based on our assessments and were identified as either marginally consistent habitat (1-2 years) or non-habitat. Thus, the refuge provides critical habitat that is increasingly isolated. Protected areas provide limited opportunity for conservation in areas such as Texas for ocelots because they constitute only 4% of total land area and are generally small, fragmented sites. A holistic approach to ocelot conservation will, by definition, require collaboration with private landowners (Kamal et al., 2015). Ocelots have been surviving on working ranchlands since at least the 1970s by coexisting with livestock (Janečka et al., 2016;Tewes, 1986;Tewes, 2019) and likely can continue to persist under a land sharing scenario for multiple-use ranchlands (Green et al., 2005). There has been significant work already to collaborate with land owners for habitat conservation and even translocation (The East Foundation, personal communication). Conservation strategies would benefit from a shared regional approach that aims to conserve the US ocelot, regardless of land ownership.
Ocelots are habitat specialists in South Texas, selecting for dense woody vegetation communities (Harveson et al., 2004;Horne et al., 2009;Jackson et al., 2005) and have been strongly linked to dense canopy cover across their geographic range (Emmons, 1988;Wang et al., 2019). For habitat specialists, higher order selection processes such as home range placement (second order; Johnson, 1980) can greatly influence environmental exposure and thus habitat use at lower orders (e.g., Holbrook et al., 2017;Morin et al., 2020;Rettie & Messier, 2000). Availability of a resource within a home range is inherently determined by home range placement, as resource selection is considered a hierarchical process (Johnson, 1980). Behavioral processes within the home range may be constrained by selection at higher orders (Rettie & Messier, 2000), indicating the importance of examining habitat use across multiple scales (Holbrook et al., 2017). Our results suggest that higher order selection processes (i.e., second order) may indeed influence and constrain ocelot habitat selection at the third order.
Despite widespread recognition of the importance of multiscale analyses of habitat relationships, a large majority of published work does not address multiple spatial or temporal scales (McGarigal et al., 2016). Our combined analysis provides a refined lens of ocelot habitat use in the context of habitat availability over multiple spatial and temporal scales as well as trends in land ownership. Previous work on threatened felids has shown that understanding habitat relationships across multiple scales is critical for protecting habitat and landscape connectivity (Macdonald et al., 2018). We observed multiscale habitat selection patterns indicative of a habitat specialist, similar to other taxa (e.g., Mauritzen et al., 2003;Rettie & Messier, 2000) including medium-sized felids (e.g., Canada lynx [Lynx canadensis]; Holbrook et al., 2017Holbrook et al., , 2019Morin et al., 2020). Our work provides context for multiscale habitat relationships for an endangered habitat specialist that is directly threatened by human impacts on the landscape over several decades. Our work represents one of the few studies that has assessed habitat selection and predicted habitat use across multiple decades for wildlife populations (e.g., Anderson et al., 2012;Harmange et al., 2019;Lone et al., 2018;Stewart et al., 2012). By coupling a long-term telemetry dataset on an endangered carnivore and an expansive time series of satellite imagery, for the first time we described ocelot habitat use and selection across its known distribution in the United States. We demonstrated the benefit of employing a long-term telemetry dataset and satellite imagery to evaluate changes in habitat use (i.e., functional responses), predict areas of consistently high-quality habitat, and evaluate trends in habitat by landownership for an endangered carnivore along the United States-Mexico border. This study leveraged the power of time series and long-term research on ocelot behavior to inform future conservation efforts. Our process can serve as a model for future habitat assessments over long time periods for threatened and endangered species.