Additive and interactive cumulative effects on boreal landbirds: winners and losers in a multi-stressor landscape.

Stressors created by multiple resource industries can result in cumulative effects over time and space. Many studies have evaluated single stressors and assumed that cumulative effects can be understood by adding stressors together. However, there is growing evidence that interactive effects are important in structuring biological communities. We evaluated whether the effects of multiple stressors in the boreal forest (linear features, energy, forestry) combine additively or interactively by testing a candidate model set of 12 cumulative effects models of abundance for 27 landbird species. We fitted paired additive and interactive Generalized Additive Models and examined model predictions in the Athabasca Oil Sands Area of Alberta, Canada, and a theoretical no-disturbance version of the study area. Of the 27 species examined, an additive disturbance model was the best for nine species, while an interactive disturbance model was the best for 11 species. In the current study area, disturbance models predicted strong increases in abundance for species associated with deciduous forest and open habitats (winning species) and moderate decreases for species associated with conifer forest (losing species). We found a 15% change in landbird community composition between the current study area, with 8.4% disturbance, and the theoretical no-disturbance study area. Complex synergistic and antagonistic interactions among stressors were observed for 39% of landbird species, with the majority of interactions observed being synergistic. Stressors with relatively small disturbance areas, such as narrow linear disturbances, frequently interacted with other stressors to affect species' responses, and energy sector stressors often had additive or interactive effects with forestry stressors. Interactive cumulative effects from multiple sectors will make it increasingly difficult for industry and land managers to manage impacts unless interactions among stressors are incorporated into cumulative effects assessments and regional land use planning processes. © Her Majesty the Queen Right of Canada 2019. Reproduced with the permission of the Minister of Environment and Climate Change.

Introduction

Human activity has altered much of the North American landscape in the last two centuries, resulting in considerable shifts in species composition (Fischer and Lindenmayer 2007). Until recently, western boreal forests have remained largely undeveloped, with the exception of agriculture at the southern edge (Hobson et al. 2002). However, over the past two decades, simultaneous development by multiple resource industries has intensified within this region (forestry, bitumen/oil sands extraction, conventional oil and gas extraction, and peat mining), particularly in Alberta (Schneider et al. 2003). This pattern of development is intensive at local scales (e.g., open pit mines) and extensive at landscape and regional scales (e.g., seismic lines). These resource development activities result in the removal of native vegetation and have led to hundreds of thousands of kilometers of edge caused by linear disturbances such as roads, railways, powerlines, seismic lines, and pipelines (Lee and Boutin 2006, Pattison et al. 2016). These activities are known to impact boreal landbird communities (Hobson and Schieck 1999, Bayne et al. 2005, 2016), and will hereafter be referred to as stressors.

The decrease and alteration of native vegetation due to resource or industrial development is known to be one of the major processes affecting the quantity of breeding habitat for landbirds in the western boreal forest (Hobson and Schieck 1999, Bayne et al. 2016), and more broadly throughout North America (Thompson et al. 2002, Faaborg et al. 2010). At a regional scale, the western boreal forest is transitioning from an intact landscape driven by natural processes to a variegated landscape driven by both natural and human processes. A variegated landscape has low to moderate landscape modification (>60–90% native vegetation cover remaining) and a high occurrence of novel ecological boundaries between native vegetation and disturbed land (McIntyre and Hobbs 1999). The increase in habitat isolation and creation of edge in areas of remaining native vegetation may affect breeding habitat quality for landbirds (e.g., altered food and prey availability, competition, predation, or parasitism), possibly resulting in changes to breeding‐season behaviors (e.g., shifted territory location, altered movement patterns or disrupted dispersal ability; McIntyre and Hobbs 1999, Fischer and Lindenmayer 2007).

Many breeding bird species exhibit negative responses to disturbance, with the most specialized species typically exhibiting the largest responses (Julliard et al. 2004, Devictor et al. 2008). Habitat disturbance and degradation negatively affect specialists through increased competition with generalists and increased extinction or extirpation risk caused by altered ecological interactions between species (Clavel et al. 2011). Specialists may be more sensitive than generalists to disturbance because of the narrow range of resources or conditions they require to maintain population viability. Community shifts, and in particular the decline of specialist species, in disturbed landscapes can lead to biotic homogenization (Devictor et al. 2008, Clavel et al. 2011). Biotic homogenization is a process by which the genetic, taxonomic, or functional similarity of two or more species assemblages increases over space and/or time (Olden 2006). It has been hypothesized that the western boreal forest may be starting to undergo a process of biotic homogenization (Mahon et al. 2016), and investigating the impacts of multiple stressors on boreal birds is key to developing effective mitigation and best management strategies for reducing cumulative effects and conserving heterogeneous bird communities.

Most research on the effects of forestry and energy activity on birds in the western boreal forest has focused on single stressors (Hobson and Schieck 1999, Bayne et al. 2005, 2016), or assumed that the cumulative effects of multiple stressors can be assessed additively (the cumulative effect is assumed to be equal to the sum of individual stressors; Mahon et al. 2014). However, it is critical to look beyond additive effects and examine the interactions of multiple stressors, which can give rise to unexpected and complex relationships (Didham et al. 2007, Crain et al. 2008, Darling and Côté 2008). There is growing evidence of interactive effects as key factors structuring terrestrial and marine biological communities (Brown et al. 2013, Cartwright et al. 2014). Interactive effects can be synergistic (the interaction of stressors increases the magnitude of the response above that seen for additive effects; A + B + A × B > A + B) or antagonistic (the interaction of stressors decreases the magnitude of the response below that seen for additive effects; A + B + A × B < A + B). Interactive effects represent major uncertainties for long‐term planning and ecosystem management (Brown et al. 2013). In particular, failure to incorporate synergistic interactions among stressors into cumulative effects assessments and land use planning processes may result in unplanned changes to ecosystems and increased rates of biodiversity change beyond those expected.

Development in the boreal forest is still in the early stages, creating opportunities for proactive conservation and strategic management to mitigate the impacts of stressors created by multiple resource sectors. However, effective and efficient management strategies require understanding additive and interactive effects on boreal bird populations. If direct effects (i.e., loss of native vegetation) are the primary influence on bird populations, then additive models of multiple stressors should be the best predictors of species’ response to resource development. Alternatively, if both direct and indirect effects (i.e., altered ecological interactions and changes in behavior associated with habitat isolation and edge effects) are influencing bird populations, then interactive models of multiple stressors may be better predictors of species’ response to development. Given the rapid development in the western boreal forest, we anticipate that landscape modification from multiple resource sectors may result in effects that cannot be predicted by loss of native vegetation alone, i.e., that both direct and indirect effects will be present. Our objectives are to (1) assess the impacts of multiple resource stressors on boreal bird density; (2) evaluate whether the effects of multiple stressors are additive or interactive, and to assess the effect size and interaction type (synergistic or antagonistic) of interactive effects; (3) compare boreal bird density and community similarity between the current study area (multiple stressors) and the same study area if human disturbances did not exist (no stressors).

