Habitat and climatic associations of climate‐sensitive species along a southern range boundary

Abstract Climate change and habitat loss are recognized as important drivers of shifts in wildlife species' geographic distributions. While often considered independently, there is considerable overlap between these drivers, and understanding how they contribute to range shifts can predict future species assemblages and inform effective management. Our objective was to evaluate the impacts of habitat, climatic, and anthropogenic effects on the distributions of climate‐sensitive vertebrates along a southern range boundary in Northern Michigan, USA. We combined multiple sources of occurrence data, including harvest and citizen‐science data, then used hierarchical Bayesian spatial models to determine habitat and climatic associations for four climate‐sensitive vertebrate species (American marten [Martes americana], snowshoe hare [Lepus americanus], ruffed grouse [Bonasa umbellus] and moose [Alces alces]). We used total basal area of at‐risk forest types to represent habitat, and temperature and winter habitat indices to represent climate. Marten associated with upland spruce‐fir and lowland riparian forest types, hares with lowland conifer and aspen‐birch, grouse with lowland riparian hardwoods, and moose with upland spruce‐fir. Species differed in climatic drivers with hares positively associated with cooler annual temperatures, moose with cooler summer temperatures and grouse with colder winter temperatures. Contrary to expectations, temperature variables outperformed winter habitat indices. Model performance varied greatly among species, as did predicted distributions along the southern edge of the Northwoods region. As multiple species were associated with lowland riparian and upland spruce‐fir habitats, these results provide potential for efficient prioritization of habitat management. Both direct and indirect effects from climate change are likely to impact the distribution of climate‐sensitive species in the future and the use of multiple data types and sources in the modelling of species distributions can result in more accurate predictions resulting in improved management at policy‐relevant scales.


