# Importance of spatio-temporal connectivity to maintain species experiencing range shifts

Jun-Long Huang, Marco Andrello, Alexandre Camargo Martensen, Santiago Saura, Dian-Feng Liu, Jian-Hua He, Marie-Josee Fortin

► To cite this version:

Jun-Long Huang, Marco Andrello, Alexandre Camargo Martensen, Santiago Saura, Dian-Feng Liu, et al.. Importance of spatio-temporal connectivity to maintain species experiencing range shifts. Ecography, 2020, 10.1111/ecog.04716 . hal-03410888

HAL Id: hal-03410888

https://hal.umontpellier.fr/hal-03410888

Submitted on 16 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ECOGRAPHY

Research

Importance of spatio-temporal connectivity to maintain species
experiencing range shifts

Jun-Long Huang, Marco Andrello, Alexandre Camargo Martensen, Santiago Saura, Dian-Feng Liu, Jian-Hua He and Marie-Josée Fortin

J.-L. Huang (https://orcid.org/0000-0002-8227-2215), D.-F. Liu (https://orcid.org/0000-0002-0929-8494) and J.-H. He (https://orcid.org/0000-0001-6706-1703) ✉ (hjianh@whu.edu.cn), School of Resources and Environment Sciences, Wuhan Univ., CN-430079 Wuhan, PR China. – JLH and M.-J. Fortin (https://orcid.org/0000-0002-9935-1366), Dept of Ecology and Evolutionary Biology, Univ. of Toronto, Toronto, ON, Canada. – M. Andrello (https://orcid.org/0000-0001-7590-2736), MARBEC, Univ. de Montpellier, CNRS, IFREMER, IRD, Sète, France. – A. C. Martensen (https://orcid.org/0000-0001-9017-7480), Centro de Ciências da Natureza, Univ. Federal de São Carlos, Brazil. – S. Saura (https://orcid.org/0000-0002-5239-0763), Depto de Sistemas y Recursos Naturales, ETSI Montes, Forestal y del Medio Natural, Univ. Politécnica de Madrid, Madrid, Spain.

Ecography

43: 591–603, 2020

doi: 10.1111/ecog.04716

Subject Editor: Timothy Keitt

Editor-in-Chief: Miguel Araújo

Accepted 22 November 2019

<figure><img src="image_2.png" /></figure>

Climate change can affect the habitat resources available to species by changing habi-
tat quantity, suitability and spatial configuration, which largely determine population
persistence in the landscape. In this context, dispersal is a central process for species
to track their niche. Assessments of the amount of reachable habitat (ARH) using
static snap-shots do not account, however, for the temporal overlap of habitat patches
that may enhance stepping-stone effects. Here, we quantified the impacts of climate
change on the ARH using a spatio-temporal connectivity model. We first explored the
importance of spatio-temporal connectivity relative to purely spatial connectivity in a
changing climate by generating virtual species distributions and analyzed the relative
effects of changes in habitat quantity, suitability and configuration. Then, we studied
the importance of spatio-temporal connectivity in three vertebrate species with diver-
gent responses to climate change in North America (grey wolf, Canadian lynx and
white-tailed deer). We found that the spatio-temporal connectivity could enhance the
stepping-stone effect for species predicted to experience range contractions, and the
relative importance of the spatio-temporal connectivity increased with the reduction
in habitat quantity and suitability. Conversely, for species that are likely to expand
their ranges, spatio-temporal connectivity had no additional contribution to improve
the ARH. We also found that changes in habitat amount (quantity and suitability)
were more influential than changes in habitat configuration in determining the relative
importance of spatio-temporal connectivity. We conclude that spatio-temporal con-
nectivity may provide less biased and more realistic estimates of habitat connectivity
than purely spatial connectivity.

Keywords: climate change, dynamic network model, landscape connectivity, Ontario, species distribution

<figure><img src="image_3.png" /></figure>

www.ecography.org

© 2020 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos
This is an open access article under the terms of the Creative Commons
Attribution License, which permits use, distribution and reproduction in any
medium, provided the original work is properly cited.
## Introduction

Climate change may cause contractions in species distributions if future conditions become unfavorable in some parts of their range, and dispersal to new habitat may be necessary for species' long-term persistence (Parmesan and Yohe 2003, Garcia et al. 2014). Furthermore, climate and land-use changes may alter the spatial composition and configuration of habitat, as well as habitat suitability, leading to declines in the availability of habitat resources (Opdam and Wascher 2004, Travis et al. 2013). Hence, species persistence will largely depend on future habitat spatial pattern (Årevall et al. 2018) and species ability to track their niche in space and time through dispersal (Schloss et al. 2012). Understanding how much habitat is reachable for species in a changing environment is, therefore, of paramount importance for conservation (Littlefield et al. 2019).

Assessing the amount of reachable habitat (ARH) requires quantifying the degree of landscape connectivity (Saura et al. 2014). To date, a large number of connectivity models and metrics have been proposed, such as those based on meta-population theory (Hanski and Ovaskainen 2000, Moilanen and Nieminen 2002), network theory (Urban and Keitt 2001, Bodin and Norberg 2006, Saura and Pascual-Hortal 2007), and on circuit theory (McRae et al. 2008). These approaches have also been employed to evaluate the impacts of climate change on habitat connectivity and species persistence (Dilts et al. 2016, Albert et al. 2017, Rehnus et al. 2018, Kanagaraj et al. 2019). Most connectivity models rely on spatial snap-shots of landscape structure to generate static connectivity estimates; however, some species predicted to experience ARH decrease or increase do not actually have immediate population declines or growth as expected (Metzger et al. 2009, Semper-Pascual et al. 2018, Lira et al. 2019). Such 'extinction debt' and 'colonization credit' phenomena (Tilman et al. 1994, Lira et al. 2019) require ecologists to consider temporal ecological processes going beyond a purely spatial perspective.

Although some recent studies highlighted the importance of accounting for temporal aspects in connectivity modeling (Alagador et al. 2014, Mui et al. 2017), only a few of them accounted for dispersal processes over generations (Saura et al. 2014, de la Fuente et al. 2018), the variation in importance of habitat patches in maintaining connectivity due to seasonal habitat change (Mui et al. 2017, Bishop-Taylor et al. 2018), or the implications of ephemeral patches and transient connectivity change in metapopulation dynamics (Reigada et al. 2015, Perry and Lee 2019). Yet, the importance of temporal dynamics of habitat patches has largely been ignored (Blonder et al. 2012). Indeed, some habitat patches may appear and/or disappear through time in a dynamic landscape. Because of these dynamics, habitat patches may coexist at certain time periods even though they do not exist together at the start or end of a given temporal period (Zeigler and Fagan 2014). Such temporal overlap of habitat patches may contribute to enhancing the stepping-stone effect. To account

for this effect, Martensen et al. (2017) developed a spatio-temporal network-based model that incorporates temporal links between habitat patches, and showed that purely spatial connectivity models underestimate the actual connectivity patterns for forest landscapes experiencing habitat loss and fragmentation. However, the importance of spatio-temporal connectivity still remains unexplored in the context of species tracking their climatic niche, especially for species displaying varying vulnerabilities to climate change.

Here, we highlight the need of spatio-temporal connectivity models (Martensen et al. 2017) to predict changes in the ARH across space and time under climate change. First, we provide a general framework to guide how to conduct spatio-temporal connectivity analyses for species experiencing range shifts. Then, we apply this framework to both theoretical simulations and case studies. The former generates a wide range of hypothetical species distributions on virtual landscapes, to illustrate under what conditions spatio-temporal connectivity can be more important than spatial-only connectivity. The latter applies ensemble species distribution models to predict potential distributions of three terrestrial mammals in North America with divergent vulnerabilities to climate change: white-tailed deer, Canada lynx and grey wolf. We expect that white-tailed deer and grey wolf will expand their range due to warmer temperatures and lower winter severity in future climates (Weiskopf et al. 2019), leading to a lower importance of spatio-temporal connectivity. We also expect that Canada lynx will contract its range due to less snowfall in winters (Hoving et al. 2005), which may produce a higher relative importance of spatio-temporal connectivity.

## Methods

### Spatio-temporal connectivity model

To assess the change in the amount of reachable habitat under varying climate conditions, we used a spatio-temporal connectivity model (Martensen et al. 2017). The spatio-temporal connectivity model is actually a specific example of the multi-layer network (Pilosof et al. 2017), in which habitat patches are considered as nodes, and two static landscape snapshots as two layers. Nodes are connected by spatio-temporal link that incorporates both intra-layer links (i.e. spatial links) and inter-layer links (i.e. temporal links) to represent their spatio-temporal interactions. Further, the spatio-temporal path is a path made of one or multiple spatio-temporal links, representing the possibility of a dispersing organism moving from a given habitat patch at time $t_1$ to another habitat patch at time $t_2$ ($t_1 < t_2$). To ensure the persistence from $t_1$ to $t_2$ in the dynamic landscape, species can either stay in a Stable patch, or disperse from a patch with habitat at $t_1$ (i.e. Loss or Stable) to another patch with habitat at $t_2$ (i.e. Gain or Stable). The schematic of spatio-temporal connectivity model is shown in Fig. 1.
<figure><img src="image_4.png" /></figure>

Figure 1. Schematic representation of the conceptual difference between spatial-only and spatio-temporal connectivity models. The spatial-only connectivity model provides assessments based on static landscape snapshots (panels (a) and (b)). At time $t_1$, the landscape has two connected habitat patches, *i* and *k*, and the interaction intensity (indicated by dispersal probability) between them is the same in both directions. At time $t_2$, patch *k* is lost while another patch *j* appears in the landscape; yet, because the distance between patches *i* and *j* is greater than species dispersal ability, there is no connection between them. Nevertheless, from the spatio-temporal perspective (panel (c)), where blue, red and green colors denote Stable, Loss and Gain types of patches, respectively; *i* and *j* may be connected if a dispersing individual first moves to *k*, which may still have habitat at an intermediate time point, then disperses again from *k* to *j* before *k* is eventually lost. Patch *k* therefore acts as a stepping-stone from the spatio-temporal perspective. The solid and dashed lines denote the essential and auxiliary links (explained in the text). In this case, the interaction between patches is directional and the intensity may not be identical in both directions.

Given the time interval from $t_1$ to $t_2$, each habitat patch in a network can be classified into one of three types through polygon change analysis: Loss (i.e. habitat in $t_1$ but not in $t_2$), Stable (i.e. habitat in both $t_1$ and $t_2$) and Gain (i.e. habitat in $t_2$ but not in $t_1$). Between any pair of habitat patches, the temporal link occurs only when two patches simultaneously exist during some intermediate time window from $t_x$ to $t_y$ ($t_1 < t_x < t_y < t_2$). If such information is unknown, the temporal link may or may not occur, and is therefore assigned a probability *w* between 0 and 1.

Links in the spatio-temporal network can also be classified into two types: 1) the essential link, which alone, is able to transport an individual from a patch with habitat at $t_1$ to another patch with habitat at $t_2$; and 2) the auxiliary link, which by itself, is unable to successfully transport an individual as the essential link does. However, auxiliary links are complementary to essential links. For example, the path ‘Stable→Loss→Stable’ is not possible if using essential links only, as no habitat will be readily available at $t_2$ when species arrives at Loss. Nevertheless, this path becomes possible after accounting for the auxiliary link: the species can first move to Loss that still has habitat in an intermediate time point through an auxiliary link, then it will need to disperse again from Loss to Stable using an essential link before Loss eventually loses habitat to ensure its persistence. Moreover, for the path ‘Stable→Loss→Gain→Stable’, the three links are auxiliary, essential and auxiliary links in sequence. The

