Staying connected: assessing the capacity of landscapes to retain biodiversity in a changing climate

Management for positive biodiversity outcomes under a changing climate requires a shift of perspective relative to traditional conservation. Here we develop a repeatable indicator for measuring the capacity of landscapes to retain biodiversity under a range of plausible climate futures, as a function of the condition and spatial configuration of native habitat. The Spatial Resilience Index extends an existing approach to assessing the potential for biodiversity associated with any location in a region to access suitable habitat in the surrounding landscape under climate change, incorporating multiple dispersal rates integrated over time and an optimised spatial structure. Derivation of the indicator is demonstrated for an Australian case study, covering the entire State of New South Wales, drawing on existing spatial datasets and models. Mapping of the Spatial Resilience Index across New South Wales suggests that different regions, and locations within these regions, vary markedly in their expected capacity to retain biodiversity, depending on the direct rate of climate change, the degree of climatic buffering (or reduction of climate velocity) afforded by landscape heterogeneity, and the degree of anthropogenic impacts on the connectedness of habitat in the landscape. The developed approach accounts for the interplay between these processes by treating them within a unified framework. The index highlights areas which can potentially benefit from adaptive management (e.g. habitat restoration) to enhance capacity to retain biodiversity under climate change, and offers an objective means of monitoring any resulting change in this capacity over time.


