An ongoing challenge in wildlife conservation is to assess how wildlife species interact with human infrastructure and how this interaction is related to changing environmental conditions and demography. Migratory tundra caribou and wild reindeer, Rangifer tarandus groenlandicus and R. t. tarandus, respectively, respond to human activities over different scales of distances, which are typically several kilometers (Vistnes and Nelleman 2008). The distance from infrastructure where the effect of disturbance on habitat selection diminishes to background levels is termed the ‘zone of influence’ (Nellemann and Cameron 1996, Boulanger et al. 2012). Aboriginal elders have also commented on caribou avoidance of mines based on their expert knowledge (Tlicho Research and Training Institute 2013). Monitoring the zone of influence (ZOI) helps verify environmental assessment impact predictions, test the effectiveness of mitigation measures, and is the basis of an approach to quantify the energetic cost of caribou avoidance of mine vicinity. Estimates of ZOI for migratory tundra caribou differ among studies (Plante et al. 2018) which is, however, not surprising given all the factors which affect caribou habitat selection such as type of development and insect harassment (Valente et al. 2020).

Zone of influence estimation is useful for barren-ground caribou and other species given that it involves estimating a ‘point of change' or threshold relative to stressors which can be anthropogenic or from other sources (Groffman et al. 2006, Ficetola and Denoel 2009, Benítez-López et al. 2010, Torres et al. 2016). One issue with estimation of ZOI is that studies use different sampling and estimation methods, some of which involve subjective interpretation of patterns in habitat selection or abundance relative to distance from infrastructure. For example, approaches that use point estimates of habitat selection at binned intervals from mines have to visually assess whether variation in habitat selection in the proximity of the mine is substantially reduced compared to background natural variation at further distance from the mine. Ideally, estimating a ‘point of change’ – the ZOI distance should be tested for statistically significance and be repeatable (statistically precise). This is why Boulanger et al. (2012) developed a segmented/piecewise regression method to estimate exact threshold (ZOI) distances where change in habitat selection is not discern-able from natural background variation in habitat selection relative to the mine. Using this approach allows statistical tests of ZOI estimates as well as confidence limits on the ZOI and the magnitude of ZOI. Boulanger et al. (2012) used this approach to estimate a 14 km ZOI (CI = 12.0–15.5 km) for 2003–2008 around the operational Ekati and Diavik mines but the complexity of estimation procedures and lower sample sizes of collared caribou available limited the inference on ZOI to multi-year periods of mine operation.

A current question is how ZOI varies with changes in environmental conditions, herd size and mine activity/development which likely occurs on an annual or sub-annual basis. This prompted our re-analyses of collar data from the Bathurst barren-ground caribou herd using updated methods for the original data (1998–2008) and nine additional years (2009–2017) to assess relative annual trends in ZOI, as well as factors influencing ZOI. The Bathurst herd has declined from an estimated herd size of about 350 000 in 1996 to about 8200 caribou in 2018 (Boulanger et al. 2011, Adamczewski et al. 2019) with corresponding reduction in range size and likely scale of habitat selection relative to mines. Between 2010 and 2018, climatic extremes in the Bathurst range have included droughts and higher levels of insect harassment raising the question whether those factors influenced the response of caribou to mine development.

The objective of our research was to estimate yearly ZOI using both aerial survey and collar data to estimate yearly trends in ZOI distance and magnitude for the Bathurst herd and explore climatic factors affecting trends in ZOI. We also analyzed the role of scale in habitat selection as overall herd distribution relative to the mines changed during our study period. To allow annual estimates, we adapted segmented regression methods (Muggeo 2003, 2008) in the R statistical package (< www.r-project.org>) to estimate ZOI. We developed a model-based framework for ZOI estimation to allow statistical assessment of the similarity of yearly ZOI estimates therefore eliminating the need to subjectively pool yearly data sets. We compared the newer estimation methods with our earlier multi-year estimates (Boulanger et al. 2012) and estimates of ZOI using different approaches (Plante et al. 2018) than segmented regression. Finally, we compared the overall utility of ZOI estimates using collar location data versus aerial survey data. We provide the base R code for our analysis and note that our approach can be applied to estimation of ZOI or related thresholds for any wildlife species.

Material and methods

Analyses conducted and data sets considered

The collar dataset from the Bathurst herd for 2009 to 2017 had increased sample sizes of collared caribou and shorter fix intervals compared to earlier data sets (1996–2008, Boulanger et al. 2012). The Supplementary information gives details on capture and collaring caribou. For this analysis, we first developed base resource selection function models to explain natural variation in habitat selection. We considered the scale of habitat selection relative to mines and how this was influenced by changes in herd size and subsequent size of the summer range. We then used this base model to test for yearly zones of mine influence using segmented (also called piecewise) regression. In addition, we re-analyzed the earlier collar (1996–2008) and aerial survey data sets (1998–2008) from Boulanger et al. (2012) using the segmented R package and using the same habitat selection models developed in Boulanger et al. (2012).

Figure 1.