essential link from Loss to Gain may or may not be possible, as whether Loss and Gain coexist in the intermediate time window is unknown; the third link is an auxiliary link, since there is no habitat readily available at $t_1$. To summarize, two nodes in the network can either be directly connected using a single essential link, or indirectly connected through a path consisting of several essential and auxiliary links, as long as there is at least one essential link. In the latter case, the stepping-stone effect is taken into account, which helps species reach more habitat than by using essential links only. The temporal dispersal probability ($p_{ij}^q$) of each link scenario is listed in Table 1 (see Martensen et al. 2017 for additional details).

Lastly, the spatio-temporal dispersal probability between habitat patches ($p_{ij}^{st}$) is the product of spatial dispersal probability ($p_{ij}^{so}$) and temporal dispersal probability ($p_{ij}^t$). The spatial dispersal probability is commonly calculated using an exponential form (Urban and Keitt 2001): $p_{ij}^{so} = e^{-kd_{ij}}$; where *k* is a species-specific coefficient indicating the species dispersal ability, and $d_{ij}$ is the spatial distance between node *i* and *j*. The distance can either be Euclidean distance or any effective distance that reflects landscape resistance. To parameterize *k* when effective distance is used, $p_{ij}^{so}$ can be fixed as 0.5, the species median dispersal distance can be multiplied by the median cost value of the cost surface, and their product can substitute $d_{ij}$ in the equation to solve for *k* (Gurrutxaga et al. 2011).
<table><caption>Table 1. Temporal dispersal probability ($p_{ij}^t$) for each type of link, adapted from Martensen et al. (2017).</caption><thead><tr><th rowspan="3">Source patch: individual location at t<sub>1</sub> for the essential links or at t<sub>x</sub> (t<sub>1</sub> &lt; t<sub>x</sub> &lt; t<sub>2</sub>) for the auxiliary links</th><th colspan="6">Target patch: individual location after t<sub>1</sub></th></tr><tr><th colspan="3">Essential link (individual location at t<sub>2</sub>)</th><th colspan="3">Auxiliary links (individual location in t<sub>γ</sub>, t<sub>x</sub> &lt; t<sub>γ</sub> &lt; t<sub>2</sub>)</th></tr><tr><th>Stable</th><th>Loss</th><th>Gain</th><th>Stable</th><th>Loss</th><th>Gain</th></tr></thead><tbody><tr><td>Stable</td><td>1</td><td>0</td><td>1</td><td>N/A</td><td>1</td><td>N/A</td></tr><tr><td>Loss</td><td>1</td><td>0</td><td>w</td><td>N/A</td><td>1</td><td>N/A</td></tr><tr><td>Gain</td><td>0</td><td>0</td><td>0</td><td>1</td><td>w</td><td>1</td></tr></tbody></table>

