Canada lynx Lynx canadensis are federally listed in the contiguous United States as a threatened species under the Endangered Species Act (U.S. Fish and Wildlife Service 2000). Therefore, determining the presence or absence of lynx across broad landscapes, such as wilderness areas or national forest units, is important to conservationists and land managers. The presence of carnivores can be especially difficult to establish as they are often cryptic, exhibit crepuscular activity patterns and occur at low densities across large spatial scales (Wilson & Delahay 2001, O'Connell et al. 2006, Kolbe & Squires 2007).

Delineations of current and historical ranges for rare or sensitive species are often based on anecdotal occurrence data (unverifiable observations or sign; McKelvey et al. 2008). Using anecdotal data to delineate species occurrence leads to errors of omission and commission that result in a misallocation of limited funding and inefficient conservation actions (McKelvey et al. 2008). Snow-tracking to detect carnivores is particularly effective because it does not require a solicited response from the animal, it documents the presence of individuals over a span of several days if snow conditions permit, and it is generally low cost compared to live-trapping or remote cameras (Reid et al. 1987, Halfpenny et al. 1997, Foresman & Pearson 1998). Snow-track surveys have been used successfully to determine the presence of a variety of carnivore species, including wolverine Gulo gulo, river otter Lutra canadensis, red fox Vulpes vulpes and short-tailed weasel Mustela erminea (Reid et al. 1987, Edelmann & Copeland 1999, Forsey & Baggs 2001, Ulizio et al. 2006, Heinemeyer et al. 2008). However, a chronic weakness of monitoring strategies based on snow-tracking is species misidentification. This shortcoming can be addressed to meet evidentiary standards of species identification by considering the snow-track a ‘collection device’ for obtaining genetic samples (Squires et al. 2004), a method demonstrated to be effective for lynx (McKelvey et al. 2006) and other forest carnivores (Ulizio et al. 2006).

Monitoring efforts based on snow-track surveys conducted on a standardized survey grid can be used to delineate local distributions of lynx (Squires et al. 2004). In the Northern Rocky Mountains, Montana, USA, delineations of local lynx distributions as determined by snow-tracks overlapped the distribution of lynx as estimated by radio-telemetry by 97% (Squires et al. 2004). Track counts can also estimate animal population density if track age and estimated daily travel rates are known (Stephens et al. 2006). However, using snow-track surveys to establish either the presence or absence of species in potential habitat may be more difficult. Although species presence can be ascertained by a single detection (assuming correct identity), absence is difficult to establish with statistical certainty. A variety of statistical concepts, including occupancy models, species distribution models and mark-recapture models, provide statistical frameworks for modeling species occurrence across landscapes (Pollock et al. 1990, McKenzie et al. 2003, Elith et al. 2006). These models are data intensive and require significant financial and time commitments. Since our primary goal was to provide managers with a relatively accessible method of estimating probability of detection and presence/absence of lynx, we chose to use a more tractable model by McArdle (1990) and Reed (1996). This model uses the estimated probability of detecting a given species to determine the number of sampling events necessary to establish presence or absence with known statistical certainty, and has proven to be effective in determining the number of sampling events necessary to detect species that occur at different densities (Kery 2002, Jackson et al. 2006).

In a previous study, we demonstrated the efficacy of snow-tracking for delineating local lynx distributions (Squires et al. 2004). In that study, we used computer simulations to model the probability of detecting lynx snow-tracks depending on the amount of search effort. Our goal in this study is to expand on this earlier work by calculating the empirical probability of detecting lynx based on snow-track surveys in areas of known presence. We focus on detectability within areas occupied by resident populations of lynx rather than in areas occupied by single dispersing individuals to ensure that underlying lynx densities are relevant for conservation. We address the following two questions: 1) What was the probability of detecting lynx with snow-track surveys conducted in areas with confirmed resident populations? 2) Given this estimated detection probability, what level of survey effort would be necessary to infer the absence of resident lynx populations with confidence? Based on these results, we offer management recommendations to optimize survey methods for monitoring lynx using snow-tracking methodology.