Methods

Study area

We assessed the impacts of human disturbance in the Athabasca Oil Sands Area (AOSA), within the larger Oil Sands Area in northern Alberta, Canada (Fig. 1). The AOSA (9.3 million ha) falls within the Boreal Plains ecozone. Summer (May, June, July, and August) mean temperatures range from 7.2°C to 20.2°C and mean total precipitation is 2.4 cm. Mesic sites in upland areas are dominated by pure and mixed stands of trembling aspen (Populus tremuloides) and white spruce (Picea glauca) mixed with balsam poplar (Populus balsamifera), white birch (Betula papyrifera), and balsam fir (Abies balsamea), while drier upland sites are dominated by jack pine (Pinus banksiana). Lowland areas are composed of wetlands in the form of marshes, swamps, and black spruce (Picea mariana) and tamarack (Larix laricina) dominated bogs and fens. The AOSA is a heterogeneous patchwork of upland and lowland habitats, where patch size for lowland and upland habitats is 13.6 ± 106.4 ha (mean ± SD), and 9.0 ± 80.3 ha, respectively.

Figure 1

Study area in the Athabasca Oil Sands Area (AOSA) in northern Alberta, Canada.

We define disturbed lands (hereafter disturbances) as all areas originating from human modification through removal or alteration of native vegetation for linear features, forestry, and oil, gas, and bitumen exploration and extraction. Areas originating from wildfire, the primary natural disturbance in western boreal forests, are included as a native vegetation type. Forest harvesting, conventional oil and gas development, and bitumen development, including open pit mines (20% of bitumen development) and steam assisted gravity drainage (80% of bitumen development), affected 8.4% of the accessible portion of the AOSA (7.4 million ha). Major stressors in the AOSA result from forest harvesting (harvest units, roads), energy exploration (seismic lines, exploration wells) and extraction (pipelines, production wells, industrial facilities), and infrastructure construction (roads, railways, powerlines). We did not include agriculture and peat mining in this analysis because they impact only small areas at the southern edge of the AOSA.

Avian data

Environment and Climate Change Canada conducted point‐count surveys in 2011–2014 within the Oil Sands Areas (for further details, see Appendix S1). We surveyed birds throughout most of the AOSA (Fig. 1), but not within the mineable oil sands area due to safety concerns. Surveys within the AOSA targeted disturbance effects and representative habitats. We grouped individual point counts into survey clusters to maximize sampling efficiency (by air or ground travel), minimize safety risks for field staff, and meet targeted or stratified sampling objectives for disturbance and habitats. Within a survey cluster, we separated point‐count sites (≥9) by a minimum of 300 m, the maximum detection radius for most boreal songbird vocalizations (Matsuoka et al. 2012). We followed the standard recommended protocols of Ralph et al. (1993) and Matsuoka et al. (2014) and conducted point‐count surveys during suitable weather conditions between official sunrise and 4–5 h after sunrise, from the last week of May to the first week of July. We conducted point counts for 10 min, using three distance bands (0–50 m, 50–100 m, and >100 m).

The majority of observations in a point count comprised zero or one individuals (unpublished data) of each bird species. Therefore, for this analysis, we combined bird observations from the nine point counts in a survey cluster in order to have a more finely discriminated response variable and ensure that we were assessing avian response at a scale that would include multiple stressors. In total, we selected 303 survey clusters from our compiled data. We ensured that survey clusters were separated by a minimum of 800 m and selected clusters that representatively captured the underlying gradients of disturbance type and intensity within the AOSA (Table 1, Fig. 2), including some survey clusters without any disturbance. The sample unit in our analysis was the survey cluster of point counts and the 300 m area around each point count. Because stressors associated with energy and forestry sectors had a non‐uniform or patchy distribution across the AOSA, we ensured that a minimum of 25% of survey clusters (and the 300 m around each point count) contained disturbance data > 0 for each variable in the statistical models (see Cumulative effects models). We restricted data analysis to 27 species that occurred in >30% of survey clusters, and for which the negative binomial model provided a good model fit.

Table 1

Mean and range of environmental, habitat, and disturbance variables tested in the 12 disturbance models across the 303 survey clusters

Predictor variable	Variable type	Prevalence (percentage of occurrence)	Mean	10th and 90th percentiles
TAREA (ha)	environmental	N/A	206	181–237
AGE (yr)	habitat	N/A	82	43–114
SWFB (percentage of area)	habitat	0.80	13.2	0–34.5
HDWD (percentage of area)	habitat	0.96	38.6	2.7–75
LWLD (percentage of area)	habitat	0.90	17.6	0–51.8
CON2 (percentage of area)	habitat	0.83	13.0	0–36.9
OPNF (percentage of area)	habitat	0.85	7.1	0–17.7
NL (percentage of area)	disturbance	0.95	1.3	0.2–2.4
WL (percentage of area)	disturbance	0.65	1.9	0–5.1
HU (percentage of area)	disturbance	0.40	5.7	0–21.3
WE (percentage of area)	disturbance	0.61	0.9	0–1.8
All Disturbance (percentage of area)	disturbance	0.95	10.0	0–25.7

The prevalence or occurrence of habitat and disturbance classes within survey clusters are calculated using the percentage of occurrence of each variable. Predictor variable definitions: TAREA, total area sampled by each survey cluster; AGE, area‐weighted stand age; SWFB, white spruce/balsam fir; HDWD, hardwood (combined trembling aspen, white birch, balsam poplar); LWLD, black spruce/larch lowlands; CON2, black spruce/pine uplands; OPNF, open non‐forested habitat (shrublands, grasslands, recent burns <20 yr old); NL, narrow, vegetated linear features; WL, wide linear features (pipelines, powerlines, all roads); HU, harvest units with harvest activity <20 yr old; WE, bitumen, oil, and gas wells; All Disturbance, total area of all disturbances.

Figure 2