| INTRODUC TI ON
Climate change is altering species distributions and shifting geographic ranges for many species along elevational and latitudinal gradients (Parmesan et al., 1999;Williams & Blois, 2018). These range shifts can result in altered community assemblages and novel, no analogue communities , particularly along range boundaries, as climate-sensitive species decline along trailing edges and other species expand ranges along leading edges (Williams & Jackson, 2007). Historically, research on the drivers of species range extents has focused on habitat as the dominant driver of species range extents; however, recent studies have indicated that particularly in climate-sensitive species and regions along southern range boundaries, climate may be equally as important and growing in importance as a driver of species range limits (Reich et al., 2022;Sultaire et al., 2016). These two drivers-habitat and climate change-do not operate independently and there is considerable overlap between them at both macro-and microscales. Climate change can alter hydrological patterns, fire regimes, and other mechanisms that directly affect biotic and abiotic components of habitat (Halofsky et al., 2020;Trenberth, 1999). Hence, species may experience direct physiological effects from altered precipitation or temperature patterns, but also indirect effects due to climate change impacts on habitats. A better understanding of the interactive effects of these drivers on species distributions will elucidate the mechanisms underlying range shifts and be important for managing biodiversity under global change.
Ecological communities and wildlife populations along ecosystem boundaries are often among the first to experience impacts of climate change (Gilman et al., 2010). The Northwoods region, spanning the Upper Great Lakes region of North America, contains one such ecological boundary that delineates where southern deciduous forests and northern coniferous ("boreal") forests meet (Andersen, 2005). While the boreal forest is one of the world's largest biomes (Gauthier et al., 2015), dynamics along the southern boundary can often differ greatly from those in the core biome, and the effects of climate change are already having a pronounced effect on community dynamics (Wilson et al., 2022). For example, subnivium habitat-the zone between fallen snow and terrain (Pauli et al., 2013)-is becoming more unstable along southern range boundaries due to the effects of climate change (Thompson et al., 2018). Because many species make use of the subnivium during winter months, and rely on the relatively stable temperatures within to persist during winter , the interaction of changing snow conditions and temperatures can negatively impact snow-adapted species. Indeed, recognition of the role of snow as habitat has prompted the creation of spatially and temporally explicit winter habitat indices to characterize dynamic snow conditions in space and time (Gudex-Cross et al., 2021).
In addition to climate-sensitive habitats, the Northwoods region is also home to multiple climate-sensitive wildlife species (Hoving et al., 2013), including important game species such as American marten (Martes americana), moose (Alces alces), snowshoe hare (Lepus americanus), and ruffed grouse (Bonasa umbellus). These species provide numerous cultural and ecosystem services yet are distinctly vulnerable to climate change. Martens rely on deep snow to avoid competition and direct predation by larger mesocarnivores (e.g.,

bobcat [Lynx rufus], coyote [Canis latrans], fisher [Pekania pennanti])
that have heavier foot loading and more difficulty moving through snow (Crête & Larivière, 2003). Additionally, loss of preferred prey species such as red-backed voles (Myodes gapperi) linked to habitat and climate change can negatively affect marten populations (Carlson et al., 2014;Scott et al., 2022). Moose experience direct effects from warming temperatures through increased heat stress, increased levels of parasitism by ecto-and endoparasites, and increased competition with white-tailed deer (Odocoileus virginianus), which are historically limited by deep snows (Weiskopf et al., 2019).
Additionally, moose may experience indirect effects due to negative impacts of climate change on boreal forest types in which they reside and forage (Weiskopf et al., 2019). Snowshoe hares have experienced northward range contractions along trailing edge boundaries (Burt et al., 2017;Diefenbach et al., 2016;Sultaire et al., 2016), partially driven by increased predation rates linked to mismatch in the timing of coat color molts with attenuated snow cover due to climate change (Wilson et al., 2019;Zimova et al., 2016). Ruffed grouse use the subnivium, the below-snow refugium that maintains a stable thermal environment during winter months, as a method to avoid predators and conserve energy during winter months (Pauli et al., 2013). The loss of consistent deep snow along southern range boundaries results in increased stress in individuals and declines in survival (Shipley et al., 2020). Despite extensive geographic overlap in distributions and general associations with colder temperatures and northern forest habitats, specific habitat requirements and mechanisms linking climate to vital rates can differ greatly among species, resulting in challenges for simultaneous management of these species in the region. For example, while early successional forests are preferred habitat for ruffed grouse (Rusch et al., 2020) and can increase survival in snowshoe hares (Wilson et al., 2019), marten prefer mature forest with structural complexity (Chapin et al., 1997). Therefore, there is a need to understand how species differentially respond to climate and habitat and identify areas of overlap where management action can have the broadest benefit for multiple species and ecosystems.
Here, we modeled the effects of climate, climate-sensitive forest cover types, and direct anthropogenic effects on the abovementioned four iconic vertebrate species in Michigan-a geographically unique state split by the Great Lakes into upper and lower peninsulas