3
Vol:. (1234567890) and ecological communities has motivated considerable reappraisal of the fundamental goals and strategies of nature conservation over recent years (Prober et al. 2019). A common theme emerging from this thinking is the concept of "managing for change" (Stein et al. 2013;van Kerkhoff et al. 2019). From a biodiversity conservation perspective, managing for change entails a shift from the traditional focus on securing and maintaining species or communities where they occur at present, towards a strategy which embraces the possibility that the distributions of at least some species will need to shift under climate change, thereby altering the composition of presentday communities. The fundamental goal of biodiversity conservation implicit in this way of thinking is about ensuring persistence of collective biological diversity, i.e. "gamma diversity", within some relatively large spatial extent (Mokany et al. 2013)-e.g. a whole landscape, a country, a continent, or ultimately the entire planet. Change in the distribution of species, and in the composition of communities at local scales, is therefore viewed as being essential to enabling retention of gamma diversity, e.g. total number of species persisting, within the whole region of interest. If couched in terms of ecological resilience theory (Holling 1973;Quinlan et al. 2016), this localscale biological reorganisation is what underpins the capacity of the broader system to retain its essential structure (in this case retention of gamma diversity) in the face of a major, yet highly uncertain, stressor (in this case climate change).
Efforts to enhance the capacity of a landscape or region to retain biodiversity under climate change focus most commonly on actions aimed at either: (1) preventing future loss in the condition of relatively intact native vegetation through habitat protection; or (2) improving the condition of relatively degraded areas through various forms of habitat restoration. Most practitioners involved in such efforts are, however, acutely aware that biodiversity persistence under climate change is likely to be a function not just of the total amount of habitat protected and/or restored in a landscape, but also of the spatial location, configuration and resultant landscape-scale function of that habitat. If the spatial distribution of fundamental climatic niches of species are shifting then the persistence of these species will depend, at least in part, on the availability of habitat within areas of suitable climate under future, not just present, climatic conditions (Hannah et al. 2007). For at least some species the ability to disperse to climatically suitable areas of habitat into the future will also depend on the degree to which present and future habitats are connected spatially (Keeley et al. 2018;Schloss et al. 2022). The role that habitat connectedness across space and time plays in enabling persistence of biodiversity under climate change therefore constitutes an important facet of the 'spatial resilience' of natural systems (Cumming 2011). In their review of the literature on conservation strategies under climate change, McLaughlin et al. (2022) identify increasing connectivity as one of the top two recommendations.
Effective consideration of these factors continues to pose significant challenges for quantitative sciencebased analyses aiming to inform large-scaled habitat protection and restoration efforts (Gillson et al. 2013;Jones et al. 2016;Donaldson et al. 2017). To help appreciate the nature of these challenges it is worth first distinguishing between two broad types of question typically addressed by such analyses. The first type of question relates to the expected capacity of a given system, as a whole (e.g. an entire landscape or region), to retain biodiversity in the face of ongoing climate change given the present condition and spatial configuration of habitat across that system. Assessing capacity at this whole-system level can help in characterising the overall need for action to protect and/or restore habitat. Repeated application of such an assessment also offers a means of monitoring the extent to which a system's expected capacity to retain biodiversity has improved or declined over time as a function of observed changes in habitat condition or configuration, including as a result of any implemented management. The second type of question relates to the prioritisation of individual management actions. Assuming that the decision has already been made to invest in a program of habitat protection or restoration within a given landscape or region, then this is the challenge of deciding precisely where best to direct this effort spatially within the system concerned.
Approaches to addressing these two types of question are often developed and applied in isolation from each other. For example, while techniques for mapping potential climate refugia (e.g. Baumgartner et al. 2018), or significant pathways of habitat connection (e.g. Nunez et al. 2013;Schloss et al. 2022), can play a valuable role in informing spatial prioritisation of habitat protection or restoration, they rarely provide an adequate foundation for assessing and monitoring the capacity of an entire system to retain biodiversity under climate change. Effectively coupling the prioritisation of individual actions and system-level monitoring requires an integrated approach to addressing these activities within an adaptive policy or management framework. A relatively straightforward strategy for achieving such integration is to view, and perform, spatial prioritisation of actions within a broader landscape as a logical extension of the analytical approach used to assess the state of that landscape at a wholesystem level (Ferrier and Drielsma 2010). This involves predicting the marginal gain in the overall state of the system as a result of the local implementation of a given action (e.g. restoration).
Here we describe a new indicator-the Spatial Resilience Index-for assessing system-level capacity to retain biodiversity in the face of climate change as a function of the condition and spatial configuration of habitat across the entire Australian State of New South Wales (NSW), covering an area of 801,150 km 2 . This work builds on a solid foundation of biodiversity-related spatial data, and associated models, already established by previous studies across NSW including existing indicators of habitat condition Love et al. (2020) and ecosystem persistence Drielsma et al. (2020) delivered under the NSW biodiversity indicator program. Also of direct relevance and value to the current work, are the data and models generated by Drielsma et al.'s (2017) analysis of potential impacts of plausible climate futures on biodiversity, and prioritisation of revegetation actions to ameliorate these impacts.
Our approach to developing the new indicator is based on, and extends, that already employed in deriving an indicator assessing capacity to retain biodiversity in the face of climate change at global scale-the Bioclimatic Ecosystem Resilience Index, BERI ). This particular approach offers an alternative to, and effectively occupies the middle ground between, the two major analytical paradigms which have dominated most other efforts to address this general challenge over recent years. The first of these paradigms focuses on modelling shifts in the distribution of individual species as a spatially explicit function of projected changes in climate, employing an extensive array of correlative and mechanistic modelling techniques (Urban et al. 2016). While advances continue to be made in extending these techniques to account for complicating factors such as population dynamics, evolutionary adaptation, and species interactions, the species-specific data and knowledge needed to undertake such modelling is rarely available for more than a small proportion of all species.
This has motivated many workers to explore the second of these two major paradigms, which focuses exclusively on analysing spatiotemporal patterns and dynamics of abiotic environmental variables alone, without any direct use of biological data or explicit modelling of biological responses (Carroll et al. 2017). Prominent examples of this paradigm include analyses of climate velocity -i.e. the speed at which hypothetical organisms associated with a given location would need to move across the landscape to track projected shifts in climate (Brito-Morales et al. 2018)-and adaptation strategies based on "conserving nature's stage" (Beier et al. 2015)-i.e. ensuring that retained areas of natural habitat encompass high levels of diversity in, and steep gradients of, climate and other relevant abiotic variables (e.g. soils).
In common with this second paradigm, the approach we employ here makes no attempt to explicitly model responses of individual species to climate change. However, unlike other applications of the paradigm, our approach works with environmental variables (including climate variables) which have been non-linearly scaled to better reflect the observed effect that changes in these variables have on the species composition of ecological communities. This is achieved by applying generalised dissimilarity modelling, GDM (Ferrier et al. 2007) to best-available occurrence records for large numbers of species to transform multidimensional environmental space such that distances within this transformed space correlate as closely as possible with observed levels of turnover in species composition between locations under present-day climatic conditions. Space-for-time substitution can then be used to predict the level of compositional turnover expected between a location under present conditions and either itself, or another location, under any future climate scenario (Fitzpatrick et al. 2011 ;Ferrier et al. 2012;Blois et al. 2013).
The existing global-scale BERI indicator (Ferrier et al. 2020) combines GDM-based modelling of compositional turnover with cost-benefit assessment of habitat connectivity (Drielsma et al. 2007a) to estimate the connectedness of a given focal gridcell to areas of supporting natural habitat in the surrounding landscape which, under a plausible range of climate futures, are predicted to support an assemblage of species similar to that of the focal cell under present climatic conditions. This indicator provides an effective means of assessing capacity to retain biodiversity from the perspective of each-and-every individual cell in a region of interest, thereby allowing variation in this capacity to be mapped across that region. Here we significantly extend this existing analytical capability to derive a truly system-level indicator of the capacity of an entire region (in this case the state of NSW) to retain biodiversity in the face of climate change. This is achieved by making further use of fitted GDM-based models to integrate consideration of spatial variation in species composition, i.e. 'beta diversity', into the estimation of the region's collective capacity to retain gamma diversity.
We demonstrate the utility of this extended approach for whole-system monitoring by generating results for the new indicator across NSW, employing high resolution habitat-condition mapping and GDM modelling of vascular-plant compositional turnover, while additionally accounting for plausible variation in future climate trajectories and in maximum velocities of biological redistribution expected for different plant-dispersal modes. Although beyond the scope of this initial study, it is also intended that the indicator will, in the future, provide a foundation for prioritising actions aimed at enhancing the capacity of this system to retain biodiversity under climate change.