Depiction of landscape disturbance gradients within the survey clusters and the larger AOSA study area for the four patterns of development tested in this study: (1) linear (model group 2); (2) linear and energy (model group 3); (3) linear and forestry (model group 4); and (4) cumulative (linear, energy, forestry; model group 5).

Our response variable was the survey‐cluster‐level density of each landbird species (number of singing individuals per ha). Raw bird counts are an incomplete measure of abundance and need to be adjusted for species‐specific detection probabilities and detection radii (Sólymos et al. 2013). To account for imperfect sampling detectability, we calculated an offset correction factor to include in our models, using a larger Environment and Climate Change Canada data set for the Oil Sands Area (>6,800 point counts). The correction factor adjusts raw counts for temporal variation in singing rate and differential detection rates in specific vegetation types, based on the effective area of detection (Matsuoka et al. 2012, Laake et al. 2015). We calculated the correction factor using multi‐covariate distance‐sampling models using the R package mrds (for further details, see Appendix S1), and adjusted the cluster‐level counts to a cluster‐level density by summing and log‐transforming the effective areas of detection for each of the nine point‐count sites within a survey cluster (Sólymos et al. 2013).

Habitat and disturbance data

We classified all land cover into one of nine classes, five natural landcover classes derived from habitat niche models developed for the same study area (Mahon et al. 2016), and four disturbance classes representing a dominant land‐use stressor: (1) white spruce/balsam fir (SWFB); (2) hardwood (combined trembling aspen, white birch, and balsam poplar; HDWD); (3) black spruce/larch lowlands (LWLD); (4) black spruce/pine uplands (CON2); (5) open non‐forested habitat (includes shrublands, grasslands, and recent burns <20 yr old; OPNF); (6) narrow, vegetated, linear features (percent cover of narrow seismic; NL); (7) wide linear features (percent cover of pipelines, powerlines, and all roads; WL); (8) harvest units (percent cover with forest harvest activity in the last 20 yr; HU); and (9) wells (percent cover of bitumen, oil, and gas wells; WE). Natural landcover classes were derived from Alberta Vegetation Inventory data (Alberta Sustainable Resource Development 2006), which is created by interpreting medium‐scale (1:60,000 or 1:40,000) aerial photographs. We updated this inventory data with forest harvest and wildfire data. After classifying vegetation by age class and dominant vegetation types, we used classification and regression tree analysis to derive the reduced natural landcover classes used in this analysis (Mahon et al. 2016). We derived disturbance landcover classes from 1:20,000 access and base feature data updated annually by the Government of Alberta. When stressors overlapped, we assigned the area to the highest priority ranking, where roads > hard‐surface industrial facilities > oil and gas wells > pipelines and powerlines > harvest units > seismic lines (Alberta Biodiversity Monitoring Institute 2012). This priority ranking of stressors accounts for the relative impact of individual stressors with respect to the degree of surface disturbance, surface substrate type, and infrastructure type. For example, we considered roads to be a higher priority disturbance because they typically have highly compacted asphalt or gravel surfaces, while seismic lines typically have light‐moderately compacted soil, wood chip, or vegetated surfaces. These categories were discrete, with no statistical collinearity among the nine landcover classes.

We quantified the amount of each natural and disturbance landcover class within the survey cluster of point counts and the 300‐m area around each point count (the sample unit in this analysis), because 300 m is the maximum detection radius for most boreal songbird vocalizations (Matsuoka et al. 2012). The total area sampled by each survey cluster (TAREA) differed because of slightly different configurations of the point‐count stations, so we rescaled the area of each landcover type to a proportion of the total area sampled. Similarly, we calculated area‐weighted stand age (AGE) within the area sampled by each survey cluster. We derived stand age for each natural landcover class using age data from the Alberta Vegetation Inventory (Alberta Sustainable Resource Development 2006), updated with recent forest harvest and wildfire data.

Cumulative effects models

We tested hypotheses for additive and interactive human‐disturbance relationships for 27 boreal landbirds using a two‐step modelling approach. First, we selected the best environmental model, which included habitat composition, spatial distribution and sampling variables. We used the environmental model to control for factors other than landscape disturbance that might influence landbird density. Second, we assessed whether the addition of disturbance variables improved model fit.

We considered 10 potential predictors within environmental models for each landbird species: percent cover of the five natural landcover classes (SWFB, HDWD, LWLD, CON2 and OPNF), AGE, latitude (LAT) and longitude (LONG) at the center of each survey cluster, and two variables that control for potential sampling differences (Julian date when survey was conducted, JDATE, and TAREA). An exploratory analysis indicated that nonlinear effects might be present for four variables (LAT, LONG, AGE, and HDWD), so we also considered polynomial terms for these variables. In addition to an environmental model with linear‐only effects, we tested models including polynomial terms for (1) LONG, (2) LAT, (3) AGE, (4) HDWD, (5) LONG +LAT, (6) LONG + AGE, (7) LONG + HDWD, (8) LAT + AGE, (9) LAT + HDWD, and (10) AGE +HDWD. All environmental models included the effective area of detection correction factor as a fixed offset.

We wanted to use a generalized additive model framework to develop environmental models, but needed to reduce the number of environmental predictors. Thus, we followed a two‐step procedure to identify the best environmental model for each species. First, in order to initially identify the most parsimonious model for each set of environmental predictors, we used a bidirectional step‐wise procedure based on Akaike's information criteria (AIC) within a generalized linear modeling framework. Generalized linear models used the negative binomial distribution with a logarithmic link function, and included quadratic terms for variables that were thought to have a nonlinear functional form. Second, once the most parsimonious model was identified for each set of variables, we incorporated this reduced set of variables into a generalized additive model, modelling any polynomial variables as a spline term with three degrees of freedom to ensure that the functional form was optimally fit. We used AIC weights to select the best generalized additive environmental model for each species.