Locations of mine areas on Bathurst caribou range considered in the analysis for the Bathurst herd in Northwest Territories, Canada. Aerial transects surveyed in 2009 and 2012, and the extent of summer/fall range as delineated by kernel 95% utilization distribution (Worton 1989) of collar points are shown for 2009 and 2017.

Base model for collar analysis

We considered daily collar location data from 15 July to 30 November, when Bathurst caribou were most likely to interact with the mine sites. We defined habitat availability for the collar dataset based on each caribou location and estimated movement rates. We determined the proportion of habitat classes in a 1 km buffer radius around collar locations, then compared each buffered point with the buffered area around six random points that were within a circle around the previous location of the collared caribou. The size of circle was the ‘availability radius’ defined by the 95th percentile of the distance moved for caribou for the interval between successive point locations (Arthur et al. 1996, Johnson et al. 2005, 2006). Locations from collared caribou that could encounter the mine sites (as indicated by including the mine in the availability radius) at least once in a given year were included in the analysis.

We first built a base resource selection function (RSF) model (Manly et al. 1993) to account for habitat selection, then added ZOI terms to test for effects of the mine. To describe habitat, first we used vegetation classes and landform features from the Land Cover Map of Northern Canada (Olthof et al. 2009), and Earth Observation for Sustainable Development of Forests (EOSD; < www.geobase.ca/geobase/en/data/landcover/index.html>) land cover classification. Esker coverage was extracted from 1:250 000 scale National Topographic Data Base maps (Natural Resources Canada; < http://geogratis.cgdi.gc.ca/geogratis/en/index.html>). We used 12 habitat classes pooled between the NLC, EOSD and eskers coverage (Table 1). We used normalized difference vegetation index (NDVI) imagery to track plant phenology and productivity within the study area (James and Kalluri 1994, Didan et al. 2015). NDVI values were generated for 10–15 day intervals for each year averaged over the entire study area to account for annual differences in seasonality (Latifovic et al. 2005) as detailed in the Supplementary information. We also considered a climate drought index based on NASA MERRA (Modern era retrospective analysis for research and applications) remote sensing data (Russell et al. 2013).

Table 1.

Covariates used in development of base resource selection function (RSF) models for the analysis of caribou collar data for the Bathurst herd, Northwest Territories.

To assess whether reduction in size and location of seasonal ranges since 2009 affected habitat selection at the edge of seasonal ranges (see next section), seasonal ranges were estimated for the Bathurst herd. Kernel 95% home ranges were estimated for each year (15 July–30 November) using the Adehabitat package in R (Calenge 2015). Smoothing parameters estimates were based on the variances of the x and y coordinates under the assumption of a bivariate normal utilization distribution. Contours from the kernel home range were overlaid on used and random locations to index locations relative to the edge of seasonal range boundaries. A binary InRange covariate was then assigned to locations that fell outside the 95% contour to model reduced habitat selection of caribou at edges of the seasonal range, which was likely due to collective behavior of the caribou rather than habitat selection.

To build a base habitat model for 2009–2017, we applied conditional logistic regression analysis tests to determine the statistical significance of individual habitat predictor variables (Hosmer and Lemeshow 2000). Yearly data were pooled for this analysis. The type of logistic regression model varied for collar and aerial survey data, however the general form of the model was:

The quadratic term (habitat variable²) was tested for when stronger associations with habitat values were likely to occur in the midpoint of the habitat variable value, as opposed to a linear relationship. We used the interactions among habitat variables and NDVI to test for the effect of phenology on habitat selection and to account for yearly differences in seasonality. We also tested for seasonal habitat selection due to insect harassment and other factors not accounted for by NDVI, by modeling the interaction of habitat type and season. Significant variables from statistical robust Z-tests were then added into a multivariate model in the same order as the univariate model (i.e. linear habitat variable, then habitat variable × NDVI, etc.).

The fit of individual terms was evaluated by robust Z-tests from the generalized estimating equation (GEE) conditional logistic regression model in R coxreg. which groups each pairing of used and available locations into strata (Prima et al. 2017). These strata centered each comparison on the habitat available to the caribou at the time when the location was taken. Individual caribou and year combinations were then defined as a cluster in the analysis. This approach avoided issues with pseudo-replication caused by pooling telemetry data from different individual caribou. All estimates from the logistic model were expressed as odds ratios given that absolute probability of presence cannot be estimated using used/availability analyses (Boyce et al. 2002). We examined the fit of habitat selection predictions for each year using k-fold cross validation (Boyce et al. 2002, Johnson et al. 2006). For the cross-validation analysis, we randomly subdivided the data into training and testing datasets based on Huberty's rule of thumb (Huberty 1994). The goodness of fit of a model developed with the training data set was then compared with the testing dataset. We estimated the Spearman correlation (Zar 1996) of successive resource selection function (RSF) score bins with the frequency of used locations in each bin (adjusted for availability area of each bin). If the model fit the data then the RSF bin score and area-adjusted frequencies should be positively correlated (Boyce et al. 2002).

Assessment of scale of habitat selection relative to mine sites and changes in herd distribution

The estimation of ZOI requires detecting a spatial gradient in habitat selection relative to the mine area with the assumption that the gradient is mainly determined by attraction or aversion to the mine area. The Bathurst herd range was reduced substantially from 1996 to 2017 with the relative position of the Ekati and Diavik mines shifting from the center to the southern part of the summer, fall and early winter range (Fig. 2). At the edge of the range, spatial gradients in habitat selection due herd-level rather than individual finer-scale selection likely occurred, therefore creating an additional gradient in habitat selection. The different scale gradients complicated segmented regression which is ideally suited to detecting single gradients of a pre-defined shape. We therefore adopted a multi-phase approach to identify gradients in the data set and parameterize the ZOI model to test for gradients closest to mine areas. The first phase of the ZOI analyses was determining larger-scale spatial gradients of habitat selection relative to mine areas using the base habitat model with generalized additive model (Hastie and Tibshirani 1990) spline terms to model the gradient in habitat selection relative to distance from mine. The range of distances from mine that were not affected by larger scale selection was assessed as indicated by minimal slope in the fitted lines. This range of distances (i.e. 0–50 km from the mines) with minimal large-scale habitat selection (termed the iteration zone) was then used for ZOI estimation. Habitat selection at distances beyond the iteration zone was modeled with a year-specific intercept term, therefore accounting for differential yearly habitat selection across distances further from the mine sites. The relative fit of models with different iteration zones was further compared as part of the model selection procedure.

Figure 2.

Spatial predictions of the base habitat model as defined by used and available locations for the Bathurst caribou herd in Northwest Territories, Canada. The seasonal range, as delineated by a 95% kernel polygon, is indicated by a brown polygon. The actual used points for each year is show in Fig. 3. Ekati and Diavik mines are indicated by red markers that are larger than the actual mine footprints.

We used two approaches to describe distance from mine areas. First we used the nearest distance to roads and then, second the nearest distance to centroids of the footprint and the perimeter of the footprint (with points occurring within a footprint being assigned a distance of 0 from the mine area) based on the yearly footprint size (Supplementary information). Using centroids allows modeling of the gradation of habitat selection including areas within the footprint, whereas using distance from the footprint perimeter assigns these areas an intercept value. The relative fit of both methods was tested to assess which method best described habitat selection near the mine footprint.

Segmented estimation of ZOI

The segmented package (Muggeo 2003, 2008) was used with the R coxreg conditional regression object (< www.r-project.org>) to estimate ZOI for each year of the analysis. The segmented regression model was constrained to emulate the underlying ZOI curve with a linear increasing slope from 0 distance from mine up to an asymptote at the estimated ZOI, after which selection should be constant as denoted by a horizontal line (Boulanger et al. 2012). We estimated two parameters: the ZOI which was the distance where selection asymptotes to background levels, and the slope of the curve up to the ZOI (β_zoi). The overall magnitude of the ZOI was the height of the asymptote which can be estimated as the product of ZOI and β_zoi. Significance of a ZOI was evaluated by whether its confidence limit overlaps 0 and β_zoi was tested using a robust statistical test from the underlying conditional logistic regression model. Model support was evaluated using Quasi information criterion (QIC) methods (Burnham and Anderson 1998, Koper and Manseau 2009) which are similar to Information theoretic methods but adapted for the GEE coxreg. Quasi information criterion methods based on the GEE model were therefore used to compare relative model fit and the model with the lowest QIC score was then used for ZOI estimates. The suite of models included the base model, year-specific ZOI estimates, models with different iteration zones with distances defined by footprints or centroids, and a single pooled-year estimates of ZOI.

We also estimated selection at successive 3 km distances (bins) from the mines as a secondary estimate of ZOI, as used in recent studies (Plante et al. 2018, Johnson et al. 2020). This analysis also tested the assumptions of the segmented regression approach that there are no directional trends in selection beyond the ZOI distance, and that the relationship between selection and distance from mine is linear prior to the ZOI asymptote. The 3 km interval was chosen based on restrictions in yearly sample sizes of caribou locations in proximity to the mines.

Obtaining year-specific ZOI estimates involved the building of custom design matrices in program R to create dummy variables that instructed the regression model to use specific portions of the data for yearly estimates. We also estimated ZOI for earlier collar data sets (1996–2008) using the base habitat model developed for 1996–2008 by Boulanger et al. (2012) with the segmented package to test for ZOI. This model had been developed using similar methods including goodness of fit tests. Use of this base model made it possible to directly compare estimates from the segmented package with previous estimates. A description of segmented R code used in this analysis is given in the Supplementary information.

Results

Collar analysis

Summary of collar data from 2009 to 2017

The number of collared caribou that encountered mines, as delineated by at least one yearly location that was within the 95th percentile movement radius (that defined available habitat), varied between 9 and 13 from 2009 to 2015 and then increased to 24–42 in 2015–2017 (x = 17.0, SD = 12.4) with 153 of 230 collared caribou encountering the mine (Supplementary information). Most encounters of collared caribou were during the summer (x = 16.7) with lower numbers of encounters in the fall/rut x = 9.2) and winter (x = 5.0) seasons until 2017 when there was a marked increase in winter encounters (Supplementary information).