Material and methods

We conducted our study in two regions of western Montana: the Purcell Mountains in the northwest corner of the state and the Seeley Lake area, including the Garnet Range, located approximately 50 km to the southwest (Fig. 1). The boreal forests in the Seeley Lake area and the Garnet Range are composed mainly of subalpine fir Abies lasiocarpa and Engelmann spruce Picea engelmannii, with a smaller component of lodgepole pine Pinus contorta and western larch Larix occidentalis; ponderosa pine Pinus ponderosa and Douglas-fir Pseudotsuga menziesii were dominant at low elevations. Forests in the Purcell Mountains are similar with the addition of western red cedar Thuja plicata and western hemlock Tsuga mertensiana in mixed stands. Elevations ranged from 800 to 2,300 m a.s.l. in the Purcell Mountains and from 1,200 to 2,500 m in the Seeley Lake area and Garnet Range.

During 2006 and 2007, we conducted winter snow-track surveys for lynx on a grid of 8 × 8 km squares following Squires et al. (2004). This grid size roughly equals an average female home range (of 72 km²) for lynx from southern boreal forests (United States and southern Canada; Aubry et al. 2000, Squires et al. 2004). We used VHF-telemetry to verify the presence of lynx in our study areas prior to surveys; however, the technicians did not know the actual location of lynx when they conducted their surveys. We selected survey routes prior to fieldwork by overlaying our 8 × 8 km grid across occupied lynx habitat (i.e. all grid cells were in occupied home ranges). Survey routes followed existing forest roads that bisected grid cells through as much dense forest as possible. Technicians used the same basic survey routes in both years and searched for lynx tracks for an average of 10 ± 2 (SD) km per grid cell from snowmobiles traveling approximately 30 km/hour. We employed the same technicians for both years of the study to address potential observer bias when detecting lynx tracks (Weckerly & Ricca 2000). Track simulations from a previous study suggested that surveys conducted three days after a track-obliterating snowfall had a 0.55 probability of successful detection (Squires et al. 2004). Thus, we conducted surveys on an average of three days (SD = 1 day, range: 0-6, N = 83) following a snow event to ensure adequate time for animals to traverse home ranges. We plotted the survey route and location of each lynx detection using Trimble GeoExplorer® 3 data logging GPS units (Trimble Navigation Ltd., Sunnyvale, California, USA).

In the Purcell Mountains, we surveyed 12 grid squares in 2006 and again in 2007. In Seeley Lake, we surveyed nine grid squares in 2006 and again in 2007. Surveys were replicated twice in each year for a total of four surveying sessions per area; all surveys were conducted between January and March. We evaluated locations of lynx track detections and survey routes with ArcGIS 9.2 (ESRI, Redlands, California, USA). We used Q-Q (quantile-quantile) plots to determine whether lynx track detections from different study areas and years warranted separate analyses or could be combined. The Q-Q plot clearly showed a difference in the distribution of lynx track detections between the two study areas (Fig. 2), while the distribution of lynx track detections between years appeared to be similar. To confirm, we performed a non-parametric Wilcoxon rank sum test to compare the distributions of lynx track detections per grid cell between years, which supported our decision to analyze years together (2006: N = 42, median = 5, 2007: N = 41, median = 4, W = 924.5, P = 0.56). We used GPS to quantify the length of survey routes and the number of lynx track detections to determine the average survey distance per detection. We also measured the distance between the start of the route and the first detection of lynx tracks to determine the minimum survey distance to the first track detection. We considered all track detections unique unless observers could clearly see that two tracks were from the same individual without leaving the survey route.

