High frequency GPS bursts and path-level analysis reveal linear feature tracking by red foxes

There is a need to quantify and better understand how wildlife interact with linear features, as these are integral elements of most landscapes. One potentially important aspect is linear feature tracking (LFT), yet studies rarely succeed in directly revealing or quantifying this behavior. In a proof-of-concept study, we employed short-term intensive GPS monitoring of red foxes (Vulpes vulpes) in a multiple-use landscape in southern Norway. Using periodic bursts of high frequency GPS position fixes, we performed modified path selection analyses to estimate the propensity of foxes to track natural and man-made linear features (roads, forest edges, and streams) once they are encountered. Foxes in our study tracked primarily forest edges and roads. Forty-three percent of bursts that encountered any linear feature resulted in LFT. LFT, although prominent, was manifested as a short-lived behavior, with overall median times to linear feature abandonment around two minutes. Movement speeds were highest along roads, perhaps due to greater ease of travel or higher perceived risk. In the highly heterogeneous habitats that characterize human-dominated landscapes, LFT may be manifested at such a fine spatio-temporal scale that it would remain hidden during telemetry studies employing conventional position fix frequencies. The approach described here may aid others studying spatial behaviors that are manifested over very short durations, yet are biologically significant.

LFT Does linear feature type influence the speed of movement during LFT? Prediction: Movement speeds will be higher along roads than other linear features due to the greater ease of travel, but also greater risk of, and exposure to, human disturbance.
The high spatio-temporal resolution yielded by the telemetry bursts, in combination with path-level analysis, allowed us to detect LFT and individual heterogeneity therein, which would have remained obscured or entirely hidden from investigations using conventional position intervals.
Methods study system. The study was carried out in southern Norway, in the municipalities Vestby (59°34.0′-59°38.3′N, 10°38.8′-10°43.5′E) and Ås (59°39.6′-59°42.7′N, 10°43.8′-10°47.7′E), at elevations 0-152 m.a.s.l (Norwegian Mapping Authority 2018; https://norgeskart.no/). In addition, one of the foxes moved south of this core study area, reaching Fredrikstad (59°13′N-10°56′E). About 50% of the study area is covered by forest managed for timber or firewood production (Statistics Norway 2017; https://www.ssb.no/), mainly mixed conifer-deciduous boreal forests. The dominating tree species are Norway spruce Picea abies and Scots pine Pinus sylvestris, whereas birch Betula spp. are the most common deciduous trees, followed by rowan Sorbus aucuparia, Salix spp. and European aspen Populus tremula. The entire landscape is influenced by human activities and land use (Fig. 2), with 31% cultivated land, 18% of other open areas and a high road density. Human population density is 159 residents/km 2 , mainly concentrated in small town centres or residential areas, but also on farms and single houses scattered throughout the so-called "cultural landscape". The result is a fragmented landscape with forest patches of varying sizes interspersed with crop fields, pastures, and human settlements. In addition to linear landscape features like forest edges and streams, purely anthropogenic linear elements like roads are common. Fox hunting is legal in the study area throughout www.nature.com/scientificreports www.nature.com/scientificreports/ the year except during the season when females are with dependent young (April 15 -July 15). On average, 0.5 foxes/km 2 /year were killed in 2013-2017 in the study area (Statistics Norway 2017; https://www.ssb.no/). The climate is Temperate-Continental with an average temperature in the coldest month (February) of −1.9 °C/−4.7 °C in 2017/2018 (Norwegian Metrological Institute 2018; https://www.met.no/).

Gps collars.
We developed and built custom GPS devices, specifically designed using off-the-shelf components together with a customized printed circuit board and a standard micro controller. The backbone of the electronics is a SIM808 GSM/GPRS module from SimCom ™ featuring built in GPS, General Packet Radio Service (GPRS) and Bluetooth functionality. This module is controlled by an 8 bit ATmega 328p micro controller running our own software, optimized for the application with main focus on robustness and reliability in harsh environments. Once active, the collar can be fully controlled by a predefined subset of commands through standard SMS messages, thus allowing any type of basic cellphone to be used as a communication device. The commands range from requesting current battery level and signal strength to changing the sampling/burst rates and the default behavior when GSM and/or GPS signals are lost.
Raw data were sent immediately over GPRS to a dedicated server and stored. GPS units and batteries (lithium polymer; 3000mAh) were housed in 3D-printed plastic cases (7 cm × 4 cm × 4 cm), which were attached to a collar made of a 2 cm wide and 1 mm thick plastic coated strap (BioThane ® ). A short cotton string was incorporated into each collar as a wear-and-tear-type drop off mechanism. Assembled GPS collars weighed on average 123 g, which was less than 2.3% of the average fox body weight (5.3 kg) in our study.
Gps position error determination. In order to quantify the true position error of our custom-built GPS devices during movement, we conducted empirical tests using a high-end APX-15 GNSS/IMU unit. The reference unit was mounted on the roof of a car together with two different GPS collars, and a test route was driven on narrow forest roads with varying canopies and at slow speed (1-5 m/s). Due to the inertial measurement unit (IMU) the reference unit is able to maintain its high accuracy at the level of a few centimeters even when the main GNSS antenna is temporarily obstructed. The raw data collected from the reference unit were stored and post processed for the whole trajectory yielding an overall precision for the reference trajectory of about 10 cm. The median true error of our GPS collars was 2.4 m (95% confidence interval, hereafter CI: 0.6 m-10.1 m), based on a sample of 300 positions.
Fox capture and handling. Foxes were captured in large wooden box traps (approximately 200 cm × 80 cm × 80 cm) with two trapdoors, which are legal and standard traps in Norway for live capture of small and medium-sized carnivores. Traps were baited with meat (primarily dead chickens Gallus gallus domesticus) and monitored daily. Captured foxes were restrained with the help of gloves, a catch pole, and/or neck tongues. In addition to GPS collar attachment, we recorded the sex and weight of each fox, and collected hair samples as voucher specimen and for future DNA analysis. Handling, between removal from the trap and release ( Fig. 1) lasted between 10-20 min. The animals' eyes were covered with a dark cloth during handling to reduce stress, as handling was performed without anesthesia. All capture and handling conformed to the current laws and regulations in Norway and the study was approved (case IDs 2016/4769 and 18/211316) by the Norwegian Animal Research Authority (FOTS) under the auspices of the Norwegian Food Safety Authority.
Gps tracking and data preparation. GPS collars were programmed to take position fixes in bursts of 20 (15 second inter-fix interval), with an inter-burst interval of 10-20 minutes. The inter-burst interval increased to 60 minutes once the battery had lost approximately 60% of its original charge. Positions were transmitted immediately after collection via GSM network to an internet server, from which they were downloaded for processing and analysis. To minimize the potential initial effects of capture and handing on movement and behavior, we excluded positions collected during the first day after each fox was released. Furthermore, since we were interested in position data collected during movement/travel, we removed those portions from the data that indicated stationary or localized behavior, e.g., time in dens and day beds or time spent investigating focal points. For each burst, we calculated sinuosity 42,43 and average movement speed between consecutive positions. We then excluded all bursts with sinuosity >0.5 and average movement speed <0.1 m/s. Data preparation and subsequent analyses were conducted using R version 3.5.1 44 .
Linear features. Forest edges, roads, and streams are the most prominent linear features in our study area (in that order of prevalence) and were therefore selected in our test for LFT. We utilized contemporary (2016) high-resolution (5 m) land cover data (AR5; https://www.nibio.no/tema/jord/arealressurser/arealressurskart-ar5) to identify edges between forest and other (primarily open) land cover types (Fig. 3). We obtained GIS data of streams and roads through the Norwegian Mapping Authority's geographic data portal (kartkatalog.geonorge.no).
Linear feature tracking (LFt). We used a modified burst-level path selection analysis to quantify the propensity of individuals to stay in close proximity to linear features once these had been encountered. We first selected bursts with at least one position within 10 m of a linear feature (i.e., road, forest edge, or stream). This proximity buffer width was chosen in order to balance the need to capture truly close association with the linear feature, i.e., the animal was either walking directly on the feature or right next to it, with the GPS position error (upper 95%CI limit of 10.1 m). We then considered all positions in a burst following and including the first position within the proximity buffer as an observed ("case") trajectory. For each case trajectory, we simulated 10 random trajectories ("control") with step lengths and relative turning angles sampled from the empirical distributions of these parameters from the respective case trajectory (Fig. 4). This represents a null model based on a correlated random walk 45 . We then fit conditional logistic regression models (R package "survival" 44,46 ) to the resulting data with observation type (case vs. control) as the response and whether a position was located within (2019) 9:8849 | https://doi.org/10.1038/s41598-019-45150-x www.nature.com/scientificreports www.nature.com/scientificreports/ 10 m of a linear feature or not as the predictor. This approach tests whether the true trajectory taken by the fox is more or less likely to exhibit LFT than alternative realizations of the true trajectory. We repeated the analysis with buffer widths of 5 m and 15 m, without qualitative changes to the main results. All analyses were performed in R version 3.5.1 44 .
Because linear features may co-occur (e.g., roads through forests also represent forest edges), all linear features were included as dummy covariates in the same model (coded 0/1 for inside and outside the distance threshold, respectively). As a result, any observation (case and control) could be placed in various combinations of overlapping linear feature buffers or outside all buffer. For simplicity, we chose to use only additive effects, but this approach could also be used to accommodate interactions. Burst ID (i.e., case trajectory ID) was included as a stratum 47 and the sequence ID of a position as a clustering variable. Separate models were fitted for each individual fox as we were interested in individual differences in the propensity of performing LFT and because models differed between individuals (not all individuals encountered all linear feature types), thus precluding the commonly used two-step approach 48 . In each individual model, we included only those linear features types that were encountered during at least 5 different bursts. Animals may perform continuous LFT by stitching together movement along multiple feature types. Therefore, in addition to the models that differentiated between linear feature types, we also ran models with roads, forest edges, and streams combined into one linear feature ("any", Fig. 3). Our approach to path-level selection analysis differs from most applications of path analysis, which generally employ random rotation and shifting of empirical paths for null model generation [49][50][51][52][53] . We chose the correlated random walk, based on empirical distributions of turning angles and distances, to generate our null model  www.nature.com/scientificreports www.nature.com/scientificreports/ (Fig. 4). In our opinion, this more closely reflects the process leading to realized path topology as a sequence of small-scale movement choices, rather than a predetermined fixed path shape. Furthermore, GPS positions were obtained every 15 seconds during bursts. At this scale, spatial autocorrelation is very high and pure step selection is not a meaningful test of LFT. Most individual steps are short and, in the case of trajectories closely tracking the line, most or all simulated control steps may fail to leave the LFT buffer around the linear feature. Association with linear features indicates possible LFT, but, even at the fine temporal scale of our data, it is not direct proof of it. Individuals may stay close to linear features by moving along them, but could also perform other behaviors near a linear feature without necessarily tracking it. In order to quantify the incidence of LFT from the GPS data and describe attributes of LFT, we defined an LFT event as a contiguous sequence of positions that a) all fell within the 10 m proximity buffer, b) consisted of at least 3 positions, c) represented a total movement distance of at least 35 m, and d) exhibited displacement-to-distance ratio of >0.5. The displacement-to-distance ratio was calculated as the average displacement of all subsequent points from the first position in a sequence divided by the average distance moved between points. Setting a threshold on this measure removed position sequences within the buffer that represented movement without apparent direction rather than actual travel.