Base RSF habitat model

Seasonal and sex-specific movement rates for males and females indicated that most caribou moved an average of approximately 10 km a day for the annual time period considered in the analysis; the 95th percentile of distances was approximately 30 km per day for both males and females. The actual range depended upon season and year (Supplementary information).

The base habitat model included main habitat terms in linear and quadratic form, many of which interacted with season (Supplementary information). In addition, selection of the sedge habitat category was influenced by phenology as indexed by NDVI. Annual drought indices influenced the selection of water bodies and moss/lichen habitats with both habitat types being more likely to be selected during drought years. The general relationship for water was quadratic with higher selection for areas of low to mean proportions of water with increasingly negative selection as proportion water increased (which would reflect random locations within lakes).

Summer to early winter range of the Bathurst caribou herd varied yearly with a sharp reduction in size of ranges occurring from 2015 to 2017 (Fig. 2). During this time, seasonal ranges were also less differentiated spatially. Ekati and Diavik were roughly in the center of annual ranges from 2000 to 2012 and then on the southern extent from 2013 to 2017. Predictions from the base habitat model do not contain ZOI terms and so do not illustrate response of caribou to the mines. The general pattern reveals selection of similar habitats around lake areas each year with increasing concentration of selected areas in later years as range size was reduced (Fig. 2, Supplementary information).