## Connectivity metrics

The probability of connectivity (PC) is the probability that two individuals randomly placed within a landscape fall into habitat patches that are reachable for each other across the habitat network (Saura and Pascual-Hortal 2007). To quantify the amount of reachable habitat, we used the equivalent connectivity (EC; Saura et al. 2011), which is calculated as the square root of the numerator of PC. Further, to distinguish the area-based and habitat configuration aspects of connectivity metrics, we partitioned the PC into three fractions following Saura et al. (2014): PC<sub>intra</sub>, PC<sub>direct</sub> and PC<sub>step</sub> (see Supplementary material Appendix 1 for their formulae). All metrics except PC<sub>intra</sub> are functions of dispersal probabilities. Spatial-only (with an 'so' superscript) and spatio-temporal (with an 'st' superscript) connectivity metrics were calculated using $p_{ij}^{so}$ and $p_{ij}^{st}$, respectively (Martensen et al. 2017). EC<sup>so</sup> indicates the habitat resource (e.g. habitat area and suitability) of a single patch that would provide the same probability of connectivity as the actual habitat pattern in the landscape. EC<sup>st</sup> indicates the resources of a single 'Stable' habitat patches that would provide the same probability of spatio-temporal connectivity as the spatio-temporal network. PC<sub>intra</sub><sup>so</sup> indicates the intra-node connectivity provided by all existing patches in the landscape, while PC<sub>intra</sub><sup>st</sup> indicates the intra-node connectivity provided by Stable patches only. PC<sub>direct</sub><sup>so</sup> is the inter-node connectivity provided by direct spatial links, while PC<sub>direct</sub><sup>st</sup> the inter-node connectivity provided by direct spatio-temporal links. PC<sub>step</sub><sup>so</sup> denotes the inter-node connectivity contributed by stepping-stone effects through indirect spatial links, while PC<sub>step</sub><sup>st</sup> denotes the inter-node

connectivity contributed by stepping-stone effects through indirect spatio-temporal links.

## General framework for the spatio-temporal network construction

Here, we demonstrate how to construct spatio-temporal network for species experiencing range shifts under climate change. First, the habitat suitability for the focal species needs to be obtained using some approach on the considered landscape (Fig. 2A). To analyze connectivity for species distributed in a continuous landscape, the landscape can be divided into multiple equal-size blocks with a coarser resolution than that of the landscape. For example, the spatial resolution of the landscape is one grid cell while the resolution of the block is five grid cells (Fig. 2B). Then, block centroids can be considered as nodes for the network analysis (Dilts et al. 2016). The mean and the sum of habitat suitability are computed for each individual block. A block is considered suitable habitat when its mean habitat suitability is >0.5. Therefore, a node is classified as Loss if the block's mean habitat suitability is >0.5 at $t_1$ but not at $t_2$, Gain if opposite, and Stable if higher than 0.5 at both $t_1$ and $t_2$. To account for the habitat suitability heterogeneity within the block, centroids can be located in a spatially weighted mean fashion (Fig. 2C). Node weights are given by the sum of habitat suitability within the block at $t_1$ for the Loss, the sum at $t_2$ for the Gain, and the mean of sums at $t_1$ and at $t_2$ for the Stable, and zero for blocks without habitat at both $t_1$ or $t_2$.

Given that landscape heterogeneity influences species dispersal, landscape resistance can be calculated as the inverse of

<figure><img src="image_5.png" /><figcaption>Figure 2. (A) An example of the generated virtual species distributions (50×50 grid cells) with probability-of-occurrence represented on each grid cell. (B) Dividing the landscape into equal-size study blocks of 5×5 grid cells. (C) The location of block centroid was weighted by the probability-of-occurrence, and the centroid was considered as the node in the subsequent network analysis. See Supplementary material Appendix 1 Fig. A2 for the example of the actual species.</figcaption></figure>
habitat suitability. This approach is suggested to be sufficient to estimate landscape resistance in the absence of empirical species movement data (Zeller et al. 2018). Next, least-cost modeling (Adriaensen et al. 2003) can be applied, using block centroids as endpoints and resistance surface with the fine resolution, to estimate effective distance among immediately adjacent blocks (van Etten 2017). For the spatio-temporal connectivity case, the effective distance can be calculated as a weighted average form to reflect spatio-temporal variations in landscape resistance:

$$d_{ij,\text{st}} = \alpha \times d_{ij,t1} + (1 - \alpha) \times d_{ij,t2}, \alpha \in (0,1), \quad (1)$$

where $d_{ij,t1}$ and $d_{ij,t2}$ denote the effective distance from $i$ to $j$ at $t_1$ and at $t_2$, respectively; and $\alpha$ is the weight given to $d_{ij,t1}$ that ranges from zero to one, since the actual time of change in landscape resistance is unknown. To keep the model simple in the general framework, $\alpha$ and $w$ (Table 1) are set to 0.5.

After obtaining node attributes and types, as well as links between nodes, spatial-only and spatio-temporal connectivity analyses can be conducted for the focal species. The relative importance of spatio-temporal connectivity is indicated by the difference proportion of spatio-temporal EC over its spatial-only counterpart (i.e. $(EC^{st} - EC^{so}_{t2}) / EC^{so}_{t2}$). Thus, positive values indicate that spatio-temporal connectivity enhances connectivity.

## Theoretical simulations using virtual species

We used the ‘virtualspecies’ package (Leroy et al. 2016) in R ver. 3.6.0 (<www.r-project.org>) to generate virtual species distributions on a two-dimensional landscape of 50×50 grid cells, based on two hypothetic environmental variables with the same spatial extent. As environmental conditions often display spatial autocorrelation (Legendre 1993), we used the ‘gstat’ package (Pebesma and Heuvelink 2016) to generate the spatially correlated random field for each variable (under current and future scenarios) using sequential Gaussian simulations with an exponential variogram model (sill number = 0.025, range parameter = 15). A total number of 120 virtual species distributions (half current, half future) were generated. Changes between current and future distributions therefore represented a wide range of species vulnerabilities to climate change. On the distribution map, each grid cell has a value ranging from 0 to 1 indicating the probability-of-occurrence (i.e. relative habitat suitability) for the virtual species. Following the framework in the last section and defining the block size as 5×5 grid cells, we conducted both spatial-only and spatio-temporal connectivity analyses for three hundred virtual species, as combinations of the sixty divergent vulnerabilities to climate change and five different dispersal abilities (5, 10, 15, 25 and 35 grid cells as the median natal dispersal distance). Connectivity analyses were performed using R scripts in conjunction with the command version of Conefor (<www.conefor.org/>).

