Integrating functional connectivity in designing networks of protected areas under climate change: A caribou case-study

Land-use change and climate change are recognized as two main drivers of the current biodiversity decline. Protected areas help safeguard the landscape from additional anthropogenic disturbances and, when properly designed, can help species cope with climate change impacts. When designed to protect the regional biodiversity rather than to conserve focal species or landscape elements, protected areas need to cover a representative sample of the regional biodiversity and be functionally connected, facilitating individual movements among protected areas in a network to maximize their effectiveness. We developed a methodology to define effective protected areas to implement in a regional network using ecological representativeness and functional connectivity as criteria. We illustrated this methodology in the Gaspésie region of Québec, Canada. We simulated movements for the endangered Atlantic-Gaspésie caribou population (Rangifer tarandus caribou), using an individual-based model, to determine functional connectivity based on this large mammal. We created multiple protected areas network scenarios and evaluated their ecological representativeness and functional connectivity for the current and future conditions. We selected a subset of the most effective network scenarios and extracted the protected areas included in them. There was a tradeoff between ecological representativeness and functional connectivity for the created networks. Only a few protected areas among those available were repeatedly chosen in the most effective networks. Protected areas maximizing both ecological representativeness and functional connectivity represented suitable areas to implement in an effective protected areas network. These areas ensured that a representative sample of the regional biodiversity was covered by the network, as well as maximizing the movement over time between and inside the protected areas for the focal population.

Atlantic-Gaspésie caribou population, so data that could help track and find these animals won't be made available publicly. However, the modeled data we used to quantify the balance between ecological representativeness and functional connectivity will be made available in a public repository (DRYAD) upon acceptance for researchers who meet the criteria for access to confidential data. and contact information or URL

Abstract
Context: Land-use change and climate change are recognized as two main drivers of the current 25 biodiversity decline. Protected areas attempt to secure the landscape from additional land modifications and likely help species cope with climate change impacts. To maximize effectiveness, protected areas need to cover a representative sample of the regional biodiversity and to be functionally connected, facilitating individual movements among protected areas in a network. 30 Objectives: We developed a methodology to define effective protected areas to implement in a regional network using ecological representativeness and functional connectivity as criteria. We illustrated this methodology in the Gaspésie region of Québec, Canada.
Methods: We simulated movements for the endangered Atlantic-Gaspésie caribou population (Rangifer tarandus caribou), using an individual-based model, to determine functional 35 connectivity based on this large mammal. We created multiple protected areas network scenarios and evaluated their ecological representativeness and functional connectivity for the current and future conditions. We selected a subset of the most effective network scenarios and extracted the protected areas included in them.
Results: There was a trade-off between ecological representativeness and functional connectivity 40 for the created networks. Only a few protected areas among those available were repeatedly chosen in the most effective networks.
Conclusions: Protected areas maximizing both ecological representativeness and functional connectivity represented suitable areas to implement in an effective protected areas network.
These areas ensured that a representative sample of the regional biodiversity was covered by the 45