We used Q-Q plots and goodness-of-fit tests to compare observed distributions of the number of tracks per grid cell with the Poisson and negative binomial distributions; two common distributions for count data. The negative binomial distribution provided a good fit to the data (Fig. 3) as expected since the observed data were overdispersed (mean < variance). The negative binomial is also sometimes useful for data with correlated observations (Bliss & Fisher 1953, Barron 1992). We then performed a bootstrap procedure to sample the observed number of tracks per grid cell 10,000 times with replacement, in order to calculate the probability of track detection as well as a rigorous estimate of variability (calculated using the ‘boot’ package in the program R 2.10.1; R Development Core Team 2008). In addition, we maximized independence between survey cell replicates by using survey cells consistent in size with lynx home ranges. For each bootstrap sample, we calculated the probability of detecting zero lynx using the negative binomial distribution, as well as the standard error and a 95% confidence interval using the adjusted bootstrap percentile interval. We then used one minus the estimated probability of detecting zero lynx as the estimated probability of detecting one or more lynx tracks in a grid cell.

We also estimated the cumulative probability of detecting lynx tracks per kilometer surveyed based on the distance between consecutive track detections along survey routes. We performed a bootstrap procedure to sample the distance from one detection to the next among all survey routes 10,000 times with replacement, and used the resulting empirical cumulative distribution function to estimate the probability of detecting one or more lynx per kilometer surveyed. We omitted the distance from the start of the survey to the first lynx track detection, since the start point was arbitrary and did not accurately reflect distance between lynx tracks. We also calculated a bootstrap estimation of standard errors and 95% confidence intervals for these values using the adjusted bootstrap percentile interval from the ‘boot’ package in R.

To determine the number of survey replicates statistically necessary to infer absence, we used the following equation, developed by McArdle (1990) and Reed (1996):

Since one of the main concerns of snow-tracking methodology is track misidentification, we performed a sensitivity analysis to determine the potential impact of track misidentification on the probability of detecting at least one lynx in a given grid cell, as well as the number of visits required to infer absence. We used the actual number of lynx track detections per grid cell, reduced by a range of percentages, to simulate a reduction in actual tracks located due to misidentification. We then used the negative binomial bootstrap procedure detailed above to calculate the new probability of detection and the associated number of survey visits necessary to infer absence.

Results

In the Purcell Mountains, we surveyed 48 routes in 12 grid cells (24 in 2006 and 24 in 2007) for a total of 489 km surveyed over all sampling events and years. We detected lynx tracks 416 times during the two years of sampling. Average survey distance per lynx detection was 2.27 km (SD = 2.29 km, range: 0.39-10.42 km, N = 48). The distance to first track detection for lynx averaged 1.8 km (SD = 1.9 km, range: 0.02-7.8 km, N = 48). We encountered an average of 8.5 lynx detections/10 km of survey distance.

In the Seeley Lake area, we sampled 35 survey routes in nine grid cells (18 in 2006 and 17 in 2007) and surveyed 371 km in total. We detected lynx 130 times during the two years of sampling. The average survey distance per lynx detection was 4.29 km (SD = 3.51 km, range: 0.7-14.9 km, N = 32). In the Seeley Lake area, we surveyed an average of 3.5 km (SD = 2.8 km, range: 0.21-10.0 km, N = 32) in sample grids before first detecting lynx. Our track encounter rate averaged 3.5 lynx detections/10 km of survey distance.

The probability of detecting lynx (p) in a given grid cell after 10 km of search effort, as estimated by bootstrap procedure, was 0.98 (SE = 0.01, 95% CI = 0.95-0.99) in the Purcell Mountains and 0.91 (SE = 0.03, 95% CI = 0.80-0.96) in the Seeley Lake area. Given these detection probabilities, the equation by McArdle (1990) and Reed (1996) estimated a 95% probability of detecting one or more lynx after 0.76 (range: 0.65-1.0) survey replicates throughout lynx habitat in the Purcell Mountains, whereas 1.24 (range: 0.93-1.86) survey replicates were necessary in the Seeley Lake area. If the desired probability of detection (α) is 99%, then 1.18 (range: 1.0-1.54) survey replicates would be needed in the Purcell Mountains and 1.91 (range:1.43-2.86) survey replicates would be needed in the Seeley Lake area.