Overview
The Spatial Resilience Index, operationalised here for reporting across NSW, builds on the cost-benefit approach to assessing habitat connectivity originally described by Drielsma et al. (2007b) and applied in a climate change context at global scale in the BERI indicator . The application within NSW further builds on existing biodiversity and environmental data layers and models already in use across the State, at 9 arc-second (250 m) resolution. Climate scenarios were drawn from the NSW/ACT Regional Climate Modelling project (NARCliM) (Evans et al. 2014) domain, which provides extensive buffering beyond the State borders as required in any study of ecological responses to climate change. Modelling of current and future patterns of turnover in species composition across space and time was based on GDM models previously fitted within the NARCLiM domain (Drielsma et al., 2015) and used to assess ecosystem persistence  under the NSW biodiversity indicator program (OEH and CSIRO 2019).
A habitat condition surface for this same domain was drawn from the Biodiversity Impacts and Adaptation Project (BIAP) undertaken by Drielsma et al. (2017), and updated within NSW to include more recent modelling of ecological condition ) also published as part of the NSW biodiversity indicator program. Here we treat habitat condition as a continuous measure which spans all land uses and ecosystem states, ranging from 0 (degraded) to 1 (intact). Each grid cell in the analysis domain therefore has a habitat condition score, and the resultant surface represents the current state of the landscape, which is tested under multiple climate scenarios. The condition surface defines the permeability of the landscape (for the calculation of least-cost paths) and is combined with pairwise compositional similarity to determine the benefit of connectedness to the test cell. In this latter role, condition is treated as the effective habitat area of the grid cell, such that a condition of 0.5 is equivalent to half an intact grid cell.
The complete workflow for derivation of the Spatial Resilience Index is illustrated in Fig. 1. Firstly, a GDM model is fitted to current climate, static landform and substrate descriptors and aggregated biological survey data. The model predicts the level of dissimilarity (or conversely similarity) in species composition expected between any pair of grid-cells if those cells were still covered by intact native habitat. This model of present-day patterns of compositional similarity between cells is then used to project, through space-for-time substitution, changes in these patterns expected in the future as a function of different climate scenarios. This is achieved by using the GDM-based nonlinear transformations of climate variables, fitted to the present-day data, to transform both the present-day and future climate grids.
A grid of current habitat condition is then used to derive the Spatial Resilience Index for each individual cell. This involves estimating the connectivity   Table S2, with an enlargement of the geometry close to the centre (b) showing the effects of pixelation on radial segments. The larger diagram (a) shows the complete radial structure, with a least-cost path (dark red) connecting the focal cell to a highlighted (red) target segment. More permeable paths span more intact condition. The value of each segment is determined by the ecological similarity of the environment at cells within this segment, relative to the current environment of the focal cell. This similarity is depicted here for current climatic conditions (green) and a single future climate scenario (turquoise) of that cell to all cells in the surrounding landscape which are predicted (based on the GDM modelling) to be able to support a similar composition of species, under a given climate scenario, to that currently supported by the cell of interest. This analysis is repeated for a representative set of climate scenarios, and a plausible range of median dispersal distances for the biological group concerned (i.e. vascular plants in the case of this study). The resulting estimates of the total amount of compositionally-similar habitat which is well connected to the cell of interest are then averaged across these analyses. The cellwise Spatial Resilience Index is derived by expressing this average as a proportion of the maximum possible value obtained if the entire landscape were covered by intact habitat, and climatic conditions are fixed at their current level.
Scaling up the cellwise calculations described above, to thereby assess the capacity of whole regions to retain biodiversity under climate change, requires consideration of spatial variation in species composition (i.e. beta diversity) both within and between these regions. To help enable this assessment, the present-day GDM model is further used to estimate the compositional uniqueness of each cell within the NARCLiM domain.
We now describe in greater depth each of these steps in deriving the Spatial Resilience Index.

