Threats to biodiversity from cumulative human impacts in one of North America's last wildlife frontiers

Land‐use change is the largest proximate threat to biodiversity yet remains one of the most complex to manage. In British Columbia (BC), where large mammals roam extensive tracts of intact habitat, continued land‐use development is of global concern. Extant mammal diversity in BC is unrivalled in North America owing, in part, to its unique position at the intersection of alpine, boreal, and temperate biomes. Despite high conservation values, understanding of cumulative ecological impacts from human development is limited. Using cumulative‐effects‐assessment (CEA) methods, we assessed the current human footprint over 16 regional ecosystems and 7 large mammal species. Using historical and current range estimates of the mammals, we investigated impacts of human land use on species’ persistence. For ecosystems, we found that bunchgrass, coastal Douglas fir, and ponderosa pine have been subjected to over 50% land‐use conversion, and over 85% of their spatial extent has undergone either direct or estimated indirect impacts. Of the mammals we considered, wolves were the least affected by land conversion, yet all species had reduced ranges compared with historical estimates. We found evidence of a hard trade‐off between development and conservation, most clearly for mammals with large distributions and ecosystems with high levels of conversion. Rather than serve as a platform to monitor species decline, we strongly advocate these data be used to inform land‐use planning and to assess current conservation efforts. More generally, CEAs offer a robust tool to inform wildlife and habitat conservation at scale.


Introduction
Conservation aims to prevent species and ecosystem loss (Soulé 1985) while still managing human uses of environmental resources (Kareiva & Marvier 2012). Yet land-use change is the single largest threat to biodiversity (MEA 2005;Turner et al. 2007) and planning efforts to manage it have failed to slow development or resulting biodiversity impacts (Butchart et al. 2010;Newbold et al. 2016).
The variety of land-use activities, from agriculture to land clearing for settlement or resource extraction (Foley et al. 2005), make tracking and managing cumulative use challenging (Raiter et al. 2014). Systematic conservation planning has emerged to address this challenge (Margules & Pressey 2000) and frameworks like cumulative effects assessment (CEA) (Halpern & Fujita 2013) are gaining traction. Cumulative effects assessments contextualize local development in a regional setting and are used to assess large-scale land-use impacts to inform small-scale planning (Baxter et al. 2001). Typically, CEAs have three primary steps focused around predefined ecological values (Spaling & Smit 1993). The first is quantifying the total regional human footprint. The chosen ecological values determine the spatial boundaries of the assessment (Therivel & Ross 2007). Thus, footprints may shift depending on species' ranges or ecosystem distributions. The second step is estimating the impact of that foot-print on ecological values. Estimating impacts is based on quantitative predictions that are refined by monitoring ecological values through time (Burton et al. 2014). The final step of a CEA is outlining future development scenarios. Using calculated footprints, estimated impacts, and future scenarios, CEAs can inform strategies that minimize risks to ecological values. Cumulative effects assessments can be used to manage at multiple scales and over many land uses, which protects conservation values while allowing sustainable development (Duinker & Greig 2006). In practice, examples of comprehensive CEAs are rare, even where they are increasingly needed.
British Columbia (BC) represents an area of high global conservation value, yet it has undergone little provinciallevel CEA and planning. Habitat diversity in BC is high; elevations range from 0 to >4000 m and climate regimes range from the very wet hypermaritime to the semiarid grasslands (Meidinger & Pojar 1991). In continental North America, range contractions of over 20% have occurred for seventeen mammals since Euro-American settlement. British Columbia plays a prominent role in habitat provision for these dwindling populations (Laliberte & Ripple 2004) because it contains large tracts of globally significant untouched habitat. Land use in BC is recent because much of the natural resource base is remote and inaccessible. Pressures on the landscape are increasing as technology opens previously inaccessible areas, and terrestrial species populations are declining across the province (BC Ministry of Environment 2014). There is significant economic reliance on natural resources, especially natural gas and lumber (BC Ministry of Finance 2016), and agriculture is prominent in the central, south, and northeast regions. This situation creates a pressing need for comprehensive land-use planning.
Cumulative effects assessments start with mapping the human footprint (Connelly 2011), often represented by the spatial extent of land-use (Toews 2016). We mapped the current footprint provincially and focused on its distribution across ecosystems and select mammal ranges. Given BC's accessibility issues and the spatial distribution of resources, we expected particular land-use types to be isolated within certain ecosystems. However, some development types such as roads are likely diffuse and thus affect all provincial ecosystem types. When narrowed to individual species' ranges, the footprint will likely shift, and wide-ranging species such as large carnivores will be the most affected.
Once mapped, the next step is estimating land-use impacts on ecological values. We used historic range estimates of mammal species to investigate local species extirpations based on individual activities, cumulative effects, and indirect effects. Range boundaries are notoriously coarse (Hurlbert & Jetz 2007), and comparing to historical estimates is difficult (Tingley & Beissinger 2009). We acknowledge these pitfalls; yet, range estimates provide a critical foundation for future species monitoring and for monitoring the effectiveness of land-use planning for conservation outcomes, of which species persistence is a key performance indicator.
By creating quantitative models of land-use relationships with range loss, we are building foundations for refined, predictive knowledge of land-use impacts on these mammals. In the final step of the CEA process, future scenarios are outlined and paired with predicted ecological impacts (Smit & Spaling 1995). From this, recommendations on sustainable development can be made. We did not extend our study into scenario predictions. Assessing cumulative impacts for multiple species is essential for understanding trade-offs and for identifying both idiosyncratic and consistent impacts. We sought to take the first major steps toward a comprehensive CEA in a globally important region by calculating the total human footprint for a set of ecological values; assessing the status of those values based on land-use cover and range loss; and creating preliminary predictive models of land-use impacts on those values.