Table 2.

Model selection for ZOI analysis of caribou collar data for Ekati and Diavik mines, 2009–2017 for the Bathurst herd in Northwest Territories, Canada. The iteration zone refers to the range of distances from mine considered in the ZOI analysis. Parameters for the base habitat model is listed in the Supplementary information. Quasi information criterion (QIC), difference between most supported and current model (ΔQIC) and the number of parameters (K) are given.

Cross validation of model fit suggested adequate fit with mean correlation of predicted and observed data of 0.98, which was significant in all re-samplings of the data set. Yearly fit was adequate as assessed by agreement of observed and expected frequencies with slight reductions of fit in 2011, 2013 and 2014, the latter potentially due to high drought conditions (Supplementary information).

ZOI of Ekati–Diavik

A generalized additive model analysis was first used to assess large scale selection gradients followed by a ZOI analysis. This analysis used the base habitat model with distance from mine modeled using spline terms. Large-scale selection gradients were relatively stable at distances of 50 km or less from the mine and therefore locations within this distance (and distances of 30 and 40 km) were used as iteration zones in the segmented analysis (Supplementary information).

A model with yearly ZOI estimates with distance from perimeter of mine footprint was most supported (model 1: Table 2). An iteration zone of 30 km was most supported indicating that this was the most appropriate scale to consider for the ZOI curve estimation. Models that used distance from centroid were less supported (model 2, 5, 6 and 8). In addition, models that assumed a constant ZOI across all years (model 7 and 8), or the base habitat model without ZOI terms (model 9), were not supported.

Overall, the most statistically significant estimates of ZOI, magnitude of ZOI (as estimated by the log odds ratio of the asymptote of the ZOI curve), and significance tests for β_zoi were for 4 of 9 years (2011, 2012, 2013 and 2015) where both the ZOI and odds ratio estimates were significant (Table 3). Of significant estimates the average ZOI was 9.5 km (min=6.2, max=12.8, SD=3.3, n=5). In 3 years (2010, 2014 and 2016) the ZOI estimate was imprecise (as indicated by a CV > 30%) but significant (as indicated by confidence intervals not overlapping 0), but the β_zoi was not significant suggesting a weaker effect of the mine. Estimates of ZOI using footprints (model 2) were 28% lower than those from centroids for 2011, 2012, 2013 and 2015.

Table 3.

Yearly estimates of ZOI derived from caribou collar data from model 1 (Table 2) for the Bathurst herd in Northwest Territories, Canada. Non-significant estimates are shaded. Significance of ZOI estimates is determined by confidence limits overlap of 0. Significance of odds ratio is determined by significances tests for β_zoi. A negative ZOI term implies attraction to mine areas.

One way to conceptualize these estimates is to view the raw predictions from the segmented model with the fitted ZOI curves, plotted as a function of distance from the mines and roads (Fig. 3). In general, used locations occur above the segmented curve because the model estimates higher selection for used locations. The spread of points is due to differential habitat selection. For example, random points that fall in the middle of lakes would have the lowest habitat selection scores. If there is no ZOI the density of used points should be similar to random points across the range of distances from the mines. Zones of influence occurs when there are higher densities of random points in the proximity of the mines compared to used points. This pattern primarily occurred in 2011, 2012, 2013 and 2015 which resulted in the segmented model estimating a significant ZOI slope and ZOI estimate (solid vertical line in Fig. 3) where the difference in densities occurred. Alternatively, there were higher densities of used points in the proximity of the mines than random points in 2014 and 2016, resulting in a negative estimates of ZOI. In 2009, 2010, 2014, 2016 and 2017 the slope of the ZOI was not significant as indicated by minimal change between habitat selection at distance from mine = 0 and the estimated ZOI.