## Statistical analysis

For the simulated virtual species, we disentangled the effect sizes of changes in 1) habitat quantity, 2) habitat suitability and 3) habitat spatial configuration on the relative importance of spatio-temporal connectivity for species with different dispersal abilities. The first explanatory variable was the difference proportion of the number of blocks with suitable habitat ($t_2$ over $t_1$). The second explanatory variable was the difference proportion of the habitat suitability averaged over the 100 blocks. The third explanatory variable was the difference proportion of the habitat aggregation index (He et al. 2000), which is independent of habitat composition and indicates the degree of habitat aggregation/isolation. As some explanatory variables did not display linear relationship with the response variable even after log-transformations (Supplementary material Appendix 1 Fig. A3), a random forests algorithm was used to investigate the variable importance given its robustness to nonlinear relationships (Breiman 2001). The number of trees was 1000, and 999 permutations for the ‘out-of-bag’ data were conducted to assess the variable importance. We used the ‘randomforest’ (Liaw and Wiener 2002) and ‘landscapemetrics’ packages (Hesselbarth et al. 2019) in R to perform these analyses.

## Case studies using real species

### Study area and species

We applied spatio-temporal connectivity to actual species in North America. The study area covers ten eco-regions in central and southern Ontario, Canada (ca 603 thousand km²; Fig. 3), with dominant climatic types of warm-summer humid continental and subarctic. The northern eco-regions were excluded as the primary land cover there is tundra, which is unlikely to become suitable habitat in the short term according to the climate change velocity map (Burrows et al. 2011). According to Species at Risk in Ontario List, there are in total 232 endangered and/or threatened species, whilst climate change is considered one of the main threats to biodiversity in the province (Ontario Biodiversity Council 2011).

Following Meurant et al. (2018), we selected three terrestrial species with different dispersal abilities, habitat needs and vulnerabilities to potential climate changes in temperate North America: white-tailed deer *Odocoileus virginianus*, Canada lynx *Lynx canadensis* and grey wolf *Canis lupus*. The latter two are carnivorous species that primarily prey on snow-shoe hare and deer, respectively. We obtained the four species presences data (year 1970–2018) from Global Biodiversity Information Facility (GBIF) and their characteristics (e.g. dispersal abilities and longevity) from the literature. See Table 2 for detailed information.

We downloaded current and projected bioclimatic variables, as well as current forest loss map from databases listed in Table 2. We defined the year 2030 as a near future time point, so that the temporal difference from current to future approximately fits the longevity of our focal species. For future climate scenarios, we adopted the four Representative
<figure><img src="image_6.png" /><figcaption>Figure 3. Species occurrence data downloaded from GBIF (time captured: 1970–2018).</figcaption></figure>

Concentration Pathway scenarios (RCP 2.6, 4.5, 6.0 and 8.5; see IPCC AR5), in which the extent of emissions and global warming increases from RCP 2.6 to RCP 8.5.

We then divided the study area into 428 equal-size 40×40 km blocks. This size was chosen to ensure that all focal species could move across blocks during one generation to facilitate network analyses. The snowshoe hare was excluded from the connectivity analysis, as its dispersal ability is smaller than the block size.

## Ensemble forecasting for species distributions

We used the 'biomod2' package (Thuiller et al. 2009) to implement ensemble species distribution projections, which are considered to increase the predictive accuracy and robustness with regard to individual models in the face of model uncertainties (Araújo and New 2007). To account for predator–prey interactions, we first predicted potential spatial distributions of snowshoe hare and white-tailed deer using the approach below, and then integrated their distributions with other environmental variables to predict potential distributions of their predators: Canada lynx and grey wolf, respectively.

As a first step to the ensemble forecasting, we added an equal number of pseudo-absences to presences data of each

focal species across study area using a surface range envelope model, which forces pseudo-absence candidates to be selected in condition that differs from a defined proportion (0.025 in this study) of presence data (Beaumont and Hughes 2002).

Then, we ran a principal component analysis (PCA) on the nineteen bioclimatic variables using the 'ade4' package (Dray and Dufour 2007) to avoid potential collinearity problems. We retained the ones potentially most relevant to the species distribution among highly correlated bioclimatic variables, based on the biological knowledge of the focal species: as the white-tailed deer is sensitive to drought events (Tosa et al. 2018), Bio12 (annual precipitation) and Bio17 (precipitation of the driest quarter) were retained; as snowshoe hare and Canada lynx are sensitive to snowfall and snow cover (Hoving et al. 2005), Bio11 (mean temperature of the coldest quarter) and Bio19 (rain precipitation of the coldest quarter) were retained; and given that the grey wolf is sensitive to prey distributional changes, the same variables as the ones retained for the white-tailed deer were selected as they also displayed great importance and low correlation for grey wolf (Supplementary material Appendix 1 Fig. A1).

