Targeting climate diversity in conservation planning to build resilience to climate change

. Climate change is raising challenging concerns for systematic conservation planning. Are methods based on the current spatial patterns of biodiversity effective given long-term climate change? Some conservation scientists argue that planning should focus on protecting the abiotic diversity in the landscape, which drives patterns of biological diversity, rather than focusing on the distribution of focal species, which shift in response to climate change. Climate is one important abiotic driver of biodiversity patterns, as different climates host different biological communities and genetic pools. We propose conservation networks that capture the full range of climatic diversity in a region will improve the resilience of biotic communities to climate change compared to networks that do not. In this study we used historical and future hydro-climate projections from the high resolution Basin Characterization Model to explore the utility of directly targeting climatic diversity in planning. Using the spatial planning tool, Marxan, we designed conservation networks to capture the diversity of climate types, at the regional and sub-regional scale, and compared them to networks we designed to capture the diversity of vegetation types. By focusing on the Conservation Lands Network (CLN) of the San Francisco Bay Area as a real-world case study, we compared the potential resilience of networks by examining two factors: the range of climate space captured, and climatic stability to 18 future climates, reflecting different emission scenarios and global climate models. We found that the climate-based network planned at the sub-regional scale captured a greater range of climate space and showed higher climatic stability than the vegetation and regional based-networks. At the same time, differences among network scenarios are small relative to the variance in climate stability across global climate models. Across different projected futures, topographically heterogeneous areas consistently show greater climate stability than homogenous areas. The analysis suggests that utilizing high-resolution climate and hydrological data in conservation planning improves the likely resilience of biodiversity to climate change. We used these analyses to suggest new conservation priorities for the San Francisco Bay Area.