Duration of LFt (t LFT ).
We used time-to-event analysis 54 to estimate the duration of LFT events on a given linear feature type. Although distance to event analysis 55 is a feasible alternative for quantifying the endurance of LFT, we were interested in relating the results to the temporal component (schedule) of the GPS position fix strategy. Only position sequences identified as LFT (see above) were included in this analysis. An event (departure from LFT) was defined as the instance when the LFT sequence was followed by a position within the same burst that exceeded the 10 m distance threshold (event = 1). LFT sequences that were not abandoned before the end of the burst were right censored, i.e., these LFT sequences were considered available for LFT abandonment during analysis, but the event was considered to have occurred after the observation period. The time variable was defined as the number of seconds from the first position of a sequence until either censoring or departure occurred. In order to further minimize bias in time-to-event estimates, we excluded LFT events for which the starting point was unknown, i.e., when the first point in the burst was already within the buffer. We then constructed Kaplan-Meier time-to-event curves for each linear feature type and estimated median time-to-event, i.e., the time by which 50% of LFT events had not yet ended, using the survfit function in the R package "survival" 46 . We used bootstrapping to derive 95% CI around median times-to-event, sampling at both the individual and burst level to account for non-independence. Although median time to event is directly linked with the proximity threshold, patterns in relative hazard (differences in the instantaneous potential of abandoning/ending an LFT event) with respect to feature type are unlikely to change drastically as long as the distances picked represent a reasonably close proximity to the linear feature. We used Cox proportional hazards models using the coxph function in the R package "survival" 46 to test for differences in the time-to-event associated with different feature types, with feature type as the predictor variable, individual ID as strata variables, and burst ID as cluster variables.