The cumulative detection probability per kilometer, as estimated by the distances between track detections, reached an asymptote at 7 km, after which track detections > 7 km apart occurred too infrequently to model. The empirical distribution of distances between lynx track detections in a grid cell was skewed, with the majority of detections occurring < 1 km apart (Fig. 4). Cumulative detection probability increased from 0.63 and 0.79 after 1 km of survey distance in the Seeley Lake and Purcell Mountain study areas, respectively, to an asymptotic high of 0.97 and 1.0 after approximately 7 km (Table 1).

The sensitivity analysis showed that the probability of detecting lynx (p) in a given grid cell was insensitive to track misidentification. In the Purcell Mountains study area, the probability of detection was approximately 0.94 for changes in track misidentification from 2 to 25%, and decreased to 0.91 for a 50% reduction in track detections (Table 2). In the Seeley Lake area, which had fewer average lynx detections per grid cell than the Purcell Mountains, the probability of detection decreased to approximately 0.80 for all levels of simulated track misidentification, from 2 to 50% (see Table 2).

Discussion

In the northern Rocky Mountains of Montana, the probability of detecting lynx in habitats that supported resident populations was between 0.80 and 0.99 for surveys conducted on an 8 × 8 km grid with 10 km of search effort per cell. Thus, the probability of detection would be sufficient to infer the presence or absence of lynx in the northern Rockies after two survey replicates with 95% probability. We speculate that the high probability of detecting lynx using snow-track surveys is due to the species' high mobility, which results in spatially well-distributed tracks across home ranges (Aubry et al. 2000, Mowat et al. 2000).

It is likely that the probability of detection will vary with lynx density (Kery 2002, Gu & Swihart 2004). Although formal density estimates are not available for our study area, we instrumented a total of 36 lynx across our two study areas from 2006 to 2007 in concurrent lynx research studies (Squires et al. 2008, 2010). From these data, we estimated an average minimum population density of 2.8 lynx/100 km² within our survey areas (J. Squires, unpubl. data). This estimate is absolute density, since the survey area was composed entirely of forested lands that included both preferred and non-preferred habitats due to natural and anthropogenic fragmentation. Thus, by surveying preferentially in dense forest, we may experience higher lynx track encounter rates compared to survey routes placed randomly within grid cells. Lynx density in our survey areas was low compared to northern populations in Canada that varied from 2.3 to 3 lynx/100 km² to 17-30 lynx/100 km² (Poole 1994, O'Donoghue et al. 1997). Few density estimates are available for southern lynx populations; however, the estimated density of lynx in Montana was similar to the 2.3 lynx/100 km² in Washington (Koehler 1990), but lower than the 9.2-13 lynx/100 km² documented in Maine (Vashon et al. 2008). The probability of lynx detection that we estimate is applicable to lynx populations of similar density and may underestimate or overestimate the detection probability at higher or lower densities, respectively.

We restricted our survey to grids that overlapped resident, occupied home ranges that included reproductive females. Therefore, the detection probabilities that we report are applicable to managers attempting to locate resident populations of lynx at similar or greater population densities to those found in Montana. However, we believe that snow-track surveys are useful to detect dispersing individuals occupying atypical habitats, such as the verified lynx records in North and South Dakota (McKelvey et al. 2000), with potentially lower detection probabilities. For example, we conducted snow-track surveys in mountain ranges of Wyoming (2000-2002) and Colorado (2010) using the same methods as in this study and successfully detected the only individuals present in these survey areas based on telemetry and local knowledge (J. Squires, unpubl. data).

Methods exist to estimate the proportion of sites occupied by a species based on repeated surveys within an occupancy modeling framework (Mackenzie et al. 2003, Royle & Nichols 2003, Mackenzie 2005). These models use repeated surveys to infer the proportion of sites occupied by a given species based on whether or not the species is detected during each survey. Our methods differ in that we estimated the number of surveys necessary to infer absence of lynx in the survey area based on detection probabilities estimated by observed track detections. We showed high detectability of lynx using snow-tracking, with two important implications: 1) survey efforts designed to detect presence/absence of lynx broadly distributed within large jurisdictions (e.g. U.S. National Forests) can infer absence with high confidence after 1-2 visits, and 2) subsequent survey efforts designed for more detailed study of local lynx distribution and density in occupancy modeling frameworks may require relatively few visits within closed population sampling windows (Mackenzie et al. 2003, Royle & Nichols 2003).