INTRODUCTION
Biological communities are already changing in response to climate change. Climate is a primary driver of species distributions; as climate warms, populations are shifting in abundance, demographic traits, and distribution (Parmesan 2006, Moritz et al. 2008, Bellard et al. 2012, Rapacciuolo et al. 2014). These shifts are raising questions about the adequacy of current protected area strategies to conserve biological diversity (Peters and Darling 1986, Noss 2001, Pressey et al. 2007, Heller and Zavaleta 2009. There is need for methods to determine where to invest in land acquisition given directional climatic changes (Hunter et al. 1988, Hannah et al. 2002, Hannah et al. 2007, Klausmeyer and Shaw 2009, Iwamura et al. 2010, Schloss et al. 2011, and, how to adapt management plans in the face of changing demography, species interactions and community composition (Millar et al. 2007, Cross et al. 2013.
The protection of diverse physical conditions on the landscape, taking into consideration the effects of topography and geologic structure on the cycling of energy and water, is one approach to climate adaptation planning (Ackerly et al. 2010, Anderson and Ferree 2010, Beier and Brost 2010, Schloss et al. 2011. This is an alternative or complementary approach to species distribution modeling (SDM), which is commonly used in climate adaptation planning (e.g., Hannah et al. 2002, Araú jo et al. 2005, Hannah et al. 2007), yet has some widely acknowledged limitations. Reliable application of SDM methods requires knowledge of species distributions, disequilibrium dynamics, fundamental niches, inter-species interactions, and other data that is key for making accurate projections, but is often not available (Araú jo et al. 2005, Gallagher et al. 2010, Veloz et al. 2012, Svenning and Sandel 2013, Gavin et al. 2014). SDM may not work well when there are many species being considered, which will show a wide diversity of responses to climatic change. SDM is also sensitive to variation in algorithms, future scenarios, the spatial resolutions of data inputs, and life history traits (e.g., Araú jo et al. 2004, Kueppers et al. 2005, Elith and Graham 2009, Franklin et al. 2013.
Collectively these issues make it difficult to interpret and apply results with confidence.
In contrast, geophysical features like soil, topography, and solar insolation interact with climate to structure the local abiotic environment in relatively predictable ways (Weiss et al. 1988, Dobrowski 2011. Gradients of abiotic features drive ecological and evolutionary processes that shape patterns of biological diversity (Lawler et al. 2015). Species composition turns over in time in response to changing climatic and biotic conditions, but the underlying physical features of the landscape, and the resulting climate gradients that emerge, are likely to endure. Thus the protection of diverse physical settings suggests a strategy to protect the ''arenas'' for biodiversity generation and maintenance that is robust in the face of uncertainty (Hunter et al. 1988, Rouget et al. 2003, Pressey et al. 2007).
There are a number of ways that local climatic diversity, which will persist even as the macroclimate changes, could affect the resilience of regional biodiversity pools. A conservation lands network that protects patches of relatively cool conditions, and patches where climate changes more slowly than in surrounding areas, would provide refugia for populations experiencing negative impacts from warming (Ashcroft et al. 2012, Gavin et al. 2014. Alternatively, relatively warm or dry patches on the landscape may harbor populations of native species that are range outliers. These outlier populations are likely foci for future expansion as the climate warms (Svenning and Sandel 2013). Protecting foci of favorably adapted native taxa will reduce dispersal distance and may facilitate the persistence of native taxa under future climates.
Furthermore, heterogeneity in climate at multiple spatial scales is expected to diminish vulnerability of the overall biodiversity within a reserve network to climate change (Ackerly et al. 2010). Topoclimatic heterogeneity buffers the impacts of rapid, global climate change, by allowing species populations to better keep pace through dispersal and re-location (Loarie et al. 2009, Willis andBhaghart 2009), and also by creating in situ refugia (Gavin et al. 2014). Heterogeneous landscapes may also support greater functional diversity within communities, v www.esajournals.org 2 April 2015 v Volume 6(4) v Article 65 increasing the likelihood that at least some species will succeed under future climates. Finally, populations spanning diverse climate space are likely to have more morphological and genetic diversity, promoting in situ adaptation (Noss 2001, Moritz 2002, Harris et al. 2006, Millar et al. 2007, Sgrò et al. 2011. In this study, we focus on the protection of coupled hydro-climatic diversity of the landscape for climate adaptation planning. We compared conservation networks designed using different conservation goals (vegetation diversity versus hydro-climate diversity) and measured at two different scales (regional versus sub-regional). We then evaluated the potential resilience of the networks based on two characteristics (climate space captured and climate stability). Climate space is measured as the range of climate conditions in the combined set of planning units in a network. Climate stability is defined as the overlap of historic and future climate conditions for the combined set of planning units in a network. While it isn't possible to evaluate resilience to climate change before it happens, it is possible to measure climate space and climate stability. High diversity of climate conditions (i.e., climate space) within a network provides the capacity for ecosystems to absorb change and re-organize on local scales, so as to retain similar function, structure, identity, and feedbacks. Higher climate stability within a network increases the chances for species to track their habitat niche. We address these questions using the San Francisco Bay Area Conservation Lands Network (CLN) as a case study, and apply our results to inform conservation priorities in this region.

Setting
The San Francisco Bay Area (California, USA) is a topographically and hydrologically complex landscape characterized by both high climate diversity and high biodiversity. The region, as defined here, spans 10 California counties with a total area of about 2 million hectares. Approximately 25% of the region is currently converted to urban, suburban or agricultural uses, and 25% is formally protected as open space, through direct acquisition or conservation easements (Walker 2007). In 2011, the Bay Area Open Space Council's (BAOSC) Upland Habitat Goals project completed a comprehensive assessment of biological diversity relative to the extent of current and potential protected open space. The result is called the Conservation Lands Network (2011 CLN), a prioritization plan for land conservation of upland, terrestrial habitat in the San Francisco Bay Area (Bay Area Open Space Council 2011). The 2011 CLN is currently used by conservation organizations throughout the region for strategic acquisition and management of protected areas.
The 2011 CLN prioritized land for conservation using the spatial planning tool Marxan (Ball et al. 2009). Marxan is a software program designed to identify a conservation network (a portfolio of planning units) that captures multiple conservation goals, while minimizing spatial fragmentation and user-defined costs. The Marxan analysis for the CLN was based on protecting sufficient habitat for biological elements, defined as coarse filter vegetation types together with fine filters of endangered taxa and unique landscape features, such as riparian corridors and serpentine soil types. The 2011 CLN prioritizes approximately 500,000 hectares for conservation, in addition to the approximately 400,000 hectares of already protected land. Climate and climatic change were not explicitly addressed as part of this prioritization plan. We used the core 2011 CLN datasets in conjunction with a set of high resolution historical and projected future hydrological-climate datasets generated from the California Basin Characterization Model (BCM) Flint 2012a, b, Flint et al. 2013) to evaluate the potential resilience of the CLN in the face of climate change.

Marxan
Marxan is a spatial conservation prioritization tool that designs a network through the identification of the best (least-cost) set of planning units needed to satisfy user-defined conservation goals and constraints. For a given study domain divided into planning units, Marxan requires data inputs about the distribution and abundance of features (such as a focal species or habitat types) within those planning units. Conservation targets for features (i.e., the amount to be protected) are specified, as well as costs for planning units that provides inforv www.esajournals.org 3 April 2015 v Volume 6(4) v Article 65 mation about their relative conservation suitability. In running Marxan, there are many parameters that can be varied to structure the final output, including the number of runs, adjustment of the boundary length modifier (BLM) to minimize fragmentation in the spatial pattern of the network, and flexibility in meeting targets (defined as the species penalty factor (SPF)). Marxan output provides a best model and a selection frequency table that shows the number of times each planning unit was selected as part of a solution from all program runs. Planning unit selection frequency indicates how important a planning unit is to achieving the conservation goals in an efficient network.

Scenarios run in Marxan
Our study domain consisted of the nine counties bordering the San Francisco Bay plus Santa Cruz County. The entire area was partitioned into a continuous mesh of 100-hectare hexagonal planning units (N ¼ 20,306). Each planning unit was assigned a cost based on its conservation suitability. Cost (the inverse of conservation suitability) increased with (1) proximity to roads (2) higher population density and (3) higher degree of parcelization. We created distribution maps within the study domain for two feature types: vegetation habitat types (v) and climate isotypes (i ) (see Methods: Model inputs). We varied the scale at which the rarity of these features was estimated, either regionally (Reg) or sub-regionally by landscape units (LU), henceforth referred to as 'planning scale'. We then used the Marxan software program to design four different conservation network scenarios, referred to as iNetLU, vNetLU, iNetReg, vNetReg.
The subregional scale was based on the 2011 CLN landscape units (LU ) (http://www. bayarealands.org/overlays/CLN_LandscapeUnits_ Map.html; see Fig. 6). Landscape units were defined as natural biogeographic units in the region, generally encompassing local mountain ranges or valley bottoms (average size 57,000 ha). They are likely to be the areas within which species dispersal is facilitated, because LUs are often separated from each other by urban or agricultural lands. Defining protection targets for conservation features within LUs spreads the conservation network across the region more evenly than if targets are defined for the region, and it helps ensure that the diversity of habitat types is conserved within the likely dispersal distances of local species pools (CLN 2011). We varied planning scale because we wanted to assess how important this step is to attaining high climate diversity at scales thought to be most relevant to species dispersal.
To design each network scenario, we ran the Marxan software program 20 times, each run with 1 million iterations, using the annealing and iterative improvement algorithm with no heuristic or cost threshold. The species penalty factor (SPF) was set to 5, which was as low as possible to allow Marxan flexibility to find an efficient solution, while ensuring targets were met in full in all four scenarios. The likelihood of any planning unit being selected in a run was only a function of its cost and its portfolio of feature targets. We did not use the boundary length modifier (BLM) constraint in order to maximize the chance for differences, driven by the spatial distribution of features, among networks to emerge. We assigned the hexagonal planning units to each of the four network scenarios: iNetLU, iNetReg, vNetLU, vNetReg based on their frequency of selection, as in the 2011 CLN. We used this method because it takes advantage of information from multiple runs revealing which planning units are highly useful in creating an efficient network. The threshold selection value for inclusion varied for each scenario such that the total number of planning units was similar to the number selected in the best run. For example, if the best run required 10,000 planning units then we selected a threshold of selection frequency (e.g., 16) that resulted in about 10,000 planning units being included in iNet or vNet. On average any planning unit selected greater than 12 times (60% of runs) was assigned to versions of iNet or vNet.

Model inputs
Vegetation habitat types.-Vegetation habitat types, and their rarity, were based on the vegetation map developed in 2011 CLN planning. This 30-m pixel map of 63 land cover types was produced using a combination of different statewide vegetation maps, and also climate, hydrological, soil, and farmland maps, as well as expert opinion (Bay Area Open Space Council v www.esajournals.org 4 April 2015 v Volume 6(4) v Article 65 2011). For example, vegetation types were stratified for serpentine bedrock using the USDA Natural Resource Conservation Service State Soil Geographic Database (STATSGO, http://water.usgs.gov/lookup/getspatial?ussoils). The approximately 400,000 ha of habitat identified as 'Annual Grasslands' by CalVeg (http:// www.fs.fed.us/r5/rsl/projects/classification/) were split into Cool, Moderate, Warm, and Hot categories based on July maximum temperature maps obtained from PRISM (800-m-scale Parameter-elevation Regressions on Independent Slopes Model; Daly et al. 2008). These steps provided a fine-grain map of vegetation types (or habitat types), many of which are defined by the dominant woody species (e.g., Redwood, Blue Oak). The rarest vegetation type 'Monterey Cypress' covers about 34 ha, and the most widespread 'Warm Grasslands' covers 212,240 ha.
Climate isotypes.-Climate isotypes, and estimates of rarity, were constructed from 270-m (;7.3 ha) resolution climate-hydrology data layers for the region Flint 2012b, Flint et al. 2013). The data set includes an ensemble of 18 different projections of global climate models (GCMs) drawn from a range of emission scenarios from CMIP3, as well as different representative concentration pathways from CMIP5 (Coupled Model Intercomparison Project, http:// www-pcmdi.llnl.gov/projects/cmip/). Downscaling to the 800-m scale was accomplished by Thrasher et al. (2013) using the Bias-Correction Spatial Disaggregation algorithm (Wood et al. 2002(Wood et al. , 2004. These data were then spatially downscaled according to Flint and Flint (2012b) and used as inputs to 270-m scale BCM. Data summaries of different climate and hydrological variables were provided as 30-year averages for the historical and recent period covering 1910-2010, and for projected futures ranging from 2010 to 2100. The 18 models were selected to cover the range of projected future climate conditions for the San Francisco Bay Area from a larger subset of 112 models (Fig. 1). Models show increased temperatures ranging from 18 to 68C by the end of the century relative to 1951-1980. Rainfall is highly variable.
Using the recent period (1981-2010), we created climate isotypes by combining the values of three climate/hydrological variables important to vegetation distributions in the region (Cornwell et al. 2012): average winter (December, January, February) minimum temperature (winter Tmin), average summer (June, July, August) maximum temperature (summer Tmax), and annual climatic water deficit (CWD). CWD is a seasonally integrated measure of how much evaporative demand exceeds available water. It is a measure of end of dry season drought stress that reflects the impacts of climate, topography, seasonality, soil porosity and depth on water balance and is directly relevant to vegetation (Stephenson 1990). Because CWD integrates the physical landscape and climate it has been show to correlate to species distributions (Flint and Flint 2012b), vegetation distributions, including in the Bay Area (Cornwell et al. 2012), changes in forest structure (McIntyre et al. 2015) and forest tree mortality (Millar et al. 2012, Das et al. 2013, van Mantgem et al. 2013. We categorized the continuous surfaces of climate/hydrological variables into discrete categories ranked 1 to 4 (breaks for CWD are 700, 800, and 900 mm y À1 , breaks for winter Tmin are 48, 58, and 68C, Breaks for summer Tmax are 248, 278, 308C). Intervals were created in order to partition the variance in the recent climate for each variable, and to create a similar number of categories as the vegetation types. Isotypes were defined by combining the ranks of variables into a 3-digit code. For example region '111' indicates low water deficit (,700 mm), cool winter temperature minimums (,48C), and cool summer temperature maximums (,248C) relative to the 10-county region. Fifty-nine of the 64 possible combinations were observed across the region, with the rarest '421', '412', '311' all occupying one pixel (;7 ha). The most widespread isotype ('323') covered 157,836 ha. To test the sensitivity of our results to other climate classifications we also experimented with k-means clustering, a non-hierarchical clustering approach, and other interval breaks. These different techniques produced networks that captured similar climate space to the isotype version reported here and are not discussed further (see also Torregrosa et al. 2013). However, the results of these alternative runs were used to suggest new priorities in the CLN (see Methods: Prioritizing the CLN ).
Defining targets.-To set targets for feature types, we quantified the abundance of each v www.esajournals.org 5 April 2015 v Volume 6(4) v Article 65 feature in intact lands (not developed) relative to the total land area. This follows the methodology of the CLN 2011. We intersected the vegetation and climate maps with land-use maps (i.e., protected, intact, cultivated, urban, and rural residential, data from Bay Area Open Space 2011) to estimate the abundance of features by land-use type for networks planned at regional scales. To plan networks at sub-regional scales, we further intersected maps with LU.
Features were then ranked from 1 to 3 based on their rarity in intact lands. For Marxan, we assigned rare features higher proportional targets for conservation than common features. For rank 1 features, found in less than 1% of intact land, targets were to attain 90% of remaining intact land with that feature. For rank 2 features, found on less than 5% of intact land, targets were to attain 75% of remaining intact land with that feature, and for rank 3, found on more than 5% of intact lands, targets were to attain 50% of intact land with that feature. There were no targets set for the 10 types of non-native vegetation habitat types (e.g., Eucalyptus forest). This scheme resulted in the identification of 53 vegetation types with a combined target of 696,852 ha for conservation, and 59 climate isotypes with a combined target of 820,572 ha for the regional network scenarios. For LU network scenarios, each vegetation or climate type in each LU is a unique target, resulting in 594 vegetation types with a combined target of 660,474 ha, and 639 climate isotypes with a combined target of 740,746 ha.
Network similarity and climate space in current climate.-To compare the similarity of networks, we calculated the Spearman rank correlation of the selection-frequency output for each planning unit between different Marxan scenarios. Scenarios that result in a similar set of planning units selected in efficient solutions will show a high correlation. To explore potential resilience, we examined the climate space captured in each network by plotting the distribution of CWD, summer Tmax, and winter Tmin using frequency counts and density plots of occurrences of each value. To assess the data in two-dimensions we made contour plots of frequency counts of occurrences of x, y pairs of values. We evaluated the climate space of each network for the region at large, and for each LU independently.
Climate stability.-To further test whether the four different network scenarios varied in potential resilience to climate change, we calculated a measure of climate stability similar to Iwamura et al. (2010) and Watson et al. (2013). Climate stability was defined as the overlap in climate space within each network between current and future climates. We assume connected networks with greater overlap in climate space between time periods would allow plants and animals to track suitable climates within the network and thus offer more potential resilience to climate change. For each model network, we calculated a convex hull volume of the climate space in three dimensions (CWD, summer Tmax, winter Tmin) for the recent climate  and then for each of 18 possible futures for three time slices, 2010-2039, 2040-2069, 2070-2099. The hull volume was calculated for the regional network as a whole and for each landscape unit independently. We then calculated climate stability as the overlap between convex hull volumes, using the Q hull program in R ( Barber et al. 1996) according to Eq. 1: where HV stands for hull volume, and c, f, and u denote current, future, and the union of current and future climate space data (further detail and code for calculating climate stability is documented in the Supplement). If current and future hulls occupy the same climate space, then complete overlap will occur and result in 100% stability. If current and future climate hulls have some but less than complete overlap, the hull volume of their union will be less than the sum of their individual hulls, because the hull volume of the union is partially overlapping, leading to stability values greater than zero but less than 100% overlap. If current and future climate hulls have no overlap, the hull volume of their union will be greater than the sum of their individual hull volumes (because it includes the intervening climate space), leading to a negative stability value. Thus stability can range from one, to less than zero, and is unbounded on the negative end, indicating a complete and potentially increasing departure from current conditions. Assuming that species can disperse within the network, those inhabiting networks with higher climate stability presumably would be more able to shift their distribution, and therefore less vulnerable to climate change than those with low climate stability.
We calculated climate stability for the region as one single domain and for each LU separately. We compared the climate stability for each future time period for the four network scenarios using pair wise Wilcoxon signed rank tests. Each of the 18 GCMs was used as a separate observation within each scenario. We also compared the stability of networks across the 18 GCMs to evaluate how well the networks preform in different futures, and among LUs to evaluate how topography affects potential resilience.
Prioritizing the CLN.-We used iNetLU models to recommend priorities within the existing 2011 CLN, and to highlight new areas for conservation from the perspective of climate diversity. We searched for the set of planning units that was consistently selected in climate-based networks, and including those that used the alternative climate categorization schemes and that varied in representation goals. For each planning unit, we summed the number of times a planning unit was selected across all runs from the three different climate characterization schemes, at three-levels of conservation area goals (full as defined in 2011 CLN, 50% smaller, 75% smaller) (n ¼ 9). Climate priority spots were defined as those planning units with a combined score of 153 or higher, meaning they were almost always selected (.80) in all nine Marxan scenarios. We intersected climate priority spots with the 2011 CLN priority classifications.

Network similarity and climate space in current climate
Different network scenarios were similar in size and distribution and all were strongly correlated (Fig. 2, Table 1). However, the different scenarios did not do a good job at meeting each other's targets (Table 2). iNet models did not thoroughly capture vegetation targets, nor did vNet models capture climate targets. Networks v www.esajournals.org 7 April 2015 v Volume 6(4) v Article 65 planned on the same feature type but at different planning scales were better at meeting each other's targets than networks planned on different feature types at the same planning scale. In particular, iNet models performed poorly at capturing the rare vegetation features and vice versa (Table 2).
A few other general patterns emerged from comparing networks. Networks planned at the sub-regional scale were more spread and frag-mented (higher boundary area) relative to regional networks. Networks planned at regional scale were more concentrated in the northern part of the study domain and toward the coast relative to the sub-regional networks indicating there is more intact land with rare climate in the northern Bay Area relative to the other parts of the region. Networks designed on climate isotypes were larger than those designed on vegetation types (iNetLU ¼ 8392 planning units, Fig. 2. Selection frequency of planning units for the four network scenarios in the San Francisco Bay Area region: the climate-based network planned at the landscape unit scale (iNetLU), vegetation-based planned at the landscape unit scale (vNetLU), climate-based planned at the regional scale (iNetReg) and vegetation-based planned at regional scale (vNetReg). Color shows the number of times a planning unit was selected out of 20 Marxan runs. Planning units shown in red and orange were included as part of each network solution. Some landscape units are white, or mostly white, because they are predominately urban.
v www.esajournals.org 8 April 2015 v Volume 6(4) v Article 65 vNetLU ¼ 7674, iNetReg ¼ 8527, vNetReg ¼ 7321). Differences in network size reflect differences in the total area targeted for conservation, which was based on the rarity of types in intact lands.
When considering the region as a single climate domain, the four network scenarios capture a similar range of climate space relative to each other, and similar to the range of available climate (Fig. 3). Contour plots of climate variables in two-dimensions showed no strong differences in shape from each other (not shown). However, when we evaluated climate for each landscape unit separately, we found that iNetLU captured a significantly greater range and variance of some climate variables than the other scenarios (Table 1). We observed no differences between iNetReg and vNetReg, and they are both less variable than LU based models. Comparing the two LU based network scenarios, we found that individual landscape units in the iNetLU were, on average, two and half times more likely to show a greater range in winter Tmin and CWD than in vNetLU, but did not differ for summer Tmax. This did not result in visual differences in the distribution of climate Network listed first has greater variance for significant comparisons than network listed second. à Wilcoxon signed rank test. * P , 0.05. ** P , 0.01, *** P , 0.001. Non-significant P-values appear in parentheses.
§ The network scenario is named by the feature type targeted (I ¼ isotype, V ¼ vegetation) and the planning scale used (LU ¼ Landscape Unit, Reg ¼ Regional). Note: The network scenarios do not meet targets in full even for the features they were designed to capture because they were defined as the subset of planning units that were most frequently selected across multiple runs (i.e., 60%). The rarity rank 3 features are especially poorly captured in the networks designed to capture those feature targets because these are very common features. As a result they can be captured in many places on the landscape, and so there is high variability across runs in Marxan in which planning units are selected to capture these common features. Thus because we designed network scenarios by selection frequency, we do not include enough planning units that have common features to met the targets in full. All individual solutions from Marxan runs met targets in full.
The network scenario is named by the feature type targeted (I ¼ isotype, V ¼ vegetation) and the planning scale used (LU ¼ Landscape Unit, Reg ¼ Regional).
v www.esajournals.org 9 April 2015 v Volume 6(4) v Article 65 variables however, as for most LUs iNetLU and vNetLU overlapped entirely, except in a handful of cases (shown in Fig. 4).

Climate stability
When climate stability was assessed for the regional network as a whole, networks showed similar climate stability over time on average (Fig. 5). In pairwise comparisons of networks within the same GCM and same time period, vNetLU and iNetLU were not significantly different, and were significantly more stable than vNetReg and iNetReg (p , 0.05, n ¼ 18). When climate stability was assessed at the scale of individual LUs, iNetLU showed consistently significantly higher climate stability than all other networks for pairwise comparisons in all time periods (p , 0.01, n ¼ 522 (18 GCM 3 29 LU); Fig. 5). vNetLU and iNetReg were not different, but were significantly more stable than vNetReg in all comparisons (p , 0.01). Climate stability is lower when networks are evaluated individually for each LU as compared to the region at large (Fig. 5). This is because networks at the LU scale are smaller, and in many cases stability values are below zero, indicating no overlap between present and future climate. While iNetLU showed the highest stability overall, network scenario type makes more of difference in some landscape units than others. iNetLU performed significantly better than vNetLU in 51% of landscape units, whereas vNetLU showed significantly higher stability in 17%. In 31% of LUs the networks were not statistically different from each other (Fig. 6). The higher stability of iNetLU models compared to vNetLU models is particularly large in North Contra Costa Valley, Solano Plains, and Tri-Valley (Fig. 6). These are the same landscape units that also showed strong differences in the distributions of some climate variables (Fig. 4).
The differences in climate stability among network scenarios are relatively small, however, compared to the differences in climate stability that are driven by the choice of GCM/concentration pathway (Fig. 7), or by the variation among landscape units (Fig. 6). More topographically complex (higher standard deviation of elevation) landscape units show higher climate stability (n ¼ 29, F ¼ 70.9, r 2 ¼ 0.72, p , 0.0001; Fig. 8). The relative climate stability of landscape units did not change over time or across GCMs, indicating that the relationship between topographic heterogeneity and climate stability is robust to uncertainty in climatic change.

Prioritizing the CLN
We found 716 priority spots of which 591 were already in the CLN, and 117 were not (Fig. 9). Priority spots emerged in all but one landscape unit (Santa Clara Valley). The highest numbers of priority spots (18) were found in the Sonoma Coast Range. Planning units in LUs farther from the coast, and with less elevational relief (e.g., Solono Plains, Conta Costa Delta), tended to be under-represented in scenarios planned at the regional scale relative to the sub-regional scale, indicating these LUs have less regionally rare climate intact habitat. These same under-represented LUs also tended to show less climate stability and fewer climate priority spots. In contrast, larger, more topographically heterogeneous landscape units, such as Mount Hamilton and Sonoma Coast Range, scored high for the number of climate priority spots, climate stability, and regional climate rarity.

DISCUSSION
Maintaining climate diversity in conservation networks is expected to improve the resilience of regional species pools to climate change. We designed a method to represent climate diversity in conservation planning by incorporating climate isotypes as conservation targets. We then tested whether conservation networks planned to represent the diversity of climate isotypes resulted in networks that captured a greater range of climate space and are more stable to projected climate changes compared to networks planned on representing the diversity of vegetation features alone. Our results suggest that at regional scales, networks based on vegetation features were successful at capturing most of the diversity of climate space, and showed climate stability similar to networks based on climate type diversity. This shows that the vegetation types used in planning the 2011 CLN effectively stratify a majority of climate across the region, as expected given the strong effects of climate on vegetation distributions (Cornwell et al. 2012).
However, at sub-regional scales (landscape unit), which are thought to be most relevant to species dispersal and meta-population dynamics, feature type made more of a difference. Networks planned on climate isotypes captured a greater range and variance in climate variables, and showed higher climatic stability to projected climate changes than networks planned on vegetation types. Our results indicate that bringing climate data into conservation planning is useful for improving potential resilience of biodiversity at spatial scales relevant to species dispersal and population persistence.

Bio-climate feature types
In this analysis, we used fine-grain maps for both climate/hydrology and vegetation features. Both feature types integrate fine-scale topographic and soil gradients on the landscape, which drive vegetation and climate heterogeneity. If the vegetation layer used in this analysis were not partially defined by climate, we would not expect it to capture climate variance as successfully. Similarly, if we used climate layers at resolutions more typical of downscaled global climate model projections (i.e., 1-12 km), and that did not incorporate hydrology, we would expect to find less overlap between vegetation and climatebased networks, both in network design and climate space representation. This is because the topoclimate gradients, which emerge at scales below 1 km, drive beta and alpha species diversity (Ackerly et al. 2010). Therefore, climate maps averaged over large grid cells do not reflect the biotic community turnover present on the landscape.
However, climate isotypes did not function as surrogates for vegetation types, because iNet models failed to sufficiently capture targets for rare vegetation types. Thus climate and vegetation based approaches may capture similar patterns at a coarse scale, but cannot substitute for each other to capture rare or localized features. Taken together these results suggest using a bio-climate approach to defining conservation features is important for climate change resilience planning, as has been suggested elsewhere (Pyke and Fischer 2005). A bio-climate approach is illustrated in the treatment of grasslands in the 2011 CLN vegetation map: temperature data were used to sub-divide annual grasslands into four categories, ranging from cool to hot. Conservation land networks formulated from targeting fine filters alone (i.e., the distribution of an endangered species), or coarse filters that are too broad (coarsely defined vegetation types, such as 'Annual Grasslands'), are far less likely to represent climate diversity. Schloss et al. (2011) reached similar conclusions when they compared conservation networks planned on abiotic features to those planned on biotic features. They found networks planned on abiotic targets revealed a pattern of prioritization that sufficiently captured coarse filter vegetation targets, but not fine filter species targets. They conclude that targeting abiotic features in conservation planning will help to create a more robust network to changing climate, but is not a substitute for planning on biotic features, especially rare features. They did Fig. 7. Climate stability at the regional level shown for iNetLU (acronym is as in Fig. 2) for each of 18 different GCM scenarios at three future time periods. Climate stability decreases over the century in all models, but the declines are most stark in the scenarios with high relative concentration pathways from CMIP5. The same result is seen for the other three conservation network scenarios. Targeting vegetation types, even fine-grain types as used in the analysis here, did not always work to capture climate. Twenty percent of the landscape units (6 out of 29) in vNetLU diverged markedly in the distribution of climate capture compared to iNetLU (Fig. 4) and four of these six LUs in vNetLU (Solano Plains, Santa Clara Valley, North Conta Costa Valley, Tri-Valley) also showed particularly large declines in climatic stability to future climate projections relative to iNetLU (Fig. 6). This suggests the distribution of vegetation types in these landscape units are strongly influenced by factors other than, or in addition to, climate; so just selecting on vegetation did not sufficiently stratify climate space. Therefore targeting climate directly in conservation planning will add more new information in some locations compared to others.
Differences between iNet and vNet emerge in some cases because rare climate types, those representing less than one percent of intact lands, are hosting common vegetation types and vice versa. For example, in Contra Costa Delta, iNetLU models identified patches of hot grasslands that showed low CWD values relative to the other hot grassland patches in the region. Those rare isotypes fell in planning units that were predominately urban and agricultural, and thus the conservation suitability was low (i.e., cost factor was high). vNetLU did not select these patches, because hot grasslands could be captured for less cost elsewhere. However, the iNetLU solution consistently selected these planning units regardless of cost because these were the only places to satisfy goals for these isotypes. It would be interesting to investigate the biological communities in these rare isotypes to determine if their species composition is different than other patches of the same vegetation type that are within a more common isotype. The rare isotypes may be important for harboring unique biological communities and warrant additional investigation, or alternatively they may become important biologically as climate change progresses, as small differences in climate conditions could be vital to the survival of individual species. Alternatively, they may be highly degraded and have minimal conservation value.

Climate stability
A robust finding in this work is that landscape units with low topographic heterogeneity do not show as much climate stability as landscape units with high topographic heterogeneity. This finding is robust to all climate change scenarios and was not particularly sensitive to conservation network design. In other words, for flat landscape units, any network design was highly exposed to climatic changes (Fig. 8). Topographic complexity provides spatial buffering to the impacts of climatic change on populations. When landscapes are topographically rough the shifts in climate occur over shorter distances (Loarie et al. 2009), and there is an increased chance of some overlap between current climate and future climate within a landscape unit, and thus higher climatic stability. This indicates that conservation investments in areas with high topographic complexity are a good bet, as these areas are inherently more resilient to climate change (see Ackerly 2012 for a similar analysis at a statewide level). At the same time, the inherent vulnerability of communities living in flat areas suggests that it will be all the more important to conserve large areas. Implementing connectivity planning across landscape units, for example from flat areas to neighboring mountainous areas, would also increase the ability of native species to track climate throughout a region. In general, creating habitat connectivity is essential to any climate adaptation strategy Zavaleta 2009, Schloss et al. 2011). In the San Francisco Bay Area opportunities to connect protected areas via conservation, especially east to west, may require establishing movement pathways through the heavily developed valley floors, which will be quite difficult. This points to the importance of greater attention to restoration and creation of habitat in developed areas (e.g., taking advantage of opportunities in backyards, urban parks, business parks, agricultural lands as an integral part of conservation practice) as part of climate smart conservation.
The climate stability analysis also illustrates that efforts to mitigate carbon dioxide emissions will have a significant impact on species exposure, more so than conservation network design. Across the 18 future scenarios considered here, the differences in climate stability between scenarios of low emissions (RCP 2.6) and high emissions (RCP 8.5) overwhelm differences between networks designed on different feature types (Figs. 5,6). Emission pathway and climate model strongly affect exposure to climate change (Fig. 7). Scenarios like GISS RCP 2.6 show high climate stability throughout the century because there is little precipitation and temperature change relative to the recent past, whereas the MIROC ESM RCP 8.5 shows negative values for stability because there is a large increase in temperature and a decrease in precipitation projected (Fig. 1). The ensemble average is most similar to GISS AOM A1B and CCSM4 RCP 8.5, which show modestly positive climate stability at the end of the century. This suggests that there will likely be some overlap in climate conditions in the future, which will help populations persist v www.esajournals.org 16 April 2015 v Volume 6(4) v Article 65 if they can move around the conservation network.

Climate targets
In this work, we targeted climate diversity as a function of rarity, without any bias toward particular zones of climate, because we theorized that higher climate diversity supports greater biodiversity. We can think of rationale for prioritizing the protection of both the cool spots in a region (refugia) and the hot spots (potential nuclei for vegetation expansion), which is why we focused on climate diversity. However, we found that the climate stability for the conservation network at the LU scale was mostly below zero. It may be possible to improve climate stability through preferential protection of cool, wet isotypes over hot, dry isotypes. This is because there is a general warming and drying trend for the region (Fig. 1) resulting in an increase in hot, dry isotypes (Torregrosa et al. 2013). For example, when we calculated isotypes for the GFDL A2 scenario for 2070-2099, half of the current isotypes were lost, and there were large gains in the frequency of observations in arid, hot codes (e.g., 444, 434). Some GCMs project increased precipitation in the future, yet, even models with increased precipitation show increased CWD values due to high evaporative demand from warmer temperatures and fundamental limitations of soil water capacity during the dry season (Micheli et al. 2012).
Thus it may be possible to increase climate stability of the conservation network by setting higher representation goals in Marxan for the cool, moist isotypes relative to the dry, hot isotypes. If we had only targeted the cool, moist isotypes in designing networks, we would expect to find less overlap between iNet and vNet in terms of the similarity of prioritization of planning units (Table 1). This is because vNet would be targeting diversity based on current vegetation rarity, and iNet would be targeting diversity based on climates projected to be rare in the future. It is hard to predict how this method would affect the representation of vegetation diversity on today's landscape as climate isotype categories span many vegetation types (on average 25 vegetation types per isotype, with 49 in the most widespread isotype '323').
The approach illustrated here, first, creating climate types, and second, using them as targets in conservation planning to prioritize acquisition could be adjusted in many different ways to fit individual systems. Identifying the extent to which variation in how conservation goals are set (e.g., targeting cool isotypes over hot isotypes) impacts the structure of priority networks or affects the climate stability of the conservation network across the next century would be an important exercise for future work.

Conclusions
Adaptation to climate change is a relatively new task for conservation science (Heller and Zavaleta 2009). It is vital to build a toolbox of approaches that will be robust to multiple futures and uncertainty. We suggest that incorporating climate diversity in conservation planning activities may increase the potential resilience of conservation plans to the effects of climate change. Climate diversity could be incorporated into targeting and goal-setting schemes along with other elements of biotic and geophysical diversity (Beier andBrost 2010, Schloss et al. 2011). Results from our study region suggest that targeting coarse vegetation types, if stratified spatially and climatically, can function to capture climate diversity, but fall short of meeting climate targets in full, and thus there is a utility to directly targeting climate in planning. We use this additional information to recommend new priorities to the existing 2011 CLN, though it is reassuring that 82% of climate priority spots are already found in the CLN. Surveys of the 18% of climate priority spots not already included in the CLN are needed to better understand their conservation value. To develop this method further, additional work is needed to determine how finely classified vegetation and climate targets need to be to function as an effective surrogate for the other, and for which areas this method is likely to fail. The utility of biotic vs. abiotic approaches, or a synthesis between them, will need to be evaluated in different regions.

ACKNOWLEDGMENTS
We thank Kirk Klausmeyer, Sam Veloz, Claudia Tebaldi, Healy Hamilton, and Adina Merenlender and other members of the Bay Area Terrestrial Biodiversity and Climate Change Collaborative (TBC3.org) for v www.esajournals.org 17 April 2015 v Volume 6(4) v Article 65 many helpful discussions in devising this analysis. We also thank two anonymous reviewers for helpful edits and feedback. We thank the Gordon and Betty Moore Foundation for financial support for this work (grant #2861), and the California Landscape Conservation Cooperative for funding awarded to Jason Kreitler. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.