Comparing alternative methods of modelling cumulative effects of oil and gas footprint on boreal bird abundance

Oil and gas activity is increasing in the western boreal forest of North America. To manage cumulative effects of this industry, a better quantification of footprint effects on wildlife is needed. We used point-count surveys to evaluate how well footprint amounts within 150 m, and proximity to seismic lines, pipelines, well sites, roads, and energy facilities predicted the abundance of 48 bird species. We developed models for each species, evaluating the best functional forms for different footprint effects, then predicted how different model structures influenced estimates of regional population size. Most species exhibited at least one nonlinear response to footprint amount (79% of species) or distance (88%). Species associated with older coniferous forests decreased more often with increased footprint amount and closer proximity to footprint, while species associated with open lands and young forests increased with greater footprint and closer proximity. In one-third of species, bird abundance versus distance changed from a positive to negative relationship (or reverse) at a threshold distance from at least one footprint. Models based on footprint proximity had better fit than those based on footprint amount for 29 of 48 species, but both model types produced similar population estimates. Both footprint amount and proximity models were useful for assessing cumulative effects on wildlife and mechanisms causing change. Models of distance to footprint can provide evidence of positive or negative edge effects for developing management buffers. Models of footprint amount provide important information on functional changes in habitat.

structures influenced estimates of regional population size. Results Most species exhibited at least one nonlinear response to footprint amount (79% of species) or distance (88%). Species associated with older coniferous forests decreased more often with increased footprint amount and closer proximity to footprint, while species associated with open lands and young forests increased with greater footprint and closer proximity. In one-third of species, bird abundance versus distance changed from a positive to negative relationship (or reverse) at a threshold distance from at least one footprint. Models based on footprint proximity had better fit than those based on footprint amount for 29 of 48 species, but both model types produced similar population estimates.