Movement speed during LFt (v )
LFT . We calculated average feature-dependent travel speeds for each fox based on all LFT events identified on or along a given linear feature type. Ninety-five percent confidence limits were derived by bootstrapping. We used linear mixed regression with R package lme4 56 to test for and quantify differences in speeds along different linear feature types, with the log-transformed average speed during an LFT event as the response, feature type (roads, forest edges, and streams) as fixed effect, and burst ID nested within individual ID as random effect.

Results
Gps tracking. Eighteen foxes (7 females; 11 males) were GPS tagged during this study, including six that could be classified as juveniles based on size and appearance (Supplementary Table S1). All foxes were GPS tracked for a short duration (5-28 days), but with high intensity (1123-4617 positions per fox). Depending on the individual, between 35% and 94% of positions (after cleaning and removal of the first 24hrs after marking) were categorized as "traveling" based on sinuosity and speed and thus included in the analysis (Supplementary Table S1). Negative coefficients, suggesting that prolonged proximity to linear features was avoided once they were encountered, emerged in four cases (roads: 1 individual; forest edges: 2 individuals; streams: 1 individual; Fig. 5A). Of bursts that encountered roads, forest edges, streams, or any linear features, 41%, 39%, 22%, and 43% resulted in LFT, respectively (Figs 5B, 6, Supplementary Table S2). LFT occurred also in fox/feature type combinations that were associated with non-significant or negative selection coefficients (Fig. 5A).