We then assessed whether adding disturbance effects, in both additive and interactive frameworks, improved the fit of the generalized additive environmental model. Disturbance models capture existing patterns of resource development within the survey cluster of point counts and the 300‐m area around each point count (the sample unit in this analysis) and the larger AOSA study area (Fig. 2, Table 2) in six model groups: (1) the null model, which includes only habitat (the set of best fitting environmental variables); (2) linear disturbance, which tests the influence of narrow and wide linear stressors (seismic lines, pipelines, powerlines, roads) because of the potential for edge effects; (3) linear and energy disturbance, which tests the influence of linear and energy (polygonal well) stressors; (4) linear and forestry disturbance, which tests the influence of linear and forestry (polygonal harvest unit) stressors; (5) cumulative disturbance (linear, energy, forestry); and (6) disturbance area only (total area of all disturbances). We combined the percent cover of narrow and wide linear disturbances to create a new variable, all linear (AL), for model groups 4–6. Within each model group, we systematically tested both additive and interactive generalized additive models to determine the strength of influence of interaction terms between disturbance variables (Table 2). We used multi‐model inference based on AIC, combined with log‐likelihood tests, to select the best of the 12 candidate disturbance models for each species and ranked models using AIC weights. We used log‐likelihood‐ratio tests to test (1) whether the best‐ranked disturbance model was significantly different from the null model (at P = 0.1, given complexity of data and the power required to detect fit) and (2) whether the best‐ranked model differed from the second‐ and third‐ranked models. When the top two models were additive and interactive models within the same model set, and the interactive model was not a better fit as determined by a log‐likelihood test, we selected the more parsimonious additive model. We used mvrnorm in the R MASS package to calculate the 10% and 90% confidence intervals, based on 5,000 simulations (Luis et al. 2015, Ripley et al. 2015).

Table 2

Candidate disturbance model set tested in the additive and interactive analysis

Model group and set	Model type	Model predictors
Habitat only (HO)
1a	N/A	Habitat and Spatial Variables Only
Linear (LI)
2a	additive	NL + WL
2b	interactive	NL + WL + NL × WL
Linear and energy (LE)
3a	additive	NL + WL + WE
3b	interactive	NL + WL + WE + WE × NL + WE × WL
3c	interactive	NL + WL + WE + WE × NL + WE × WL + WE × NL × WL
Linear and forestry (LF)
4a	additive	AL + HU
4b	interactive	AL + HU + HU × AL
Cumulative (CU)
5a	additive	AL + WE + HU
5b	interactive	AL + WE + HU + HU × AL + HU × WE
5c	interactive	AL + WE + HU + HU × AL + HU × WE + WE × AL + HU × WE × AL
Disturbance area (DA)
6a	additive	all disturbances summed (combined AL, WE, HU)

All models include habitat and spatial predictors selected for each landbird species. Stressor definitions: NL, narrow linear (percent cover of seismic lines); WL, wide linear (percent cover of pipelines, powerlines, all roads); AL, all linear (percent cover of narrow and wide linear combined); HU, harvest unit (percent cover harvest units with harvest activity <20 yr old); WE, wells (percent cover with bitumen, oil, and gas wells).

All models were robust and stable, based on examination of residual diagnostic plots and stability of individual predictor terms (ratios of predictor coefficients to their standard error and odds ratios). Chi‐square tests of deviance confirmed that the negative binomial distribution was suitable for the data. All models fit well based on proportion of deviance explained.

Estimating impacts to bird populations

To evaluate how disturbance and habitat loss have impacted landbird densities in the AOSA, we compared predicted densities for each of the 27 species for the current landcover layer based on the observed human disturbance within the study area and for a no‐disturbance landcover layer that theoretically represents the current study area if there was no human disturbance (sensu Nielsen et al. 2007). Our theoretical no‐disturbance landscape was primarily based on forest inventory data from 1990 onward, the period of time in which most resource development within this region has occurred. Because natural and disturbance landcover types were derived from separate data layers, in the vast majority of habitat polygons we were able to “reset” the disturbed habitat polygon to its condition in the original Alberta Vegetation Inventory data (Alberta Sustainable Resource Development 2006). In polygons where the original stand composition was not available, the disturbance was reassigned to a stand composition type based on adjacent land cover. Given the high proportion of the study area with known “historical” condition, we expect this approach to be far more accurate than simulation models, which require numerous assumptions. We used the best disturbance model for each species to calculate mean predicted densities and standard errors within each landscape scenario (current, no‐disturbance), using the 5,000 mvrnorm simulations.

To compare the similarity of the bird communities between the two study areas (current, no‐disturbance), we used multivariate Procrustes analysis (following Jackson 1995). We first created data matrices containing the predicted densities for the 27 species in 303 survey cluster sample units. We then created nonmetric dimensional scaling (NMDS) ordination plots for the current and no‐disturbance study areas (based on the Bray–Curtis dissimilarity measure) and used the Procrustes test to assess the degree of concordance between the two study areas. The Procrustes test compares the two community ordinations and produces a numerical measure of concordance scaled between 0 and 1 (1 is complete similarity). We bootstrapped the Procrustes analysis 5,000 times to calculate mean values and 90% confidence intervals. We conducted all analyses and data summaries in program R (Version 3.2.1), using the packages MASS, mgcv, lmtest, vegan, and ggplot2.

Results

Of the 27 species analyzed with the candidate disturbance model set, the null model (model 1a) was selected as the best model for seven species: Black‐and‐white Warbler (Mniotilta varia), Cape May Warbler (Setophaga tigrina), Least Flycatcher (Empidonax minimus), Magnolia Warbler (S. magnolia), Palm Warbler (S. palmarum), Ruby‐crowned Kinglet (Regulus calendula), and Winter Wren (Troglodytes hiemalis) (Table 3). The inclusion of disturbance effects improved models for the remaining 20 species; additive models were the best for nine species, while interactive models were the best for 11 species (Table 3). A single best additive disturbance model could not be selected for Boreal Chickadee (Poecile hudsonicus) or Chipping Sparrow (Spizella passerina), where the two best models were 4a (linear and forestry‐additive) and 6a (disturbance area only, additive). These were the only species to demonstrate strong selection for the disturbance area model.

Table 3

Model results for the 12 candidate disturbance models tested for 27 landbird species (sorted by habitat grouping)