Introduction
Wildlife conservation in landscapes undergoing industrial development requires understanding species' responses to combined or cumulative effects of multiple human footprints (Johnson and St.-Laurent 2011). Many wildlife studies involving human footprint place wildlife surveys in different footprint types (i.e., control-impact field designs) (e.g., MacDonald 2000; Manel et al. 2000; Bart and Earnst 2002;Bayne et al. 2016). Most control-impact studies have smaller spatial extents and focus on statistical significance of one footprint type rather than quantifying the effect size of different footprints (Johnson and St.-Laurent 2011). Quantitative cumulative effects assessment (CEA) along gradients of industrial development is needed to: (1) assess if there are nonlinear effects or threshold values of environmental stressors at which wildlife responses change (Toms and Villard 2015); (2) determine if socially accepted thresholds or criteria have been surpassed; and (3) determine relative impacts of different land uses (Bender et al. 1996;Roloff and Kernohan 1999;Burgman et al. 2011;Mahon and Pelech 2021). How to model cumulative effects is not well established, however.
One way to quantitatively assess cumulative effects of human footprint on wildlife is a "dose-response" approach (Karr and Chu 1997), in which a response (metric of habitat suitability or quality for wildlife like abundance, occupancy, nest success, or survival) is compared to the numeric amount or intensity ("dose") of a particular stressor per sample unit. This is a relatively easy approach for quantifying how species react to loss or change of forest due to different types of human disturbance happening within a particular area, after accounting for other variables affecting the response.
In contrast, zone-of-impact studies measure wildlife responses at various distances from human disturbance(s). Such studies can be used to identify buffer zones around human disturbances and evidence for edge effects within habitat patches, by identifying threshold distances from the edge of the patch at which environmental variables (e.g., species abundance within the patch) and distance from human footprint changes sign (e.g., from positive to not positive) (Laurance and Yensen 1991;Ewers and Didham 2007;Johnson and St.-Laurent 2011). Edge effects (separate from edge-related habitat loss or creation) may cause species abundance within habitat fragments to change at distances large enough for survey areas to contain no actual edge. Understanding the depth of an edge effect helps to separate direct habitat loss from effects due to habitat fragmentation. However, whether the distance to the nearest disturbance fully quantifies the magnitude of cumulative effects in areas with multiple overlapping disturbances is not well established.
Due to proximity of multiple footprint types to each other, boreal wildlife may be affected by multiple types of human footprint simultaneously (Johnson and St.-Laurent 2011). A large body of studies has tested for effects of forestry (e.g., Schieck and Song 2006;Huggard et al. 2014) and other human footprint types (e.g., oil and gas) on boreal forest wildlife [e.g., ungulates (Hebblewhite 2011); small mammals (Darling et al. 2019); birds (Bayne and Dale 2011;Bayne et al. 2016;Mahon et al. 2019)]. Energy development effects on wildlife include increased direct and indirect mortality and reductions in individual fitness and reproductive success (Bayne and Dale 2011;Hebblewhite 2011;Johnson and St.-Laurent 2011); temporary or permanent loss of older forests accompanied by increases in early successional, open, and non-vegetated lands (Bayne and Dale 2011); edge effects (e.g., changes in microclimate, noise, predation) extending from the footprint into the surrounding forest (Bayne and Dale 2011); and persistent disturbance and forest fragmentation patterns very different from both natural disturbances and timber and pulpwood harvest (Pickell et al. 2015). Oil and gas footprint types in boreal forests vary with respect to each other and to forest harvest activity in terms of individual footprint size, intensity of disturbance, cumulative extent, and the relative amount of edge bordering boreal forest wildlife habitats. Processing facilities and other infrastructure have the highest disturbance intensity and largest individual footprint sizes but are limited in number and cumulative extent. In contrast, linear footprint like seismic lines, pipelines, and roads have small footprint sizes locally, but are numerous and extensive, creating considerable edge habitat relative to processing facilities (Jordaan et al. 2009). Well sites are intermediate in size, wider than most linear features but smaller than most processing facilities (Boxall et al. 2005). Well sites result in lower amounts of edge than linear footprint (Jordaan et al. 2009) but are numerous and may result in considerable habitat loss (Wells et al. 2020).
In the sedimentary basin of western Canada, there have been > 30 years of studies on boreal bird habitat associations (Cumming et al. 2010). Conclusions about energy sector effects from these different study methods broadly agree that bird species associated with young forests and open lands likely benefit from oil and gas disturbances (Bayne and Dale 2011;Bayne et al. 2016). These studies have either compared effects of footprint within different treatments or examined effects of footprint amount rather than distance to footprint. We do not know how well dose-response and zone-of-impact methods of measuring footprint effects agree quantitatively, if one method has greater predictive accuracy, or how different regional population predictions can be using one method over another. Thus, we used a standardized data set collated from past boreal bird studies to directly compare dose-response and zone-ofimpact methods of modeling cumulative effects of oil and gas footprint on individual boreal bird species. Dose-response methods have been used to analyze effects of the amount of different oil and gas sector footprints on boreal bird abundance within spot mapping plots or point counts as the response (Bayne et al. 2005;Lankau et al. 2013;Mahon et al. 2019;Wilson and Bayne 2019;Sólymos et al. 2020).
Our objective was to compare dose-response and zone-of-impact methods of modelling energy development impacts on boreal bird abundance. First, we examined if both methods produced similar direction of effects. Second, we determined if both methods produced similar regional population estimates of species. Finally, we compared the predictive ability of the best dose-response and zone-of-impact models in each of 48 boreal bird species.