Generalised dissimilarity modelling
A Generalised Dissimilarity Model of compositional turnover in vascular plants (Drielsma et al. 2015) was used to map the distribution of cells predicted to have a similar ecological environment in the future to that currently associated with a given focal cell. The GDM scales the similarity in environment between two locations in terms of the predicted Sørenson-index similarity in their species composition. The model was fitted within the NARCLiM climate domain (full extent shown in Fig. S2) which includes all of NSW, the Australian Capital Territory, and Victoria, and is further buffered into Queensland and South Australia in order to adequately represent climate change processes in NSW.
Two categories of environmental variables were used in the modelling, those that change over time as a function of projected climate change, and landform and substrate descriptors which are in this case assumed to be constant over the time period.
NARCLiM provides a full suite of 20-year climate scenarios for a baseline 'present' period 1990-2009, and two future periods, of which we selected one, 2060-2079 for this analysis. The dataset covers four Global Climate Models (GCMs) spanning a range of possible futures: CSIRO mk 3.0 (warm, dry); MIROC 3.2 (warm, wet); MPI ECHAM 5 (hot, dry); and CCCMA3 (hot, wet). For each of these GCMs, three Regional Climate Models (RCMs) are used to downscale the data to a 250 m grid for the entire domain. A single socio-economic scenario, SRES A2 was modelled in NARCliM. This scenario remains consistent with the world's current trajectory and has equivalent forcing to Representative Concentration Pathway 8.5. Summary climate metrics, comprising standard BIO-CLIM variables, were derived using ANUCLIM 6.1 (Hutchinson and Xu 2015) and supplemented by radiation-adjusted statistics derived from monthly summaries (Reside et al. 2013;Harwood et al. 2016).
A suite of environmental layers were provided as candidate predictor variables for the GDM model described below, and Table S1 lists the 25 variables included in the final model and their sources.
Subsequent spatio-temporal analysis and projection of the fitted GDM model is achieved by transformation of the environmental layers for each climate scenario. The predicted pairwise similarity (s) between any two grid cells (i,j) separated by space and/or time can then be calculated as a function of the sum across all layers of the absolute pairwise difference between grid cells (|x i -x j |) within each transformed layer (x) following (Ferrier et al. 2007) as This variable is scaled from 0 (compositionally distinct) to 1 (compositionally identical). For any grid cell (i), it is also possible to calculate a continuous measure of the relative area of similar ecological environments in an intact pre-industrial landscape as the sum of pairwise similarities ( ∑ j=N all j=1 s ij ) to all (N all ) other cells (j). The inverse of this sum provides a measure of the original compositional uniqueness of each cell (Supplementary Figure S2a). Those ecological environments that are relatively unique, such as the high alpine areas will have a lower ∑ s ij than more common environments. This surface was generated for the baseline climate and was then saved for �xi−xj� subsequent use in summarising the Spatial Resilience Index for regional reporting units.

Habitat condition
Habitat condition is here used to describe the cellwise intactness of native vegetation relative to its preindustrial state. This is represented as an index from 0 (totally transformed) to 1 (totally intact) for each grid-cell and is the key test variable being assessed by the indicator. At the resolution of an individual cell, within which we assume homogeneity, habitat condition can be viewed as a proxy for local system resilience (sensu Cumming (2011)). The impact of changes in landscape-scale habitat condition, through ongoing degradation or restoration are reflected in the Spatial Resilience Index. However, it is worth noting that the assessment of condition beyond the site scale is challenging, and that issues such as inaccuracies or resolution issues in input data and the requirement to average condition across all states within grid cells will always be limiting factors. Here we used a condition surface previously generated for reporting in NSW. This surface synthesises land cover, land use, fractional cover dynamics and other landscape features ) (Drielsma et al. 2015) and represents a habitat-condition baseline (h) centred on 2013. Key assumptions underpinning the extrapolation of the surface beyond NSW (following Drielsma et al. 2015) are that locations which have a similar structural type and annual variability will have similar condition under similar management, that particular land-use or management practices have the same consequences in the same vegetation type, and that these relationships (calibrated for NSW) are constant across the study area. Whilst these areas beyond NSW have reduced resolution, they are used only as a buffer region to avoid edge effects in the analysis.

Calculation of the spatial resilience index per grid cell
The Spatial Resilience Index was calculated for the full set of NARCLiM futures, taking into account a range of effective dispersal distances. The cellwise calculation of the index for each dispersal distance was based on the Ferrier et al. (2020) extensions of the Drielsma et al. (2007b) cost-benefit approach, thereby enabling assessment of habitat connectivity in a climate-change context. Our application of this approach employs a new formalised geometric structure to maximise the rigour with which habitat connectivity is assessed within the neighbourhood around each focal cell (Fig. 2).

