The zone of influence (ZOI) is the area in the vicinity of industrial development where avoidance by caribou Rangifer tarandus or other wildlife species is observed. Here we examine ZOI around two diamond mines in the Northwest Territories (NWT), Canada from 1998 to 2017. In this paper, we further develop segmented/piecewise regression methods to analyze collar location and aerial survey data with a focus on yearly trend in ZOI for the Bathurst caribou herd. A base habitat model was initially formulated to account for habitat selection around mines followed by estimation of ZOI distance and magnitude. Seasonal ranges of the herd contracted from 2009 to 2017 due to decline in herd size, which influenced the distribution of caribou relative to the mines as well as larger scale habitat selection. Models with year-specific estimates of ZOI were more supported than models assuming a constant ZOI across years. Significant ZOI's occurred for aerial survey and/or collar data in 9 of 15 years from 2003 to 2017 when both mines were in full operation, with ZOI distances ranging from 6.1 to 18.7 km. Non-significant ZOI's occurred from 1998 to 2002 before both mines were fully operational. Caribou were attracted to lakes in drought years which significantly influenced distribution relative to mines as well as the magnitude of the ZOI detected. The ZOI extent averaged 7.2 km (CI = 3.8–10.5) when standardized for mean levels of drought. Our analysis suggests that ZOI varies both annually and spatially because of the location of mines relative to habitat selection and seasonal range size. Therefore, exacting analysis methods that account for these sources of variation are required for robust ZOI estimates. Segmented regression methods have become available in the R statistical program that allow flexible ZOI estimation for implementation of the methods in this study for caribou or other wildlife species.
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
The study area is approximately 300 km northeast of Yellowknife, NWT, Canada (Fig. 1) and lies within tundra of the Southern Arctic ecozone (Ecological Stratification Working Group 1996). The landscape has esker complexes, boulder moraines and numerous lakes. Shrub communities of willow Salix spp., shrub birch Betula spp. and Labrador tea Ledum decumbens dominate areas with adequate soil development. Mats of lichens, mosses and low shrubs are found across exposed rocky and gravel sites.
The Bathurst herd of caribou annually migrates from wintering ranges near or below tree-line, to calving and summer ranges on the open tundra. The winter ranges have contracted and shifted north as the herd declined which has reduced the extent of pre-calving migrations. Caribou were previously most likely to be in the vicinity of the mines during the post-calving through fall seasons, however in recent years caribou were also in the vicinity of the mine in early winter (Fig. 1, Supplementary information).
The two operational diamond mines are the main Ekati and Diavik mines which are approximately 30 km apart. Both mines are open pit mines (Ekati also had three underground mines during the study period) with accommodation complexes, ore-processing buildings and airstrips. Ekati has a separate camp and open pit (Misery) which is connected by the 29 km all-weather Misery road to the main Ekati site. The Misery camp and pit are 7 km north of Diavik mine, which is restricted to a large island in Lac de Gras (Supplementary information). Because of the juxtaposition of the Ekati and Diavik operations, we modeled these mines as a combined unit. Caribou data were available for the mine sites from construction periods (Ekati: 1996–1998, Diavik: 2000–2002) through to full operational periods (2003–2017) which allowed description of caribou distribution relative to mine areas across a range of mine footprints and activities. The footprints of the Ekati and Diavik mines in 2017 were 35 km2 and 12 km2, respectively.
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).
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).
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 variable2) 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.
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.
Aerial survey analysis
Logistic regression analyses based on presence/non-detection of caribou groups in transect cells, for surveys conducted from 1998 to 2009 and 2012, were used for ZOI analysis (Supplementary information). The base habitat model developed by Boulanger et al. (2012) was used in the aerial survey analysis given that the caribou observation data set used in the analysis was the same as that in Boulanger et al. (2012), except for the addition of survey data from 2009 to 2012. This model included terms for relative occupancy (relative population size of caribou on the grid during the survey), sedge/wetlands, water, low shrub, tundra and phenology as indexed by NDVI (Table 1). We used presence/non-detected as the response rather than group size given the gregarious nature of caribou and subsequent non-independence of individuals within a group. We assumed sightability of caribou was constant across the study areas. Both of these topics are discussed in detail in Boulanger et al. (2012).
Base RSF habitat models were run using the glm (< www.r-project.org>) or geepack (Yan 2002) packages in R. Once a base model was developed, it was run through package segmented (Muggeo 2003, 2008) to estimate ZOI distance and associated magnitude (as determined by the odds ratio, OR) for each survey year. We examined whether ZOI models adequately fit the data by testing the final segmented models for goodness of fit using ROC tests (pROC package (Xavier et al. 2011). If a model adequately fit the data then the ROC area under the curve should be ≥ 0.7 (Boyce et al. 2002, Boulanger et al. 2012). Analyses were conducted mainly in R with additional plotting of maps and charts using the ggplot2 R (Wickham 2009) and sf (simple features) (Pebesma 2018) R packages and QGIS software (QGIS Foundation 2020).
Analysis of trends in the distance and magnitude of ZOI
We conducted a weighted least squares regression analysis of trends in ZOI distance and magnitude with a yearly trend term, data type (collars or aerial surveys), as well as a term for yearly drought index for time periods when significant ZOIs were detected. In addition to analysis of ZOI distance and magnitude, the effect size of ZOI, which is the product of the estimated ZOI and the log odds ratio divided by 2, was used to assess trends in ZOI. This effect size is equivalent to the area under the ZOI curve up to the estimated ZOI and provides an estimate of the overall strength of the ZOI effect. Variance for the effect size were estimated using the delta method (Thompson 1992) based on the combined variance of ZOI distance and the ZOI odds ratio terms. Estimates were weighted by the inverse of their variance in the regression analysis to account for different levels of precision of estimates as well as meeting the assumption of equal variances of response variables across the range of predictor values (Kleinbaum and Kupper 1978). All estimates were included in the analysis if the slope term for the ZOI curve (βzoi) was significant at α = 0.2 to allow a full range of effect sizes to be considered in the analysis. A Huber–White sandwich estimator was used to estimate p-values for regression parameters (Kleinbaum and Kupper 1978, Milliken and Johnson 2002).
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).
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.
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.
Aerial survey analysis
Sample size summary
Aerial surveys were carried out annually from 1998 to 2009, then again in 2012 (Table 4). The number of surveys per year varied from 17 to 9 with the most surveys in 2008. The proportion of 1 km transect cells with caribou varied from 1.6 to 10.3% (Supplementary information). Details of the aerial survey analysis are presented in Boulanger et al. (2012) and in the Supplementary information.
Year-specific estimates of ZOI
A significant ZOI was not detected until 2003 (when both mines were in full operation) with both ZOI and the odds ratio overlapping 0 for years prior to 2003 (Table 4). Beginning in 2003, ZOIs with a mean at approximately 14 km were estimated for all years except 2007 where a ZOI of 7 km was estimated. A significant odds ratio was detected (at α = 0.05) in all years after 2002 except for 2004, 2009 and 2012 (Fig. 5). ROC scores suggested reasonable fit for yearly models.
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.
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.
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.
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.
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).
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.
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.
Dr. Vito Muggeo (Inst. of Social Statistics, Univ. of Palermo, Palermo, Italy) provided assistance in the use of program segmented including adaption of the R code for generalized estimating equation models. The mine aerial survey data sets were supplied by Dominion Diamond Ekati Corporation (Harry O'Keefe) and Diavik Diamonds Mines Inc. Don Russell (Yukon College, Whitehorse, Yukon, Canada) supplied drought index data. We thank Andrea Patenaude (Environment and Natural Resources (ENR), Government of Northwest Territories, Yellowknife, NWT) for support of this project and ENR for providing the collar data.
Funding – This analysis was funded by Environment and Natural Resources, Government of Northwest Territories, Yellowknife, NWT.
Author contributions – JB conducted primary analyses and writing of paper. KP advised on analyses and wrote portions of the paper. AG and JA wrote portions of paper including interpretation of results. JW conducted GIS and spatial analyses and reviewed the paper.
Supplementary information (available online as Appendix wlb-00719 at < www.wildlifebiology.org/appendix/wlb-00719>).