Next, we modeled future forest cover change. The future (2012–2030) forest loss amount across the study area was predicted based on the assumption that the future period
<table><caption>Table 2. Data used in this study. All raster data were transformed to a uniform projection (Canada Lambert Conformal Conic) and spatial resolution (1 km) prior to modeling. The mismatch between the time range of species presence data and climate variables was considered unimportant, given the large overlap between the two time periods (1970–2018 and 1970–2000; Mestre et al. 2015).</caption><thead><tr><th>Data category</th><th>Name</th><th>Data format</th><th>Spatial resolution</th><th>Time captured</th><th>Source</th><th>Description</th></tr></thead><tbody><tr><td>Species occurrence data</td><td>Snowshoe hare <i>Lepus americanus</i><br>White-tailed deer<br><i>Odocoileus virginianus</i><br>Canada lynx <i>Lynx canadensis</i><br>Grey wolf <i>Canis lupus</i></td><td>csv (exported to point files)</td><td>–</td><td>1970–2018</td><td>Global Biodiversity Information Facility (GBIF; available at: &lt;www.gbif.org/&gt;)</td><td>Order: <i>Lagomorpha</i>.<br>Habitat types: Forest, shrubland.<br>Longevity (years): 5.<br>Median dispersal distance (km) based on natal dispersal: 6.0 (Albert et al. 2017).<br>Number of presences data from GBIF: 143.<br>Order: <i>Artiodactyla</i>.<br>Habitat types: Forest, Savannah, grassland.<br>Longevity (years): 6–14.<br>Median dispersal distance (km) based on natal dispersal: 36.0 (Nixon 1994).<br>Number of presences data from GBIF: 803.<br>Order: <i>Carnivora</i>.<br>Habitat types: Forest, shrubland, grassland.<br>Longevity (years): 15–20.<br>Median dispersal distance (km) based on natal dispersal: 88.0 (Poole 1997).<br>Number of presences data from GBIF: 398.<br>Order: <i>Carnivora</i>.<br>Habitat types: Forest, shrubland, grassland, rocky areas.<br>Longevity (years): 6–14.<br>Median dispersal distance (km) based on natal dispersal: 129.7 (Sutherland et al. 2000).<br>Number of presences data from GBIF: 61.<br>The nineteen bioclimatic variables used as input data for species distribution models (SDMs).<br>The future climate projections downscaled to North America. The general circulation model (GCM) of the used dataset is MOHC_HADGEM2_ES (Collins et al. 2011).<br>Forest loss data were used to predict future forest cover. Forest cover data were used as the input data for SDMs.<br>Used to delineate the boundary of the study area. Five northern eco-regions were excluded, as there are few presences of our focal species, and the primary land-cover type is tundra.</td></tr><tr><td>Climate data</td><td>Current bioclimatic variables<br>Future bioclimatic variables</td><td>Raster datasets (tif)</td><td>30°</td><td>1970–2000</td><td>WorldClim (&lt;http://worldclim.org/&gt;; Fick and Hijmans (2017))</td><td></td></tr><tr><td>Land cover data</td><td>Forest canopy cover<br>Forest loss</td><td>Polygon</td><td>30 m</td><td>2000</td><td>Hansen et al. (2013)</td><td></td></tr><tr><td>Ecoregion data</td><td>Ecoregion</td><td>Polygon</td><td>–</td><td>2000–2012</td><td>–</td><td>Land Information Ontario (&lt;www.ontario.ca/data/ecoregions&gt;)</td></tr></tbody></table>
would have the same average annual forest loss rate as in the past (2000–2012). The spatial allocation of future forest loss was achieved by assuming that the block that lost more forest in the past would also have higher probability of losing forest in the future. To do so, we computed forest loss amount in the past period for each block and created a probability-of-loss surface. We then generated future forest loss points based on the total loss amount and probability-of-loss surface using the ‘Create Spatially Balanced Points’ tool in ArcGIS 10.2, and converted spatial points to the raster map. By overlapping raster layers of forest loss (actual 2000–2012 and predicted 2012–2030) with the forest cover in the year 2000, we obtained the potential forest cover in 2030.

Further, we built the ensemble models for each species by combining eight algorithms available in the ‘biomod2’ package: GAM, GBM, GLM, CT, ANN, FDA, MARS and RF. To evaluate the predictive accuracies of these models, we conducted three-fold cross-validations by randomly splitting the dataset into two parts: 80% for calibration and 20% for evaluation. Then, we used the true skill statistics (TSS) as the model evaluation metrics following Allouche et al. (2006). We selected models with TSS >0.7 for the final ensemble forecasting.

Finally, the ensemble models were projected to the current and future climate scenarios based on four ensemble algorithms: mean possibility, median possibility, weighted mean possibility and committee averaging, and only the one with the highest TSS was retained for further analysis.

We generated potential 1) habitat suitability (range: 0–1000) maps and 2) binary range maps for our focal species. We defined the habitat threshold as the minimum habitat suitability within the species range, rather than a uniform threshold in the general framework. The thresholds for white-tailed deer, Canada lynx and grey wolf were: 781, 286 and 349, respectively.

We then conducted connectivity analyses for the three focal species following approaches elaborated in the section 2.3. We also conducted a sensitivity analysis to test whether spatio-temporal connectivity assessments are robust to parameterizations of *w* (Table 1) and *α* (Eq. 1). Nine assessments (as combinations of *w*=0.25, 0.5, 0.75 and *α*=0.25, 0.5, 0.75) were generated and standard deviations of values of spatio-temporal connectivity metrics were calculated for each focal species.

# Results

## Virtual species

In the 60 climate change scenarios, 27 of them (45%) led to net increments in habitat quantity from $t_1$ to $t_2$, and 36 of them (60%) had net increases in habitat suitability (Supplementary material Appendix 1 Fig. A4). Out of the 300 virtual species, 44% of them displayed higher (positive) importance of spatio-temporal connectivity relative to spatial-only connectivity in the amount of reachable habitat

(EC; Fig. 4). Cases in which spatio-temporal connectivity had higher relative importance occurred mostly when habitat quantity, suitability and aggregation were predicted to decrease in $t_2$. Conversely, when habitat quantity, suitability and aggregation were predicted to increase, spatio-temporal connectivity had lower relative importance.

The relationships between the three components of connectivity and changes in landscape were similar to the relationship observed for EC. In particular, $PC_{direct}$ and $PC_{step}$ were more affected by changes in landscape than $PC_{intra}$, since only 8% of the virtual species showed higher importance in the $PC_{intra}^{st}$ fraction, while 32% of them showed higher importance in $PC_{direct}^{st}$ and 45% in $PC_{step}^{st}$.

The relationships between the relative importance of spatio-temporal connectivity and the three landscape variables did not display much variation among species dispersal abilities. However, species dispersal ability changed their relative effect sizes in explaining the importance of spatio-temporal connectivity. Specifically, habitat aggregation was important only for short-distance dispersers, while its importance decreased with increasing dispersal ability; yet, changes in habitat quantity and suitability always had larger effect sizes than habitat configuration across all dispersal abilities (Supplementary material Appendix 1 Fig. A6).

## Focal species in Ontario

The goodness-of-fit (TSS values) of the ensemble forecasting of the species distributions were 0.93, 0.94, 0.95 and 0.94 for snowshoe hare, white-tailed deer, Canada lynx and grey wolf, respectively; indicating that our models had considerably high predictive accuracy. In current conditions, white-tailed deer, Canada lynx and grey wolf are likely to occur across southern and central Ontario, respectively; and in all future scenarios, white-tailed deer and grey wolf were predicted to expand northwards, whereas Canada lynx's range was predicted to shrink (Fig. 5).

The number of ‘Gain’ blocks ranged from 68 (under RCP 6.0) to 168 (under RCP 8.5) for white-tailed deer, and from 120 (under RCP 6.0) to 210 (under RCP 8.5) for grey wolf (Supplementary material Appendix 1 Fig. A5A). In contrast, habitat quantity of Canada lynx declined considerably under all predicted future climate scenarios, and RCP 8.5 produced the most severe range contraction. In terms of changes in habitat suitability, RCP 8.5 produced the largest improvement for white-tailed deer and grey wolf, while causing the largest reduction for Canada lynx (Supplementary material Appendix 1 Fig. A5B).