T A X O N O M Y C L A S S I F I C A T I O N
Global change ecology, Landscape ecology, Spatial ecology ( Figure 1a). Our models utilized the Integrated Nested Laplacian Approximation (INLA; Rue et al., 2009) with Stochastic Partially Differential Equation (Lindgren et al., 2011), which is a computationally efficient and accurate method for analyzing spatial data and accounting for the effects of spatial autocorrelation. For our models, we integrated multiple forms of animal occurrence data, including radio telemetry, harvest records, and observations recorded in conjunction with citizen science programs. While data collected through formal field survey has inherent advantages, namely standardized sampling methods allow for more accurate modeling of detection probabilities, it is often expensive and yields a relatively small amount of data compared with community science data, which lack repeat sampling but can yield large amounts of data rather cheaply.
Focusing solely on one form of data can often yield misleading results (Kamp et al., 2016), whereas incorporating both can yield more robust estimates. To capture changing dynamics of snow as habitat, we also incorporated recently developed winter habitat indices, which are intended to more accurately describe the mechanisms affecting individual fitness (Gudex-Cross et al., 2021). Finally, we spatially predicted the distributions of our four focal species given model results in order to visualize the geographic heterogeneity of habitat suitability.
We predicted that multiple species would be associated with the same forest cover types, potentially allowing for mutually beneficial F I G U R E 1 Maps of study areas with occurrence data of focal species. (a) Map of Northwoods region of Michigan, with ecoregions labelled (per Albert, 1995), the greater Northwoods area is outlined in orange. (b) Maps of occurrence data collected and used to construct distribution models for Ruffed Grouse (Bonasa umbellus), moose (Alces alces), American marten (Martes americana), and snowshoe hare (Lepus americanus). management of forest habitat. We predicted that we would observe species-specific relationships with climate variables, but that measures of snow conditions (i.e., winter habitat indices) would outperform more general measures of climate, such as temperature, particularly for marten, grouse, and hares-species with recognized interactions with snow. Finally, we predicted that the southern extent and contiguity of distribution would vary between species as species respond to different habitat and climatic factors. These distributional models provide insight into effects of climate, habitat, and direct anthropogenic activity on climate-sensitive species distributions and allow for future projections of the direct and indirect effects of climate change on climate-vulnerable species in this region. Knowledge of how climate affects wildlife species directly and indirectly, via habitat alteration, can facilitate efficient use of resources to manage habitat at policy-relevant scales (e.g., local and state level) and buffer species from the negative impacts of climate change resulting in observed range shifts.

| Study area
Our study area comprised the Northern Lower Peninsula and Upper Peninsula of Michigan (Ecological section VII per Albert, 1995; Figure 1a), which combined is considered part of the Northwoods ecosystem. While these two peninsulas are unconnected and connectivity between populations in the lower and upper peninsulas is unlikely, they do share similar regulatory and conservation history. This region is highly seasonal with mean annual temperatures of 6.2°C, while mean summer temperatures and mean winter temperatures were 18.3°C and −6.4°C, respectively. Annual precipitation ranges from 71 to 86 cm across the study area, with lake effect snow being a defining feature across much of the coastal region and annual average snowfall ranging from 101 to 356 cm. Elevation changes in the region were minimal with a mean elevation of 296 m (SD = 88.6) and minimum and maximum elevations of 173 m in the southeast and 613 m in the western upper peninsula, respectively.

| Predictors
Seven forest cover types were considered to be under moderate-tohigh vulnerability from climate change (Handler et al., 2014). To represent them, we used a spatial layer containing the sum of the basal area (m 2 /ha) for two to five tree species representative of each cover type (Table 1; Dickmann & Leefers, 2003;Handler et al., 2014). Basal area is sum of cross-sectional surface areas of each live tree within a plot measured at breast height (Bettinger et al., 2017). Basal area for tree species were obtained from a US Forest Service data product containing raster maps of live tree basal area for tree species at a 250-m resolution (Wilson et al., 2013), obtained by integrating vegetation phenology from MODIS imagery and field plot data from the Forest Inventory and Analysis database between 2000 and 2009 (Wilson et al., 2013). This dataset was validated for spatial accuracy by comparing modelled data to observed data from Forest Inventory and Analysis (FIA) plots and comparing accuracy metrics (Riemann & Wilson, 2014).We created separate layers for each forest cover type by summing the basal areas for all representative species within each cell ( Figure S1). Basal area is a common measure of tree density used in forestry to represent aboveground biomass (Bettinger et al., 2017) as it contains information on number of trees and size (Balderas Torres & Lovett, 2013), and is often used to represent forest structure in resource selection models (Irwin et al., 2020;Parsons et al., 2019).
We used two types of climate predictors to model species occurrence. First, we used three different metrics of temperature: mean annual temperature, mean temperature of the warmest quarter (i.e., summer), and mean temperature of the coldest quarter (i.e., winter).
Mean values for temperature metrics from 1970 to 2000 were obtained at a 30-s resolution from WorldClim (Fick & Hijmans, 2017).
Next, we used three winter habitat indices developed to represent more functional aspects of winter ecology, specifically snow season length, snow variability, and duration of frozen ground without We used four different metrics to represent direct effects of human activity and urbanization: housing density, road density, distance to major roads, and distance to conservation lands. Both housing and road density are well-established indices of human density (Forman et al., 2003;Lewis et al., 2011), and represented levels of human activity and development within the study area.
Highways and other major roads can also serve as a barrier to dispersal and animal movements (Forman et al., 2003). Additionally, citizen science and other unstructured data sources frequently have detection biases in relation to roads due to increased accessibility (Cretois et al., 2021). Conservation and recreation lands comprised both private and public lands managed for either conservation (e.g., state wildlife areas and national forests) or recreation (e.g., hunting clubs and golf courses). Both protected areas and working lands (i.e., rangeland, agriculture, and forested areas used in commercial enterprises) can provide suitable habitat for a variety of species within human dominated landscapes and may act as refugia for climatesensitive species (Pacifici et al., 2020; but see Parks et al., 2023).
Housing density was represented by the 2010 household density (households/km 2 ) value obtained from a dataset detailing census block-level housing change between 1990 and 2010 (Martinuzzi et al., 2015). Road density (m/km 2 ) was represented by the cumulative length of road features obtained from Michigan's Open GIS Data portal (Allroads v17a; Michigan DOT, 2015) in each 1-km resolution grid cell. Distance to major roads (m) was represented by the distance to nearest road with a National Functional Classification

