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

may provide less biased and more realistic estimates of habitat connectivity than purely spatial connectivity.


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 andYohe 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 andWascher 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 metapopulation theory (Hanski andOvaskainen 2000, Moilanen andNieminen 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 andLee 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 spatiotemporal 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.

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.
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 t ) 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) ; 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). Figure 1. Schematic representation of the conceptual difference between spatial-only and spatio-temporal connectivity models. The spatialonly 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.
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 intra , PC direct and PC step (see Supplementary material Appendix 1 for their formulae). All metrics except PC intra 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 so 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 st 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 intra so indicates the intra-node connectivity provided by all existing patches in the landscape, while PC intra st indicates the intranode connectivity provided by Stable patches only. PC direct so is the inter-node connectivity provided by direct spatial links, while PC direct st the inter-node connectivity provided by direct spatio-temporal links. PC step so denotes the inter-node connectivity contributed by stepping-stone effects through indirect spatial links, while PC step st 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 Table 1. Temporal dispersal probability ( p ij t ) for each type of link, adapted from Martensen et al. (2017).
Source patch: individual location at t 1 for the essential links or at t x (t 1 < t x < t 2 ) for the auxiliary links Target patch: individual location after t 1 Essential link (individual location at t 2 ) Auxiliary links (individual location in t y , t x < t y < t 2 ) Stable Loss Gain Stable Loss Gain 595 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: 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 α 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, α 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 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.

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 2 ; 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 snowshoe 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 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 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 (2017)) The nineteen bioclimatic variables used as input data for species distribution models (SDMs would have the same average annual forest loss rate as in the past (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(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-ofloss 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 whitetailed 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.

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, whitetailed 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  . 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. than EC so and the relative increase ranged from 1% to 106%, with RCP 8.5 producing the largest increment. The three PC st fractions of Canada lynx also exceeded their spatial-only counterparts in most cases, and showed up to a 250% relative increase for PC direct st and 550% for PC step st 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 step st fraction is nearly five times higher than PC step so 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 6. (A) Amount of reachable habitat from the spatial-only and the spatio-temporal (w = 0.5, α = 0.5) perspectives (indicated by EC so and EC st , 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 1 . For the current condition, EC st is missing as no spatiotemporal 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. 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 habitat quantity, suitability and spatial configuration (Villard and Metzger 2014). Numerous studies have debated the relative effect sizes of these factors on population persistence for different 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 spatiotemporal connectivity models to evaluate the impact of climate 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 lessbiased results for species with expanding ranges. Nevertheless, it does not mean that spatio-temporal connectivity can guarantee the long-term persistence of species in the face of drastic 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, primarily increasing habitat quantity and suitability, to prevent eventual species extinctions.