Conceptual framework and uncertainty analysis for large-scale, species-agnostic modelling of landscape connectivity across Alberta, Canada.

Sustainable land-use planning should consider large-scale landscape connectivity. Commonly-used species-specific connectivity models are difficult to generalize for a wide range of taxa. In the context of multi-functional land-use planning, there is growing interest in species-agnostic approaches, modelling connectivity as a function of human landscape modification. We propose a conceptual framework, apply it to model connectivity as current density across Alberta, Canada, and assess map sensitivity to modelling decisions. We directly compared the uncertainty related to (1) the definition of the degree of human modification, (2) the decision whether water bodies are considered barriers to movement, and (3) the scaling function used to translate degree of human modification into resistance values. Connectivity maps were most sensitive to the consideration of water as barrier to movement, followed by the choice of scaling function, whereas maps were more robust to different conceptualizations of the degree of human modification. We observed higher concordance among cells with high (standardized) current density values than among cells with low values, which supports the identification of cells contributing to larger-scale connectivity based on a cut-off value. We conclude that every parameter in species-agnostic connectivity modelling requires attention, not only the definition of often-criticized expert-based degrees of human modification.

www.nature.com/scientificreports www.nature.com/scientificreports/ the contribution of an area to larger-scale connectivity, so as to compare the impact of alternative development scenarios and prioritize areas for conservation and restoration.
After deriving a main wall-to-wall map of current density as a measure of species-agnostic connectivity for Alberta based on the degree of human modification, we use this example to compare the contributions of different factors to uncertainty. We thus assess the sensitivity of this map to (i) the definition of the degree of human modification (i.e., the relative weight given to degree of physical human footprint vs. intensity of human use), (ii) the decision whether water bodies are considered barriers to movement, and (iii) the scaling function used to translate degree of human modification into resistance values. While the influence of most of these parameters on connectivity patterns has been previously studied independently, we aim to directly compare their contributions to uncertainty. We further (iv) evaluate the effect of such uncertainty on the classification of cells as being important or less important for maintaining larger-scale connectivity using different cut-off values. This study prepares the ground for future research to evaluate model validity with Alberta's extensive biodiversity monitoring data. conceptual framework. We argue that a species-agnostic approach is conceptually better suited for integrating landscape connectivity into land-use planning, whereas a focal-species approach is better suited for conservation management (Fig. 1). Understanding the difference between these perspectives can help clarify conceptual differences, guide researchers in making decisions about how to model connectivity, and inform practitioners about the potential and limitations of resulting maps. Conservation is often based on focal species, e.g., in species-at-risk management, where it is paramount to adopt an organism perspective 60 to define critical habitat and to consider the organism's ability to move between habitat patches. This will result in a model of potential functional connectivity for the specific organism of interest, and a multi-species model can be derived by overlaying models for a representative suite of species 10,22 . In contrast, land-use planning focuses on the sustainable development of multi-functional landscapes. In this human-centred perspective, land parcel ownership and administrative boundaries define the relevant spatial scale and the degree to which landscape development can be influenced by policy, which in the case of Alberta includes the introduction of ecosystem services and biodiversity markets that play a major role in balancing environmental considerations with socio-economic drivers 59 . Regarding landscape connectivity, the focus thus lies on how human landscape alteration affects the connectivity of the remaining natural heritage system and how to compare the expected effect of local development alternatives on larger-scale connectivity. This focus is highly compatible with the modelling of connectivity based on human modification, which will result in a model of structural landscape connectivity.
A coarse-filter approach to biodiversity conservation management falls somewhere in-between 61,62 . The focus in this approach lies on preserving a community, such as grassland, wetland, or forest interior species, by maintaining a sufficient amount, quality, and connectivity of the respective ecosystem in the landscape. This implies that the species within a community have similar needs and characteristics, including their ability to move between habitat patches. Regarding landscape connectivity, we suggest that species-level data, such as mark-recapture, radio-telemetry, and molecular genetics data, which provide a direct measure of actual functional connectivity, may be used to test and compare species-agnostic (top-down) or multi-species (bottom-up) connectivity models. At the community-level, however, data such as biodiversity monitoring that is systematically Conceptual framework: different perspectives on the landscape (a) and corresponding goals and approaches to connectivity modelling (b). The horizontal axis (black) illustrates how the landscape definition shifts from a multi-functional landscape in the context of land-use planning to an organism perspective in the context of species conservation. The vertical axis (grey) illustrates how the focus shifts from natural processes in basic ecology to human values and actions in a policy context. These values are reflected in policy that aims to preserve ecosystem services, biodiversity in general, or designated species of concern. Species-agnostic modelling of landscape connectivity may be most suitable in a land-use planning context, whereas focal species approaches may be most suitable in a conservation management context. When biodiversity data are available for many taxa, multi-species overlay can also be an efficient option to model connectivity for land-use planning and conservation programs for a representative suite of species. Validation of structural connectivity or potential functional connectivity based models could also be done with appropriate biological data, to assess actual functional connectivity.
collected across a large geographic area can be used as indirect measures of functional connectivity. If a connectivity model can explain, e.g., variation in species composition that remains unexplained by local site characteristics, this may indicate that it successfully captured the shared response of a broad range of species to human landscape modification. Note that there is a lack of empirical studies that compare a general-use connectivity assessment based on the overlay of many species-specific models across taxonomic groups to a species-agnostic model.
Another distinction occurs along the vertical axis in Fig. 1a: while ecology primarily considers natural processes, policy addresses human values and actions. The two perspectives meet in the consideration of the human impact on ecosystems. Based on this framework, a species-agnostic approach to modelling landscape connectivity quantifies how human actions and their manifestation as non-natural landscape features constrain ecological flow across the underlying natural fabric of the landscape. The main goals of this approach are to identify critical linkages and to support the evaluation of alternative development scenarios, with the general aim of maintaining ecosystem services and biodiversity in the context of sustainable development. This cannot replace fine-filter approaches for specific conservation goals, such as species-at-risk management. Note also that ecosystems can have an impact on humans as well, such as the spread of animal-transmitted diseases (e.g., West-Nile virus or Lyme disease). Here, we focus on human impacts on natural ecosystems.
To implement a species-agnostic approach, resistance values are assigned to non-natural landscape features based on expert opinion 33 . Conceptually, this involves two steps: (1) quantifying the degree of human modification, and (2) using a scaling function to assign resistance values based on the degree of human modification. We argue that beyond an inevitable element of subjectivity in expert-defined values, different interpretations of what constitutes degree of human modification can result in considerable differences between values. For instance, two types of roads may have the same physical footprint in terms of deviation from natural conditions (impermeable surface) but differ vastly in their intensity of use (traffic volume), which will affect e.g. road mortality rates. The relative scaling of resistance values is known to have a large effect on connectivity modelling, and a sharp contrast between low, medium or high resistance values has previously been suggested and used 34,38,39 .
In addition to human modification, some authors including Dickson et al. 34 also considered natural barriers to the movement of terrestrial organisms and applied resistance values to features like water bodies and topography. Conceptually, this introduced an element of potential functional connectivity as it makes assumptions about the movement ability of organisms. In practice, this raises the concern that, e.g., water bodies may be traversable for some terrestrial species but not for others. Specifically, water bodies are a special case in such a framework as they are not a terrestrial feature, and while they are mostly natural, some authors have treated them as a barrier. However, in the context of species-agnostic models, much less thought has been given to the decision how to conceptualize water bodies than to other factors.
Here, we define a conceptually justified range of variation for the three factors identified above, with the goal of directly comparing their relative importance, as they have been addressed independently in previous studies on species-agnostic connectivity modelling: (i) the assignment of degree of human modification values by degree of deviation from natural conditions (footprint) or by the intensity of use; (ii) the assignment of resistance value to water bodies between minimum (natural) and a value close to the maximum of all human features combined (near-complete barrier); and (iii) the choice of scaling function resulting in a weak or strong contrast between non-natural and natural landscape features. We use intermediate values 33 to create the main map, and we assess uncertainty in a factorial design using the extremes of each factor.
Circuit theory models connectivity as current density, where areas with high current density are interpreted as contributing to larger-scale landscape connectivity and thus 'ecological flow' across the study area 32,42 . We thus assess the contribution of each factor to variation among maps in the quantification of current density as the contribution of a cell to large-scale connectivity. We argue that visual interpretation of current density maps is affected by the choice of colour ramp, which should reflect a conscious decision about how cells with important contribution to large-scale connectivity (high current density values) are identified. J. Bowman (Trent University, pers. comm.) suggested that visualizing and interpreting current density values that have values higher than one standard deviation above the map mean (i.e., z scores> 1) as contributing areas provides a good balance for application. Here, we further vary this cut-off between the mean (i.e., z scores> 0) and two standard deviations above the mean (i.e., z scores> 2) to assess whether the uncertainty related to the factors above varies with the cut-off level used.