| Analysis
We implemented a Bayesian hierarchical approach in INLA, using the R package INLA (Rue et al., 2009) to model occurrence of each focal species using a generalized linear mixed model. INLA uses the Stochastic Partial Differential Equations (SPDE) approach for the spatial effect, approximating a Gaussian Random Field where the correlation between locations is Matèrn (i.e., covariance between two points is related to the distance between points). The random fields served as error terms to measure spatial autocorrelation and other uncertainty not explained by fixed effects included in the models. We constructed a random field (RF) by creating a two-dimensional triangulated mesh using guidelines provided by Lindgren and Rue (2015). We used mesh vertices to represent background locations, and supplemented these points with a regularly spaced grid over terrestrial regions.
We extracted values from each predictor layer for occurrence points and all background points. All continuous predictors were scaled by subtracting the mean and dividing by standard deviation.
The dataset for each species was then divided randomly into training (80%) and testing (20%) subsets. Observations of our focal species are modelled as a Bernoulli point process such that Z(s) indicates the presence (1) or absence (0) of the species at location s, with the probability of presence given as s . This relationship can be expressed as, where 0 is the intercept term, X s is a vector of values for each predictor, habitat , climate , and anthropogenic are the coefficients for each habitat, climate, and anthropogenic predictors respectively, RF s is the spatially structured random effect and ∈ NN(s) is a predictor included to capture spatially structured effects of aggregated individuals (i.e., clustering) (Humphreys et al., 2017).
We used a multistage process for model selection. First, we checked for issues with collinearity using the R package corrplot, and excluded any potential models that contained predictors with Bayesian information criterion, and while issues remain regarding the use of WAIC in structured datasets such as spatial models (Gelman et al., 2014), it is frequently used in hierarchical spatial modelling (see Leach et al., 2016;Sultaire et al., 2022;Williamson et al., 2022) and is preferred to other criteria such as DIC (Doser et al., 2022;Duncan & Mengersen, 2020). Coefficients of all top-ranked models were examined, and statistical significance was determined by comparing 95% credible intervals of the effect coefficients to zero.
To facilitate the interpretation of parameters, we converted beta coefficients for fixed effects to relative selection strength values by applying an exponential function to each coefficient value (Avgar et al., 2017). Relative selection strengths (RSS) can be interpreted as the relative intensity of use of the fixed effect by the species between locations that differ by one unit of the fixed effect if all other effects are equal (Fieberg et al., 2021). To visualize these relationships, we predicted relative probability of presence values at all available locations using the fixed effects and random field values for final top-ranked model of each species and plotted these values against the variable values for each predictor at the available location (Avgar et al., 2017). Relative probability of presence was scaled from zero to one by dividing each value by the maximum value. We then used the gam function in the R package mgcv (Wood & Wood, 2015) to fit a function to the available data, while using smoother parameters <5 to avoid overparameterization (Perrig et al., 2020).
To compare fitted and predicted spatial models (i.e., with and without inclusion of a spatially explicit random field), we created and spatial-effect model, respectively (Table 3). Additionally, the fitted model had higher values for TSS, sensitivity, and specificity (Table 3).
Fitted and predicted spatial models showed strong similarities and both predictions showed a mid-to-high probability of presence throughout the Upper Peninsula, with lower habitat suitability across the Lower Peninsula (Figure 4). Postmodel fitting, the random field showed only slight variation across the study region indicating that our model predictors explained much of the spatial variation present in the region ( Figure S3). Examination of ecoregions indicated that probability of presence was high across the Upper Peninsula, but lower in the Lower Peninsula ( Figure S4). Model validation indicated only minor improvements in model performance due to addition of the spatial effect, with only negligible improvements in AUC, TSS, sensitivity, and specificity values (Table 3).