Radial geometric structure
In common with other spatial analyses, the geometric configuration of spatial units employed directly affects the results. The Spatial Resilience Index is calculated using a radial grid of connected "segments"-defined here as the intersection between sectors and radial bands-increasing in size with distance from the focal cell. Habitat condition is then averaged across all 250 m cells falling within each segment. Drielsma et al. (2007b) do not define a particular geometry in their original description of the cost-benefit approach, but show a figure with irregular 'petals'. More recent iterations such as Drielsma et al. (2017) have defined a radial structure with approximate doubling of petal area in concentric rings. The BERI index for tropical forests , working at 1 km grid-resolution globally, uses a grid split into 8 radial sectors (N, NE, E, SE, S, SW, W, NW), with concentric circles at radii of 2, 5, 10, 25, 50, 100, 200, 300, 400, and 500 km selecting by the cell centres. Whilst these regular approaches seem straightforward, and are well behaved at band widths of greater than 5 or 6 cells, their behaviour breaks down close to the focal cell (where connectivity is critical) as the size of the radial segments approaches that of the underlying raster's resolution and the Cartesian and polar-coordinate systems interact. Here there is a tendency for adjacent segments to change drastically in shape and orientation, making it difficult to treat them as equivalent units, and with some configurations resulting in non-concentric geometries. Since connectivity in the immediate vicinity of the test cell is critical, a more reliable spatial definition is required.
In order to formalise the geometric structure relative to the raster cell size, we established some simple guidelines. Firstly, the 8 radial sectors centred on the cardinal and ordinal directions from BERI were used, distributing segments equally in all directions. Secondly, the principle (Drielsma et al. 2017) of a regular stepwise (near doubling) increase in area in concentric bands was set as a target. Lastly, visual inspection of the geometric structure close to the test cell was carried out to examine the extent to which each radial band of segments conformed to the shape of the radial segments at greater radii. An exploratory environment was set up, and different approaches were iteratively tested against these criteria. Whilst there are only a limited number of possible geometries close to the focal cell, a generalised formulation that is equally applicable at all scales was defined relative to the Cartesian grid. Bands were defined (Supplementary Table S2) by including all cells with cell centres within a specified radius and specific radii were generated up to a maximum radius of 200 grid cells.

Dispersal rates
The potential rate of movement of plant populations under climate change will be determined in the first instance by the rate of dispersal of seed. Only once a seed germinates and establishes at a new location does the dispersal of pollen and the subsequent persistence of a population become relevant. Consequently, whilst pollen can move long distances the reproductive viability of individuals is determined largely by seed dispersal. Differences in seed dispersal strategies result in a range of median dispersal distances, and this means that the connectivity of the landscape will vary markedly for different species. In order to address this, we calculated the Spatial Resilience Index for multiple dispersal distances representing a range of plant dispersal rates, and averaged the final index scores across this range. Corlett and Westcott (2013) describe annual dispersal rates of up to 1500 m year −1 although most species disperse more slowly. Median dispersal distances (λ) were calculated for a 50-year period and four were used: 3.125, 6.25, 12.5 and 25 km, representing annual dispersal rates of 62.5, 125, 250 and 500 m year −1 . These were supplied to a negative exponential dispersal kernelp = e − d ∕ , where d is the habitat-condition weighted path length between locations. Annual exponentially-distributed dispersal rates were scaled up to take into account the implications of overlapping distributions over time. The exponential distribution defines a two-dimensional probability field (p t ) around each grid cell i for a single year's dispersal to other cells j. The probability of dispersal from this original source i in the following year can then be calculated from all cells j with a probability p t j using the same exponential distribution. A new expanded two-dimensional probability field (p t+1 ) is then calculated as the product of each cellwise probability from time t (p t j ) and the exponential distribution. This was explored through simulation for 50 years and, in line with the Central Limit Theorem, the resultant probability density functions rapidly converged to a Gaussian distribution as P = e −d 2 ∕ Λ , where Λ increases linearly with the number of iterations (g) asΛ = 5.79g . The probability (P) of reaching a location at distance d from the source, based on the resultant 50-year dispersal kernel (g = 50) can therefore be expressed in terms of the original exponential median distance λ as To derive habitat-weighted path lengths (d) across graph edges between segments, the Euclidian distances between the centres of adjacent segments were multiplicatively weighted by the mean habitat condition (h) of their respective segments inversely scaled between 2.5 (h = 0) and 1 (h = 1). This results in the traversal of a completely degraded segment (h = 0) costing 2.5 times that of a completely intact petal (h = 1) of equal size, in line with generic cost factors adopted by Love et al. (2020) and Drielsma et al. (2022). The habitat-weighted length (D) of the least cost path from the source (i) to segment (ψ) through multiple segments (q) was then calculated as the sum of the habitat-weighted lengths across all n traversed segments as