Results
The main map (Fig. 2) shows prominent connectivity pathways extending mainly from the southwest to the northeast of Alberta. More specifically, areas of high current density are concentrated in the foothill region, connecting the Rocky Mountain range to the forest-dominated northern part of the province. Patterns of high current density were also observed in the south-eastern grassland and in the north-eastern Canadian Shield region of the province.
The uncertainty analysis showed large variations in absolute current density values among maps, both in their ranges (Table 1) and distributions (Fig. 3). Differences in the shape of distribution of current density values were most pronounced for values below the mean. The proportion of cells with z > 1 (or z > 2) was quite constant, but more variable for cells with z > 0 (Fig. 3). These variations resulted in different patterns of high current flow pathways between maps (Fig. 4).
Cell-by-cell correlations in current density values between maps indicated a wide range in concordance between maps varying from a minimum correlation of r = 0.34 to 0.94 (Table 2). Redundancy analysis showed that map correlations were most affected by differences in water resistance and scaling function ( www.nature.com/scientificreports www.nature.com/scientificreports/ explained 12.5% of the variation and 16.7% of variation among maps could not be explained by marginal effects. Such unexplained variation could be related to specific combinations of factor levels (i.e., interactions) or random variation, e.g., due to random selection of pairs of nodes.
When using a cut-off (z ≥ 0, 1, or 2) to classify cells as important or unimportant for connectivity based on the z-score of their current density value, a similar pattern emerged. However, consistencies between map classifications were more affected by scaling function (39.4-47.8%) than water resistance (28.6-33.1%). The effect of variation in the H index was considerably smaller, increasing from 1.8% for z ≥ 0 to 3.7% for z ≥ 2, whereas the proportion of variation unexplained by marginal effects was larger (21.8-27.8%).
Qualitative interpretation of the correlation matrix in Table 2 suggests considerable interactions between factors (which we did not formally include in the db-RDA model to avoid over fitting). For pairs of maps with the same water resistance level, the differences between H F and H U maps were more pronounced with the high-contrast scaling (correlation r = 0.78 and 0.79) than for the low-contrast scaling (r = 0.94). With high-contrast scaling, the effect of water was more pronounced when using the H U index (r = 0.54) than when using the H F index (r = 0.73). With low-contrast scaling, water resistance had a smaller impact on the maps (r = 0.85), irrespective of the choice of H index.

Discussion
We presented a conceptual framework for species-agnostic connectivity modelling (Fig. 1) and applied it to model connectivity as current density across the entire province of Alberta, Canada. We assessed the effect of conceptual and computational decisions involved in implementing the approach, and provided a much-needed direct comparison of the degree of uncertainty related to these decisions. First, we examined the effect of the conceptualization of human modification (H index) represented by the degree of physical footprint H F and the intensity of human use H U on current density maps. Second, we analysed the impact of the scaling function that was used to convert H values to resistance, and therefore the contrast between cells with low, intermediate, or high degrees of human modification in their ability to constrain current flow. Thirdly, we assessed the effect of the conceptualization of water bodies as a near-complete barrier on the resulting current density maps. Lastly, we evaluated whether the sensitivity of current density maps to these three parameters is affected by the cut-off  ; water resistance (either 0 or 1000; dark grey); and scaling function (either linear or derived power; light grey). Analysis was done separately for each z-score cut-off used to identify important cells, (either z = 0, 1, or 2) and for the cell-by-cell correlations in current density values between maps (Cor). Variances are presented as percentages (%). (2020) 10:6798 | https://doi.org/10.1038/s41598-020-63545-z www.nature.com/scientificreports www.nature.com/scientificreports/ threshold used to identify cells that contribute to larger-scale connectivity. We discuss the importance of these factors in an applied land-use planning context.
Since resistance values commonly rely on arbitrary definition by expert opinion, connectivity modelling is often criticized for human decision bias that might affect the resulting maps 23 . In a species-specific approach, expert opinion can be based on a thorough understanding of the life history and ecology of the focal species (e.g. 21 ), but key parameters may be unknown or incorrectly assessed. While new optimization methods have  www.nature.com/scientificreports www.nature.com/scientificreports/ been proposed for fitting resistance values based on measures of actual functional connectivity, such as molecular genetic data for a focal species 41 , inference of resistance values has been shown to vary between landscapes even for a single species 63 . The species-agnostic approach relies on the assumption that ecological flow is negatively related to degree of human modification, though this concept is often vaguely defined and may include many dimensions of human modifications. In our study, we conceptualized the degree of human modification for each human-footprint category using two independent parameters, the degree of physical human footprint (H F ) and the intensity of its use (H U ). An expert panel defined both sets of values, separately for the forest-and agriculture-dominated areas of Alberta. In line with the land-use planning perspective of a species-agnostic modelling of landscape connectivity (Fig. 1b), the assigned values reflect expert consensus about degree of human modification, and not resistance to organismal movement. We believe that such conceptual clarity is important as it makes species-agnostic models more transparent and better justifiable. For the main map, we then used the mean of the two indices H F and H U . Despite differences in quantitative and relative values of human modifications attributed to human-footprint categories, we found that the resulting current density maps showed limited sensitivity to these variations. Pairwise correlations between maps (based on cell-by-cell current density values) showed a moderate degree of variation between H F and H U maps with high-contrast scaling. When a cut-off was used to identify cells that contribute to larger-scale connectivity, however, there was a high degree of consistency between maps. This suggests that the differences may be largely due to variation in the lower end of the distribution of current density values, where we found considerable differences in the shape of distributions (Fig. 3). From an applied perspective, consistency in the upper tail (identification of areas important for larger-scale connectivity) arguably is more important than minor distinctions among areas that contribute relatively little to larger-scale connectivity.
As expected from the literature 38,39 , the scaling method had a major impact on current density maps. We found that a high-contrast scaling amplified the effect of other parameter settings (water resistance, H index). High-contrast scaling function increased sensitivity to conceptual decisions regarding the constraining effect of human landscape alterations and topographic features. Such sensitivity is not necessarily bad; rather, large contrasts may be required to effectively model constraints on connectivity. For instance, Koen et al. 28 assigned resistance values of {10, 100, 1000} to natural, unnatural but permeable to movement, and unnatural impermeable land cover types, and successfully validated their current density map for South-eastern Ontario with road mortality data for 20 reptile and amphibian species and with radio-telemetry data for fishers (Pekania [Martes] pennanti). As a preliminary step in their study, Dickson et al. 34  Connectivity modelling requires a careful consideration of not only humanmodification type and degree of intensity, but also of natural landscape features that could form potential barriers to movement. We largely followed Dickson et al. 34 in assigning resistance values for water bodies and slope, and Koen et al. 28 likewise treated the resistance of water bodies as barriers equivalent to impermeable non-natural landscape features. However, this is a move away from the conceptually 'clean' modelling of landscape connectivity as a function of human modification. Assigning resistance values to natural landscape features implies assumptions about organism movement ability and thus a taxonomic reference. Inevitably, examples can be found where a natural topographic feature may act as a barrier for one species but not for another, thus crossing over from a pure species-agnostic approach to an implicitly multi-species approach (Fig. 1).
Our sensitivity analysis showed a very large effect of the resistance value attributed to water on the resulting current density maps. Based on pairwise correlations between maps, the variation due to resistance to water was considerably larger than the variation due to scaling function. (Note that we did not vary the parameters for resistance due to slope, but its range is more limited, where a slope of 100% results in an increase of resistance values by 25 only). Given the magnitude of the effect of water resistance, the lack of awareness of its importance in previous studies is surprising. For the main map, we thus deviated from Dickson et al. 34 and used an intermediate resistance value for water.
Studies that apply circuit theory to model landscape connectivity are lacking an agreement on a cut-off to distinguish cells that are contributing to the large-scale connectivity. In our study we explored the sensitivity of our maps to the change in three cut-offs represented as z-scores for standardized current density values. The shape of the lower tail of the distribution of current density values varied considerably between maps. The proportion of cells with standardized current density values of z > 1 (or z > 2), however, was relatively constant, but more variable for cells with z > 0. Map comparison with db-RDA showed that overall, cell classification was surprisingly robust to the definition of human modification (H F vs. H U ), which becomes slightly more important at higher cut-offs. Taken together, these results indicate that a cut-off of z > 1, as suggested by Bowman (pers. comm.), may be a robust threshold for identifying cells that are contributing to connectivity.
Species-agnostic approaches for modelling landscape connectivity have been criticized mainly for not being based on species traits. We argue that the power of such species-agnostic connectivity models relies on their generality and direct link to multifunctional land-use planning, as H values are assigned to policy-relevant human-footprint categories 28 . As we indicated in our conceptual framework, species-agnostic models are not meant to replace species-specific models designed explicitly to serve specific conservational goals, such as Scientific RepoRtS | (2020) 10:6798 | https://doi.org/10.1038/s41598-020-63545-z www.nature.com/scientificreports www.nature.com/scientificreports/ species-at-risk management. It is important to acknowledge that even the use of a species-agnostic model implies assumptions about organismal movement, albeit in a generalized way. The decision of how to treat natural elements like water or slope has been discussed above. An implicit assumption is also that human landscape features affect most species in a similar way. However, examples from the literature indicate, e.g., that some species may actively use roads as corridors, whereas for other species, forest may impede movement 64 . Acceptance of species-agnostic models of landscape connectivity will likely depend on validation with independent ecological data, as successfully performed by Koen et al. 28 for a limited set of species and at a smaller spatial scale. While testing the different models is beyond the scope of the present paper, we plan to use the unique biodiversity monitoring data (with 1656 sites spaced 20 km across Alberta; ABMI) and derived species-level distribution models available through the ABMI to test our connectivity maps. This could be done by assessing which proportion of species turnover across large taxonomic and functional groups can be explained by the current density (i.e., value of being near a "biodiversity highway") that cannot be explained by current species distribution models based on local site conditions and composition-based landscape metrics (including habitat amount). More research is needed to empirically compare the performance of top-down vs. bottom-up connectivity modelling for land-use planning purposes, i.e., to ascertain whether overlaying many species-specific models would result in a substantially different general-use connectivity assessment than species agnostic models.
The main goal of a large-scale connectivity map is to identify and visualize 'biodiversity highways' , or important routes of ecological flow. Note that a large-scale connectivity map is not designed to evaluate connectivity in a network of patches (e.g., a natural heritage system). To evaluate connectivity between two specific patches, considering all possible movement paths by using current density, nodes would have to be placed in the two patches 42 . In a network model, it is possible to quantify the contribution of each patch or link to overall probability of connectivity 65,66 . However, this should not be done based on current density, as such network analysis aims to distinguish between alternative paths, and hence least cost path distances are better suited for this goal. Linkages between pairs of patches can then be modelled probabilistically based on assumptions about organism's dispersal ability 50 .
Large-scale, high-resolution modelling of current density remains a computational challenge, despite the new, computationally efficient implementation of Circuitscape in Julia 52 or the alternative implementation in GFlow 49 . Using GFlow on the Niagara supercomputing facility (SciNet HPC Consortium), we still encountered technical limitations that prevented us from modelling current density at a spatial resolution below 100 m across the province. We calculated resistance values at a high-resolution of 10 m to minimize artefacts related to rasterisation of linear features 67,68 . However, the barrier effect of linear features likely was dampened by aggregating resistance within 100 m cells prior to current density modelling.
Land-cover changes due to human activities and climate change are main drivers of land-use planning, and their study is the basis of land system science (as defined in Verburg et al. 69 ). This discipline focuses on both the drivers and impacts of land-use change as part of global change. Land systems also offer solutions to global change through adaptation and mitigation and can play a key role in achieving a sustainable future for Earth. Recent studies propose to incorporate landscape ecology issues when designing land systems to understand how and why governance impacts human land change decisions and how those land change patterns influence, and are affected by, the underlying ecological processes (e.g., 70 ). Species-agnostic connectivity modelling based on human modification provides a powerful tool in evaluating the effect of human decisions and actions on multi-level ecological processes and dynamics and simulating the effect of alternative future scenarios.

Methods
Study area. Our study area, Alberta, is a western Canadian province bordered by British Columbia to the west, Saskatchewan to the east, the Northwest Territories to the north, and Montana (USA) to the south (Fig. 6a). We divided Alberta into two regions (Fig. 6). The forest-dominated area (61% of area), which includes the northern half of the province and the Rocky Mountain area, is mostly covered by native and managed boreal forest and contains most Crown lands and public lands. The agriculture-dominated area (39% of area), which includes the southern half and the central-western part of the province, is dominated by agricultural lands and grasslands and includes most of the urban areas.
Spatial human footprint data. To map landscape structural connectivity in Alberta as current density, we used ABMI's Human Footprint Inventory (HFI) derived from the detailed 2014 land-use dataset (for complete information on data collection and computation of the land-use geodatabase, see Geospatial Centre ABMI, 2018). This dataset was derived from 2014 SPOT6 Satellite Imagery (1.5 m Color SPOT 6 Mosaic) provided by Alberta Environment and Sustainable Resource Development, Informatics Branch (https://www. alberta.ca/environment-and-parks.aspx). The original geodatabase consisted of 21 polygon-based layers including 115 land-uses grouped into 6 main categories: agricultural, forestry, transportation, urban, energy, and human-created water bodies 71 .
To avoid boundary effects 72 , we mapped human footprints in a buffer zone surrounding Alberta, using several datasets collected at the federal levels of Canada and the USA. Land cover in the Canadian buffer zone was mapped using raster thematic data originating from classified Landsat 5 and Landsat 7 ortho-images circa 2010 (resolution: 30 meters) 73 . The physical footprint of roads and railways in the Canadian part of the buffer zone was mapped using the National Road and Railways Network databases 74,75 . Land cover in the USA buffer zone was characterized using the National Land Cover Database 2011 for the conterminous United States (NLCD 76 ; www. mrlc.gov). The physical footprint of roads and railways in the USA was mapped using TIGER 2014 data (http:// www.census.gov/geo/maps-data/data/tiger-line.html). After homogenizing and grouping variables for each dataset, we retained 84 human-footprint categories in total, with 66 categories for Alberta, 8 for the Canadian buffer zone, and 10 for the USA (Appendix 1). (2020) 10:6798 | https://doi.org/10.1038/s41598-020-63545-z www.nature.com/scientificreports www.nature.com/scientificreports/ Mapping landscape resistance. We modified the methodology by Theobald 33 and Dickson et al. 34 to map landscape resistance and connectivity across Alberta based on the degree of human modification, H. At a given location, H is the sum of the individual degrees of human modifications h caused by all human-footprints at this location 33 . As a major modification of the method, the value of H attributed to each human-footprint was conceptualized as the combination of two parameters: degree of physical footprint and intensity of human use (Appendix 1). Values of degree of footprint H F were based on the degree of change from the native state and the degree of modification. Values of intensity of human use (H U ) were based on the assumed average of the amount of human passage, presence, and density throughout the year. H F and H U values were assigned summarily to each human-footprint category (not, e.g., for individual roads) and were allowed to differ between the forest-and agriculture-dominated areas of the province (Appendix 1). Both the H F and H U indices were calculated on a scale from 0.0 (low footprint, low human use, respectively) to 1.0 (high footprint, high human use, respectively). Values of H F and H U indices were assigned based on expert opinion of the Alberta Human Footprint Technical Committee (Alberta Human Footprint Mapping Program). We explicitly asked the experts to focus on the degree of physical footprint and intensity of human use, respectively, without any consideration of particular taxa and their relation with considered human-footprints. The main, average map was calculated as the mean between H F and H U indices, hereafter referred to as H FU , whereas for the sensitivity analysis, we used maps with either H F or H U but not both.
For each human-footprint j (n = 84) separately, we attributed an individual value of either 1 or 0 to each 10 × 10 m cell of a raster grid covering the study area, depending on whether or not the cell overlapped with polygons of the considered human-footprint category. We then applied human modification values h j according to Table 1, separately for H F , H U , or H FU , to cells of value 1. Finally, we combined the 84 layers h j into a single, overall metric of the degree of human modification H, separately for H F , H U , and H FU . In the absence of human footprint, we used a value of h = 0. When a cell was impacted by multiple human footprints j, we assumed that it should have a higher degree of human modification than a cell with a single human-footprint 33 . We used the fuzzy algebraic sum 77 for which the result is always at least as large as the largest contributing factor, so that the effect is additive, but never exceeds 1 33 . The degree of human modification H i at each cell i was then calculated, separately for each H F , H U , and H FU , as: where h j is the degree of modification for an individual human-footprint layer j (for j = 1…k).
Dickson et al. 34 derived resistance values R by rescaling H values with a scaling function and adding resistance values for areas with steep slopes s and for water bodies w (in their case, major rivers), to account for their possible constraining effect on organismal movement: 10 In this equation, the maximum resistance due to H is (1 + 1) 10 = 1024, a slope percentage of s = 100 carries a penalty of 100/4 = 25, and a water body (w = 1) carries a penalty of 1000, which corresponds to a value of H = 0.999 (Note that in the rare case of a cell with intensive human modification (H = 1) and water body (w = 1), the www.nature.com/scientificreports www.nature.com/scientificreports/ total resistance could technically reach a value of up to 2024). Percent slope s was derived from the US National Elevation Dataset for the entire study area 78 .
We modified Eq. 2a for the main map to assign water bodies a lower resistance value of 57.7, which corresponds to an intermediate value of H FU = 0.5 (Eq. 2b): For the sensitivity analysis, we independently varied the H index (H F or H U ), the scaling function, and resistance of water bodies in Eq. 2a as shown in Table 3, resulting in eight combinations. Raising (H + 1) to the power of 10 results in a very strong contrast between natural and non-natural landscape features 34 . As a low-contrast alternative, we used a linear scaling function, (1 + H * 1000), which preserves the proportional importance of human modification values. This means that for instance, resistance values R based on H F with low-contrast scaling of H and low resistance of water bodies w were calculated as: Because of computational limitations, we aggregated the 10 × 10 m resolution resistance layers R to 100 × 100 m resolution before assessing current density. Cells were aggregated using the mean resistance value of combined 10 × 10 m cells (aggregate function in R, raster package 79 ).
General procedure for current density computation. We derived a current density map from each resistance map R at the 100 × 100 m resolution using GFlow 49 . GFlow is a new software which parallelizes the computation of circuit theory 42 and allows for simultaneous large-extent and fine-grained connectivity modelling.
We simulated current density between pairs of points (nodes) randomly chosen at each iteration among 50 points that were evenly distributed (approx. every 100 km) along the outer margin of the buffer zone (Fig. 6). This methodology reduces node location bias compared to randomly selecting nodes within the study area and requires fewer pairwise computations 28 . A buffer zone with a linear width of at least 20% of the width of the study area has been shown to be sufficient to remove the effects of node placement on current density 28 . We used the "Buffer by Percentage" plugin in QuantumGIS 80 with a setting of 250% to ensure a sufficient buffer size around Alberta's polygon (Fig. 6).
To ensure that we considered a sufficient number of random pairs of nodes to converge upon a solution, we set the convergence factor implemented in GFlow (correlation function) to 3 N (=0.999). The convergence value is calculated for each pairwise solve (i.e., iteration) and ranges from 0 to 1. With an increasing number of iterations, the current density maps converge, so that additional connections have a marginally small effect on the result. When the convergence threshold is reached, the computation stops and the current density map is completed 49 . As a final step, we removed the buffer zone to only keep Alberta's extent.
Map comparison and evaluation of the uncertainty. We standardized current density values as z-scores (z = x /sd) by subtracting the mean value x , and dividing by the standard deviation sd. We used three different thresholds (z ≥ 0, 1, or 2) for identifying cells with important contribution to connectivity. We assessed the relative contribution of the three factors listed in Appendix 1 to the variance between current density maps in a factorial design with variation partitioning using distance-based redundancy analysis (db-RDA 81 ). First, we derived an 8 × 8 correlation matrix by calculating the Pearson correlation r between each pair of maps. For each of the three cut-off levels, we derived an 8 × 8 similarity matrix by counting, for each pair of maps, the proportion of cells that were classified the same way (either as 'important' or as 'unimportant' in both maps). All four matrices were converted to dissimilarity matrices. For each dissimilarity matrix, we performed a db-RDA 82 where the three explanatory factors (H index, scaling function, and water resistance) were coded as dummy variables. We used the function varpart in the R package vegan to obtain adjusted R 2 values for the contribution of each factor.
Unless otherwise stated, all computations and analyses were performed in R 3.4.2 83 on Niagara supercomputer at the SciNet HPC Consortium.

Data availability
The full human-footprint dataset of Alberta can be found at https://abmi.ca/home/data-analytics. Dataset and R script used for the ANOVA / statistical analysis is available at https://doi.org/10.5683/SP2/KJVW8Y.

Factor
Level Parameter values H index H F See Table 1 H U See Table 1 Scaling function Low contrast (1 + H * 1000) High contrast (H + 1) 10 Water resistance Low 0 High 1000 Table 3. Independent factors and levels considered in the sensitivity analysis and their levels, and their associated parameter values.