Discussion
Our study revealed pronounced linear feature tracking by red foxes in a highly fragmented landscape. High frequency GPS bursts and a path-level null model approach allowed us to detect fine-scale behavior that would have remained hidden had we employed conventional fix frequencies in the range of several minutes or hours. The incidence and strength of LFT varied among foxes and feature types. LFT was most pronounced in association with forest edges and roads (Fig. 5). The duration of LFT (t LFT ) was shorter along streams than along roads and forest edges. Travel speeds (v ) LFT were on average higher along roads than along forest edges and streams. We found strong evidence of linear feature tracking in foxes, confirming our first prediction (Figs 5 and 6). The presence of linear features such as roads can have contrasting effects at different spatial scales. Individuals can avoid areas with high road density at large spatial scale (e.g., when establishing a home range 57 ), but may select for them during daily activities at a smaller scale (e.g., moving 58 ). Here we found that foxes selected roads at a very fine scale by performing LFT. Between 60% and 72% of foxes exhibited a propensity to remain in close proximity to linear features once they were encountered. Linear feature tracking events occurred even in fox/feature combinations that did not emerge with significant positive coefficients from the path analysis (Fig. 5). Overall, apparent linear feature tracking events occurred in 43% (95%CI: 37-48%) of bursts that came in proximity of either of the linear feature types used in our study.
While all foxes apparently tracked linear features occasionally, time-to-event analysis indicated that the median time until a linear feature was abandoned once LFT started was relatively short, 120 seconds (95%CI: 90s-165s) along the combination of all three linear feature types. Although longer LFT events also occurred, these were rare. The brief duration of LFT is likely attributable to a combination of species-specific and environmental characteristics. Our study area is highly fragmented, with an abundance of intersecting and tortuous linear features (Fig. 2). High levels of human disturbance/activity, including intensive hunting, the risk of intra-guild  Supplementary  Table S1. www.nature.com/scientificreports www.nature.com/scientificreports/ predation posed by Eurasian lynx Lynx lynx 59 and a generalist foraging strategy of foxes, may further explain their variable and tortuous movement trajectories.
Using linear features, which are generally characterized by open habitat, allows foxes to move faster, as shown for wolves 10 . However, LFT may also make foxes more detectable to predators and hunters, and predation risk www.nature.com/scientificreports www.nature.com/scientificreports/ influences habitat selection 60,61 . The recurrent, but short-lived LFT exhibited by foxes may thus illustrate the fine scale behavioral decisions of a medium carnivore in a landscape of fear 62 , where foxes trade off ease of travel with predation and hunting risk.
Given the short duration of LFT by foxes in our study area, the phenomenon would likely remain hidden from telemetry studies employing conventional position fix rates in the order of hours or even several minutes. Our study was located in a fragmented, human-dominated landscape; other study systems with less tortuous linear features and more homogeneously structured landscapes 14 may readily detect LFT using much longer inter and within burst intervals, perhaps even omitting bursts all together. We also expect that longer monitoring durations (months instead of days), accompanied by an increase in the total number of GPS positions, will improve detection of selection for linear features. Regardless of the study system, evidence of association within the individual home range (e.g., through conventional step selection analysis) is not a direct evidence of LFT, which requires that an animal's path coincides with, not only periodically intersects, an underlying geographic feature. To detect true LFT, the temporal grain of the data must be fine enough to allow detection of intermittent deviations from a particular feature.
Developments in GPS tracking technology 63 and the increasing availability of high-resolution spatial environmental data, have motivated a growing number of studies to employ position fix frequencies in the order of minutes, rather than hours or days, or a combination of both 64,65 . This provides opportunities to better understand the interactions between linear features and wildlife. For example, our approach allowed us to detect and quantify LFT and not only the selection for linear features 30,58 . In addition, our framework could be extended to help distinguishing the intrinsic and extrinsic drivers that prompt individuals to perform LFT. This would involve the use of individual and/or site specific covariates in the conditional logistic regression when quantifying LFT, or in the time to event analysis when quantifying time to departure from the linear feature.  www.nature.com/scientificreports www.nature.com/scientificreports/ LFT is not the only fine-scale behavior that may be missed or its prevalence underestimated even at moderately high telemetry fix frequencies. Clusters of positions are often used to locate and subsequently investigate specific behaviors such as resting 32 , foraging 66 , and to study kill rates 67 , inter-specific interactions 68 , or responses to human disturbance 32 . Some key behaviors may only take minutes or even seconds to execute, making them virtually undetectable with conventional telemetry schedules. For example, studies quantifying predation and its spatial determinants may be biased towards larger prey items that require longer handling times and thus are more likely to result in observable position clusters, whereas small prey may go undetected 69 . Similarly, high frequency position data and a null model approach could help quantify how infrastructure mitigation structures (such as wildlife crossings on roads) are used by wildlife and, assuming appropriate scope of the study, how this use may scale up to the level of landscapes and populations.
Being able to obtain fine scale data, "may give us the potentially false impression that fine spatial and temporal scale dynamics are relevant to ecology or, most critically, conservation" 70 . We argue that LFT, or other behavior that may be so short-lived that GPS intervals lasting seconds or minutes are needed to detect and properly describe it, are relevant because they can have key impacts on life history and populations in general. For example, such behaviors can be associated with significant spikes in foraging opportunities, which in turn can be important to maintain population density 71 . Even if linear features are not selected for or are actually avoided within the home range, animals may utilize them once they are encountered. Overall, in highly fragmented human-impacted landscapes such as our study area, even random walk-type movement or moderate avoidance would lead to frequent encounters of linear features and thus opportunities for interactions with them at fine scales. In such fragmented scenarios, fine temporal and spatial scales can help us better understand the mechanisms involved in wildlife behavioral responses, which in turn can help predict the effects of human-induced environmental change on species and communities 72 . The increasing availability of fine scale environmental information, such as the detailed land cover maps used in this investigation, also boosts the utility of and opportunity to use fine scale movement data 73,74 .
Managing the tradeoff between resolution and longevity is an ubiquitous challenge in wildlife telemetry 33,75 . Configuring GPS fix schedules into a sequence of bursts provides a practical way to balance the need for fine-scale information about animal movements with the simultaneous desire for home-range level inferences. Technical characteristics of wildlife GPS telemetry devices and the goals of a given study will guide the choice of inter-burst interval and the number and frequency of position fixes within bursts. For increased utility and power efficiency, the initiation of bursts could be made conditional on times, positions, or inferred activity, with the help of schedules, proximity rules (e.g., geo-fencing), and integrated accelerometer data. Even in cases where scheduling is not constrained by other considerations, GPS position accuracy will provide a lower limit for the scale at which inferences can be drawn 34 . We concede that, although bursts provide a reasonable tradeoff between the desire to obtain position information at both fine and coarser scales, investigators should be cautious when designing telemetry schedules. Study designs well-suited for targeting one or a few questions are generally preferable over designs that attempt to answer many questions at once. The main goal of the study should determine the most appropriate position fix interval 34,35 , within the overall constraint of battery life.
Finally, technical and analytical developments in peripheral devices integrated into wildlife GPS collars may allow for inferences about behavior at fine scales, without the need for very high GPS fix rates. Accelerometers can be used to prompt periods of higher fix rates 33 or to reconstruct detailed movement paths 76 . Similarly, telemetry collars with integrated photo and video cameras could help detect direct evidence of fine-scale behavior such as LFT. For instance, intensive fieldwork is necessary to study prey selection and kill rates of carnivores 77 , and neck-mounted cameras help estimate kill rates accurately 78 .

Conclusions
Linear feature tracking, especially in carnivores, is frequently mentioned, but has rarely been tested for or quantified as such. In our study, the use of GPS bursts, in combination with a path-level null model approach, allowed us to detect and quantify linear feature tracking occurring at a fine spatial and temporal scale in a highly fragmented landscape. We also detected individual variation in the propensity to perform LFT; we encourage further studies to determine potential individual and environmental drivers of this variation. While the importance of fine-scale data may at times be overemphasized 70 , increasing fragmentation in human-dominated environments may require finer-scale movement decisions, and, in turn, observations of sufficiently high spatio-temporal resolution. Another requirement is the availability of analytical procedures that do not only account for the autocorrelation inherent in fine scale data, but make use of the signal about underlying structure contained therein 52 . Further developments in telemetry technology and its application, together with flexible analytical approaches, will help study behaviors that are manifested over very short durations, yet are biologically significant.

Data Availability
Data used in this analysis are available from the corresponding author upon request.