Fit of the segmented model was cross validated by comparing estimated selection at 3 km intervals from the mine as predicted by the base habitat model (Fig. 3). Selection estimates at 3 km intervals (green line in Fig. 3) were unconstrained and therefore allowed a test of variation in habitat selection prior to and after the estimated ZOI. In general, mean predictions of the binned model were included within the confidence limits of the segmented model. In some years (2012, 2013 and 2015) the binned model predicted lower selection relative to the mines than the segmented model. This was most likely due to low sample sizes of used points near the mine and subsequent low precision of estimates. In most years, there were no directional trends in selection indicated by the binned model after the estimated ZOI. In 2011 and 2014 there was random variation in selection, however, the variation was contained within the confidence limits of the segmented model. Firm estimation of ZOI from the binned models was problematic given the 3 km intervals and variation among mean scores at each interval.

Spatial predictions of the segmented model illustrate variation in distribution of caribou relative to the mine site at the scale displayed (approximately 50 km on either side of the mine area). In general, used points (collar locations) occurred in the proximity of areas predicted to have higher selection (Fig. 4). Areas of lower selection in the proximity of the mines were most pronounced in 2015 but also evident in 2012. In 2017, caribou were mainly concentrated to the north of mine areas (Fig. 2) and therefore potential ZOIs were possibly confounded with the edge of the seasonal range.

Collar data from 1996 to 2008 (Boulanger et al. 2012) were re-analyzed using the segmented package (Supplementary information). Low sample sizes prevented yearly estimates of ZOI. Instead, estimates were derived by mine phase with non-significant estimates of ZOI for the Ekati construction (1996–2009) and Ekati operation/Diavik construction phase (2000–2002). For the Ekati and Diavik operation phase (2003–2008), a ZOI estimate of 15.5 km (CI = 6.6–24.3) was estimated, however, the corresponding odds ratio was not significant (Supplementary information). The ZOI distance estimate was similar to that obtained by Boulanger et al. (2012; 11 km, CI = 1–17 km) for the 2003–2008 period.

Overall trends in estimates from collar and aerial survey data

ZOI estimates from collars suggested high yearly variation (Fig. 5) and similar to the aerial survey ZOI, the ZOI estimates were not significant prior to 2003 when mines were not in full operation. Statistically, the ZOI estimates were stronger in 2011, 2013 and 2015 where both ZOI and β_zoi terms were significant. ZOI estimates from collars in years previous to 2009 were pooled by mine phase (Supplementary information) and are displayed by the midpoint of each period. There were two years of overlap (2009 and 2012) of the aerial survey-based estimates with the more substantive (2009–2017) collar data. In 2009, both ZOI's were not significant and in 2012, the estimates of ZOI were 11.4 and 7.1 km from aerial survey and collar data, respectively. Two years of collar data (2014 and 2016) had negative ZOIs but the β_zoi terms for these years were not significant and therefore a reliable estimate of ZOI was not possible.

Figure 3.

Predicted selection (log odds ratio) for segmented model (blue line) with confidence limits (hashed lines) as a function of distance from Ekati and Diavik (footprints and roads) by year for the Bathurst caribou herd in Northwest Territories, Canada. Note the different scales on the y-axis. Estimated zones of influence are shown as vertical lines with shaded areas indicating confidence limits. ZOI estimates in which the slope of ZOI was not significant (at α = 0.05) are shown as dashed line. Raw predictions from a model that estimated selection at 3 km intervals are shown with random (grey dots) and used points (red dots) delineated. A green line summarizes the mean predictions from the midpoint of each 3 km interval for the binned selection model.

A weighted regression estimate of overall trend in ZOI from 2003 to 2017 was –0.53 km per year (CI = 0.31–1.40, t = –1.36, p = 0.20) suggesting a minimally significant negative trend (at α = 0.2), however, drought index as an additive term to year was significant (β = –0.39, CI = –0.693 to 0.09, t = –2.5, p = 0.025). Odds ratios did not exhibit directional yearly trends (t = 0.31, p = 0.75) but also showed a negative association with drought index (β = –0.06, CI = –0.008 to –1.11, t = –2.5, p = 0.026; Fig. 6). When added to the trend and drought models, data type (collars or aerial surveys) was not a significant predictor of ZOI (β_aerial = –1.42, t = –0.47, p = 0.65) or odds ratio (β_aerial = –0.92, t = –0.90, p = 0.38). Predictions of ZOI and log (odds ratio) from the regression models at mean levels of drought were 7.2 km (CI = 3.8–10.5) and 1.22 (CI = 0.48–1.96), respectively, with ZOI and odds ratio decreasing as drought index levels increased.

Figure 4.

Spatial predictions of segmented model with used locations indicated in blue for the Bathurst caribou herd in Northwest Territories, Canada. Red areas are areas of avoidance with green areas being selected. A figure with all years in the analysis is given in the Supplementary information.