Introduction
Habitat change is recognized as the main driver of the current declines of terrestrial species [1-3]. Securing habitats by creating or expanding protected areas networks is part of the solution to the challenge of biodiversity loss [4]. A regional protected areas network could be considered 55 ultimately effective insofar as it can sustain the region's biodiversity into some reasonably foreseeable future. Such effectiveness is not guaranteed [5,6] and could be limited by many factors [7]. We consider three such factors: ecological representativeness [8], functional connectivity among protected areas within the network [9] and resilience to climate change [10].
Ecological representativeness measures the degree to which the various non-60 anthropogenic habitats within a focal region are available within a protected areas network in proportion to their regional abundance [11,12]. When protected area locations are skewed towards certain types of habitats [13,14], usually for economic or social reasons, ecological representativeness will be low and high valued habitats for some species may be underrepresented [8]. Conversely, when ecological representativeness is high, it is reasonable to 65 assume that the habitat requirements for most species will be satisfied within the protected areas network. This assumption is usual in conservation and representativeness is one of the core concepts in systematic conservation planning [12,15].
A high degree of ecological representativeness is a necessary condition for an effective protected areas network, but it is not sufficient. There are species whose requirements are not 70 automatically satisfied by ecological representativeness alone. One example would be wideranging species with habitat requirements varying among seasons or life history stages.
Functional connectivity is "the degree to which the landscape facilitates or impedes movement among resource patches" [16] or, as here, among protected areas within a protected areas 6 network. Functional connectivity is species-or population-specific [17,18]. A protected areas 75 network with high functional connectivity facilitates the movement between different protected areas for individuals of a given species, increasing their access to resources and ultimately the rates of recruitment or survivorship [19]. Increasing functional connectivity may thus increase population size and decrease extinction risk for the vulnerable species [19]. These effects would increase the effectiveness of a protected areas network. 80 Climate change is among the main current drivers of ecosystem change and its negative impacts on biodiversity have increased rapidly over the past century [1,3,20,21]. Climate change disrupts environmental patterns and species' habitats globally [22,23]. As a result, species distributions and individual movement patterns are impacted [24,25]. The effectiveness of a fixed protected areas network designed for current conditions may decrease in the future, if 85 either ecological representativeness or functional connectivity decline. Enhancing or maintaining functional connectivity inside protected areas networks is thus one approach to helping species cope with climate change [10,26]; individuals would be better able to access resources available in distant protected areas, or even to migrate from less to more favorable areas. Because of the "cost of waiting", managers should proactively account for future climate change effects when 90 implementing new protected areas networks [27] or expanding existing networks. One important effect would be a loss of functional connectivity. Therefore, the ability of a protected areas network to sustain functional connectivity under is another dimension of effectiveness.
Many methods exist for designing protected areas networks to achieve ecological representativeness. In the systematic conservation planning literature [12], variations of the site 95 selection problem are posed, where one seeks a subset of available sites that, in aggregate, achieve some measure of ecological representativeness at a near-minimum of total area or cost.

7
Variations of these approaches exist that can also partially satisfy other conservation objectives such as topological connectivity, where selected sites are spatially aggregated in some degree [e.g. 28]. 100 The aim of our paper is to present a new method to include functional connectivity in the design process, while accounting also for its persistence under climate change. Our method applies individual-based movement models to simulate maps of habitat use by a population, under present conditions and under hypothetical future conditions reflecting climate change scenarios. From these maps, we derive numerical indices of functional connectivity that allow 105 alternate protected areas network designs to be compared. These are coupled with a novel variant of the site-selection algorithm that constructs protected areas networks of specified total area that can achieve a high degree of ecological representativeness while also satisfying secondary constraints on connectivity and intactness; the variant is presented in detail in Saucier et Our goal is to identify protected areas network designs that simultaneously achieve high degrees of ecological representation and functional connectivity under present conditions and future climates. Achieving this is a multi-objective optimization problem [29,30] that does not admit a unique solution. Instead, one can define a tradeoff surface by those points where no single objective can be increased without decreasing at least another. Such points are said to be 115 Pareto optimal; points on the interior of this surface are sub-optimal in that other solutions exist that are better in terms of all the dimensions considered. Our design methodology uses randomization to generate a large sample of feasible solutions for a given protected areas network design problem. Using these samples, we can approximately define the tradeoff-surface and identify a set of feasible, near-equivalent solutions in the vicinity of any specified point on 120 that surface.
We illustrated our method in the Gaspésie region (Québec, Canada) using the endangered Atlantic-Gaspésie population of woodland caribou (Rangifer tarandus caribou) as our focal population. This isolated herd, that has been identified as one of 12 Designatable Caribou Units considered irreplaceable components of Canada's biodiversity [31], is declining since the late 125 19 th century and reached a critical abundance level [32]. To measure functional connectivity with respect to this population, we applied an individual-based model of animal movement previously developed for this population [33]. The government of Québec is currently engaged in expanding the existing network of protected areas in this region to attain a proportional area target of 12% [34]. To place our work in the context of this ongoing conservation planning exercise, we 130 employed a variant of our design methodology that constructs protected areas networks by adding new protected areas to the existing network. We identified areas that were of high relative importance in achieving ecological representation in the Gaspésie region and functional connectivity for the Atlantic-Gaspésie caribou under present conditions and future climates.