Study Area
British Columbia covers 945,000 km 2 and is the westernmost province of Canada. Vegetation communities have been comprehensively described and classified by the Ministry of Forests (Pojar et al. 1987) using a system known as the Biogeoclimatic Ecosystem Classification (BEC). The BEC zones are determined primarily based on climate, vegetation, and soil data (Pojar et al. 1987). Originally established to map forest types and commercial tree occurrence, BEC zones contain different levels of extractable resources and biodiversity. There are 16 BEC zones in BC, and each zone was treated as an individual ecosystem. Wildlife use of zones tends to have high overlap; only 12% of terrestrial vertebrate species are thought to be zone specific (Bunnell 1995). Thus, ecosystem analysis based on the BEC zones captures climate, soils, and vegetation rather than habitat and range size of individual species.

Data Collection
We collected spatial data on land use and mammal range estimates (historic and current). For land-use information, we used data publicly available through GeoBC (http://geobc.gov.bc.ca/), a subset of the BC Integrated Resource Operations Division that oversees baseline spatial data and Provincial Crown Registries on land development. Details of all impact shapefiles is in Supporting Information. Species ranges were mapped previously (Laliberte & Ripple 2004) and details on historic and current range-estimate creation is in Supporting Information.
We transformed all shapefiles into rasters. Each raster was approximately 35 million cells at 250 × 250 m. This resolution is significantly finer than the typical management scale (Halpern & Fujita 2013) and allowed detailed analysis over the study area.

Land Use
We separated land use into categories within which impacts will generally be of a similar type but vary in intensity. Infrastructure covered urban and residential development and a small area for mining. These impacts are diffuse and require large amounts of clearing. Roads were analyzed independently because they dissect landscapes on large scales and are linked to compositional changes, abiotic shifts, and vertebrate mortality (Forman & Alexander 1998;Coffin 2007). Oil and gas development in BC consists of isolated physical structures associated with extraction (e.g., drills) connected by a network of access roads and pipelines. Agriculture and rangeland were grouped together.
Our last category was logging, a primary source of BC economic revenue (BC Ministry of Finance 2016). In the first 10 years after logging, vegetation tends to be open and contain high levels of forage. Through time, vegetation thickens and dense stands offer grazers protection from predators (Fisher & Wilkinson 2005). Gradually stands thin and achieve old-growth designation according to their species. The minimum age for old-growth designation is 120 years (Ministry of Forests, Lands, and Natural Resources 2003). We thus considered different times since last logging as different impact levels (0-10 years as of December 2016, 11-60 years, 61-120 years) and then all areas that have been logged plus all land under logging tenure. Tenured land in BC is licensed to private companies for active management and exploration of timber resources (Zhang & Pearse 1996). Though tenured land may also represent future human impacts, human presence and activities are heightened compared with untenured forested land. Thus, we considered it current land use.
Human development has effects beyond the physical footprint. We estimated indirect impacts around a subset of land uses but focused on roads and oil and gas development. Both are structures on the landscape that have known impacts beyond their physical presence. Because the impact to wildlife varies between and within species (Toews 2016), we estimated indirect impacts for increasing distances based on mammal avoidance data. Data collected in BC on avoidance behavior is focused on caribou, whose avoidance patterns range from 250 m (Dyer et al. 2001) to 2,000 m (Polfus et al. 2011). We used these 250-m and 2,000-m endpoints and intermediate distances of 500 m and 1,000 m to represent a range of discrete areas around direct impacts.
We did not estimate indirect impacts around other land uses. Activities such as urban and residential development or agriculture are diffuse rather than singular structures. Reflecting this, data on these land uses were available as large polygons with smoothed boundaries likely incorporating the indirect area distances we chose. For logging we chose the time-since-logged classification to represent different levels of direct to indirect impact. In addition to the temporal component, the tenured land boundaries included the majority of indirect-effect areas that would have been calculated around recently logged land had we used the same buffering protocol as outlined above.

Ecosystem-Level Indicators
Within a BEC zone, we calculated proportional area lost to each of the impacts and calculated fragmentation patterns based on direct land use: infrastructure, roads, oil and gas, agriculture, and logging within the last 120 years. In each zone, we calculated average mean patch area (square kilometers), total edge of all patches (meters), average perimeter-to-area ratio (PAR), aggregation (number of like adjacencies divided by maximum possible number), and the number of patches. We calculated metrics in Fragstats (McGarigal et al. 2002).
We calculated the range extent affected by land-use and indirect effects and assessed fragmentation patterns due to cumulative direct effects.

Population-Level Impacts
We compared current ranges with historic range estimates. Distribution maps can overestimate habitat by smoothing edges and excluding small-scale extirpation (Hurlbert & Jetz 2007), leading to optimistic estimates of species occurrence (Rondinini et al. 2005). In our case, the scale of a provincial-study on multiple species makes finding smaller-scale range data difficult. The distribution maps we used were at a provincial scale and likely provide the appropriate level of detail.
Scale also drove our choice of presence-absence rather than abundance data. Abundance data across the province are available only for caribou and grizzly bears. The smaller-bodied species are far more difficult to count and consequently have less available data. The temporal scale of our study is large. We compared historic range estimates from the 18th century with those of the 20th century. Presence-absence data captured the current outcome of that extended period for each species. In the case of delayed response to recent development (e.g., oil and gas), presence-absence data may not be adequate. Each individual species result was assessed in that light.
To enable statistical modeling of species persistence, we divided the landscape into 25 × 25 km nonoverlapping landscape parcels. Each of the resulting 1,714 parcels had a present, absent, or extirpated designation for each species and an amount of habitat loss and fragmentation. We classified each parcel as either having had local extirpation of any species or having maintained all known populations. We modeled persistence with a series of generalized linear models with a binomial response (0, extirpated; 1, persistent). To incorporate spatial dependency between parcels, we used an autocovariate regression (Dormann et al. 2007). We calculated an autocovariation metric (Augustin et al. 1996) with the spdep package (Bivand & Piras 2015) that was a weighted average of the successes (here, persistence) among all parcel neighbors. The autocovariation was included as a fixed effect in models.
The candidate set had 19 models for each species: a null model with only the spatial covariation, proportion of habitat loss to individual categories, proportion total loss, proportion total loss plus differing levels of indirect effects, and each fragmentation metric. Total loss did not differentiate among landuses. Overlap and spatial distribution of land uses led to high correlation between each use, which left them inappropriate for separate predictors in a single model. We recorded coefficient estimates and the Akaike information criterion (AIC) to assess support for each model. All analyses were completed in R (R Core Team 2014).
We used caution in interpreting model results. Provincial data lack details of land uses such as mining or nongovernmental development (e.g., private logging). Thus, our collected land-use data and coarse range estimates are conservative and likely underestimate total impacts and range contraction. We did not have access to historical land-use data, which limited our analysis to a temporal snapshot potentially underestimating historical impacts.

Land Use
Approximately 13% of BC has been directly modified by humans (Fig. 1). When indirect effects within 2 km of a land use were included, 35% of the landscape was affected. The most widespread use was logging; 7% of the total landscape was under logging tenure. Agricultural development occurred over 5% of the province and oil and gas development over 2.5%. The spatial distribution of some impacts was limited (e.g., oil and gas), although collectively the impacts were widespread. Human land use reduced the average intact patch size in each zone by 62%.

Ecosystem-Level Indicators
The distribution of land-use was as we predicted, with some land uses, such as logging, diffuse across the province. Logging tenures were present in every zone, and active logging occurred in 12 out of 16 of the zones within the last 10 years. In contrast, other land uses were concentrated in relatively small areas. Agriculture and infrastructure occurred largely in a subset of the BEC zones, thus, 3 zones (Bunchgrass, Ponderosa Pine, and Coastal Douglas fir) had 58-74% land conversion (Fig. 2). Some impacts were even more localized; almost 90% of oil and gas development was in the Boreal White and Black Spruce zone. Indirect effects further emphasized these patterns; the three zones under pressure from agriculture and infrastructure all exceeded 85% of their area under direct and indirect effects, and 24% of the landscape covered by oil and gas development had direct and indirect effects in the Boreal White and Black Spruce zone.
Given the spatial distribution of land uses, it was unsurprising that a subset of zones was relatively intact. The eight zones with the least amount of total impacts were all remote or difficult to access: alpine and high-elevation zones. Coastal Western Hemlock was one exception. It had a relatively low footprint and was at relatively low elevation. Large portions are away from population centers and generally surrounded by either ocean to the west or mountains to the east, making access difficult for most activities.
Despite the often localized distribution of land use, there were large levels of habitat fragmentation over most of the province. On average land use increased the number of patches within a single zone by 13 times (Table 1). Increased patch numbers were strongly correlated with increased road cover. Total edge either stayed the same or decreased once impacts were considered; the largest decreases occurred in zones with the highest impacts. In contrast, the PAR increased for all zones due to consistent decreases in mean patch area. For 7 of 16 zones, the mean  patch area was reduced by 90%. Aggregation remained relatively unchanged.

Population-Level Indicators
Bighorn sheep had the smallest range and the highest impact, with 18% direct impact and 45% direct and indirect impact within its range (details in Supporting Information). Direct land use was linked to considerable spatial change within species' ranges. Each distribution was heavily fragmented, ranging from 1,001 individual patches within the bighorn sheep range to 15,299 patches within the wolf range (Fig. 3). Mean intact patch area was reduced by an average of 97% across species.

Population-Level lmpacts
Estimated range losses were from 1% of the historic range for wolves to 42% for fisher (Fig. 3). In all cases, models of persistence ignoring land use-the null models-were not within the top models (Table 2). When we modeled the probability of all species persisting, the null model was the worst and all relationships with land use were negative. Oil and gas development was the only land-use with a consistently neutral or positive relationship with species persistence. It was also the most recent form of development, so effects on measured species may be forthcoming. All other impacts tended to have neutral or negative relationships with persistence, often at great improvement to the null model. Only one species, fisher, showed consistently positive relationship between persistence and land use. Models including total cumulative use were worse than the null for bighorn sheep, caribou, and elk and better than the null for mountain goat, grizzly, and wolf.
Beyond these larger patterns, species varied in their apparent response to land use. Bighorn sheep persistence was related only to recent logging. In contrast, caribou persistence related negatively to all but agriculture and oil and gas on an individual basis but was not related to total effects; models with the lowest AIC values all involved logging. Mountain goat and grizzly persistence related negatively with all human land uses except oil and gas development. Wolf persistence was negatively related to all human uses except oil and gas and logging that occurred <60 years ago.
Including indirect effects did not improve models for any species. In general, species persistence was positively related to mean patch area and negatively related to PAR. These variables were the only improvements on the null for elk. For those species linked to cumulative habitat loss (mountain goat, grizzly, wolf, and all species), we also found positive relationships with aggregation and total edge. Fisher persistence was the opposite-negatively related to mean patch area, total edge, and aggregation and positively related to PAR.

Discussion
The scale and extent of global land-use change is staggering (Wilcove et al. 1998;Pimm et al. 2014), and managers at all scales are struggling to plan for it. In areas such as BC, where rapid and relatively recent land use threatens larger continental-scale values, regional-scale CEAs can inform land-use and conservation policies. Our results support recent claims (Auditor General of British Columbia 2015) that current land-use planning has not prevented substantial losses to ecological values in BC. The ongoing impact of habitat conversion on conservation is consistent worldwide (Foley et al. 2005), and more effective methods must be implemented to achieve global conservation goals. Our study lays the groundwork for a full CEA for a region critical to North American mammal and habitat diversity. Future steps require refinement of the predictive models presented here, scenario creation across the province, and application to land-use decisions at all scales.
As predicted, some land uses were clustered in particular zones, leading to high losses in individual ecosystems. Lower-elevation zones have undergone the largest changes. All three of the most affected zones occur from 0 to 1,000 m. In contrast, five of the six least affected zones occur at 1,000 m and higher. This is unsurprising, given that lower-elevation sites are more accessible for agriculture, resource extraction, and urbanization. Importantly in BC, these are often the zones along the southern border, where some high-and low-latitude species' range limits intersect (Swenson & Howard 2005). These zones are generally the most diverse in the province and have the highest numbers of threatened or rare species (Gibson et al. 2009;Fraser et al. 2011). Continuing on the same development trajectory in the worst affected zones may lead to substantial losses in provincial-level diversity.
Despite only 13% total direct impact, high levels of fragmentation have occurred. Large spatial changes were not always linked to total land-use cover; for example, oil and gas development typically covers a small total area but is composed of scattered linear features. Boreal White and Black Spruce, the zone most developed for oil and gas, has around 20% total direct effects of land use but a 98% decrease in average patch area and a 4,900% increase in the total number of patches. In general, larger ecosystem patches are expected to contain greater species richness than smaller patches and to be exposed to fewer edge effects such as microclimate shifts and altered nutrient cycles (Saunders et al. 1991;Haddad et al. 2015). For individual species, it is likely that fragmentation interacts with habitat loss, potentially compounding the singular   effect of habitat loss (Andrén 1994). Here, the strong correlation between fragmentation and habitat loss made statistically teasing apart the role of each, and their interactions, impossible. Although experienced levels of fragmentation differ among species (e.g., some species may experience roads as corridors rather than barriers [Toews 2016]), such dramatic changes in spatial configuration are broadly concerning for ecological diversity and ecosystem functioning. Total human footprint was largest in the distributions of bighorn sheep, elk, fisher, and wolf. Yet wolves have lost the least amount of range, and models of species response showed the fisher was positively related to many land uses, elk were neutrally related to all land uses, and bighorn sheep were negatively related only to logging. All 3 species were associated with extirpations, but their extirpation areas did not align with intense land use. Thus, they may be more able to adapt to land use than other species. Alternatively, the effects of land use on these species may be better explained by impacts not captured here or at different scales than those we considered. For example, fisher populations have declined historically due to trapping (Weir 2003) and are still subject to trapping pressure (Weir & Corbould 2006). It is unknown whether pressure from human behavior such as trapping is intense enough to cause recorded range losses. Fisher also function in smaller-scale home ranges than species such as caribou and grizzly bears, so fisher responses to land use may be reflected in smaller-scale models. The type of footprint calculated and the relevant scale at which it is modeled is likely to shift on a speciesby-species basis.
Caribou offer an excellent case study for interpreting these data because they have been extensively researched and there are well-supported hypotheses on local extinction drivers. The indirect effects of land use on caribou are linked to interspecific competition. Moose (Alces alces) often move into early serial logging sites near caribou herd ranges (Potvin et al. 2005). Their presence provides prey year-round for wolves, whose populations spike due to winter resource-limitation release (Seip 1992). When summer caribou ranges overlap with moose ranges, predator-induced mortality reaches up to 30% of adult females and 100% of herd calves (Seip 1992). In our models, wolves responded negatively to all human land uses except forest logging within the last 60 years (Table 2), whereas in the best caribou-extirpation model caribou responded negatively to forest logging within the last 60 years and exhibited the largest effect size for forest logging within the last 10 years. Our models effectively highlighted these complex ecological relationships, which implies that other models presented here may accurately reflect natural processes. For instance, we did not find a relationship between elk or fisher and land use; managing these species may need to focus on other threats such as overharvest that may play a larger role in range contraction.
There are also known links between caribou decline and oil and gas development (Hebblewhite 2017) that were not reflected in our results. Rather, oil and gas development was positively linked to caribou persistence because much of the remaining caribou territory is being developed for petroleum. Although herd numbers across Canada have been reduced dramatically (Wittmer et al. 2005; Johnson et al. 2015), population management in BC has been active and extirpation has not occurred, despite precipitous declines. All quantitative models must be complemented by detailed ecological knowledge and monitoring-based refinement of hypotheses (Burton et al. 2014), particularly given that observational models such as this do not test causality. For application in future scenario planning, our models would benefit from data such as abundance, historical land-use, and other forms of human impact.
These results emphasize the presence of threatened ecosystems and species. Ecosystems such as the BEC zone Coastal Douglas Fir, Bunchgrass, and Ponderosa Pine have been largely developed for human use and are vulnerable to further change. Beyond the heavily managed caribou, there is evidence that large carnivores are particularly sensitive to land use. Grizzly bears lost an estimated 14% of their historic range, and persistence was negatively related to all land uses. Wolf persistence had similar trends, despite low range loss. It is well known that large predators are sensitive to land-use change and fragmentation (Crooks 2002;Crooks et al. 2011) because they require large tracts of habitat, large-bodied prey, massive quantities of forage (grizzlies), and protection from conflict with humans (Prugh et al. 2009;Ripple et al. 2014).
Particular types of impact are featured in our results. Logging tenures were found in every BEC zone. Five of the seven species' persistence patterns had negative relationships with logging, including bighorn sheep, a species not found to have relationships with any other land-use type. Although we found that logging recently occurred on <10% of BC, total impacts are likely underreported (Fig. 4). The difference between public information on logging and the on-ground reality underscores the amount of data that may be missing from our calculations. Thus, our footprint estimations may be conservative to varying degrees. A more accurate assessment would likely shift model results and provide stronger hypotheses for management planning.
More generally, our results emphasize that the resources and tools are available for comprehensive CEAs. Land-use data such as those we used are readily available for many regions, and the methods applied are accessible to any manager. Yet CEAs that create future forecasts and assess cumulative impacts regionally, rather than make decisions on a project-by-project basis, are extremely rare (Baxter et al. 2001;Duinker & Greig 2006;Halpern & Fujita 2013). Cumulative effects analysis offers a robust landuse planning tool in changing landscapes but only if they inform decisions at early stages of project development, a process that involves establishing stakeholder-driven ecological values, highlighting areas for conservation or sustainable use, and bridging decision making across regulatory agencies (Johnson 2011). History shows that species in decline have an elevated probability of extinction (Woinarski et al. 2017), and methodically applying available tools is critical for preventing other species from sharing the same fate. This study begins that process for an area of high conservation concern in North America by generating quantitative relationships between land use and probability of extinction that can be applied to future cumulative-effects assessments. and J. Rosenfeld provided insightful feedback that shaped many of the final concepts. We acknowledge the financial support of the Hakai Institute, Mitacs Institution, Pacific Institute for Climate Solutions, Canada Foundation for Innovation, the Natural Sciences and Engineering Research Council of Canada, and The Ian McTaggart Cowan Professorship at the University of Victoria.

Supporting Information
Land-use data sets and descriptions (Appendix S1) and Species range estimates (Appendix S2) are available online. The authors are solely responsible for the content and functionality of these materials. Queries (other than absence of the material) should be directed to the corresponding author.