Study area
Data were collected in boreal forest landscapes north of Edmonton, Alberta (Fig. 1). Dominant vegetation consisted of: lowland bogs and fens where black spruce (Picea mariana) and tamarack (Larix laricina) are the primary tree species; mesic uplands dominated by deciduous and mixed-wood forests of trembling aspen (Populus tremuloides), balsam poplar (P. balsamifera), and white spruce (Picea glauca) with smaller amounts of balsam fir (Abies balsamea) and paper birch (Betula papyrifera); and xeric uplands dominated by jack pine (Pinus banksiana). Open vegetated lands in the study areas were associated with non-treed bogs and fens, wetlands including beaver ponds, and shrublands and grass growing up after forest fires, harvest, or in disturbances associated with oil and gas footprint (Alberta Biodiversity Monitoring Institute 2017).

Point counts
We used point counts (N = 20,055 visits, 11,569 stations) from multiple bird studies completed between 2000 and 2019. Point count methods differed among studies in point count duration (3-10 min) and whether distancesampling (0-50 m, 0-100 m, and 0-unlimited distance) was used (Boreal Avian Modelling Project 2015). Overall, 71% of the point counts came from control-impact studies where stations were centered within different forest stands and ages or energy sector disturbances by design (e.g., Bayne et al. 2016;Mahon et al. 2016). Only 6% came from 3-min roadside Breeding Bird Survey counts from 2000 to 2005, with 1-2 randomly selected visits from different years per station. More recent BBS data are available for Alberta but we lacked temporally accurate vegetation and footprint data for surveys after 2005. The remaining 23% of point counts were in 10 × 10 square arrays of point counts spaced 600 m apart: these arrays were located near steam-assisted gravity drainage facilities in northern Alberta and varied in the density of oil and gas footprint. Point counts in these arrays were conducted using autonomous recording units (Shonfield and Bayne 2017) and were not necessarily centered on or away from any energy sector disturbances (Fig. 1).
To account for varying detection probabilities for each species due to differences in survey methods, time of year and day of survey, and varying sound attenuation with distance in different vegetation types, we used the QPAD approach (Solymos et al. 2013). In brief, QPAD creates a statistical offset that adjusts the mean count based on known relationships between bird abundance and the area sampled (i.e., point count radius) and duration of point count. The QPAD approach also enabled us to predict bird abundance in a common unit of density (#males/ha) (Sólymos et al. 2013). While the QPAD approach can be used to estimate detection offsets for non-songbirds, we limited our analyses to detections of singing birds for 48 songbird species. Code for QPAD is available at https:// github. com/ psoly mos/ QPAD/.

Quantifying vegetation and human footprint
Forest age classes (including harvested forests) and nonforest habitats [shrubland, grassland, wetland, seismic lines, pipelines, well sites, roads, and facilities (e.g., power plants, compressor stations, processing facilities)] were quantified within 150 m of each point count. We chose this spatial scale for analysis as it is the approximate radius over which most passerine songbirds are detected using autonomous recording units (Shonfield and Bayne 2017). We buffered locations in ArcGIS and intersected vegetation and footprint data from the 30-m wall-to-wall Alberta Biodiversity Monitoring Institute Vegetation Inventory (v. 6.1 vegetation/back-filled layer), aligning the point count data with land cover data within 150 m from the year of a particular survey (Alberta Biodiversity Monitoring Institute 2017, 2018). We summarized forest structure and composition using four variables: proportion of all non-forested lands excluding human footprint apart from harvest areas and including any stands < 10 years old (proportion open); area-weighted mean forest age within 150 m (forest age), proportion of all stands within 150 m that were conifer-dominated (proportion conifer), and proportion of all lands within 150 m that were forested or non-forested wetlands (proportion wetland). Shrubland and grassland were included within proportion of open land; other types of open land that were human footprint (e.g., cropland, surface mines) occupied negligible cover (< < < 1% of land within 150 m) on very few point counts. We did not model effects of forest harvest per se in our analysis because our primary focus was on comparing alternative methods of measuring oil and gas footprint effects, rather than building the best possible cumulative effects models for birds. We used the "Near" tool in ArcGIS 8.3 to calculate the nearest distance to each of these footprint types from each point count (Environmental Systems Research Institute 2002).

Modelling each cumulative effects assessment approach
We ran models for 48 species of songbirds that were detected on at least 1% of point count visits ; Table 1). These species varied in their preferences for open land versus complex older forests, coniferous versus deciduous forest, and upland versus wetland habitats  and their responses to different types of energy sector development in previous studies Mahon et al. 2019).
We used generalized linear models ["glm" function in R (R Core Team 2020)] with a Poisson error distribution to create both dose-response and zone-ofimpact models. We first built a set of habitat models for each species that also included a linear term for year. While some sites had multiple visits within or across years and/or were clustered spatially, the large number of random effects required, the difficulty of assigning random effects in an objective manner, and the large number of single survey point counts within our data made the use of mixed-effects models prohibitive. Instead, as latitude and longitude were not multicollinear with other predictors at point counts, we used latitude and longitude to partially account for spatial correlation in the data. Different habitat models for each species included proportion open, forest age, proportion conifer, and proportion wetland within 150 m of point counts. Both linear and quadratic functions of each variable were assessed in case species were more abundant at intermediate values of those variables except for proportion of wetland, since including a quadratic term for that predictor increased multicollinearity. Starting with a global model, we compared models where all four habitat variables were present along with models where one or more habitat variables or terms of habitat variables were excluded. A null model (year + latitude + longitude, but no habitat predictors) was also evaluated in the habitat modelling stage. We selected the model with the lowest Akaike's Information Criterion (AIC) (Burnham and Anderson 2002) as the best habitat model for each species, and ran models using the dredge function in the "MuMIn" package in R (Barton 2009). Since model selection via lowest AIC could favor more complex models, we checked whether habitat effects and terms were significantly positive or negative (i.e., 95% confidence intervals for effect size excluded zero).
In dose-response models, using the best habitat model for each species as a null model, we used the amount of each footprint (seismic line, pipeline, road, well, facility) within 150 m of each point count to measure disturbance effects. We did not model amounts of railways, gravel pits, towns, cropland, mines, or transmission lines because there were negligible amounts of these footprint types within 150 m of few point counts. We used AIC to identify the best functional form of each footprint amount, evaluating linear and quadratic functions to see if species varied linearly or with intermediate amounts of footprint, inverse functions to see if species strongly varied with even small amounts of footprint, and squareroot functions to see if rates of species increase or decrease tapered off beyond some value of footprint amount. To model an inverse functional relationship, all footprint proportions < 0.01 were set to 0.01. Once we determined the best functional form for seismic line, pipeline, well site, road, and facility amount, we developed a global dose-response model containing these functional forms.
For zone-of-impact models, we measured the shortest distance to each of the five footprint types instead of the amount of each footprint. In contrast to dose-response models, negative effects of a given footprint are associated with positive coefficients for distance to that footprint. We evaluated linear, inverse and square-root functions of footprint distance, setting the minimum distance to each footprint to 1 m when modelling inverse functions. We did not consider quadratic functions of distance amount since we did not expect for bird species to reach maximum or minimum abundance at intermediate distances from footprint. We then developed a global zone-of-impact model containing the best functional forms for each footprint. We identified possible evidence of edge effects of footprint when species were best explained by an inverse function of distance to that footprint, but only if a sharp change in edge response occurred at a distance > 150 m. At that distance, there was none of that footprint within 150 m of a point count and therefore no habitat lost to that footprint. Although habitat could be lost to footprint at larger scales, 150 m was a large enough distance to contain multiple territories of the songbird species we analyzed.
Predictive ability and population consequences of each modelling approach We compared how well different approaches predicted boreal bird responses to cumulative energy sector development. We used AIC and pseudo-R 2 as measures of explanatory power and relative goodness To model inverse functional response to footprint seen in these plots, proportions of footprint within 150 m < 0.01 were set to 0.01 and minimum distance to footprints was set to 1 before obtaining predictions  Mahon et al. (2016) but were assigned to those habitat associations based on their ecology. Within each habitat association, species are listed in order of decreasing prevalence on point counts in this study  Using the final dose-response and zone-of-impact model for each species, we predicted total abundance in a landscape of adjacent 300-m cells in a 6-millionhectare area of NE Alberta that aligns with the Alberta Pacific Forest Industries Ltd. Forest Management Area (Al-Pac FMA). We generated 300-m resolution raster layers (i.e., ~ equal to the width of the 150-m radius point counts that we modelled) of each model predictor term that potentially could be included in the final models generated for each species. We created raster layers of forest vegetation variables from national 250-m layers clipped to Alberta, resampling via nearest-neighbours methods to 300-m resolution (Beaudoin et al. 2014). We created raster layers from human footprint feature layers, including but not limited to the footprint types within our models (Alberta Biodiversity Monitoring Institute 2018). Although forest harvest was among these human footprint types and is the largest source of human disturbance within the Al-Pac FMA, we assumed that its main effects on forest structure were to reduce mean forest age and increase the proportion of open vegetated lands within the raster layers of forest vegetation. Forest harvest activity also changes the temporal and spatial patterns of forest stand ages and types across the landscape, which will influence the total availability of suitable habitat for boreal birds associated with different forest types, ages, and structural stages. However, these spatial patterns of forest stand ages and types would be captured at least at a coarse spatial scale by the raster layers we used. We also estimated the shortest distance from the centroid of each cell to each of the five footprint types analyzed (seismic lines, pipelines, well sites, roads, facilities), using the "Near" proximity tool in ArcGIS, and created five 300-m distance raster layers, where the value in each cell was the nearest distance to one of the five footprints.