Results of weighted regression analysis of ZOI effect size (area under the ZOI curve) suggested significant yearly trends (β = –0.73, CI = –0.35 to 1.10, t = –4.24, p = 0.001) as well as an additive effect of drought index (β = –0.31, CI = –0.13 to 0.49, t = –3.71, p = 0.003). When added to the trend and drought models, data type was not a significant predictor of effect size (β_aerial = –9.58, t = –1.79, p = 0.10). These results suggest an overall reduction in ZOI effect size which was negatively influenced by drought level. In years of higher drought, effect size was reduced. A plot of yearly predictions from the year + drought model suggests negligible ZOI estimates when drought level was higher compared to years with lower levels of drought (Fig. 6). The model did not predict higher ZOI levels at lower levels of drought suggesting that other factors likely influenced ZOI levels.

Discussion

Analysis of caribou collar and aerial survey data from 2003 to 2017 showed that the ZOI extent averaged 7.2 km (CI=3.8–10.5) when standardized for mean levels of drought observed during this time period. The ZOI effect size significantly declined during 2003–2017 through yearly trends and an additive negative effect of drought. The lack of support of models with a single constant ZOI illustrates the need for estimation of ZOI across yearly scales rather than pooling all the data into a single period. The high degree of annual variation is to be expected (Fig. 5) given that the ZOI is an interplay between habitat selection, including annual variation in forage quality, and the perceived level of disturbances (vehicles, aircraft, blasting, etc.). Additionally, habitat selection was affected as the Bathurst herd sharply declined and its summer–early winter ranges contracted and shifted relative to the mine site (Fig. 2) and yearly drought conditions further affected habitat selection (Fig. 6). All of these factors influenced the magnitude of ZOI (Fig. 5, 6). Given the dynamic nature of ZOI, monitoring the trend in ZOI will require also monitoring the covariates that influence caribou distribution relative to mines.

Table 4.

Estimate of ZOI from aerial surveys for the Ekati and Diavik mines, 1998–2012 for the Bathurst caribou herd in Northwest Territories, Canada. The magnitude of ZOI (log odds ratio), slope of ZOI curve and goodness of fit of each yearly model is displayed. Details on the aerial survey analysis are given in the Supplementary information. Non-significant parameters are shaded in grey.

We used the segmented regression approach as it is a model-based framework to evaluate anthropogenic impacts with unbiased detection of thresholds (Ficetola and Denoel 2009, Boulanger et al. 2012). More exactly, the segmented regression constrains the relationship between habitat selection and distance from mine to test when changes in habitat selection cannot be distinguished from background levels. By using the segmented regression, a statistical estimate of ZOI, with confidence limits, can be generated. Alternative methods, such as assessment of habitat selection at successive distance from mine intervals (Plante et al. 2018, Johnson et al. 2020) or polynomial regression (Flydal et al. 2019), test the same hypothesis but do not allow an exact statistical test of both ZOI extent/distance and the magnitude of ZOI and do not estimate standard error of ZOI which are useful in determining the statistical reliability of estimates. For example, we were able to use ZOI variance as weighting terms in trend estimation therefore accounting for varying levels of precision of estimates. Comparison of segmented regression to binned selection intervals (Fig. 3) in our study suggests relatively similar estimates, however, the resolution of the binned estimates was based on 3 km intervals which make direct comparison to the segmented estimate problematic. Use of both approaches does allow testing of assumptions of the segmented curve, namely if directional selection is occurring after the estimated ZOI, and the relative linearity of the relationship between predicted selection prior to the ZOI.

We note that our approach could be applied to other species to estimate change in abundance or other response metrics relative to infrastructure. Other studies have used binned estimates of abundance as a function of distance to human infrastructure to estimate an ‘Area of influence’ (which is similar to ZOI) for a variety of wildlife species in Europe (Torres et al. 2016). Ficetola and Denoel (2009) assessed ecological threshold analyses using simulations and concluded that piecewise or segmented regression as well as generalized additive modelling provided the most robust estimate of thresholds. The segmented R package can be used with output from most regression packages in R, thus it should be possible for researchers to easily estimate ZOI once base habitat models are developed (Supplementary information). The program R approach is more flexible than a previous method which relied on assessment of peaks in likelihood scores of regression models (Boulanger et al. 2012).

Our analysis illustrates the strengths and weaknesses of using collared caribou to estimate ZOI. A strength of collar data is that at the individual scale they provide additional inference on habitat selection relative to mines at range-wide scales. Temporally, collar data extend across seasons therefore providing an estimate of ‘average' ZOI in comparison to aerial survey data that are based more upon a ‘population-level snapshot’ of distribution. Spatial scale can also confound collar analyses. The generalized additive model approach illustrated larger scale gradients in habitat selection at the edge of seasonal range which were most prevalent in later years of the study as the range size contracted. The second gradient to consider is the scale of zone of iteration where ZOI is estimated. The ZOI analysis is basically testing if a caribou will select habitat equally in the proximity of a mine relative to other areas. The actual effect would most likely occur within the distances from mine that are within the availability buffer, which was in the range of 30 km in the analyses (Supplementary information). As distances increase from the mine it is less likely that habitat within the proximity of the mine would be available and therefore the magnitude of selection based on ZOI would be diluted. The 30 km scale is also the approximate buffer extent of aerial survey transects around the mine sites (Fig. 1). Determining the iteration zone is subjective and therefore the data-driven method we used (considering models across a range of iteration zones) provides a more objective and transparent means to determine the appropriate scale. The final scale to consider is whether distance from mine should be measured from the footprint perimeter or centroid. Our results suggest distance from footprint is more supported. One reason for this is that when footprint is used any locations within the mine are used to estimate the intercept term (distance from mine = 0) which assumes similar habitat selection within the mine footprint. In both the case of iteration zone and method of estimating distance from mine the data are determining the zone that best described potential effects of the mines (ZOI) and appropriate areas of comparison away from the ZOI. This is one of the main advantages of the model-based ZOI segmented analysis compared to more subjective methods such as binned-intervals. The bottom line in terms of use of collars is that while they provide the potential to estimate ZOI, issues with scale and sample size of caribou relative to mines can limit application unless collar sample sizes are large.