Cellwise calculation
The Spatial Resilience Index was calculated for each grid cell (i), using the formalised petal framework, which was laid over the 9″ grid surrounding that cell. The condition (h ψ ) of each segment (ψ) was calculated as a weighted sum of the habitat condition (h j ) values of all cells falling within this segment, with the contribution of each cell weighted by its pairwise similarity (s ij ) in predicted species composition to the P = e d 2 (5.79g ) ∑ n j=1 s ij h j ). The least-cost path connecting each segment (ψ) to the focal cell (i) was then determined, and the connectedness (χ), given the condition surface being tested, relative to that for an intact condition surface calculated as χ i = D i ∑ n j=1 s ij h j by estimating the shortest habitat-weighted path distance as described in Drielsma et al. (2007b). Overall connectedness (c) of cell i to all segments under a given scenario k was then calculated as c ik = ∑ n =1 χ i For each dispersal distance, the index was derived for all 13 (m) climate scenarios (current climate and 12 future scenarios) and the Limited Degree of Confidence statistic (McInerney et al. 2012;Bryan et al. 2014;Ferrier et al. 2020) applied to provide a robust summary metric of connectedness across all scenarios. The subsequent four results were then averaged to provide the cellwise estimate ( i ) of projected connectedness conferred by the surrounding landscape configuration: where c i0 is idealised connectedness of an entirely intact landscape under current climate. Consequently the index will be 1 for all cells under current climate within an intact landscape.

Calculation of spatial resilience index by summary regions
The global BERI indicator  aggregates the cellwise calculations (ρ i ) to country and regional scales by taking the arithmetic mean of all cells within a region (N region ). However, this approach does not account for the effects of beta diversity on the collective persistence of biodiversity at a regional level (i.e. gamma diversity), as outlined in Ferrier et al. (2004). The Spatial Resilience Index addresses this problem by explicitly accounting for varying levels of compositional similarity, through the sharing of species, between grid cells both within and between reporting regions. This is achieved by weighting the averaging of cellwise index values scaling aggregation according to the compositional uniqueness layer, i.e. the inverse of ( ∑ j=N all j=1 s ij ), generated under base- line climatic conditions (Supplementary Figure S2a). This builds on the approach previously described by Ferrier et al. (2004) and Allnutt et al. (2008) for weighting the contribution that individual cells make to regional indices of gamma diversity: Initial reporting for NSW was conducted using the Interim Biogeographic Regionalisation for Australia (IBRA) bioregions (Department of Sustainability Environment Water Populations and Communities 2012). The analysis was conducted within the NSW State boundary. This reflects both higher confidence in the mapping of habitat condition within NSW and the requirements of state-level reporting. Where parts of a given IBRA bioregion lay outside NSW, these cells were excluded from the analysis, and N region was therefore defined by the intersection of NSW and a given region. However, the contribution of areas outside NSW was taken into account through ∑ j=N all j=1 s ij . A further output to aid interpretation was generated in the form of summary histograms for each region, recording the frequency of the cellwise index in bins of width 0.02. The distribution of values within a region provides useful information in addition to the single summary statistic.

Results
Key interactions between components shaping cellwise results for the Spatial Resilience Index are illustrated for a selected portion of NSW in Fig. 3. The results in panels b-d were generated for all dispersal distances following the procedure outlined in the Cellwise Calculation section above but varying the inputs. This figure shows how the spatial configuration of habitat condition (panel a) determines the current connectedness of habitat in the absence of climate change (panel b: condition as panel a, current climate only), and how this combines with the expected effect of climate velocity under plausible climate futures (panel c, future climates but all cells have intact condition) to shape results for the full index (panel d, condition as panel a, all climate Fig. 3 Interaction between components of the Spatial Resilience Index for the area inland from Coffs Harbour in northern NSW: a current habitat condition; b the Spatial Resilience Index derived using current condition but assuming, hypothetically, no future change in climate-this highlights the effect of habitat connectedness on spatial resilience; c the Spatial Resilience Index derived using future climate scenarios but assuming, hypothetically, that NSW is covered by intact native habitat-this highlights the effect of climate velocity on spatial resilience; and d the full Spatial Resilience Index combining the effects of both habitat connectedness and climate velocity. All four maps share the same colour ramp. Four geographical areas are highlighted: 1. The Pilliga Forest, 2. Mount Kaputar National Park, 3. Nambucca River Catchment 4. Walgett to Lightning Ridge scenarios). Four areas are highlighted to explore these interactions. Two largely intact areas, the Pilliga Forest (area 2) and Mount Kaputar National Park (area 3) show how the moderating effects of terrain can enhance the capacity of landscapes to retain biodiversity under climate change. Despite the Pilliga Forest being a larger area in good condition, and therefore providing good habitat connectedness under current climate ( Fig. 3a and b), its flatter terrain means that climate velocity will be higher for this landscape (Fig. 3c), resulting in lower overall spatial resilience (Fig. 3d) than for the complex dissected terrain of Mount Kaputar which should allow many species to move a relatively short distance (e.g. upslope) to track shifts in suitable habitat under climate change.
A similar situation is demonstrated by the other two landscapes highlighted in Fig. 3, where land-use change has been greater-i.e. Walgett to Lightning Ridge (area 4) and Nambucca River Valley (area 3). The fragmentation of intact habitats in the Walgett to Lightning Ridge area ( Fig. 3a and b), combined with the high climate velocity expected in this flat terrain (Fig. 3c), results in very low spatial resilience (Fig. 3d). While habitats in the Nambucca River Valley are also quite fragmented, the effect this has on spatial resilience is offset, to some extent, by the topographic complexity and therefore lower climate velocity, of this landscape.
State-wide results for the Spatial Resilience Index are presented in Fig. 4, both at raw grid-cell resolution and summarised by Interim Biogeographic Regionalisation for Australia (IBRA) bioregions (Department of Sustainability Environment Water Populations and Communities 2012). The patterns highlighted in Fig. 3 hold true across the whole state (Fig. 4a), with highest spatial resilience in the topographically buffered and more intact areas near the coast, and reduced resilience in the more heavily utilised areas, and/or in areas of flatter terrain, and therefore higher climate velocity, further inland (i.e. the western two-thirds of the state) even where the habitat of these areas remains relatively intact. The mapping of results summarised by IBRA bioregions (Fig. 4b) following the aggregation calculation for region broadly reflects the patterns exhibited by the cellwise results. However, it should be noted that most of the bioregions near the Fig. 4 Spatial Resilience Index at raw grid resolution (a), summarised by IBRA region (b) and as frequency histograms (c). Two example bioregions are highlighted with similar sum-mary scores for the index, but exhibiting very different distributions of scores at the resolution of individual grid-cells coast have a lower aggregated score than might be expected. This is because the consideration given to beta-diversity patterns, when summarising the index by regions, is accounting for compositional differences between the more intact steeper-terrain portions of these regions, exhibiting high spatial resilience, and the less intact portions on flatter terrain, exhibiting much lower resilience.
Individual components of the index, as illustrated in Fig. 3, were calculated for all NSW to aid interpretation of the final index. A numerical breakdown of these components by IBRA bioregion using the aggregation calculation for region is shown in Fig. 5, where the committed impact of projected climate change is separated from the added impact of past habitat loss on connectedness, which can potentially be addressed through habitat restoration. Fig. 5 Breakdown of the Spatial Resilience Index for each reporting unit (IBRA bioregion) in NSW. The final index is a product of the combined effects of projected climate change (grey) and habitat loss (orange). The sum of the blue and orange bars indicates the maximum level of spatial resilience achievable in the absence of climate-change mitigation. The currently realised proportion of this potential is shown in the pie charts to the right