Study area
The Gaspésie natural region (latitude extent: 47.98 to 49.20º N, longitude extent: 64.11 to 67.57º W; Fig. 1) is a physiographically defined area of approximately 25,000 km² at the eastern end of the Gaspésie peninsula, in eastern Québec, Canada [34]. Except for narrow coastal bands, it 140 belongs to the balsam fir -white birch bioclimatic domain [35]. The climate is maritime with 9 abundant precipitation. Wildfire is infrequent; the main natural disturbance is spruce budworm (Choristoneura fumiferana) outbreaks [35]. Approximately 90% of the region is covered by forests and 80% of these are on public lands [34]. Forest harvesting and the associated extensive road network are the main proximate human disturbances to date. Only 34% of the Gaspésie 145 natural region is free from measurable human footprint [34].
An existing protected areas network covers 5.5% (1371 km²) of the region (Fig. 1). This falls short of the target of an ecologically representative 12%, set by the government of Québec The Gaspésie National Park (802 km²) is currently the largest protected area in the region [37]. This park helps conserve 42 endangered and vulnerable species of plants and animals [34], including the Atlantic-Gaspésie caribou population. Because of its designated status as endangered [38], its capacity to move widely [39] and the vulnerability of its habitat [40], we 155 chose it as our focal population to define functional connectivity. The population is designated because of its small sizecurrently ~70 individuals [32] and the observed long-term decline [41]. These caribou rely on alpine tundra [42,43] as a predator-free refuge against coyotes and black bears [44], a habitat possibly threatened by climate change [45,46]. Caribou also use midelevation, old fir forests during winter [43,44]. The fir forests around the park are affected by 160 forest harvesting, and also possibly by climate change [47].

Building protected areas networks
We adapted the BEACONs (Boreal Ecosystems Analysis for Conservation Networks) approach 165 to build protected areas networks [48,49]. This approach uses mapped hydrological catchments, rather than the cells of an arbitrary grid or other tessellation, as the basic spatial units. This allows the design of protected areas that simultaneously achieve terrestrial and hydrological connectivity [50,51]. Although functional connectivity for aquatic species was not a major concern in our study area [34], it is increasingly recognized as a desirable goal in general [52]. 170 Of more importance to the present application is that the construction does impose some limits on the shape and size of protected areas because of the size of catchments and the method of stream network traversal. However, the size of the catchments was small, implying high spatial resolution to the design. Candidate protected areas are assembled as contiguous, hydrologically connected sets of catchments satisfying minimum size and intactness criteria. The construction 175 proceeds from an initial "seed" catchment, by a modified breadth-first traversal of the stream network. Catchments are added as they are encountered, provided they satisfy a catchment-level intactness criteria. The process continues until the target size is achieved, or all pathways are blocked by headwater or non-intact catchments.
We used a custom catchment layer for Gaspésie defined from the 1:50,000 scale National 180 Hydrological Network (Version 14 NRCAN 2011-2014) and the Canada Digital Elevation Model (NRCAN 2001-2012). The catchment's average size was 2.3 km² (SD = 1.3 km²). We defined catchment intactness using existing 250 x 250m raster maps of the human footprint in Gaspésie [34]. Six disturbance types were defined: forestry activities, roads and trails, agriculture, power-line rights-of-way, urban areas and "other". We gave each type equal weight. 185 As only the first two types are widespread in our study area, we are primarily basing intactness on the absence of recent forest harvesting, roads and trails, We defined a binary 250 x 250m intactness raster, taking values of 0 where no disturbance of any type was present, and 1 otherwise. Catchment-level intactness was calculated as the mean intactness over all raster cells within a catchment, or equivalently, as the proportion 190 of intact cells. We defined our catchment-level intactness criteria as the median catchment intactness over all catchments on public lands, outside of existing protected areas (= 0.026). In order to create feasible designs given land ownership in Gaspésie, we set the intactness of catchments completely overlapping private lands [34] to zero so that they were excluded from the construction. 195 We prioritize the inclusion of headwater catchments in protected areas, because of their hydrological importance and to reduce the potential for upstream river contamination inside protected areas [53]. Accordingly, we used all intact headwater catchments as seeds. We used the mean size of the new protected areas proposed by the MELCC (Fig. 1; mean = 85.5 km², SD = 61. 4 km²; [34]) as the target size. The construction adds entire catchments one at a time, so the 200 final size may exceed this target. Because the mean catchment area is small (2.3 km²) relative to the target size, we consider such discrepancies negligible. We then applied the candidate protected area construction to each seed catchment, with target size and catchment intactness criteria as defined above. The procedure returns a list of the catchments selected, their total area, and their area-weighted mean intactness. We did not apply any protected area-level intactness 205 criteria to sample candidate protected areas as only intact catchments were used to construct them. To make our size criteria more comparable to that used by the MELCC [34], we selected candidate protected areas with a total area above 31.2 km², the size of their smallest proposed new protected area.
We used the set of selected candidate protected areas to generate a sample of 500,000 210 candidate protected areas networks. Each candidate network was constructed by adding a random sequence of candidate protected areass to the existing network, until the summed area of unique catchments exceeded the MELCC area target of 3080 km². The candidate protected areas included within a given candidate network may overlap with each other or with elements of the existing network, but this does not affect the candidate network area. It is easy to show that the 215 properties of hydrological connectivity and intactness are conserved by aggregating overlapping candidate protected areas. Thus, the only material consequence of such overlap is that spatially disjunct network elements may exceed the minimum size criteria, and may be fewer than expected based on the number of candidate protected areas that were added. However, other things being equal, larger protected areas are preferred [54], and the number of spatially 220 separated components is not one of our network design criteria, so we do not consider this overlap of relevance to the present study.