$EC^{so}$ of white-tailed deer and grey wolf was predicted to increase in 2030, with RCP 8.5 resulting in the largest increases, while the opposite was true for Canada lynx (Fig. 6A). On the other hand, $EC^{st}$ of white-tailed deer and grey wolf were smaller than $EC^{so}$ at 2030 (Fig. 6B), and the relative reductions in $EC^{st}$ ranged from ca –40% to –20%. Overall, larger relative reductions in spatio-temporal connectivity metrics were observed in scenarios predicting larger range expansions. Conversely, Canada lynx had larger $EC^{st}$
<figure><img src="image_7.png" /><figcaption>Figure 4. Relative importance of spatio-temporal connectivity for virtual species under a range of changes in (A) habitat quantity; (B) habitat suitability; and (C) habitat configuration (indicated by aggregation index; AI). Negative values on the x-axis correspond to decreases in habitat quantity, habitat suitability and habitat aggregation/connectedness, respectively. Positive values on the y-axis correspond to cases where spatio-temporal connectivity is more important than spatial-only connectivity. The number above each panel indicates species dispersal ability, corresponding to the product of median natal dispersal distance (i.e. 5, 10, 15, 25 and 35 grid cells) and the median cost value (i.e. 2, as the reciprocal of the median habitat suitability of 0.5) on the resistance surface. Colors indicate the connectivity metric: EC, equivalent connectivity; $PC_{direct}$, probability of connectivity provided by direct connections; $PC_{intra}$, probability of connectivity provided by habitat amount; $PC_{step}$, probability of connectivity provided by indirect connections that use stepping-stones.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5. Predicted distributions for the three focal species under different climate scenarios, with relative habitat suitability (0–1000) represented on each grid cell.</figcaption></figure>
than EC⁶⁰ and the relative increase ranged from 1% to 106%, with RCP 8.5 producing the largest increment. The three PC⁵⁴ fractions of Canada lynx also exceeded their spatial-only counterparts in most cases, and showed up to a 250% relative increase for PC⁵⁴ and 550% for PC⁵⁴ under the RCP 8.5 scenario. Meanwhile, all connectivity metrics did not vary significantly under a range of *w* and α values (Fig. 6B), indicating the robustness of the spatio-temporal connectivity model to the choice of these parameters.

## Discussion

Previous research has shown that purely spatial connectivity models may substantially underestimate the actual amount of reachable habitat (ARH) and colonization-extinction rates in landscapes with high levels of land-use change (Martensen et al. 2017). Here, we showed that this is also true in the context of climate change. However, the relative importance of spatio-temporal connectivity depends on species responses to climate change, and on the magnitude of climate change. Moreover, we found evidence that changes in habitat quantity and suitability have greater effects on the relative importance of spatio-temporal connectivity, compared with changes in habitat configuration and species dispersal ability.

## Importance of the spatio-temporal connectivity given different vulnerabilities of species to climate change

Our results showed that spatio-temporal connectivity was important for those species predicted to experience range contractions under climate change. For example, the shift of winter precipitation to less snowfall but more rainfall (Supplementary material Appendix 1 Fig. A7) in Canada lynx's current range is likely to have negative implications for the species (Hoving et al. 2005). Consequentially, the ARH of Canada lynx was predicted to decrease by 21–81% in 2030 from the spatial-only perspective, in agreement with other studies documenting the declining trend of available habitat for Canada lynx in this region (Hornseth et al. 2014). Yet, the spatio-temporal connectivity estimation indicated that the decline in connectivity might not be so dramatic, due to the stepping-stone effects enhanced by spatio-temporal connectivity. The evidence is shown in Fig. 6B, where the PC⁵⁴ fraction is nearly five times higher than PC⁵⁴ under the most severe scenario. As such, spatio-temporal connectivity may partly explain the delayed extinctions given habitat reductions and fragmentation (Tilman et al. 1994, Martensen et al. 2017). Indeed, Hooftman et al. (2016) found that connectivity was one factor strongly correlated to the likelihood of delayed local extinction, and the level

<figure><img src="image_9.png" /><figcaption>Figure 6. (A) Amount of reachable habitat from the spatial-only and the spatio-temporal (*w* = 0.5, α = 0.5) perspectives (indicated by EC⁶⁰ and EC⁵⁴, respectively) of each focal species under different climate scenarios. Parameter *w* denotes the temporal dispersal probability that may or may not be possible, and α the weight given to landscape resistance at *t*₁. For the current condition, EC⁵⁴ is missing as no spatio-temporal connectivity was calculated from past to current. (B) The relative importance of spatio-temporal connectivity compared with spatial-only connectivity for each focal species under different climate scenarios, with the error bar indicating the standard deviation of nine relative differences under a range of *w* and α values.</figcaption></figure>
of connectivity could be enhanced by the existence of more
temporally dynamic populations acting stepping-stones,
which is in congruence with our findings.

Nevertheless, the relative importance of spatio-temporal connectivity decreased for species expected to increase their future habitat amount. For example, warmer temperatures and less snowfall in the winter (Supplementary material Appendix 1 Fig. A7) favor the survival and population spreading of white-tailed deer (Weiskopf et al. 2019), which may also lead to population expansion of its predator, grey wolf (Fig. 5). As a result, the spatio-temporal ARH of these two species was smaller than the spatial-only one (Fig. 6A). This finding agrees with Martensen et al. (2017), which showed that the positive contribution of spatio-temporal connectivity decreased with increasing habitat gained. From the model perspective, when there are additional habitat patches gained in the landscape, none of them can provide intra-node connectivity in the spatio-temporal case, as only ‘Stable’ patches are considered as places where the intra-patch connectivity occurs. Yet, these patches can provide intra-node connectivity in the spatial-only case as they indeed exist in the landscape at time $t_2$. From the ecological standpoint, such reductions in spatio-temporal connectivity relative to its spatial-only counterpart may also partly explain the ‘colonization credit’ (Lira et al. 2019), which is an inverse concept to ‘extinction debt’ and refers to the number of individuals or species yet to colonize a focal habitat due to positive landscape changes. Hence, there should be many instances where spatio-temporal connectivity may provide less biased and more realistic estimates of habitat connectivity, given the presence of extinction debts and colonization credits (Kuussaari et al. 2009, Lira et al. 2019) and the stepping-stone effect across space and time (Saura et al. 2014).

Likewise, the intensity of climate change also changed the relative importance of spatio-temporal connectivity, as it affected habitat suitability within the patch. Previous studies have suggested that habitat suitability is an important factor affecting metapopulation dynamics (Franken and Hik 2004) and should be incorporated into connectivity models (Visconti and Elkin 2009). Here, we quantified habitat suitability based on species ecological niches, habitat requirements and food resources (i.e. prey-predator interactions), and used it as node attributes. Therefore, habitat quantity and quality jointly determined the amount of habitat resources. Also, cases where habitat quantity increased and/or decreased often corresponded to the overall improvement and/or decline in habitat suitability, so these two factors affected the importance of the spatio-temporal connectivity in the same direction (Fig. 4).

Effects of habitat amount versus effects of habitat
isolation and species dispersal ability