Discussion
Advances made in assessing spatial resilience Here we have successfully advanced an approach to assessing, at high spatial resolution, the capacity of landscapes to retain biodiversity in the face of climate change through three major refinements and extensions of the BERI indicator previously applied at global scale. Firstly, formalisation of the geometry of the radial grid employed in this approach has enhanced the rigour and consistency with which the connectedness of habitat close to each focal cell is assessed. Secondly, the incorporation of multiple dispersal-distance settings has allowed differences in dispersal strategies to be accounted for more effectively in assessing impacts of habitat connectedness on biodiversity persistence. Thirdly, incorporation of beta diversity into the aggregation of cellwise results to report the Spatial Resilience Index at whole-region level has ensured that such reporting better reflects expected changes in the collective persistence of gamma diversity for any given region, as opposed to a simple average of local persistence levels within that region.

Caveats and limitations of the approach
The technique we have used to calculate path lengths across the landscape necessarily assumes that agents will pass through all segments along the least cost path to the destination. This may be a reasonable assumption for many animals, but it is worth noting that some dispersal mechanisms involving, for example, wind-borne propagules or flying vectors may be able to ignore areas of locally degraded habitat, and dispersal in such cases may be better represented by a simple dispersal kernel (Shaw et al. 2006). A further implicit assumption is that dispersal distances remain constant across the whole analysis region, and there is no effect of local landscape heterogeneity on the requirement or ability of species to disperse.
In its current form, the Spatial Resilience Index does not take the environmental suitability of intervening habitat into account. The permeability of habitat is estimated purely as a function of its intactness or condition, not as a function of the type of habitat relative to that of the focal cell-e.g. an intact area of rainforest is assumed to offer the same permeability as an intact area of open woodland, regardless of whether the focal cell is itself rainforest or woodland. Whilst it would be possible to address this in the current framework for a static climate, doing so for a changing climate would require adoption of a more mechanistic approach such as the metacommunity modelling of Mokany et al. (2012) to fully account for the complex interactions between dispersal rates, available habitat and environmental suitability.
Another critical aspect of dispersal is time. The index is here calculated for a 2070 future, allowing 50 years of dispersal, but could readily be applied to other future periods. Importantly, however, appropriate kernels for multiple years of dispersal are not simple accumulations of annual dispersal kernels. The spatio-temporal analysis, and associated equation for multi-year dispersal, developed here in relation to the negative exponential kernel have potential for broader application in landscape ecology.
The results presented here for NSW are based on a GDM model for vascular plants. While compositional turnover in plants is well suited to environmental modelling due to a good coupling to local growing conditions and can be a good surrogate for other biological groups given the role plants play in defining habitat for other species, there are limitations to this assumption. Many mobile species will utilise specific landscape features when moving across the landscape, and may therefore be more affected by changes in the structural, rather than compositional, attributes of vegetation.
The threat of invasive species (Finch et al. 2021) under climate change is not directly addressed. Whilst migrating native species are theoretically considered in the modelling of compositional turnover, non-native species better adapted to the emerging climate, and taking advantage of weakened ecological communities, are likely to pose a complex management problem.
A final caveat is that the Spatial Resilience Index does not consider the possibility that the fundamental climatic niche of many species might be broader than that suggested by the currently realised distribution of the species concerned. Nor does it consider potential for species to undergo evolutionary adaptation in response to climatic conditions shifting beyond previous limits of tolerance (Bush et al. 2019).
Alternative approaches to assessing climate change connectivity Our approach addresses a particular dimension of connectivity under climate change, and in particular does not address the significance of the pathways themselves. Drielsma et al. (2022) show how a similar approach can be combined with an analysis of spatial links to provide more comprehensive support for adaptive management under current climate by highlighting the parts of the landscape which are important to maintain connectivity. A recent application of the Omniscape instance of the Circuitscape tool (Schloss et al. 2022) illustrates how such an approach can be extended to address connectivity under climate change. In common with climate velocity approaches (Brito-Morales et al. 2018) this work uses unscaled climate to define analogues, but addresses an aspect of the climatic quality of intermediate habitat using topoclimatic heterogeneity. The resultant maps show the parts of the landscape which under current configuration are important for the adaptive movement of species, and therefore important to preserve. In common with other prioritisation approaches, such patterns and priorities can be expected to change with the addition or removal of patches of habitat.