Species name	Species code	Habitat	Best group and model	AIC weight (w i)	Explained deviance	Best model type	Key disturbance predictors
Bay‐breasted Warbler	BBWA	CO	LE—3a	0.32	0.18	additive	WL(+), WE(−)
Cape May Warbler	CMWA	CO	HO—1a	0.12	0.42	null	N/A
Red‐breasted Nuthatch	RBNU	CO	LI—2b	0.38	0.27	interactive	NL(−), WL(+)
Winter Wren	WIWR	CO	HO—1a	0.40	0.36	null	–
Boreal Chickadee	BOCH	SB	LF/DA—4a/6c	0.32/0.28	0.17	additive	HU(−), AL(−)/AD(−)
Chipping Sparrow	CHSP	SB	LF/DA—4a/6c	0.29/0.29	0.18	additive	HU(−), AL(−)/AD(−)
Canada Jay	CAJA	SB	LE—3a	0.66	0.37	additive	NL(−), WE(+)
Dark‐eyed Junco	DEJU	SB	LF—4b	0.45	0.42	interactive	HU(−), AL(+)
Palm Warbler	PAWA	SB	HO—1a	0.35	0.57	null	N/A
Ruby‐crowned Kinglet	RCKI	SB	HO—1a	0.46	0.62	null	N/A
Black‐and‐white Warbler	BAWW	DE	HO—1a	0.34	0.31	null	N/A
Black‐capped Chickadee	BCCH	DE	LE—3c	0.57	0.32	interactive	NL(−), WE(+), WL(+)
Least Flycatcher	LEFL	DE	HO—1a	0.23	0.26	null	N/A
Mourning Warbler	MOWA	DE	LI—2a	0.25	0.40	additive	NL(−), WL(+)
Northern Flicker	NOFL	DE	LF—4a	0.29	0.22	additive	HU(+), AL(+)
Ovenbird	OVEN	DE	LE—3b	0.77	0.67	interactive	NL(−), WL(+)
Red‐eyed Vireo	REVI	DE	CU—5c	0.77	0.52	interactive	WE(+), HU(+), AL(+)
Rose‐breasted Grosbeak	RBGR	DE	LE—3b	0.73	0.42	interactive	WE(+), NL(−)
Alder Flycatcher	ALFL	OP	LF—4a	0.38	0.33	additive	HU(+)
American Robin	AMRO	OP	LF—4b	0.44	0.21	interactive	HU(+), AL(+)
Hermit Thrush	HETH	OP	CU—5b	0.36	0.36	interactive	HU(−), AL(+)
Lincoln's Sparrow	LISP	OP	CU—5b	0.25	0.24	interactive	HU(+), AL(−), WE(+)
Yellow‐bellied Sapsucker	YBSA	OP	LE—3b	0.69	0.28	interactive	NL(−), WE(+)
Common Raven	CORA	GE	LE—3c	0.91	0.23	interactive	WE(+)
Magnolia Warbler	MAWA	GE	HO—1a	0.10	0.12	null	N/A
White‐throated Sparrow	WTSP	GE	LF—4a	0.53	0.26	additive	HU(+), AL(−)
Yellow‐rumped Warbler	YRWA	GE	LF—4a	0.43	0.37	additive	HU(−), AL(+)

The best model(s) and model type are shown for each species where a best disturbance model was selected. The most important predictors (those with the strongest impact on landbird density) in each model are shown, where + indicates a positive (increasing) response of landbird density to stressors, and − indicates a negative (decreasing) response (N/A indicates habitat‐only models lacking disturbance predictors). Species shown in boldface type are niche specialists (following Mahon et al. 2016). Habitat groupings (following Mahon et al. 2016): CO, conifer (white spruce, balsam fir); DE, deciduous (trembling aspen, white birch, balsam poplar); GE, generalist; OP, open; SB, lowland and upland black spruce, Model group definitions: HO, habitat only; LI, linear; LE, linear and energy; LF, linear and forestry; CU, cumulative; DA, disturbance area. Stressor definitions: NL, narrow linear (percent cover of seismic lines); WL, wide linear (percent cover of pipelines, powerlines, all roads); AL, all linear (percent cover of narrow and wide linear combined); HU, harvest unit (percent cover forest harvest activity in the last 20 yr); WE, wells (percent cover with bitumen, oil, and gas wells). AIC, Akaike information criterion.

Interactive models were best supported for: species associated with deciduous forest (Black‐capped Chickadee [Poecile atricapillus], Ovenbird [Seiurus aurocapilla], Red‐eyed Vireo [Vireo olivaceus], Rose‐breasted Grosbeak [Pheucticus ludovicianus]); species associated with open habitat (American Robin [Turdus migratorius], Hermit Thrush [Catharus guttatus], Lincoln's Sparrow [Melospiza lincolnii], Yellow‐bellied Sapsucker [Sphyrapicus varius]); species associated with conifer forest (Red‐breasted Nuthatch [Sitta canadensis]); species typically found in upland and lowland black spruce (Dark‐eyed Junco [Junco hyemalis]); and generalist species (Common Raven [Corvus corax]). The linear and energy model group contained the best interactive model for five species (Black‐capped Chickadee, Common Raven, Ovenbird, Rose‐breasted Grosbeak, Yellow‐bellied Sapsucker). The cumulative model group contained the best interactive model for three species (Hermit Thrush, Lincoln's Sparrow, Red‐eyed Vireo), the linear and forestry model group contained the best interactive model for two species (American Robin, Dark‐eyed Junco), and the linear model group contained the best interactive model for one species (Red‐breasted Nuthatch).

Landbird species demonstrated complex, and often differing, responses to individual stressors (Figs. 3 and 4, Table 3; see Appendix S2: Table S1 for full model coefficients). Conifer‐associated species (spruce, fir, pine, black spruce) exhibited consistently negative responses to three stressors (narrow linear, harvest units, all disturbance) and positive responses to only one stressor (wide linear). Deciduous‐associated species exhibited negative responses to only one stressor (narrow linear) and positive responses to four stressors (wells, harvest units, wide linear, all linear). Open and generalist‐associated species exhibited positive responses only to wells. Conifer species exhibited mixed responses (e.g., positive and negative) to wells and all linear, while open and generalist species exhibited mixed responses to harvest units and all linear.

Figure 3

Predicted densities (±90% confidence envelopes) for each stressor based on the best disturbance model(s) for species that had additive cumulative effects models: (A) Alder Flycatcher (ALFL), (B) Bay‐breasted Warbler (BBWA), (C) Boreal Chickadee (BOCH), (D) Canada Jay (CAJA), (E) Chipping Sparrow (CHSP), (F) Mourning Warbler (MOWA), (G) Northern Flicker (NOFL), (H) White‐throated Sparrow (WTSP), and (G) Yellow‐rumped warbler (YRWA). Results for the two best models are presented simultaneously in the same graph for Boreal Chickadee and Chipping Sparrow. The x‐axis is the range of stressor percent cover values observed in this study. The predicted densities for each stressor are calculated while holding all other model predictors constant at their mean values.

Figure 4

