Using multi-scale distribution and movement effects along a montane highway to identify optimal crossing locations for a large-bodied mammal community

Study area

Our study was conducted along Southern Trans-Provincial Highway 3 (hereafter Hwy 3) between the towns of Creston and Cranbrook, in south-eastern British Columbia, Canada (Fig. 1). The study area is located in the Purcell Mountain Range, which ranges from 620 m to 2,087 m in elevation, and is comprised of Interior Cedar Hemlock and Interior Douglas Fir Biogeoclimatic zones (Meidinger & Pojar, 1991). We chose this study area for its ecological importance as a trans-boundary priority area (Yellowstone to Yukon Conservation Initiative, 2013) that connects small populations of carnivores such as grizzly bears (Ursus arctos horribilis) and Canada lynx (Lynx canadensis) along the Canada – USA border. Hwy 3 bisects this important corridor, possibly leading to negative effects on the connectivity of this movement corridor for mammal populations. The average annual traffic volume (AADT) for this highway section was 3050 cars/day in 2007, with a seasonal (December to March) average of 2020 vehicles/day (British Columbia Ministry of Transportation and Infrastructure, 2010).

Data collection

We monitored species movement through the study area by recording tracks in the snow where animals attempted to cross Hwy 3, as well as along ten transects approaching the highway, set back from any highway right-of-way (distance from transects to highway ranged from 10 to 900 m, mean 175 m). We pre-defined our transects as survey lines marked with flagging tape, roughly parallel to the highway. Highway and transect tracks were recorded over two winter seasons, January to March 2007 and December 2007 to February 2008 (all observations recorded by RS).

Highway and transect track surveys were conducted using methods similar to Van Dyke, Brocke & Shaw (1986), and Alexander, Waters & Paquet (2005). Briefly, we conducted highway crossing attempt surveys along a 95 km length of Hwy 3, at least 12 h after the last snowfall. Each survey was conducted from a moving vehicle with a speed of approximately 10–15 km/h. When a track was observed, the investigator stopped the vehicle and conducted an on-foot inspection to identify the track. In total, we investigated tracks for 12 mammal species: coyote (Canis latrans), fox (Vulpes vulpes), wolf (Canis lupus), cougar (Puma concolor), bobcat (Lynx rufus), lynx, marten (Martes americana), wolverine (Gulo gulo), elk (Cervus canadensis), moose (Alces alces), white-tailed and mule deer (Odocoileus virginianus and Odocoileus hemionus, respectively). When we were uncertain of the identity of a track, we recorded track pattern measurements, took photos and later consulted field guides (Sheldon, 1997; Elbroch, 2003) for identification. Data at a total of 463 crossing sites were georeferenced with a handheld, Garmin eTrex Summit GPS receiver (WGS 1984, ±10–40 m). If multiple tracks were found for one species at a single crossing area, we recorded the total track count. We also recorded the success of a crossing attempt, here defined as the presence of a continuing set of tracks on the opposite side of the road. When tracks of the same species were found within 300 m of a crossing site, it was not recorded as an individual crossing attempt, but rather as a possible repeat crossing of the same individual (Alexander, Waters & Paquet, 2005). Surveys were suspended when continuous heavy snowfall covered tracks during data collection.

Transects were established off-road in suitable areas close to the highway. Suitability was contingent upon minimal disturbance from residential areas, and no barriers to observer access (i.e., lakes, steep terrain, fences or private property). Seven transects had a linear distance of 1 km, while one was 2 km (Transect 6) and one was 5.4 km in length (Transect 10). Only the first kilometer of transect 10 was surveyed during the second season of data collection, and this was classified as Transect 9 for ease of data handling. We recorded tracks of the same species according to the protocol of the crossing attempt surveys, and georeferenced a total of 308 individual track locations along the transects. We surveyed transects between 12 and 96 h after snowfall, usually starting the day following a road survey, with 5 to 7 km of transect being surveyed per day. Due to the limited number of tracks recorded for carnivores (coyote, bobcat, cougar, wolf, fox, lynx, marten, wolverine, see Results) we grouped all the above species into one category ‘carnivores’, while evaluating the remaining species of ‘ungulates’ (moose, elk, deer) separately for landscape variable preference models and predictive mapping, and both separately/combined for estimates of permeability across the highway.