Figure 5.

Estimates of ZOI and odds ratios for Ekati and Diavik analyses for collared caribou and aerial surveys for the Bathurst herd in Northwest Territories, Canada. ZOI estimates are classified as significant (both ZOI and OR estimates were significant), partially significant (ZOI estimate was significant but slope term of ZOI was not significant suggesting a weak ZOI) and not significant (both ZOI and β_zoi were not significant). Estimates of ZOI from collars prior to 2009 (Supplementary information) are displayed by the midpoint of pooled years in the analysis The predicted linear yearly trends from weighted least squares regression (year + drought index) for ZOI and Log odds ratio are shown as a black line (with shaded confidence limits).

Aerial surveys are a more precise measure and finer scale sampling of caribou distribution in the proximity of mine areas. A difficulty with aerial surveys is obtaining adequate sample sizes of cells with caribou to estimate base model and ZOI parameters. Our analyses showed that it is possible to estimate ZOI on a yearly basis if within-year sample sizes are adequate. Analysis of the 2008 data suggest that at least six surveys (7865 km of transect with caribou present in 140 of the 7865 1-km cells) are needed to estimate the parameters of the underlying habitat model and the ZOI terms (Boulanger unpubl.). After this sample size is achieved, the main effect of increasing sample size is improved precision of the ZOI estimates. This result is similar to the ‘rule of 10's’ (Hosmer and Lemeshow 2000) which indicates that at least 100 cells with caribou detected are needed to support the 10 parameter base habitat/ZOI model. When sample size was below 140 the GEE model did not converge. There is also some stress to the caribou associated with aerial surveys although we acknowledge the capture of caribou for fitting collars is also stressful. Aerial surveys are expensive, particularly if flown by helicopter, which leads us to suggest that aerial surveys should be used intermittently to cross-validate collar-based ZOI estimates.

Figure 6.

Estimates of ZOI effect size for Ekati and Diavik analyses for collared caribou and aerial surveys as a function of drought index for 2013 (the midpoint year in the analysis). ZOI estimates are classified as significant (both ZOI and OR estimates were significant at α=0.2) or partially significant (ZOI estimate was significant but slope term of ZOI was not significant suggesting a weak ZOI). Non-significant estimates were not used in the analysis. The mean drought index is displayed as a hashed line in the figure.

An assumption of our trend analysis as well as many other ZOI studies of caribou (Flydal et al. 2019) is that ZOI estimates derived from aerial surveys will produce similar ZOI estimates as from collar studies. There were only two years that had adequate collar sample sizes and aerial surveys (2009 and 2012) and the majority of aerial survey data was for earlier periods of the analysis (1998–2008, 2009 and 2012) which precluded direct statistical testing of estimates from the two data types. We found no statistical difference between the data types when added as an additional term in the drought and trend regression models. However, there are differences in the type of data and corresponding response of the two data types. Namely, aerial survey utilizes a used/not detected design providing a response metric similar to occupancy at the scale of transect cells, whereas collars employ a used-available design providing a selection ratio at the scale of daily caribou movement. In both cases it is assumed that the distribution of caribou on transect cells or selection of individual caribou relative to available habitat represent the overall population of caribou (Manly et al. 1993) with additional assumptions needed to relate each selection metric to population density (Boyce and McDonald 1999). We note that the scale of distances from mine being considered for ZOI estimation for each data type is roughly 30 km which corresponds to the extent of aerial survey transect cells relative to mine areas and the iteration zone used to search for ZOI's in the collar analysis. Also, ZOI estimation assesses the change in habitat selection relative to mine rather than the direct estimates of the habitat selection from each data type. Given this, the ZOI estimates should be relatively robust to type of method used to estimate ZOI as long as analyses are conducted at similar scales relative to the mines. We suggest that future projects intersperse aerial surveys with collar efforts to allow a direct cross-validation of ZOI estimates from both methodologies.