In a climate change context, variations in the ARH mainly
come from changes in landscape structures, including habi-
tat quantity, suitability and spatial configuration (Villard and
Metzger 2014). Numerous studies have debated the relative
effect sizes of these factors on population persistence for dif-
ferent taxa (Heller and Zavaleta 2009, Hodgson et al. 2011,

Fahrig 2017, Fletcher et al. 2018). Their effects on the relative importance of spatio-temporal connectivity, however, has yet to be fully understood. This issue has significance for conservation science, as it helps to identify conditions under which extinction debts are likely to occur in the landscape (Semper-Pascual et al. 2018), and guides conservation practices to avoid eventual local extinctions. We found that the effect sizes of these three factors varied depending on species dispersal ability (Supplementary material Appendix 1 Fig. A6). However, the effect size of habitat configuration was always smaller than that of the other two factors, and the effect size continuously decreased with increasing dispersal ability. This response of the effect size of habitat configuration to dispersal ability is logical because, for long-distance dispersers, their dispersal ability is already high enough to reach even the most isolated habitat patches. The implication of this finding for conservation science is that, given limited conservation budgets, it may be more effective to focus on increasing habitat amount (quantity and suitability) than improving habitat connections (Hodgson et al. 2011). Conversely, for species with smaller dispersal ability, enhancing habitat connections (e.g. via movement corridors; Gilbert-Norton et al. 2010) could increase species movement and the amount of reachable habitat to cope with climate change.

In conclusion, our study highlights the need of spatio-
temporal connectivity models to evaluate the impact of cli-
mate change on species' amount of reachable habitat across
space and time. Specifically, our results suggest an increasing
importance of spatio-temporal connectivity for species with
contracting ranges by enhancing the stepping-stone effect.
Moreover, spatio-temporal connectivity can also provide less-
biased results for species with expanding ranges. Nevertheless,
it does not mean that spatio-temporal connectivity can guar-
antee the long-term persistence of species in the face of dras-
tic environmental change: it may only increase the 'relaxation
time' (i.e. the time lag to extinction; Kuussaari et al. 2009).
We still need to put a lot of effort in habitat restoration, pri-
marily increasing habitat quantity and suitability, to prevent
eventual species extinctions.

Data availability statement

Data are available from the Dryad Digital Repository: <http://dx.doi.org/10.5061/dryad.gb5mkkwk8> (Huang et al. 2019).

*Acknowledgements* – We thank Louis Donelle for the help during data analysis.

*Funding* – This study is funded by CSC scholarship to J-LH, National Natural Science Foundation of China (grant no.: 41771429 and 41871301) to J-HH and D-FL, and NSERC Discovery Grant to M-JF.

*Conflicts of interests* – None.

*Author contributions* – JLH designed the study, performed data processing and analysis, and led the writing. MA, ACM and MJF contributed to study design and writing. SS, DFL and JHH helped edit the manuscript. All authors conceived ideas to this study.
## References

Adriaensen, F. et al. 2003. The application of 'least-cost' modelling as a functional landscape model. - Landscape Urban Plan. 64: 233-247.

Alagador, D. et al. 2014. Shifting protected areas: scheduling spatial priorities under climate change. - J. Appl. Ecol. 51: 703-713.

Albert, C. H. et al. 2017. Applying network theory to prioritize multispecies habitat networks that are robust to climate and land-use change. - Conserv. Biol. 31: 1383-1396.

Allouche, O. et al. 2006. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). - J. Appl. Ecol. 43: 1223-1232.

Araújo, M. B. and New, M. 2007. Ensemble forecasting of species distributions. - Trends Ecol. Evol. 22: 42.

Årevall, J. et al. 2018. Conditions for successful range shifts under climate change: the role of species dispersal and landscape configuration. - Divers. Distrib. 24: 1598-1611.

Beaumont, L. J. and Hughes, L. 2002. Potential changes in the distributions of latitudinally restricted Australian butterfly species in response to climate change. - Global Change Biol. 8: 954-971.

Bishop-Taylor, R. et al. 2018. Evaluating static and dynamic landscape connectivity modelling using a 25-year remote sensing time series. - Landscape Ecol. 33: 625-640.

Blonder, B. et al. 2012. Temporal dynamics and network analysis. - Methods Ecol. Evol. 3: 958-972.

Bodin, Ö. and Norberg, J. 2006. A network approach for analyzing spatially structured populations in fragmented landscape. - Landscape Ecol. 22: 31-44.

Breiman, L. 2001. Random forests. - Mach. Learn. 45: 5-32.

Burrows, M. T. et al. 2011. The pace of shifting climate in marine and terrestrial ecosystems. - Science 334: 652.

Collins, W. J. et al. 2011. Development and evaluation of an earth-system model - HadGEM2. - Geosci. Model Dev. 4: 1051-1075.

de la Fuente, B. et al. 2018. Predicting the spread of an invasive tree pest: the pine wood nematode in southern Europe. - J. Appl. Ecol. 55: 2374-2385.

Dilts, T. E. et al. 2016. Multiscale connectivity and graph theory highlight critical areas for conservation under climate change. - Ecol. Appl. 26: 1223-1237.

Dray, S. and Dufour, A. B. 2007. The ade4 package: implementing the duality diagram for ecologists. - J. Stat. Softw. 22: 1-20.

Fahrig, L. 2017. Ecological responses to habitat fragmentation per se. - Annu. Rev. Ecol. Evol. Syst. 48: 1-23.

Fick, S. E. and Hijmans, R. J. 2017. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. - Int. J. Climatol. 37: 4302-4315.

Fletcher, R. J. et al. 2018. Is habitat fragmentation good for biodiversity? - Biol. Conserv. 226: 9-15.

Franken, R. J. and Hik, D. S. 2004. Influence of habitat quality, patch size and connectivity on colonization and extinction dynamics of collared pikas *Ochotona collaris*. - J. Anim. Ecol. 73: 889-896.

Garcia, R. A. et al. 2014. Multiple dimensions of climate change and their implications for biodiversity. - Science 344: 486-.

Gilbert-Norton, L. et al. 2010. A meta-analytic review of corridor effectiveness. - Conserv. Biol. 24: 660-668.

Gurrutxaga, M. et al. 2011. Key connectors in protected forest area networks and the impact of highways: a transnational case study

from the Cantabrian Range to the western Alps (SW Europe). - Landscape Urban Plan. 101: 310-320.

Hansen, M. C. et al. 2013. High-resolution global maps of 21st-century forest cover change. - Science 342: 850-3.

Hanski, I. and Ovaskainen, O. 2000. The metapopulation capacity of a fragmented landscape. - Nature 404: 755-758.

He, H. S. et al. 2000. An aggregation index (AI) to quantify spatial patterns of landscapes. - Landscape Ecol. 15: 591-601.

Heller, N. E. and Zavaleta, E. S. 2009. Biodiversity management in the face of climate change: a review of 22 years of recommendations. - Biol. Conserv. 142: 14-32.

Hesselbarth, M. H. K. et al. 2019. landscapemetrics: an open-source R tool to calculate landscape metrics. - Ecography 42: 1648-1657.

Hodgson, J. A. et al. 2011. Habitat area, quality and connectivity: striking the balance for efficient conservation. - J. Appl. Ecol. 48: 148-152.

Hooftman, D. A. P. et al. 2016. Reductions in connectivity and habitat quality drive local extinctions in a plant diversity hotspot. - Ecography 39: 583-592.