We advise managers to search all grid squares that overlay forested areas where lynx are expected to be present to maximize the overall detection probability. If complete coverage is infeasible, then adaptive cluster sampling may provide a promising method of achieving a representative sample across very large landscapes (Thompson 1992). It is also important to conduct snow-track surveys at least 48 hours, and preferably 72 hours, after a track-obliterating snowfall (Squires et al. 2004). While we surveyed 10 km per cell, the results from the bootstrap estimation on the distance between tracks indicate that the probability of detection asymptotes after approximately 7 km. Thus, additional survey distance after 7-8 km may not greatly increase the probability of detecting lynx in a given area. We conducted surveys on roads and trails because of deep snow and difficult mountainous terrain. Concerns that this approach may bias lynx detections could be addressed through true random sampling (on skis or snowshoes), but at greatly increased cost and difficulty. In a previous study, we found no evidence that lynx avoided forest roads in the northern Rocky Mountains, based on 577 km of winter backtracking (Squires et al. 2010). Thus, we believe the increased efficiency warranted conducting track surveys from snowmobiles.

We recognize that genetic verification is necessary for snow-track based surveys to meet evidentiary standards for species identification (McKelvey et al. 2006, Ulizio et al. 2006, McKelvey et al. 2008). We did not collect genetic samples for track confirmation in this study because a concurrent radio-collaring study allowed us to confirm that lynx were present. Also, the technicians that conducted our surveys had multiple years of experience in trapping and snow-tracking lynx and were confident in their track identifications. However, correct track identification cannot be assumed when using snow-tracking methods to survey and monitor as part of a management protocol. This is particularly true for areas where lynx presence is unknown and therefore personnel have little direct experience. DNA verification can conclusively confirm or disprove species identification from questionable or low quality tracks (McKelvey et al. 2006). McKelvey et al. (2006) estimated a 0.55 probability of obtaining DNA verification of a lynx track after backtracking for 1 km. Given this probability, and our estimated probability of initially detecting a lynx track in 1 km of survey effort (0.55-0.83), the probability of obtaining verifiable DNA evidence after surveying 1 km of survey route and backtracking a lynx track detected on that route for 1 km is 0.30-0.46. This probability can be increased to 0.86-0.91 if a lynx track is backtracked 3 km (probability of obtaining DNA verification = 0.91) with a survey distance per cell of 7 km. After genetic verification of lynx presence in a given management area, the survey could be expedited by collecting DNA from a subsample of subsequent track detections. We recognize that DNA collection and verification may be infeasible in some situations due to cost and other logistical constraints. Although we recommend collecting genetic samples from tracks, we found that detection probabilities were robust to track misidentification in surveys used to detect resident populations. However, when track surveys are used to document new and unverified lynx populations, we stress that collecting genetic samples is a necessary component to track-based surveys.

Management implications

Our results indicate that snow-track surveys are effective for documenting the presence or absence of lynx within broad landscapes. We suggest conducting at least two replicates of the survey to infer presence or absence of a resident lynx population with 95-99% confidence. We also recommend a survey effort of at least 7 km per 8 × 8 km survey cell to maximize detection of resident lynx at low densities, and a back-tracking distance of at least 3 km per track to obtain adequate genetic evidence for species confirmation based on verifiable DNA evidence (McKelvey et al. 2006).