Estimates of ZOI from elsewhere on migratory tundra caribou ranges vary. The Leaf River herd in northern Quebec is exposed to a large but mostly underground mine on the calving and post-calving ranges and the ZOI based on satellite collar data was measured at 19–23 km (Plante et al. 2018). On the Leaf River herd's winter range the ZOI averaged 9 km for a road with hunting compared to 2–3 km for a road without hunting (Plante et al. 2018). In Alaska, despite 35 years of exposure to the large Prudhoe Bay oilfield complex of drill sites, above-ground pipelines and roads, caribou from the Central Arctic herd showed ZOIs of 1, 2 and 5 km for mosquito season, post-calving and calving, respectively (Johnson et al. 2020). Analysis using satellite collared data from the Porcupine herd in Yukon during winter revealed ZOIs between 6 and 18.5 km for roads and communities (Johnson and Russell 2014).

Climatic variability influenced yearly habitat selection as well as ZOI, with caribou more likely to be associated with water bodies and wet moss/lichen areas in drought years (Table 4). The effect size of ZOI decreased in drought years, possibly because caribou may have been attracted toward large water bodies in the proximity of mines which through micro-climate effects may buffer the effect of drought on vegetation. In addition, insect abundance indices were higher during drought years (Russell et al. 2013) which would lead to higher aggregation and behavior oriented towards relief of insect harassment including a cooler microclimate closer to large lakes (Bergerud et al. 2008) which may have reduced any avoidance of mine areas. Plante et al. (2018) described that during the particularly hot summer in 2010, caribou from the Leaf River herd did not avoid the mine. Likewise, ZOI relative to the oilfield for the Central Arctic herd were reduced during high mosquito activity periods (Johnson et al. 2020).

Our study does not attempt to estimate the demographic impacts of ZOI or how demography influences ZOI. We note that range reduction did occur during much of the study period (Fig. 2) as related to population decline with potential depensatory demographic effects that may have influenced the response of caribou to the mine. Also, it is likely that behavior during drought years when caribou were attracted to water bodies near mines represented a tradeoff between disturbance and potential reduction of forage quality near mines. All these factors could influence demography with mines creating a potential ecological trap if caribou are forced to tradeoff habitat selection for avoidance of mine areas. Ongoing demographic analyses suggest that droughts and related insect harassment do reduce calf survival as well as pregnancy rates (Boulanger unpubl.). Therefore, it is possible that ZOI covaries with demography as well as climatic factors.

The exact mechanisms that cause caribou to avoid mines, roads and oilfields has not been clearly identified. The effects of mine activities that may affect caribou habitat selection include visual disturbance, noise and dust-fall (Boulanger et al. 2012). Dust-fall may include either the larger particles of fugitive dust or the very fine particles that are carried by the wind for 10s of km. Dust is caused by blasting, piling and sorting rock and to a large degree from road traffic, and affects vegetation including lichens (Chen et al. 2017). Accumulation of dust on vegetation will often depend on wind patterns as well as rainfall events. In addition, noise propagated by mining surface vehicles can potentially be detected by caribou at 6 km, helicopters at 15 km, blasts at 40 km and larger aircraft at 100 km distance (Leblanc et al. 2018). These types of covariates were not available for our study and we suggest that studies should plan systematic collection of mine-related covariates as part of the monitoring program. Our study also did not address the effects of seasonal winter roads to access mines which may also create response and disturbance patterns (Paton et al. 2017) as well as demographic effects due to harvest near winter roads (Boulanger et al. 2011).

Other possible factors driving habitat selection and avoidance of particular areas are learning, memory and social behavior (Foss-Grant et al. 2018, Martinez-Garcia et al. 2019). Caribou show high fidelity to calving and summer ranges and older cows likely pass on their accumulated knowledge to younger caribou. This likely includes knowledge of key water crossings and migratory routes. Avoidance of areas such as mines as a learnt behavior could also be passed from older to younger caribou. Fagan et al. (2013) reviewed recent understanding of memory and movement and emphasized how cognitive memory is a key toward efficient movement through complex landscapes. Correlated movements among individuals may result from similar responses to a habitat attribute (Calabrese et al. 2018).

Analysis of collared caribou data cannot necessarily address which of these factors – dust, noise or behavioural learning, or a combination of these – is responsible for the ZOI around human infrastructure. A more detailed understanding of what exactly caribou are avoiding would be important in designing adaptive mitigations that could reduce ZOI detected for caribou. As well as not yet capturing the role of memory, our dependence on collaring individuals has not yet quantified how movements are correlated among individuals, which is essential in describing how individuals respond.

Our analyses demonstrate that ZOI estimation reveals yearly and spatial variation in how caribou react towards mine sites. In this context, the effect size of ZOI may be the best indicator of trends in response to mines given that it accounts for both spatial displacement as well as the magnitude of habitat selection related to ZOI. We suggest that ZOI should be estimated yearly given variation in climate, demographic and mine-related covariates. A regression-based estimate of ZOI standardized for mean covariate values provides a way to assess overall ZOI while accounting for climatic or other covariate variation. The analyses in this paper demonstrate how satellite collars can be used for this objective; however, analyses need to account for spatial scale and environmental variability that may also show temporal trends. The added advantages of collars are that they provide information on demographics such as adult survival as well as other useful information for management and research. Therefore, estimation of ZOI from collared individuals can be a complementary part of an integrated system of caribou management.