Network ecological representativeness
Ecological representativeness is measured at the candidate network level [48], in terms of four 225 environmental attributes: elevation, surficial deposit, drainage class and potential vegetation type [34]. Potential vegetation was defined by the Ministère des Forêts, de la Faune et des Parcs du Québec (MFFP) as the vegetation present on a site or potentially present, in the absence of disturbance [36]. Surficial deposit, drainage class and potential vegetation type had been used for a gap analysis conducted to inform the proposed expansion of the existing protected areas 230 13 network [34]. The four attributes define habitats independent of human activities such as harvesting history. All four attributes were available for the entire study area at adequate spatial resolution (source: MELCC). Catchment-level attributes were calculated using standard spatial data operations.
We measure candidate network representativeness using non-parametric, two-sample 235 univariate dissimilarity measures. For each attribute, we obtain its distributions over all catchments within the candidate network, and over all other catchments within the study region.
For continuous attributes (e.g., elevation), we calculate a two-sample Kolmogorov-Smirnov statistic. For categorical attributes (e.g. surficial deposit, drainage class, and potential vegetation type), we calculate the Bray-Curtis statistic. These statistics range from 0 to 1, and are 0 only 240 when the two distributions are identical. That is, they measure a candidate network's deviation from perfect representation with respect to an attribute. The univariate statistics for multiple attributes are combined into a univariate distance metric by calculating a Euclidean norm [48].
We took the inverse of this distance metric as the ecological representativeness score for each candidate network, so that a larger score indicates a greater degree of representativeness. 245 We also calculated the representativeness score of the proposed MELCC network. To do this, the protected areas within the MELCC network were approximated to the resolution of hydrological catchments. Because these catchments were relatively small, we assumed approximation errors to be negligible. We note that this measure of ecological representativeness does not depend on the specific choice of ecological attributes or on the specifics of our design 250 methodology; it can be applied to any existing protected areas network using whatever ecological attributes are or interest.

Network functional connectivity
We defined candidate network functional connectivity for the Atlantic-Gaspésie caribou 255 population by simulating, using a spatially explicit individual-based movement model adapted climate scenarios were relatively small. We considered evaluating differences between them to 15 be unimportant relative to our core message concerning the inclusion of functional connectivity in a protected areas network design. Accordingly, and because no climate change scenario was defined as more likely than another, we used the mean over the four climate change scenarios to represent functional connectivity under future climate conditions. We calculated the functional 280 connectivity of the MELCC network under current and future conditions, in the same way.