I(+) S(+) Q(+) I(+) Alder Flycatcher YP I(−) S(+) S(+) S(+) S(+)
Predicted density layers for each species (#males/ ha) was multiplied by 9 (# hectares per 300-m cell), then abundance in each cell was summed for the whole study area, excluding the cells containing human footprint that were not included in our models (e.g., railways, gravel pits, towns, cropland, mines, transmission lines: 116,554 out of 789,259 cells or 15% of the cells within the Al-Pac FMA, mostly along major highways, and within lakes, settlement areas, and mined areas north of Fort McMurray, Alberta). While these human footprints occurred in the study area (northern Alberta), they occupied negligible cover (< < < 1% of land within 150 m) on few point counts. We produced maps and plotted the relationship between predicted population estimates from dose-response and zone-of-impact models for each species, using each species as an observation.

Habitat models
In the best habitat models for the 48 species, those predicted to be associated with younger forests and open productive ecosystems generally responded positively to the proportion of open land and declined with mean forest age within 150 m of point counts. Species associated with older upland deciduous forests decreased with increasing coniferous cover. Species associated with older upland coniferous forest age increased with coniferous cover and with forest age and decreased with open land within 150 m. Species associated with relatively unproductive lowland  Mahon et al. (2016) but were assigned to those habitat associations based on their ecology. Within each habitat association, species are listed in order of decreasing prevalence on point counts in this study Table 3 The "best" model function for describing the zone-of-impact effect of the nearest distance to different footprint types from each point count, from the global zone-of-impact models predicting 48 species of boreal songbirds in northern Alberta

Direction and magnitude of oil and gas footprint effects
Out of 240 possible effects of each of footprint amount and distance to footprint (48 species, 5 footprint variables per species), more species associated with young, productive ecosystems or with older deciduous forests ) responded positively to footprint amount within 150 m (34%) and closer footprints (42%) than responded negatively to footprint amount (13%) or proximity (22%). In contrast, more species associated with older coniferous upland forests responded negatively to footprint amount (34%) and proximity (36%) than responded positively to amount (8%) or proximity (12%) (e. g., Black-throated Green Warbler (Setophaga virens), Fig. 2). For species associated with either older deciduous upland forests or with lowland coniferous forests, there were slightly more positive than negative species responses to footprint amount and slightly more negative than positive species responses to proximity to footprint (Tables 2 and 3, model coefficients and plots in Online Appendices 2-3).
Most species exhibited at least one non-linear response to footprint in both dose-response models (79% or 38 of 48 species) and zone-of-impact models (87.5% or 42 of 48 species). One-third of species varied with an inverse function of distance to at least one footprint type, indicating that abundance of these species strongly changed at a threshold distance from footprint. However, we found few species significantly changed in abundance (i.e., 95% confidence intervals of effect size excluding zero) at threshold distances greater than 150 m from footprint, suggesting a lack of strong edge effects separate from habitat loss (Tables 3 and 4, Online Appendix 3).