Predicted densities for each stressor based on the best disturbance model(s) for species that had interactive cumulative effects models: (A) American Robin (AMRO), (B) Black‐capped Chickadee (BCCH), (C) Common Raven (CORA), (D) Dark‐eyed Junco (DEJU), (E) Hermit thrush (HETH), (F) Lincoln's Sparrow (LISP), (G) Ovenbird (OVEN), (H) Rose‐breasted Grosbeak (RBGR), (I) Red‐breasted Nuthatch (RBNU), (J) Red‐eyed Vireo (REVI), and (K) Yellow‐bellied Sapsucker (YBSA). The top row shows the best interactive cumulative effects models with ±90% confidence envelopes, while the bottom row shows a comparison of interactive models (solid line) with complementary additive models (dashed line). The x‐axis is the range of stressor percent cover values observed in this study. The predicted densities for each stressor are calculated while holding all other model predictors constant at their mean values.

Species associated with open habitat or deciduous forest such as Alder Flycatcher (Empidonax alnorum), American Robin and Red‐eyed Vireo, and generalists such as Common Raven and Northern Flicker (Colaptes auratus), had higher densities where disturbance feature area was higher. In contrast, two conifer‐associated species, Boreal Chickadee and Chipping Sparrow, had lower densities where disturbance feature area was higher. However, many species showed opposite responses to different types of stressors. For example, Yellow‐rumped Warbler (Setophaga coronata), Dark‐eyed Junco and Hermit Thrush densities showed negative responses to forest harvesting and positive responses to linear features, in contrast to the responses of the early‐seral associated species, White‐throated Sparrow (Zonotrichia albicolis) and Lincoln's Sparrow. For species where the best‐fitting model included wells, the responses to wells were generally positive or neutral for most species, with the exception of Bay‐breasted Warbler (Dendroica castanea).

For species where interactive models provided a better fit, the interpretation of stressor effects depends on the magnitude of the other stressors (Fig. 4). As an example, Common Raven (Fig. 4C) increased in density in response to wells to a substantially greater degree when linear features were present. In most cases, interactive effects changed only the magnitude of effects observed, not the direction of response; however, the response of Ovenbirds to wells changed from negative to neutral when we included the interactive effects of linear features, and the response of Hermit Thrush to wells changed from negative to positive when we included the interactive effects of linear features and harvest units. We observed both synergistic and antagonistic responses to individual stressors. Six species exhibited all synergistic interactions: American Robin, Black‐capped Chickadee, Common Raven, Ovenbird, Red‐eyed Vireo, and Rose‐breasted Grosbeak. For these species, the combination of individual stressors within interactive models resulted in more extreme responses to stressors than would have been estimated in additive models. Antagonistic interactions were observed for Dark‐eyed Junco and Lincoln's Sparrow. Hermit Thrush showed mild antagonistic responses for two stressors (linear features and harvest units), but a synergistic response for wells. The remaining two species, Red‐breasted Nuthatch and Yellow‐bellied Sapsucker, showed only minor differences between interactive and additive cumulative‐effects models.

Finally, the mean model predictions comparing the current and theoretical no‐disturbance study areas predicted substantial (>20%) changes in density for 12 of 27 species (Fig. 5). Several species associated with deciduous habitats (Black‐capped Chickadee [+150%], Red‐eyed Vireo [+219%]) and open habitats (Alder Flycatcher [+110%], American Robin [+183%]) showed large increases in density in the current study area relative to the no‐disturbance study area (Fig. 5A). In contrast, species associated with conifer habitats (Red‐breasted Nuthatch [−41%]) and black spruce habitats (Boreal Chickadee [−24%], Canada Jay [Perisoreus canadensis; −27%]) showed decreases in the current study area relative to the no‐disturbance study area (Fig. 5B). An additional four species showed moderate increases in the current study area relative to the no‐disturbance study area, Common Raven (+24%), Dark‐eyed Junco (+25%), Hermit Thrush (+37%), and Rose‐breasted Grosbeak (+34%; Fig. 5). When comparing the entire community of 27 landbird species, the observed similarity between the current and theoretical no‐disturbance study area was 0.85 (Procrustes similarity index; 0.845–0.890 90% CI).

Figure 5

Predicted densities (mean ± SE) for species that demonstrated >20% change between the theoretical no‐disturbance study area and actual current study area where a best disturbance model(s) was selected. The trends in density are identified as + for species that are increasing in the current study area relative to the theoretical no‐disturbance study area and − for species that are decreasing in the current study area relative to the theoretical no‐disturbance study area. Deciduous and open habitat species are shown in panel A, while conifer, black spruce, and generalist species are shown in panel B. Species codes are defined in Table 3.

Discussion

This is the first study to evaluate both additive and interactive effects of human disturbances within the western boreal forest, a region undergoing rapid intensification of development by multiple resource industries. Our results provide rigorous quantitative evidence of the presence of complex, interactive cumulative effects on landbirds. Although the combinations of habitats and stressors analyzed are specific to the western boreal forest, coexisting resource sectors occur throughout many geographic regions of North America and Europe. Understanding how multiple stressors cumulatively impact wildlife, including consideration of interactions among stressors, is needed to develop effective mitigation and management strategies. Our work represents a repeatable framework to identify, model, and quantify both additive and interactive cumulative effects on terrestrial wildlife.

Past research has focused on examining changes in species abundance in response to individual stressors, including harvest units (Schmiegelow et al. 1997, Schieck and Song 2006), seismic lines (Bayne et al. 2005, Machtans 2006, Lankau et al. 2013), wells (Thomas et al. 2014, Bayne et al. 2016) and linear features (Ball et al. 2009, Bayne et al. 2016). During early stages of resource development in the western boreal forest, there was relatively little overlap among these stressors in space and time. However, in the current multi‐use landscape, industrial forestry, energy exploration and extraction, and infrastructure development have created an extensive network of overlapping and adjacent stressors, which necessitates study designs and modelling approaches that can assess the cumulative effects of combinations of multiple stressors. Cumulative‐effects studies must be designed, implemented, and analyzed at spatial scales that appropriately capture the non‐uniform distribution of multiple stressors across the gradient of development within the study area, because the scale at which analyses are conducted can influence our ability to detect landbird response to stressors (Howe et al. 2014).