Identifying priority conservation areas
By plotting the distributions of indicators (i.e., ecological representativeness, current and future functional connectivity), the outlines of the tradeoff surface may be visualized. To identify one 285 location on this surface, we incrementally decreased the quantile Q from 1 until at least 500 networks had all three indicators above the corresponding sample quantiles. This identified a subsample of 0.1% of the candidate networks that were highly ranked under all three criteria.
This represents a tradeoff which values all three attributes equally. We then extracted the individual candidate protected areas contained within the subsample of candidate networks and 290 calculated the proportional frequency of their inclusion in the subsample. We kept the candidate protected areas the most frequently selected by using as the frequency threshold the inflexion point in the proportional frequency curve when ranking them from most to least selected in the subsamples networks. Then we mapped the candidate protected areas within the study region color-coded by this proportion. Frequently selected candidate protected areas may be interpreted 295 as suggesting priority conservation areas for inclusion in a new or expanded protected areas network [28].

Protected areas networks 300
We constructed 690 unique candidate protected areas satisfying our size and intactness criteria.
The proposed MELCC network (Fig 1) had an ecological representativeness score of 4.65 ( Fig. 2a), well below our sample mean. Its current and future functional connectivity were 4.59e +05 and 4.76e +05 , respectively, both well above our sample means (Fig. 2a and 2b). 315 However, our methodology yields many candidate networks that surpass MELCC in both representativeness and functional connectivity (Figs. 2a and 2b). We conclude that the MELCC design is not optimal with respect to these three indicators of network effectiveness. # Figure 2 approximately here # 320

Identifying priority conservation areas
The quantile Q=0.925 (92.5 %-ile) yielded a subsample of 501 candidate protected areas scoring above the corresponding sample quantiles of all three indicators simultaneously (Figs 2a and 2b).
The ecological representativeness quantile was 5.68. The quantiles for current and future functional connectivity were 4.53e +05 and 4.73e +05 , respectively. All 501 selected designs had 325 higher ecological representativeness scores than did the MELCC design. Most of these 501 designs did not exceed MELCCs functional connectivity scores, but some did (Figs. 2a and 2b; 78/501 selected designs had both current and future functional connectivity higher than did the MELCC design). That is to say, we identified some designs for the given tradeoff that outperformed MELCC's in all three indicators. 330 All but one of the 690 candidate protected areas were included in at least one of the 501 selected networks. However, candidate protected area selection frequency within these networks was highly skewed (Fig 3). For example, only 32/689 candidate protected areas were included in more than 35 of the 501 selected networks (Fig 3), where 35/501 represents the inflexion point in the curve of selection frequency for the candidate protected areas (Fig. 3). These 32 frequently 335 selected candidate protected areas were fairly widely distributed over the study region (Fig. 4), but with some evidence of spatial clustering. Notably, several of the most frequently selected are adjacent to the important caribou breeding habitats in the GNP (Fig 4). with tourism if protected areas act as parks [60], or any of the constraints included by the MELCC in their scenario. Target features and constraints can be defined by local managers to help meet local biodiversity goals, and could easily be implemented in our method to improve the design of regional protected areas networks.

Network ecological representativeness
Low ecological representativeness may compromise the effectiveness or efficiency of protected areas networks. Ideally, managers would work towards achieving a better representation of the regional biodiversity when implementing new protected areas. However, to achieve biodiversity protection, protected areas networks should aim to represent the landscape in its pre-industrial 365 form and not to protect habitat types resulting from human disturbances. Gap analyses are a common tool to quantify ecological representativeness of extant networks and to identify underrepresented features [e.g. 8,61]. Here, we adapted an alternate two-stage approach that was recently developed in Canada [48,49] to support systematic conservation planning in the boreal region. The first stage of that method used geospatial analysis tools to construct potential or 370 candidate protected areas. Sites were then assembled into networks and a multivariate distribution matching methodology indicated networks which minimized gaps with respect to a chosen set of environmental covariates. We attempted to address some of the potential shortcomings of such a coarse filter approach by ammeding a fine filter approach that accommodates very different criteria, which is builds in an almost orthogonal dimension to the 375 conservation problem. By doing this, we can identify tradeoffs and understand where such potential shortcomings of representativeness can be augmented to current conservation reality.