Predictive ability and population consequences
The global zone of impact model for each species had a better fit (lower AIC, higher pseudo-R 2 ) than the habitat or global dose-response model in 29 of 48 species (Table 5). For most species, the relative difference in population size predicted by each type of global model was small [relative ratio ~ 1 (Fig. 3), Table 6]. However, there were also some dramatic differences in population estimates and density distribution maps for some species (e.g., Boreal Chickadee) (Fig. 4, Table 6, Online Appendix 4). Nevertheless, there were strong positive correlations between population estimates from the habitat, dose-response and zone-of-impact models for each species overall (Spearman r > 0.9 for all pairwise comparisons).

Discussion
We compared regional-scale predictions of boreal bird abundance as a function of oil and gas footprint based on dose-response versus zone-of-impact

(−) S(+) S(+) L(−) L(−) Palm Warbler LC S(+) 0 L(−) S(−) S(−)
Type of function: I inverse; L linear; S square root; 0 = no effect [footprint not included in the best model] or the overall effect was non-significant (i.e., 95% confidence intervals for effect size overlapped with zero) even if it was included in the best dose-response model; (−) = negative effect of footprint; (+) = positive effect of footprint on species abundance. In contrast to doseresponse models, positive effects result when bird abundance decreases with footprint distance and negative effects result when bird abundance increases with distance from footprint. *Implausible model coefficient estimates from global model, based on the "best" functional forms of each footprint. Habitat associations (from Mahon et al. 2016): YP young productive ecosystems (young deciduous forests, shrublands, open lands, marsh); OD older, structurally complex deciduous forests; OC older, structurally complex upland coniferous forests; LC lowland coniferous forests. Species in italics were not mentioned in Mahon et al. (2016) but were assigned to those habitat associations based on their ecology. Within each habitat association, species are listed in order of decreasing prevalence on point counts in this study models and evaluated the predictive accuracy of both sets of models. Zone-of-impact models more often had higher explanatory power than dose-response models for a given species. As well, footprint impacts were detected for more bird species using zone-of-impact than dose-response models. Our dose-response models considered amounts of footprint only within 150 m of point counts rather than at larger, potentially influential spatial scales (e.g., Mahon et al. 2019). In contrast, effects based on distance to footprint were not linked to a particular spatial scale and could be a function of both local and landscape-scale influences of footprint that were not present in our dose-response models. Whether this characteristic of zone-of-impact models helps explain their greater model fit for some species could be assessed by modelling effects of footprint amount within larger radii around point counts. That dose-response and zone-of-impact models generally predicted similar regional population outcomes was not surprising given the reasonably strong negative correlations (Spearman r < − 0.5) between amount of footprint within 150 m and nearest distance to footprint for most stressor types (seismic: − 0.70; pipeline: − 0.51; well site: − 0.51; roadside: − 0.50). There was a weaker negative correlation between the proportion of and distance to facility footprint (− 0.35), which may be conditional on how our study was designed and analyzed. Our collated data set included data from different studies designed to test specific gradients of energy sector footprint, distances to edge, and other objectives. Thus, we did not have a balanced experimental design with similar numbers of surveys across different footprint types, systematically located across our study area. Further, modelling dose-response and zone-of-impact effects at local spatial scales probably weakened the correlation between facility amount and distance to facility. Facility footprint is the least extensive at large spatial scales, being less numerous and more concentrated in smaller areas than other footprint types (Boxall et al. 2005). Further, closer to a facility, there is no room for additional footprints to be present except the lines feeding the facility. As a result, there is a much larger range of distances to the nearest facility than any other footprint type. Facilities can be very large and are more likely to cover the entire area sampled for birds (footprint proportions range from 0 to 1 for facilities) than the other footprint types (i.e., footprint proportions range from 0 to 0.25 for seismic). Since inherently different ranges are possible for footprint types, different approaches to modelling footprint types might be required when creating large-scale regional models of cumulative effects (e.g., Sólymos et al. 2020). Dose-response models may be more useful for examining additive and interactive effects of multiple widely distributed footprint/disturbance types in complex multi-stressor, multi-sector landscapes where footprint types are not spatially or temporally independent (Mahon et al. 2019). In areas with lower amounts of footprint types, greater dispersion of footprint types, or fewer footprint types, zoneof-impact models might be useful because distance effects can be examined independently (spatially independent) and accurately.
Strong correlations between dose-response and zone-of-impact measurements for most footprint types suggests that either model could be used in regulatory decisions. We caution against relying entirely  Canada Jay  LC  -----Palm Warbler  LC  -----Thresholds indicate where abundance of a species strongly changes with increasing distance. Cells are emphasized according to the overall effect of each footprint type on abundance of each species (bold italics = positive [abundance decreases with increasing distance from footprint], bold = negative [abundance increases with increasing distance from footprint], "−" = inverse function for a particular footprint was not included in the best zone of impact model for a given type of footprint for a particular bird species, so there was no threshold. Within each habitat association, species are listed in order of decreasing prevalence on point counts in this study on one framework over the other. Instead, creating and presenting both kinds of models can provide important complementary insights into the processes generating impacts of human footprint on wildlife. Dose-response models may be more appropriate for modelling the effect of total area disturbed by different combinations of footprint at varying spatial scales, which cannot be done using zone-of-impact models. However, the mechanisms (i.e., habitat loss versus fragmentation) are more challenging to separate with dose-response models. Non-forest vegetation and unvegetated footprint types can serve as barriers for some species or introduce edge effects (Murcia 1995) but dose-response designs do not allow us to explicitly test for the magnitude of habitat loss versus habitat fragmentation.
Zone-of-impact models may be used to detect edge effects of footprint on species of wildlife (Laurance and Yensen 1991;Ewers and Didham 2007). There is stronger evidence for an edge effect of footprint if a species' abundance strongly changes at distances large enough that there is no footprint within species' territories within survey areas. For songbirds, evidence of edge effects could be associated with threshold distances greater than the radius of a point count: a 150-m point count is large enough to contain multiple territories of most forest songbird species. Further, many forest songbird species have song characteristics that make it unlikely to detect those species within forests beyond 150 m using point counts, or detect conspecifics using habitat including footprint outside   of the point count radius (Matsuoka et al. 2012). Cumulative edge effects could be significant given the amount of habitat edge created by well-sites, linear footprint, and forest harvest in boreal forests (Wells et al. 2020) and could be positive or negative for different species (Murcia 1995). Although we did find strong changes in abundance of species at threshold distances, these distances were not large enough to exclude the possibility of habitat loss or creation by footprint. Even if we had identified large thresholds, mechanisms of edge effects (e.g., noise, microclimate changes, resources, predators) for different species in this paper have not been identified. Another explanation is that as there are multiple human footprint types (including harvest) that are close to or overlapping with each other, i.e., impacts which may not be independent. In our study, point count distance from pipelines was positively correlated with distance from roads (r = 0.72); otherwise, correlations between nearest distances to different footprint types were weak (r < 0.5). However, correlations among nearest distances to different footprint types vary spatially, so that potential edge effects of one footprint could vary with the proximity of other footprint types. Further, the collection of survey locations we used did not have equal numbers at varying distances in all vegetation types. The way a species responds to distance-to-edge will depend on the suitability of the vegetation type for that species. Future analysis should evaluate whether different thresholds are observed in different vegetation types, or with vegetation age or recovery within footprint (Lankau et al. 2013;Leston et al. 2018). Understanding edge effects is essential as it has major implications for conservation or management of wildlife habitats as estimates of suitable habitat and Fig. 4 Predicted density distribution or difference in predicted densities (#males/300-m cell) for Boreal Chickadee (Poecile hudsonicus) in #males/ha in the Al-Pac FMA based on models for this species. Cells containing footprints other than the types included in models were excluded from these maps and when calculating population estimates in study area. a Doseresponse model predictions; b zone-of-impact model predictions; c dose-response minus zone-of-impact model; d habitat minus zone-of-impact model population size change considerably when negative edge effects are found (Laurance and Yensen 1991;Ewers and Didham 2007), in addition to whatever vegetation has been lost due to oil and gas development (Johnson and St.-Laurent 2011). We found more negative effects of footprint on species associated with older coniferous upland forests. This difference could be caused by faster vegetation regrowth on abandoned disturbances in deciduous forest (Nijland et al. 2015). Importantly, old coniferous and mixed-wood forest are some of the least common habitats in our study region and at greatest risk of loss/ disturbance from other extensive and intensive industrial activities such as industrial forestry that targets both deciduous and mixedwood forests for pulp and paper, but also conifer forests for saw timber (Hobson and Bayne 2000). Forestry remains the largest human activity in the boreal forest in terms of area disturbed (Wells et al. 2020) but is usually a temporary disturbance. While timber and pulpwood harvest size and distribution can differ strongly from natural disturbance like forest fires, some harvest plans attempt to emulate natural disturbance patterns in long-term effects on wildlife communities (Schieck and Song 2006;Huggard et al. 2014). Nevertheless, effects of forestry could compound energy sector effects on wildlife, especially if both industries are operating in the preferred stands and forest ages used by a species.
Future cumulative effects studies can be improved in a few ways. First, such studies can consider a wider variety of footprints or whether different footprint types can be combined in larger categories (Mahon et al. 2019). Second, future studies could consider effects of footprint amount at larger spatial or multiple spatial scales with co-occurring footprints including harvest [e.g., Bayne et al. (2005), Mahon et al. (2019)]. Third, abundance or occupancy studies will not necessarily assess how human footprint affects habitat quality (Van Horne 1983). Abundance or occupancy studies should be accompanied by studies measuring individual fitness (e.g., nest success, productivity, survival) as indicators of habitat quality (Ball et al. 2009), although the latter types of studies are typically labor-intensive and limited in extent. Fourth, future studies can consider interactive effects of different footprint types in proximity to each other in the same landscapes (Johnson and St.-Laurent 2011;Mahon et al. 2019). Finally, future studies can combine dose-response and zone-of-impact predictors within the same models. author contributions LL processed and analyzed data, organized results presented in this manuscript, and contributed to writing. EB procured funds, provided data, proposed the analytical design, and contributed to writing. JT and CLM provided data, proposed the analytical design, and contributed to writing. AC, PS, JB, SJS, FKAS, DS, and TD contributed to writing. Data availability Data used in modelling and scripts used in this analysis are stored at: https:// github. com/ Lione lLest on/ Cumul ative-effec ts-models-Alber ta. Model outputs can be recreated from this data. Additional data for predicting regional populations is available from the Boreal Avian Modelling Project at http:// borea lbirds. ca.

Declarations
Competing interests There are no competing interests that would have influenced the results of analysis in this manuscript. Co-authors have no financial stakes in the outcomes from this study and results from this study are independent from the position of the funding agency.