Road barrier effect

We standardized the highway crossing attempt and transect survey data by the number of 12 h periods that had elapsed since the time of the last snowfall to correct for time effects (Thompson et al., 1989). For calculation of the road barrier effect, we standardized survey data for the highway and transects by kilometers surveyed: ${\displaystyle \text{Crossings per km}=\frac{\text{Total number of tracks}}{\text{Total length of surveys}}}$ We then calculated the permeability of the highway by standardizing the crossings per km of highway with the crossings per km of transect: ${\displaystyle \text{Permeability}=\frac{\text{Highway crossings per km}}{\text{Transect crossings per km}}}$ We also constructed track accumulation curves along the 95 km of highway for all four species groups to identify areas of the highway with greater crossing intensity for each mammal group.

Multi-scale landscape variables

To develop our micro-habitat assessments, we imported the collected GPS data into ArcGIS 9.3 (ESRI, 2009). The GPS points from the highway surveys and the transect surveys were set on top of a georeferenced (WGS 1984, UTM Zone 11N) orthophotograph layer from 2004, with a spatial resolution of 1 m, provided by a Web Map Service (WMS) of GeoBC (http://www.geobc.gov.bc.ca). For each GPS point, we created a circular buffer of 200 m to represent the perceptual area of the animal directly influenced by the surrounding landscape predictors (e.g., Lingle & Wilson, 2001), which we define as ‘perceptual area polygon’. For each buffer area, we digitalized polygons for predefined landscape predictors and georeferenced them using the orthophotograph layer and Google Earth, as the latter provided more recent images of the research area. We used the following landscape predictors, adapted from Dickson, Jenness & Beier (2005): forested (forest + woodland), shrub, herbaceous (grassland + agriculture), riparian, water, non-vegetated (gravel, rock + dirt), highway ( + shoulder), road/path, railroad, residential, developed, disturbed and wetland (Table 1). We then calculated the percentage of each buffer area overlapped by each landscape predictor.

Because large mammals might respond to both fine and coarse scale habitat features (e.g., Mayor et al., 2007), we developed a series of variables describing macro-habitat landscape features at three spatial scales: 200 m, 500 m, and 1 km. For modeling species abundance along the highway and transects, we chose candidate predictor variables based on their ability to predict species abundance at site and landscape levels in similar studies (e.g., Malo, Suarez & Diez, 2004; Guisan & Thuiller, 2005). All remotely sensed predictors (Table 1) were derived from the following sources: Terrain Resource Information Management (TRIM, Province of BC 1992) and Earth Observation for Sustainable Development Landcover (EOSD LC 2000, Wulder et al., 2008). Our dataset comprised 12 predictor variables from the perceptual area polygons and 11 from remote sensing on 3 scales (200 m, 500 m, 1 km; Table 1), derived at each of 463 highway locations and 308 transect locations. All remote sensing predictors were created using Geospatial Modelling Environment (Beyer, 2012) in conjunction with ArcGIS 10 (ESRI, 2010) and R v. 2.15.2 (R Development Core Team, 2012). Due to their widely varying scales, all predictors were standardized to mean = 0, sd = 1 to ensure that their importance was not driven by measurement scale (White & Burnham, 1999).

Landscape variable preference models

Since predetermining the appropriate data distribution for our count data from ecological knowledge alone was not possible, we modeled abundance incorporating both a Poisson distribution (P) and negative binomial (NB) distribution to account for potential overdispersion (e.g., Zeileis, Kleiber & Jackman, 2008). Because of the large proportion of zero values included in our data-set, we also applied zero-inflated models (ZIP, ZINB; Lambert, 1992), which are mixture models that combine both count data and a binomial model. To determine which of these distributions best represented our species data, we visually inspected the data and compared the log likelihood, AIC, and number of correctly predicted zeros for each distribution model fits using intercept-only models. To test for differences among distribution functions, we used likelihood ratio tests to compare the Poisson and negative binomial distributions, since the Poisson distribution is a restriction of the more general negative binomial distribution (Hilbe, 2008). We tested H₀ for no difference between the two and H₁ that the negative binomial was a better fit to the data. We tested the same hypothesis using the zero-inflated Poisson and zero-inflated negative binomial. Next, we used a Vuong test (Vuong, 1989; Greene, 1994) to evaluate whether the zero-inflated models were a statistically better fit to the data than their base model (Hilbe, 2008). The Vuong test is generally formulated as: ${\displaystyle V=\frac{\mathrm{sqrt}\left(N\right)*\mathrm{mean}\left(m\right)}{s_m*m}m=ln\left(\frac{{\mu }_1}{{\mu }_2}\right)}$ where μ₁ = predicted probability of y for the zero-inflated model, μ₂ = predicted probability of y for the base model, s_m = standard deviation of m, and N = number of observations in each model, where both must use the same observations. The test statistic V is asymptotically normal. If V > 1.96, the zero-inflated model is preferred; if V < −1.96, the base model is preferred; and if the value of V is between −1.96 and 1.96 neither model is preferred (Hilbe, 2008). To perform these tests we used the function vuong in R package pscl v.1.04.4 (Zeileis, Kleiber & Jackman, 2008; Jackman, 2012).

We compared our measures of model selection (AIC, LogLik, predicted zeros, Vuong) for all four distributions (P, NB, ZIP, ZINB) throughout each of the model building stages of this study to avoid bias in predetermining the distribution with intercept-only models. For each of the four mammal groups (deer, elk, moose, carnivores), we compared models of predicted habitat preference using each distribution (P, NB, ZIP, ZINB) and data set (Highway, Transect) across six separate spatial approaches: (i) 200 m scale, (ii) 500 m scale, (iii) 1 km scale, (iv) all 3 scales combined; (v) perceptual area polygons; (vi) all scales and perceptual area polygons combined. For spatial approaches (iv) and (vi), we created an iterative model fitting procedure that starts with an intercept only base model, individually adds predictor variables, records the results of each fitted model, and retains all top-ranked models (ΔAICc ≤ 2) at each iteration as base models for subsequent iterations as long as there is reduction in AICc (R function in R-Code S1). The fitting procedure constitutes an extension of a more restrictive routine that only included the top ranked model for each iteration in subsequent iterations (Schuster & Arcese, 2013). We opted for the iterative approach because creating all possible models for approach (iv) and (vi) would have resulted in 2⁴⁵ models each. For the remaining approaches, we created models for each possible combination of predictors. For both ZIP and ZINB we further expanded our model lists using (a) intercept only models for the zero-inflation component, while using predictors in the count component and (b) including the same predictors for both zero-inflation and count components. In post-processing we reduced the candidate set of models from each approach based on the statistical significance of all predictors, using p-values as a general and liberal criterion for retaining models. We selected a cutoff value of p = 0.15 as it serves as the default for many stepwise model selection approaches (e.g., Rawlings, Pantula & Dickey, 1998). We chose this approach to reduce the probability of including non-informative models (i.e., those stuck at local maxima for parameter estimation) in subsequent model averaging. Additionally we checked for and removed models with unrealistically high (>50) parameter estimates and/or Standard Errors, which would indicate lack of model fit. For model selection, we ranked all remaining candidate models by AICc and averaged those with ΔAICc ≤ 2 from the top ranked model (Burnham & Anderson, 2002). All analysis were conducted using R v.2.15.2 (R Development Core Team, 2012); package MuMIn 1.8.0 (Barton, 2012) was used for AICc calculations; package MASS v. 7.3-22 (Venables & Ripley, 2002) for NB models; and packages pscl and Formula v.1.1-0 (Zeileis & Croissant, 2010) for ZIP and ZINB models.

To determine the distribution and scale that provided the best relative fit to the data for each mammal group we compared and ranked the models with the lowest AICc of each approach and determined the approach resulting in the overall lowest AICc value. We further contrasted these results with the initial results from the intercept only models to determine whether initial models were sufficient to identify the error distributions that were most appropriate for a given data set or whether predictor variables had to be included first.

Predictive maps of preferred crossing sites

Using the previous model results we created predictive maps for each species/group for both the abundance of animals approaching the highway (transect models) and the abundance of animals, having reached the highway, crossing it (highway models). For each map we chose the approach with the lowest AICc values out of the 4 remote sensing-derived frameworks (200 m, 500 m, 1 km, 3 scales), as landscape level data was not available for the perceptual area polygons. For predictive polygons, we used 30 × 30 m polygons to follow the EOSD resolution. For the highway predictions we created polygons around the highway (line feature buffered 15 m on each side of highway) resulting in predictions for 7374 polygons. For transect predictions we expanded the buffer around the highway to 1 km, resulting in 213442 polygons. Next, we generated a predictor set for each polygon centroid that was identical to those used for survey points, and then estimated abundance based on our averaged models for each of the focal species groups. We combined individual mammal group abundance estimates into 10 quantiles to consolidate focal species maps into an index of site preference, multiplied those scores for each polygon and standardized them by dividing by 1000, resulting in community site preference scores between 0 (being the lowest preference) and 10 (highest preference).

Road barrier effect

Highway permeability values for the majority of groups were extremely low (where a value of 1 indicates full permeability across the highway and 0 represents no permeability), indicating that Hwy 3 likely acts as a barrier to mammal movement (Table 2). The permeability values for carnivores were only one third those of ungulates (moose, elk, deer, combined) on the investigated section of Hwy 3 (Table 2), indicating that deer, elk, and moose were much less affected by the highway in terms of movement than carnivores. Track accumulation curves for all mammal groups indicate that in all four cases there were areas of the highway where the focal group rarely or almost never crosses the highway (Fig. 2).

Landscape variable preference models

Based on intercept only model comparisons using likelihood ratio tests, Vuong test and AIC ranking, the best supported distributions for each mammal group were: ZINB (Deer-Hwy, Elk-Hwy, Deer-Trans), NB (Moose-Hwy, Elk-Trans, Carnivora-Trans), ZIP (Carnivora-Hwy), P (Moose-Trans), indicating that ZINB and NB were the most commonly supported distributions (Table 3). In 6 out of 8 cases, the modeling approach which included predictors from all three scales and the digitized polygons was selected as the top model based on AICc (Table 4). In only two cases for the transect data did other approaches result in lower AICc values: Deer (1 km scale) and Moose (500 m scale). In direct comparison between predictors derived from remotely sensed data and hand digitized data, the remotely sensed model framework resulted in lower AICc values in all eight cases. When comparing the remotely sensed data approach at the same scale as the digitized (200 m) data, digitized predictors resulted in lower AICc values in 7/8 cases. A comparison of the extended modeling results (Table 4) with the initial distribution tests (Table 3) indicates that in 4/8 cases the results from the initial tests were rejected and different distributions formed the basis of models with the lowest AICc values. For Highway and Transect data from each mammal group, model coefficients describing preferred (positive values) and avoided (negative values) habitat variables from the best supported model (above) are depicted in Tables S1 and S2. Summed values of preferred and avoided landscape variables for the entire mammal community are presented in Table 5.

Predictive maps of preferred crossing sites

To map preferred crossing sites, we used averaged model results for each mammal group based on the framework with the lowest AICc value out of the 4 remotely sensed model sets (Cells marked with an asterisk in Table 4; maps and shapefiles in Figs. S1–S4 and Data S2 respectively). The combined predictions of preferred (green) and avoided (red) crossing sites for all investigated species within the mammal community are illustrated in Fig. 3. Based on predictions generated from landscape variables (above), certain regions of the study area exhibited high preference scores from both approach (transect) and crossing (highway) models (e.g., Fig. 3 insert A), indicating that these locations likely represent areas of high priority when implementing mitigation measures for all species considered in our study. Conversely, certain regions of the study area exhibited high preference scores for one of the model sets (crossing vs. approach), but not the other (e.g., Fig. 2 insert B), indicating that these may represent less-ideal locations to implement mitigation measures such as crossing structures. Areas of unambiguous preference for particular crossing sites (i.e., those where crossing and approach preference scores overlap) differ for each mammal group considered in our study (Figs. S1–S4), indicating that mitigation strategies aimed at mammal communities may differ substantially from those aimed at a target species.