Our empirical models provide evidence of both additive and interactive effects across the suite of boreal landbird species examined. The addition of disturbance effects significantly improved models for 20 of the 27 landbird species considered (74%). We used habitat and niche relationships defined for landbirds within the same study area (Mahon et al. 2016) to group species into regionally relevant habitat guilds and identify specialist and generalist species. Our results suggest that conifer species like Bay‐breasted Warbler (a niche specialist), Red‐breasted Nuthatch, Boreal Chickadee, Chipping Sparrow, and Canada Jay are being negatively impacted by stressors within this landscape. Harvest units remove suitable breeding habitat, while narrow and wide linear features may affect the placement or location of breeding territories. In contrast, deciduous species (niche generalists with the exception of Ovenbird) appear to benefit from novel human stressors (wells, harvest units, wide linear, all linear), most likely due to patterns of post‐disturbance regeneration in many upland habitat types. Pioneer species of graminoids, forbs, and shrubs, and subsequent colonizers such as shade‐intolerant shrubs and regenerating alder (Alnus sp.), willow (Salix sp.), and poplar species dominate upland sites after disturbance. In some cases, specific stressors may provide new high‐quality breeding season habitats. Exploration wells may provide high‐quality foraging habitat for some boreal landbirds, although quality may depend on well size, moisture level, ground substrate, vegetation growth, and the adjacent habitat. In this region, we have observed wells being used as communal or shared foraging areas by species including thrushes (American Robin, Hermit Thrush), warblers (Yellow‐rumped Warbler), and sparrows (Chipping Sparrow, Dark‐eyed Junco, White‐crowned Sparrow) (C. L. Mahon, unpublished data). These findings support our prediction that both direct and indirect effects of resource development are affecting landbird populations in the region. These are noteworthy results, given that the analysis included species with diverse habitat associations (upland and lowland forests, bog, fen, swamp, marsh, shrubland, grassland, and burn; Mahon et al. 2016) and represented a landscape where the average total area of resource stressors is low (i.e., <10%). We detected significant changes in the bird community in a study area with 8.4% human disturbance; however, a 2015 estimate found 18.1% total human disturbance in Alberta's boreal forest (Alberta Biodiversity Monitoring Institute 2017), so the total impacts on the boreal bird community may be greater than demonstrated here.

We also found support for complex interactive effects (synergistic and antagonistic) in 11 landbird species (40% of the species assessed) including deciduous, conifer, open habitat and generalist species (for a complete list see Table 3 and Appendix S2: Table S1). It is noteworthy that most interactive effects detected were synergistic, rather than antagonistic. In some cases, the responses we observed were somewhat non‐intuitive and additional research will be required to confirm the validity of the interactions observed and investigate the mechanisms underlying these species‐specific responses; for example, some species responded more negatively to narrow linear features than to wide linear features (e.g., Bay‐breasted Warbler, Black‐capped Chickadee, Mourning Warbler Oporornis philadelphia). We note several limitations to our approach to modelling and quantifying additive and interactive cumulative effects on boreal landbirds. First, at the spatial scale of this study, we were unable to account for factors that could potentially interact with the variables assessed. For example, we know that the state of vegetative regeneration influences use of these features by landbirds (Lankau et al. 2013, Wilson 2017), but this information was unavailable for the vegetated features in our study, such as exploratory wells, pipelines, and seismic lines. Similarly, the response of landbirds to specific stressors may vary among habitat types (Bayne et al. 2016). Second, we were unable to examine species‐specific mechanisms that could explain species responses to complex synergistic and antagonistic interactions among stressors. Third, we were unable to include all resource stressors within the study area in our modelling framework. Industrial sites and oil sands mines were not surveyed, primarily due to health and safety risks to field staff, and as a result, our understanding of these stressors and how they interact with other stressors on the landscape to impact bird populations is unknown.

Despite these caveats, our work suggests that boreal landbirds in the AOSA are responding to the removal and degradation of native vegetation in ways that cannot be predicted by simple additive models. Accounting for these responses will require consideration of all disturbance types that act as stressors, not just those stressors that have large disturbance areas. We show quantitatively that stressors such as linear features, which have relatively small disturbance areas, interact with other resource stressors to affect the stressor response. Further, interactions that include wells or linear features demonstrate that the continued perforation and dissection of the boreal forest by small polygonal disturbances (e.g., wells, small industrial sites) and narrow linear features (e.g., seismic lines, pipelines) is changing habitat quantity and quality for breeding landbirds either directly (through removal of native vegetation) or indirectly (through changes in food and prey, nest and foraging sites, predators and parasites). Additional studies of both generalist and specialist species using diverse and innovative methods (e.g., spot mapping, radio telemetry, behavioral observations) and metrics associated with habitat quality (e.g., reproductive measures; space use; territory, nest and forage area selection; foraging rate and distance; prey availability and selection; predator density; dispersal distance) may help explain changes in ecological and behavioral processes such as reproduction, habitat selection (i.e., territory, nest site and foraging areas), and predator and parasite dynamics.

We predicted an overall 15% change in the landbird community between the current study area and one with no disturbances, nearly twice the area of disturbance (8.4%) in the accessible portion of the AOSA (7.4 million ha). The predicted community change includes substantial increases in density for many deciduous (Black‐capped Chickadee, Red‐eyed Vireo) and open habitat (Alder Flycatcher, American Robin) landbird species. In contrast, declines in density are predicted for conifer (Red‐breasted Nuthatch) and black spruce‐associated (Boreal Chickadee, Canada Jay) species. These community changes reflect the general patterns of response we detected, whereby deciduous species exhibited positive responses to both polygonal and linear stressors and conifer species exhibited generally negative responses to stressors. Our results suggest that the “winning species” (increasing in the current study area, sensu McKinney and Lockwood 1999) are primarily niche generalists within the western boreal (Mahon et al. 2016). These results match current land‐use patterns in the western boreal forest, where mixedwood and conifer forests are being replaced with younger deciduous forests as a result of industrial forestry operations (Hobson and Bayne 2000), and large areas of black spruce‐associated habitats are being subdivided by oil, gas, and bitumen development (Schneider et al. 2003, Lee and Boutin 2006). In this type of resource‐affected landscape, open or early seral species are predicted to increase in abundance (Hobson and Schieck 1999, Davidson and Knight 2001).

In other disturbed areas of the world, the most specialized species typically exhibit the strongest response to landscape disturbance (Julliard et al. 2004, Devictor et al. 2008). Here, we found instead that the strongest responses to disturbance were for generalist deciduous and open habitat species (Black‐capped Chickadee, Red‐eyed Vireo, Alder Flycatcher, American Robin), consistent with a previous study in this region that found that the largest changes are for species positively associated with human disturbances (Sólymos et al. 2014). Within the western boreal forest, disturbances create habitats for these species in areas where they were previously uncommon or absent. Our results suggest that, even at relatively low levels of human disturbance, we are starting to see evidence of generalists and early seral associates increasing, and specialists and conifer associates decreasing (McKinney and Lockwood 1999, Clavel et al. 2011). These changes in community composition are consistent with those expected under the hypothesis that the western boreal forest is undergoing biotic homogenization with increased resource development (Mahon et al. 2016), but do not provide direct evidence for the hypothesis.