Hornseth, M. L. et al. 2014. Habitat loss, not fragmentation, drives occurrence patterns of Canada lynx at the southern range periphery. - PLoS One 9: e113511.

Huang, J.-L. et al. 2019. Data from: Importance of spatio-temporal connectivity to maintain species experiencing range shifts. - Dryad Digital Repository, <http://dx.doi.org/10.5061/dryad.gb5mkkwk8>.

Hoving, C. L. et al. 2005. Broad-scale predictors of Canada lynx occurrence in eastern North America. - J. Wildl. Manage. 69: 739-751.

Kanagaraj, R. et al. 2019. Predicting range shifts of Asian elephants under global change. - Divers. Distrib. 25: 822-838.

Kuussaari, M. et al. 2009. Extinction debt: a challenge for biodiversity conservation. - Trends Ecol. Evol. 24: 564-71.

Legendre, P. 1993. Spatial autocorrelation - trouble or new paradigm. - Ecology 74: 1659-1673.

Leroy, B. et al. 2016. virtualspecies, an R package to generate virtual species distributions. - Ecography 39: 599-607.

Liaw, A. and Wiener, M. 2002. Classification and regression by randomForest. - R News 2: 18-22.

Lira, P. K. et al. 2019. Temporal lags in ecological responses to landscape change: where are we now? - Curr. Landscape Ecol. Rep. 4: 70-82.

Littlefield, C. E. et al. 2019. Connectivity for species on the move: supporting climate-driven range shifts. - Front. Ecol. Environ. 17: 270-278.

Martensen, A. C. et al. 2017. Spatio-temporal connectivity: assessing the amount of reachable habitat in dynamic landscapes. - Methods Ecol. Evol. 8: 1253-1264.

McRae, B. H. et al. 2008. Using circuit theory to model connectivity in ecology, evolution and conservation. - Ecology 89: 2712-2724.

Mestre, F. et al. 2015. Combining distribution modelling and non-invasive genetics to improve range shift forecasting. - Ecol. Model. 297: 171-179.

Metzger, J. P. et al. 2009. Time-lag in biological responses to landscape changes in a highly dynamic Atlantic forest region. - Biol. Conserv. 142: 1166-1177.

Meurant, M. et al. 2018. Selecting surrogate species for connectivity conservation. - Biol. Conserv. 227: 326-334.

Moilanen, A. and Nieminen, M. 2002. Simple connectivity measures in spatial ecology. - Ecology 83: 1131-1145.
Mui, A. B. et al. 2017. Using multiple metrics to estimate seasonal landscape connectivity for Blanding's turtles (*Emydoidea blandingii*) in a fragmented landscape. - Landscape Ecol. 32: 531-546.

Nixon, C. M. 1994. Behavior, dispersal and survival of male white-tailed deer in Illinois. - Biological Notes no. 139.

Ontario Biodiversity Council. 2011. Ontario's biodiversity strategy, 2011: renewing our commitment to protecting what sustains us.

Opdam, P. and Wascher, D. 2004. Climate change meets habitat fragmentation: linking landscape and biogeographical scale levels in research and conservation. - Biol. Conserv. 117: 285-297.

Parmesan, C. and Yohe, G. 2003. A globally coherent fingerprint of climate change impacts across natural systems. - Nature 421: 37-42.

Pebesma, E. and Heuvelink, G. 2016. Spatio-temporal interpolation using gstat. - R J. 8: 204-218.

Perry, G. L. W. and Lee, F. 2019. How does temporal variation in habitat connectivity influence metapopulation dynamics? - Oikos 128: 1277-1286.

Pilosof, S. et al. 2017. The multilayer nature of ecological networks. - Nat. Ecol. Evol. 1: 101.

Poole, K. G. 1997. Dispersal patterns of lynx in the Northwest Territories. - J. Wildl. Manage. 61: 497-505.

Rehnus, M. et al. 2018. Alpine glacial relict species losing out to climate change: the case of the fragmented mountain hare population (*Lepus timidus*) in the Alps. - Global Change Biol. 24: 3236-3253.

Reigada, C. et al. 2015. Metapopulation dynamics on ephemeral patches. - Am. Nat. 185: 183-195.

Saura, S. and Pascual-Hortal, L. 2007. A new habitat availability index to integrate connectivity in landscape conservation planning: comparison with existing indices and application to a case study. - Landscape Urban Plan. 83: 91-103.

Saura, S. et al. 2011. Network analysis to assess landscape connectivity trends: application to European forests (1990-2000). - Ecol. Indic. 11: 407-416.

Saura, S. et al. 2014. Stepping stones are crucial for species' long-distance dispersal and range expansion through habitat networks. - J. Appl. Ecol. 51: 171-182.

Schloss, C. A. et al. 2012. Dispersal will limit ability of mammals to track climate change in the Western Hemisphere. - Proc. Natl Acad. Sci. USA 109: 8606-8611.

Semper-Pascual, A. et al. 2018. Mapping extinction debt highlights conservation opportunities for birds and mammals in the South American Chaco. - J. Appl. Ecol. 55: 1218-1229.

Sutherland, G. D. et al. 2000. Scaling of natal dispersal distances in terrestrial birds and mammals. - Conserv. Ecol. 4: 16.

Thuiller, W. et al. 2009. BIOMOD - a platform for ensemble forecasting of species distributions. - Ecography 32: 369-373.

Tilman, D. et al. 1994. Habitat destruction and the extinction debt. - Nature 371: 65-66.

Tosa, M. I. et al. 2018. Increased overwinter mortalities of white-tailed deer (*Odocoileus virginianus*) fawns during a drought year. - Can. J. Zool. 96: 55-61.

Travis, J. M. J. et al. 2013. Dispersal and species' responses to climate change. - Oikos 122: 1532-1540.

Urban, D. and Keitt, T. 2001. Landscape connectivity: a graph-theoretic perspective. - Ecology 82: 1205-1218.

van Etten, J. 2017. R Package gdistance: distances and routes on geographical grids. - J. Stat. Softw. 76: 1-21.

Villard, M.-A. and Metzger, J. P. 2014. REVIEW: Beyond the fragmentation debate: a conceptual model to predict when habitat configuration really matters. - J. Appl. Ecol. 51: 309-318.

Visconti, P. and Elkin, C. 2009. Using connectivity metrics in conservation planning - when does habitat quality matter? - Divers. Distrib. 15: 602-612.

Weiskopf, S. R. et al. 2019. Climate change effects on deer and moose in the midwest. - J. Wildl. Manage. 83: 769-781.

Zeigler, S. L. and Fagan, W. F. 2014. Transient windows for connectivity in a changing world. - Mov. Ecol. 2: 1.

Zeller, K. A. et al. 2018. Are all data types and connectivity models created equal? Validating common connectivity approaches with dispersal data. - Divers. Distrib. 24: 868-879.

Supplementary material (available online as Appendix ecog-04716 at <www.ecography.org/appendix/ecog-04716>). Appendix 1.