Network functional connectivity
Many studies have included connectivity in conservation planning using measures of distance or 380 the cost of movement between protected areas [62] or applying graph and circuit theory techniques [63,64]. A novel feature of this study was that we derived network functional connectivity measures from movement simulations using an individual-based model [33,55] that reproduced known characteristics of the focal population within the study region. The model accounted for process-based complex movement behaviors (e.g., seasonal site fidelity, foray loop 385 movement) that would be difficult to represent in static habitat-based models. Being processbased, the models reflect behavioral responses to environmental conditions and should be 20 preferred to make predictions under future conditions compared to models based on only observed habitat where the suppositions underlying current empirical relationships may no longer hold [65]. 390

Accounting for climate change
It is more efficient to take future climate change impacts into account now instead of waiting for the changes to occur and then reacting [27]. To account for climate change, we defined the future functional connectivity criteria using simulated movements on hypothetical landscapes resulting 395 from different climate change scenarios. Network effectiveness will likely depend on the realized landscape outcome which is presently unknown. Therefore, networks selected using an average functional connectivity over a wide range of climate change impacts would be sub-optimal for any particular climate change scenario. However, they might perform reasonably well under a large range of possible future environmental conditions. The differences in movement predicted 400 for the different climate change scenarios resulted from the assumptions made about the future environmental conditions. We used this relatively simple approach as the functional connectivity differences among the scenarios were not particularly large. This suggests that using coupled climate change and vegetation change dynamic models may not dramatically improve the choice of protected areas, in this region. Furthermore, since no climate sensitive model of vegetation 405 dynamics was available for Gaspésie natural region, using the scenario approach with averaging of results seemed to us a reasonable compromise to account for uncertainty in future conditions. Our approach could readily be adapted to cases with much greater divergence among expected future conditions, at the price of increasing the dimensionality of the tradeoff surface.

21
Management decisions in that case would ideally be based on weighting the relative costs of 410 possible solutions and the probabilities associated with each alternate outcome.

Case study
The lack of functional connectivity during the last 20 years between the different major mountain summits of the Gaspésie National Park has led to a division of the caribou population into two 415 sub-populations genetically distinct and is now jeopardizing the persistence of this isolated population [66]. Consequently, using an approach that could identify key-elements to preserve in order to maintain functional connectivity could benefit to the conservation of an endangered caribou population, especially under a changing climate.
In the Gaspésie region, there is only one caribou population occurring primarily inside the 420 Gaspésie National Park [32,43]. In consequence, in our analyses, networks having high functional connectivity tended to include many protected areas adjacent or near to this park. On the other hand, to achieve high ecological representativeness, protected areas needed to be more or less evenly distributed over the entire region to capture habitat diversity. There was therefore a negative relationship between representation and connectivity in our case study, as others have 425 found [e.g. 67]. It is, in general, not possible to simultaneously optimize for multiple design criteria. Our sampling-based design methodology allows to approximate a multidimensional tradeoff surface. Given the management decision about the choice of tradeoff, we can identify a large number of solutions that are close to this point. Our methodology showed how to design new protected areas for a regional network that satisfy multiple, divergent criteria to (nearly) the 430 highest degree possible.

22
The MELCC network expansion proposed by the Québec government to achieve the 12% coverage target that we used as a comparison was different than the sample of candidate protected areas suggested from the selected best networks identified here (Fig. 4). The MELCC network achieved a lower ecological representativeness than our selected networks. However, 435 they had to respect design criteria and constraints, like socio-economic issues or the inclusion of rare ecosystems that we did not consider. This could explain the sub-optimal ecological representativeness achieved by their network. Regarding functional connectivity, their network seems well connected from the Atlantic-Gaspésie caribou point of view. The performance of their scenario is surprisingly quite good for this feature, considering it was not explicitly part of 440 their design. The networks we selected may not achieve such a high functional connectivity because our specific choice of intactness criteria resulted in some areas being excluded from any of our network solutions. In particular, we excluded some non-intact areas close to the Gaspésie National Park that were included in the MELCC design.