If “losing species” are specialists within the boreal landbird community, there could be implications for ecosystem function and resilience across boreal regions of Canada and Europe. Two examples using specialist species found in western boreal forests highlight the risk to ecosystem biodiversity. First, strong and weak cavity nesters (e.g., Black‐backed Woodpecker, Red‐breasted Nuthatch) are key excavators within interconnected nest webs (Martin and Eadie 1999). The decline or loss of species that excavate nest cavities for a variety of secondary cavity nesters like bluebirds, swallows, and squirrels will disrupt the composition (Martin et al. 2004) and function (Blanc and Walters 2007, Cockle and Martin 2015) of the nest‐web community. Second, multiple species are recognized for their demographic and functional responses to insect outbreaks that provide a short‐term resource pulse for birds that can prey on larvae or emerging adults in hemiboreal and boreal forests. Wood warblers like Bay‐breasted and Cape May Warbler respond to spruce budworm (Crawford and Jennings 1989, Patten and Burger 1998, Venier et al. 2009), weak cavity nesters like Red‐breasted Nuthatch respond to mountain pine beetle (Norris and Martin 2008, 2010, 2014), and strong cavity nesters like American Three‐toed Woodpecker respond to spruce beetle (Fayt et al. 2005). These specialized members of the landbird community can contribute to the resilience and productivity of the boreal forest by acting as highly mobile predators that converge on areas of increasing insect prey abundance, increasing overall diversity and food‐web complexity (Eveleigh et al. 2007).

Historically, wildfire was the primary disturbance modifying boreal habitats and landscapes, but the present cumulative disturbance patterns share few similarities with the disturbance and regeneration patterns created by large and frequent wildfires (Stocks et al. 2002, Parisien et al. 2011). Currently, all habitat types are being perforated and dissected by stressors associated with multiple sectors that target both above‐ and belowground resources. These stressors are diverse in terms of size (1‐ha wells to 1,000‐ha aggregated harvest units), shape (polygonal and linear), frequency (single or repeated visits), and degree of site disturbance (surface to bedrock). It may be productive to quantitatively compare current human disturbances with natural disturbances, to determine the implications of these new disturbances on the ecosystem (Roberts 2004, 2007). Moreover, simulation models suggest that the impacts of climate change could further alter forests in this region over the coming decades (Stralberg et al. 2018). The magnitude and rate at which increases in climate‐driven disturbances (primarily wildfire, but also disease and insects) will interact with human disturbances to alter the landbird community is currently unknown, but climate change is projected to reduce the availability of upland conifer and mixedwood forests, and to increase the availability of grasslands, in the coming decades (Stralberg et al. 2018). As such, climate change may further increase the habitat available for generalist deciduous and open habitat species, and decrease habitat availability for conifer species, thereby interacting synergistically with the impacts of human disturbances observed in this study. The interacting cumulative effects of diverse local‐scale stressors (e.g., resource stressors) and large‐scale climate‐driven disturbances will make it increasingly difficult for land managers to meet ecological objectives.

Management Implications

Formulating effective management decisions within an ecological context is frequently difficult due to the inherent complexity of ecological interactions; this study empirically demonstrates that challenge. At this time, it is premature to provide prescriptive advice on how to optimize habitat management and land‐use planning to minimize impacts on boreal bird species, in part because specific recommendations cannot be adequately formulated in the absence of defined regional or sub‐regional management objectives. Moreover, in light of the complexity we demonstrated in this study, prescriptive approaches risk recommending solutions that are inappropriate for some situations. Nevertheless, some broad recommendations are well supported by this study and the growing body of evidence that human disturbances impact bird communities. The time required for vegetation on seismic lines, wells, and other energy‐sector stressors to regenerate to the point where there are no longer significant impacts on bird communities is not well understood, although available results suggest that moderate‐tall shrub growth is required for impacts to begin to ameliorate (Lankau et al. 2013, Wilson 2017; J. D. Toms, unpublished data). Vegetative recovery from stressors to structural stages where pre‐disturbance bird communities return can take many decades (Lee and Boutin 2006, Schieck and Song 2006, van Rensen et al. 2015, Wilson 2017) and, currently, it appears that the rate of stressor creation exceeds stressor recovery within the Oil Sands Area. If management goals are to maintain all boreal bird species within the historical range of natural variation, resource development may need to be limited to some degree until regeneration and restoration of existing human disturbances are sufficient to prevent further impacts on species of concern. Similar recommendations have been made for management of woodland caribou (Rangifer tarandus caribou) in Alberta (Athabasca Landscape Team 2009), suggesting that implementation of these recommendations will have broader biodiversity benefits (as also seen in other studies; e.g., Bichet et al. 2016). The most effective approach would appear to be one that is holistic, ecosystem‐based, and focused on management of cumulative effects at regional scales.

Our study indicates that regional land‐use planning and environmental assessments need to consider interactions among stressors to better predict the cumulative effects of multiple industries on a single land base. Furthermore, we need to understand whether local‐scale stressors and large‐scale climate‐driven disturbances act additively, synergistically, or antagonistically; local mitigation and management efforts can be highly effective when local impacts interact synergistically with large‐scale impacts, but may be counter‐productive if they interact antagonistically with large‐scale impacts (Brown et al. 2013). For example, the targeted protection of conifer forests across the landscape, particularly those upland conifer forests most likely to persist under climate change (Stralberg et al. 2015), may be warranted, given that species associated with these forest types appear to be disproportionally negatively affected by both resource stressors (this study) and climate‐driven disturbances (Stralberg et al. 2018). The results of this study are preliminary and correlative, not conclusive, but suggest that synergistic and antagonistic interactions are common among the stressors present in this region. Further studies will be needed before we can fully account for this complexity in land use planning and environmental assessments.

Most importantly, this study indicates that impact assessments, and associated land‐use planning decisions, should be precautionary due to the potential for unanticipated synergistic interactions. Effective, integrated resource management, e.g., increased coordination of planning and operations among sectors, may be able to reduce the impacts of resource stressors (Canter and Ross 2010). In general, the effective mitigation and management of cumulative effects will require coordinated regional land‐use plans that incorporate evaluation of complex cumulative effects and climate change to ensure that development is consistent with management objectives.

Click here for additional data file.

Click here for additional data file.