Implications of initial results for policy and management
The results presented in Fig. 5 indicate that climate change over the next 50 years is likely to have a strong effect on the availability and accessibility of habitat in climatic environments similar to those currently occupied by species throughout NSW. These results suggest that, even if NSW were still covered entirely by intact native habitat, the State's capacity to retain biodiversity under climate change would be only 52% of that expected if climatic conditions remained unchanged. In other words, future climate change alone is responsible for a 48% reduction in spatial resilience, as shown by the grey bar for NSW in Fig. 5. This loss can be reduced only through global climate-change mitigation and is therefore dependent on factors largely beyond the control of land-management actions within NSW. Given that the maximum level of spatial resilience achievable through land-management actions is limited to 52% it makes sense to assess the need and potential for such actions relative to this constraint. This is the approach adopted in the pie charts depicted in Fig. 5, where the spatial resilience associated with the current condition and configuration of native habitat across NSW, and within each of the State's bioregions, is expressed as a percentage of the maximum possible under the climate scenarios analysed here. These results indicate considerable potential, across all bioregions, to increase capacity to retain biodiversity under climate change through on-ground actions designed to restore and enhance the area, condition and connectedness of native habitat.

Extending application of the indicator to prioritising actions
The results presented in Fig. 5 help to highlight regions where the impact of climate change is likely to be felt more strongly than others, and where most potential exists to ameliorate this impact through land-management actions. However, it is important not to treat the cellwise results for the Spatial Resilience Index presented in this paper as a guide to prioritising where best to locate such actions within any given region. In its raw form, the Spatial Resilience Index for a grid-cell is purely a measure of the capacity of the surrounding landscape to support biodiversity currently associated with that cell, in the face of climate change. Considerable potential exists, however, to use the analytical capability established here as a foundation for building analyses aimed at prioritising on-ground management actions. For example, to assess the relative benefit of restoring habitat within a given parcel of land, the Spatial Resilience Index would be derived twice-first using the current condition layer for NSW (as for the results presented here), and then using a modified version of this layer in which all grid cells within the proposed restoration area are set to best possible condition. The difference between the two resulting values of the index then provides a measure of the expected benefit of the proposed action. If this marginal benefit is calculated for each individual grid-cell, in turn, across NSW then this analysis could be used to generate fine-resolution gridded surfaces indicating the relative priority for applying a given action across the State, as per the general approach to priority mapping described by Ferrier and Drielsma (2010).

Conclusions
Development of the Spatial Resilience Index provides a repeatable indicator for measuring the capacity of landscapes to retain biodiversity under climate change, as a function of the condition and spatial configuration of native habitat. Initial application of this indicator across NSW using the biodiversity indicator program's mapping of current habitat condition has revealed considerable variation in spatial resilience both between and within the State's bioregions. This variation results largely from the interplay between variation in levels of habitat loss and fragmentation, and variation in the extent to which climate velocity is likely to be moderated by topographic complexity. However, considerable potential exists across all regions to improve current levels of spatial resilience through on-ground actions aimed at restoring or enhancing habitat condition and connectedness. The analytical capability established by this study now offers a solid foundation for future mapping of priority areas for action, at fine spatial resolution across the State.