| Moose
Relative habitat suitability for moose was positively related to basal area of Upland Spruce-Fir (β = 0.21, 95% CI = 0.015, 0.40; Figure 3).  without spatial effects, with AUC, TSS, sensitivity, and specificity values being negligibly improved by inclusion of the spatial effect variable (Table 3).

| Ruffed grouse
Ruffed grouse exhibited a positive relationship between habitat suitability and basal area of Lowland Riparian cover types (β = 0.13, 95% CI = 0.064, 0.20; Figure 3). An increase of 4.43 m 2 /ha in the basal area of Lowland Riparian species resulted in a 1.1× (95% CI = 1.1, 1.2; Figure 2) increase relative habitat suitability. The variable for Jack Pine was also included in the top-ranked habitat model ( distance to conservation lands (β = −1.2, 95% CI = −1.6, −0.72). Both spatial predictions with and without the spatial effect indicated grouse distribution across the Upper and Lower Peninsulas, but that probability of presence was closely linked to linear features, namely roads ( Figure 4). Additionally, the fitted model indicated a band of high relative probability of presence across the Lower Peninsula, while this band was absent in the predicted model ( Figure 4). Despite these differences in visualization of predictive models, the random field showed a relatively homogenous surface indicating that the variables included in the model were likely sufficient for explaining present spatial structure ( Figure S3). Examination of ecoregions indicated few areas with a high probability of presence ( Figure S4), though the Southern Superior Uplands and Northern Highlands still maintained relatively high levels of probability of presence. We observed notable differences in model goodness-of-fit between our predictive models with and without the random spatial effect.
Inclusion of the random field improved the model fit metrics AUC, TSS and specificity, though there was a slight decline in specificity (Table 3). Despite this improvement, both models had a relatively poor goodness-of-fit as compared to the other modelled species (Table 3).

| DISCUSS ION
As predicted, distributions of four climate sensitive wildlife species in the Northwoods region were driven by the effects of climate, habitat, and human activity. While the relative importance of climatic variables varied among species, the effect of habitat availability (portrayed as amounts of climate sensitive forest cover types) on occupancy probability affected all species. Indeed, the final model for each species contained multiple habitat predictors, and similarities among species were identified. Grouse and marten presence were positively associated with lowland riparian cover types, and moose and marten positively associated with upland spruce-fir cover type.
Additionally, snowshoe hares were positively associated with lowland conifer cover type, which was strongly correlated with lowland riparian cover type because Northern white-cedar occurred in both.
Finally, while jack pine-red pine appeared in the final models for all species, coefficient posterior means overlapped zero indicating that it was only weakly informative. However, preservation and conservation of jack pine-red pine stands benefit wildlife species, particularly if applied at optimal locations. For example, American marten in the northern Lower Peninsula select for den sites in large DBH trees located in mature red pine plantations, which is unique to this population (Sanders et al., 2017). As Jack Pine-Red Pine forests are vulnerable to the effects of climate change (Table 1), loss of these habitats will likely result in demographic consequences for this isolated population. Additionally, areas of high basal area likely represent older forests, while some species may preferentially use young forests (Pietz & Tester, 1983), which may not be apparent in our results. Important wildlife cover types in this study represent a range of vulnerabilities to climate change. Upland spruce-fir is considered highly vulnerable to climate change due to sensitivity to temperature and precipitation changes, while lowland riparian is considered moderately vulnerable (Handler et al., 2014;Reich et al., 2022). While specific mechanisms threatening each of these cover types differ, generally each is threatened by changes to the hydrological cycle, insect or disease outbreaks and increased herbivory by white-tailed deer (Handler et al., 2014). Due to the multitude of threats to their distinct cover types, it is likely that each of these vertebrate species will experience negative indirect effects from habitat loss due to climate change.  (Evans & Mortelliti, 2022), yet we observed no relationship between marten habitat suitability and any climatic variables, including those describing snow.
While we had predicted that more mechanistic measurements of local climate, such as winter habitat indices, would outperform broad temperature metrics, there are several factors potentially affecting this result. First, at the spatial extent and resolution at which we examined these relationships all climate predictors were strongly correlated and were not included in the same candidate model. It is likely had we looked at finer resolutions (i.e., closer to the 500-m resolution which the habitat indices were created at), and across a broader geographic range we may have seen more variability in the winter habitat indices and potentially stronger relationships.
Furthermore, while these indices have been shown to be strong predictors of species richness (Gudex-Cross et al., 2022), they may be less suited for constructing distribution models from presence only data. Although our analysis was unable to determine the more specific, mechanistic drivers of range contraction, we do show that in some circumstances simpler metrics such as temperature may serve as a valuable proxy for these mechanistic indices. This can be valuable as models predicting temperature under future climate change scenarios have higher confidence than those predicting precipitation (Kapnick & Delworth, 2013). Finally, recent advances in the modelling of climate-sensitive species have shown that spatial non-stationarity may be particularly important for populations along range boundaries and that the effect of climate variables on habitat suitability may vary spatially (Sultaire et al., 2022). Indeed, studies of snowshoe hares incorporating nonstationarity have indicated that at broad scales, snow cover defines distributions, but temperature modulates the strength of the relationship across space (Sultaire et al., 2022). Incorporating nonstationarity into species distribution models may result in more precise estimates of the effects of climatic variables, and improve our understanding of the interactions between space, climate, and habitat that define species' range boundaries . Nonstationarity is likely to be particularly important along range boundaries, as dynamics can often differ greatly from those in the core of distributional ranges (Sexton et al., 2009).
We observed strong anthropogenic effects on predicted occurrence of all focal species. Many of these relationships, such as the negative relationships with housing density observed by martens, hares and grouse, likely reflect the tendency of wildlife to avoid areas dominated by humans (Lewis et al., 2021); however, some effects may reflect bias in the sampling method in which location data were collected, a well-known problem with unstructured citizen science data (Dickinson et al., 2010). Many wildlife species show spatial and temporal avoidance of areas of high human activity, typified by high housing and road density and accompanying anthropogenic noise and light pollution (Gaynor et al., 2018;Wilson et al., 2021). Martens, hares, and grouse exhibited evidence of this avoidance, with increased probability of presence near conservation areas and in areas of low housing density. These results match the ecology of these species, particularly in relation to habitat selection, with martens reliant on large tracts of forest with vertical structure, hares reliant on forest patches with high stem densities or complex understories(particularly young aspen and late-seral conifers), and grouse reliant on early successional forests during spring and fall. These results also underscore importance of connectivity among conservation areas, as each of these species are expected to or are currently undergoing range contractions (Burt et al., 2017;Pomara & Zuckerberg, 2017). Increased connectivity may allow for prolonged persistence of populations along southern range edges by allowing for immigration and emigration among isolated patches within climate or habitat refugia and populations within the core of the species' range. We observed other relationships that more likely resulted from how occurrence data were assembled. For instance, the increased probability of grouse presence with increased road density likely stems from use of eBird citizen science data as the only data source for grouse presence. Similar relationships were observed for marten and moose, likely due to the proximity of trapping locations to accessible roads and the use of roadkill data for marten and moose, respectively. Indeed, moose typically exhibit scale and context-dependent relationships with roads, showing avoidance of landscapes with high road density and avoiding major roads at finer scales, but potentially showing preference for smaller roads and trails when they facilitate movements in deep snow or allow access to forage and sodium (Beyer et al., 2013;Laurian et al., 2008;Wattles et al., 2018). Spatial bias is a well-known problem within citizen science data collection (Tiago et al., 2017), and inclusion of additional data sources where surveys were conducted away from roads may have improved this model. It is also important to note that we only examined the presence of these species rather than population densities, and effects of climate and habitat on abundance rather than mere occurrence may reveal otherwise hidden effects on population dynamics. For instance, further examination of grouse abundance and population cycles are likely to reveal declines across much of the southern range boundary as declines have been observed in other regions and linked to climate variability (Pomara & Zuckerberg, 2017).
Our comparison between spatial models with and without use of a random spatial effect revealed some advantages of using random fields in fitting species distributions. Use of the spatial-effect nominally improved model goodness-of-fit for most species, but was particularly important for American marten. This was likely due to the abundance of telemetry locations at a number of focal areas used in the marten dataset, which made accounting for spatial autocorrelation critical. Examination of the random field post-model fitting revealed spatial regions where our predictors did not adequately predict the spatial structure observed from the data and provided further insight into each species' history in the study region that would have been unavailable otherwise. For example, fitted distributions of marten showed two discontinuous distributional patches in the Lower Peninsula, while these patches were absent in the pre- Peninsula. This is due to suitable climatic and habitat conditions in this area, rather than actual detections of moose, but indicate that at present this region could sustain a population of moose. However, suitable habitat/climate space along this southern edge is unlikely to persist given future projections of climate change in the region.
The measurement of the random field and resulting large values in this region likely indicate that climate and habitat alone cannot fully explain the presence or absence of these populations, and management action (i.e., translocations and regulatory protection) must be accounted for to fully understand the drivers of species distribution.
As community science data continue to become more commonly used and the use of presence-only models continues to proliferate, random spatial effects may be an effective diagnostic tool to identify areas where unexplained spatial variation occurs and to evaluate alternative drivers of spatial distribution. Despite the labelling of each of our focal species as climatesensitive, contemporary responses to climate and habitat varied greatly between species indicating that future responses to climate change may differ as well. Evidence of these differential responses may already be observable when considering the variation of suitable habitat and climatic space along the southern edge of these species' ranges. Indeed, differences in the precise mechanisms that negatively influence fitness and survival can influence the velocity at which the distributional ranges of these species shift northward.
Additionally, while we have looked at each of these species individually, climate change is also altering biotic interactions within these ecological communities (Blois et al., 2013), which will undoubtedly affect the population dynamics of these species as habitat and climatic conditions change. Nonetheless, our analyses have emphasized the importance of climate-sensitive forest cover types for climate sensitive vertebrate species, highlighting the likelihood of both indirect and direct effects of climate change on these species. While climate change and its direct effects may be most actionable on global or regional scales, habitat management can occur at local scales and the indirect costs incurred by the effects of climate change on habitat can potentially be mitigated by effective management, slowing the velocity of climate change and potentially creating refugia for species unable to persist given the negative effects of climate change.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data (Wilson et al., 2023a)