Forecasting species distributions: Correlation does not equal causation

Identifying the mechanisms influencing species' distributions is critical for accurate climate change forecasts. However, current approaches are limited by correlative models that cannot distinguish between direct and indirect effects.


| INTRODUC TI ON
A variety of extrinsic and intrinsic factors influence a species' distribution, including climate, biotic constraints, demographic limitations, biogeographic history, and evolutionary adaptations (Connallon & Sgrò, 2018;Sexton et al., 2009). These factors have varying degrees of influence that can be attributed to the spatial and temporal scale (Wiens & Graham, 2005) and range position (Louthan et al., 2015;Sirén & Morelli, 2017). Prior research has focused more on understanding correlative relationships, especially the role of abiotic factors on the formation of geographic ranges (i.e., a species fundamental niche; Sexton et al., 2009), due to the increasing threat of global change. This approach is logical as climate is often considered to be the ultimate driver of species distributions (Wiens & Graham, 2005). Correlative characterizations of distribution dynamics, however, risk overlooking potentially important causal mechanisms that determine species distributions and range limits (Filazzola et al., 2020;Lyons & Kozak, 2020). Given the current unprecedented rate of climate and habitat change, understanding the determinants, rather than correlates, of current range limits is a critical component for conserving biodiversity.
The concurrent development of new theory (Godsoe et al., 2017;Sirén & Morelli, 2017), statistical approaches (Lefcheck, 2016), and field methods (Sirén et al., 2018;Steenweg et al., 2017) provides an opportunity to evaluate the causal mechanisms that form species range limits. A long-standing hypothesis, dating back to at least Charles Darwin, is that negative biotic interactions shape the lower latitudinal limits of a species' range (Louthan et al., 2015). However, support for this idea is mixed (Cahill et al., 2014), potentially due to the paucity of data on biotic variables at broad spatial scales (Wisz et al., 2013) or because biotic interactions are only prominent at local scales (Soberón, 2007). Another explanation is that abiotic factors can be correlated with biotic interactions and mask the effects of competition and predation (Godsoe et al., 2017). Indeed, many trailing edge populations have positive associations with abiotic stress, supporting the hypothesis that climate mediates biotic interactions along lower range limits (Sirén & Morelli, 2017). This theoretical framework, along with quantitative approaches that allow for the inclusion of direct and indirect effects, provides the opportunity to disentangle the factors that form range limits.
In a previous article, we showed that a combination of direct and indirect abiotic and biotic factors shape the current regional range limits of North American mammalian carnivore assemblages (Sirén et al., 2021). In particular, snow had a strong indirect effect on the distribution of Canada lynx (Lynx canadensis) by mediating the negative effects of its competitors, primarily the closely related bobcat (Lynx rufus). However, this pattern was not consistent for all species.
For example, snow had a stronger direct effect on American martens (Martes americana) than its indirect effect, which was manifested through coyotes (Canis latrans). Furthermore, biotic factors such as habitat and prey availability imparted strong bottom-up effects on range limits (Sirén et al., 2021). Although previous work has identified the importance of abiotic and biotic factors for many species along trailing edges (Aubry et al., 2007;Jensen & Humphries, 2019;Peers et al., 2013;Sultaire et al., 2016), correlative approaches may limit our understanding of the underlying mechanisms. Given that quantifying abiotic and biotic relationships is central to predicting future species distributions under climate change, an obvious and important question is as follows: are predictions of current and future patterns consistent when contrasting correlative and causative analytical frameworks?
This study focused on data from a large-scale, 6-year camera trapping effort at 257 sites across the U.S. states of Vermont and New Hampshire (Figure 1). Although other carnivores, such as coyotes and fishers, (Pekania pennanti) negatively influence lynx (Bayne et al., 2008;McLellan et al., 2018), our previous work, and that from other regions, indicates that bobcats have the greatest impact on lynx occurrence (Peers et al., 2013;Scully et al., 2018;Sirén et al., 2021). Furthermore, we only found coyotes, rather than fishers, to have a negative effect on marten occurrence (Sirén et al., 2021), supporting recent work in the region (Jensen & Humphries, 2019).
However, the indirect effect of snow on martens through coyotes was weaker than that of lynx vs. bobcats (Sirén et al., 2021). Therefore, we focused our analysis on species pairs (Canada lynx vs. bobcat; American marten vs. coyotes) that had the strongest total effect on each other, yet differing relationships with snow.
We used causal and correlative models to predict current (2014-2019) and future (2080s) distributions of two pairs of species with abiotic and biotic factors prominently influencing their spatial distributions (Canada lynx vs. bobcats; American martens vs. coyotes).
We hypothesized that including both direct and indirect effects via a structural equation modeling (SEM) approach (Lefcheck, 2016), that is, a causal modeling framework (Grace & Irvine, 2020), would predict lynx distributions better than a regularly applied correlative projection approach (occupancy modeling; MacKenzie et al., 2017).
Specifically, we expected that a failure to jointly account for the mediating effect of snow depth on biotic interactions (e.g., competition) and bottom-up effects from habitat availability (modeled forest biomass) would result in poorly resolved predictions of both current and future distributions. Moreover, we expected these differences to vary between species pairs depending on the relative importance of abiotic versus biotic determinants of species range limits. Our inclusion of martens vs. coyotes provided a comparison with lynx vs. bobcats to evaluate these expectations. To consider the range of uncertainty introduced by climate modeling, we compared changes in future distributions using two global climate models (GCMs), F I G U R E 1 Location of 257 remote camera sites (black dots) in the northeastern US (pink dot within inset map of USA and Canada) used to evaluate differences between causal and correlative species distribution models. The camera design (upper left inset) included a snow stake with visual and olfactory lures placed 3-5 m from a camera. Elevation within the region ranged from sea level (green) to 1906 m (white) and that of the cameras ranged from 3-1487 m | 759 SIRÉN et al.
HadGEM2-ES and GISS-E2-R, that projected a range of snow depths under the RCP8.5 scenario. Our formal evaluation of how projections of species distributions vary depending on assumptions and strengths of abiotic and biotic relationships offers the much needed empirical context to advance the discussion about the importance of considering biological mechanisms in projections of ecological systems (Urban et al., 2016).

| Study system and data collection
We used data from 257 camera-trap sites operating from 9 January 2014-8 May 2019 (Figure 1). Cameras were established in nonoverlapping 2 × 2 km grids and placed along game trails to increase detection probability. Sites were established along altitudinal transects and a latitudinal gradient to capture the climate and landcover/land-use gradients of the region (Sirén et al., 2021). Cameras were positioned facing north on a tree, 1-2 m above the ground, and pointed downward at a 2-5° angle toward a snow measurement stake positioned 3-5 m away. We used skunk (Mephitis mephitis) lure and wild turkey (Meleagris gallopavo) feathers as olfactory and visual attractants, respectively, to increase detection probability for occupancy modeling. Previous work has shown that this approach outweighs any positive bias imparted by lures that may occur at local scales (see Stewart et al., 2019). We set cameras to take 1-3 consecutive pictures every 1-10 s when triggered and checked sites on an average of three (1-9) times each season to refresh attractants, download data, and ensure that cameras were operating properly.
To generate species occurrence data for the causal and correlative models, we restricted camera data from autumn to spring (16 October-15 May) for each year (2014-2019). This seasonal range was chosen as it approximates demographic and geographic closure and is based on species behavioral responses to leaf and snowpack phenology of the region (Sirén et al., 2016;Vashon et al., 2008). We identified species in photographs and used consensus from multiple observers when identification was ambiguous (Thornton et al., 2019). We organized camera data into weekly occasions using CPW Photo Warehouse (Ivan & Newkirk, 2016) and output detection/ non-detection data for each species for use in causal and correlative models. In total, we tallied 69, 487, 188, and 740 weekly detections of lynx, martens, bobcats, and coyotes, respectively. A complete summary of the detection history and number of sites occupied by each species is provided in Table S1.

| Snow and forest biomass data
We used gridded snow depth (Wang et al., 2016) and forest biomass (Duveneck & Thompson, 2019) datasets that included current and future projections of these variables under a high greenhouse gas emission scenario (representative concentration pathway [RCP]; hereafter, RCP8.5; Meinshausen et al., 2011) to model species distributions. We chose RCP8.5 because the forest biomass data were limited to this scenario, and the purpose of the study was to understand the sensitivity of species distributions to substantial changes in snow with future warming. We used ClimateNA software (Wang et al., 2016) to generate gridded outputs of snow depth (precipitation as snow; mm), using a resampled digital elevation model from the National Elevation Dataset (http://ned.usgs.gov/) that matched the resolution of the forest biomass dataset (250 m). We generated snow data for the current period (2014-2019) and for the future period (2080s: average for years 2071-2100; Figure 2). To confirm that our choice of the 'current' period was not affected by large year-to-year variability in precipitation, we compared the snow depth values over 6 years (2014-2019) sampled in this study with a 30-year average . This comparison confirmed that the 6-year period we used to calculate projections reflects baseline/ normal snow-depth conditions for the current period ( Figure S1).
We chose the HadGEM2-ES and GISS-E2-R GCMs for the future period because, compared to other GCMs, they best reproduce winter precipitation patterns for northeastern US, yet have differing projections for temperature (Karmalkar et al., 2019). That is, HadGEM2-ES has about 2°C (3.6°F) additional warming compared to GISS-E2-R. Thus, these GCMs provide a plausible range of future snow depths in the region with GISS-E2-R projecting a minimal decrease from the current period and HadGEM2-ES projecting an extreme decrease ( Figure S2).
For forest biomass predictions, we leveraged previously published forest change simulations that provided spatial layers of forest biomass density over time (Duveneck & Thompson, 2019). These layers included estimates of aboveground forest biomass density (g/m 2 ) simulated with the LANDIS-II forest landscape model (Scheller et al., 2007). Aboveground biomass included successional interactions of seed dispersal, competition, cohort growth and mortality, and water balance. We included these biomass estimates for the current and

| Modeling approach
We analyzed camera trap data using correlative and causal frameworks to predict current and future distributions of focal species (lynx, marten) and competitors (bobcat, coyote) across their range in northeastern US.
For correlative predictions, we fit a static occupancy model (MacKenzie et al., 2017) for each species with snow depth and forest biomass as occupancy predictors and a combination of observational and site covariates to evaluate detection probability (Table S2). 'Year' was fit as a fixed effect in the occupancy and detection sub-models to account for multiple years of sampling (Table S2). Prior to modeling, all detection covariates were screened for multicollinearity using a variance inflation factor (VIF) test (Zuur et al., 2010); all covariates had VIF scores <2, indicating no multicollinearity and were included in the modeling process. We evaluated all possible combinations of detection covariates, resulting in a total of 48 models, while fitting a constant occupancy model (snow + forest biomass + year). Model performance was evaluated using Akaike Information Criterion (AIC) with the 'AICcmodavg' package (Mazerolle, 2017), and the most parsimonious model with the lowest AIC score was considered the best model (Burnham & Anderson, 2002). All occupancy analyses were performed in R (R Core Team, 2019) using the 'unmarked' package (Fiske & Chandler, 2011). We also conducted goodness-of-fit tests for each species to evaluate assumptions of closure and determine how well the models fit the data; tests were performed using the 'parboot' function in the 'unmarked' package, running 1,000 bootstrapped iterations on the top performing model. For the causal model, we used structural equation modeling (SEM) (Lefcheck, 2016) and a modification of the causal system described by Sirén et al. (2021) that included only snow depth and forest biomass as exogenous variables and bobcats and coyotes as competitors/predators of lynx and martens, respectively. SEM is generally described as a series of univariate regressions within a causal graph or network of paths that provides for the evaluation of complex and competing hypotheses (Grace, 2008). Correlative approaches are well-suited for studying single processes or responses and, of course, only measure associations. SEM, on the other hand, allows one to study the direct and indirect effects that influence processes within systems. For an overview of SEM and how it can be used for evaluating ecological hypotheses, see Grace (2008) and Grace and Irvine (2020). We chose a simplified causal system, F I G U R E 2 Current (2014-2019) and projected snow depth (mm) for the 2080s and forest biomass (g/m 2 ) for 2060 given a high-emission scenario (RCP8.5) in northeastern US. The projected snow depth for the 2080s was generated using the HadGEM2-ES global climate model without prey species, to provide a more similar comparison with the correlative models and because it is not yet known how to obtain model predictions and standard errors for complex piecewise SEMs (J. Lefcheck, personal communication). It should be noted that for the SEM, we used 'detection corrected' occupancy estimates from static occupancy models with a null state model and maximal detection model (Duclos et al., 2019;Mills et al., 2020); estimates for each species were obtained using the 'ranef' and 'bup' functions in 'unmarked' (Fiske & Chandler, 2011).
We performed a SEM for each species pair (lynx vs. bobcat; marten vs. coyote) using the 'piecewiseSEM' package in R (Lefcheck, 2016). Each SEM contained a series of univariate generalized linear mixed-effects model (binomial distribution with logit-link function) fit using the 'lme4' package (Bates et al., 2015); camera sites were specified as a random effect due to repeated measurements and variability in effort across years. The goodness-of-fit of each SEM was evaluated using Pearson's χ 2 statistic of a Fisher's C test, where a p > .05 indicates adequate fit of the observed data and conditional independence (Shipley, 2009). If SEMs fit the data, the path coefficients (i.e., relationships between nodes) from the univariate regressions were used to calculate direct and indirect effects using a standard approach. In brief, direct effects are interpreted as path coefficients between connected nodes, and indirect effects are those separated by one node and calculated by taking the product of two direct path coefficients (Shipley, 2009).
Predictions and uncertainty (standard errors) under the correlative occupancy modeling framework were predicted directly using standard prediction routines (Kéry et al., 2013). The conditional structure of the SEM, on the other hand, required a sequential bootstrapping approach, where realizations of each model in the conditional hierarchy of the SEM were simulated and included in the next sub-model. Operationally, this approach required simulating spatially explicit occupancy states from the competitor models (bobcat and coyote) and carrying those over into the focal species models (lynx and marten), and repeating 1000 times. We calculated standard errors from these simulations to map uncertainty in realized occupancy. We then used each model to make two sets of landscape predictions of occupancy for the competitor and focal species: one for the current period (2014-2019) and one for the future (2080s), using snow depth and forest biomass as predictor variables.

| Model performance
To evaluate the predictive ability of causal and correlational models from the current period, we compared spatially explicit predictions of occupancy probability with independent data on lynx and martens.
We used snow track survey data collected during the same timeframe and region as our primary independent data source ( Figure S4). Briefly, we conducted snow track surveys along established routes 1-6 times/winter following a protocol designed for detecting mesocarnivores (Squires et al., 2012). All track intercepts of lynx and martens were recorded with a GPS (Garmin GPS 62S, Garmin International), and detections/non-detections were spatially assigned to a 2 × 2 km grid to match our camera survey data ( Figure S4). In total, 91 surveys (X = 2; range = 1-6 surveys/year) were conducted over a 6-year period, covering 476 unique grids and 2536 km (mean = 2.27 km/survey). We used 50 detections of lynx and 889 detections of martens from these surveys to validate the causal and correlative occupancy models. Because snow track surveys encompassed only a portion of the area we sampled with cameras (see Figure S4), we assigned six lynx and 15 marten occurrence records collected by state agencies (New Hampshire, Vermont) during the same timeframe to the grids. Given that lynx and martens are listed as Species of Greatest Conservation Need and both states have been surveying them extensively over the past decade, we assumed lynx and martens were absent for any surveyed grid without an occurrence record. We did not have this information for bobcats and coyotes, so we could not make this assumption; therefore, we did not evaluate model performance of these species. To evaluate the performance of each model to predict occurrence, we used the area under the receiver operator characteristic curve (AUC). Following previous guidelines, we deemed models with AUC values (1) between 0.5 and 0.7 to have low accuracy, (2) between 0.7 and 0.9 to be adequate at predicting occurrence, and (3) >0.9 to be highly accurate (Manel et al., 2001).
Finally, to evaluate differences in spatial predictions (and uncertainty of the predictions), we compared differences in predictions of both periods between modeling approaches and, for future projections, choice of GCMs. Our sampling was within the range of variation for the region with focus more on forested areas with deeper snow ( Figure S3); as such, we assumed our comparisons between methods, and periods were valid because this pattern was consistent across space and time.

| Correlative models
The top-performing detection models included several covariates that explained detection probability for each species (Table S3); covariates from top models were held constant (along with 'Year' fit on the detection and occupancy submodels) to evaluate the effect of snow depth and forest biomass on occurrence. For all species, the summed square of residuals (SSEs) of the top models were well within the distribution of the bootstrapped SSEs, indicating they fit the data well (Table S4).
For correlative models, lynx occupancy was positively correlated with snow depth and negatively correlated with forest biomass (Table 1, Figure 3a). Bobcat occupancy, on the other hand, was negatively correlated with snow depth, whereas forest biomass had a weak positive effect (Table 1, Figure 3a). These patterns in occupancy were similar for martens and coyotes (Table 1, Figure S5a).

| Causal models
The causal model fit the data well for lynx and bobcats (Fisher's C = 0.44, p = .80) and martens and coyotes (Fisher's C = 4.71, p = .10), indicating conditional separation, which is a core assumption of piecewise SEM (see Shipley, 2009). Lynx occupancy, evaluated using a causal framework, was influenced by a combination of abiotic and biotic factors. Snow depth had a direct and indirect positive effect on lynx occupancy, the latter realized via a path through bobcats (Table 1, Figure 3b). Specifically, snow depth had a direct negative effect on bobcat occupancy, and bobcats had a direct negative effect on lynx occupancy (Figure 3b).
Forest biomass also had a negative effect on lynx occupancy; lynx were more likely to occur in areas with lower forest biomass (  Figure S5b).

| Model performance
We evaluated the predictive ability of lynx and marten causal and correlative models using independent data collected in the same region and timeframe as our camera data. The predictive performance, as measured using the area under the receiver operating characteristic curve, revealed that the causal model was better at predicting occurrence (AUC = 0.89) than the occupancy model (AUC = 0.80; Figure S6) for lynx. However, there was no difference in predictive ability between the causal (AUC = 0.91) and correlative (AUC = 0.91) models for martens ( Figure S7). Asterisks indicate significance level (*p < .05, **p < .01, and ***p < .001).

| Model differences for the current period
Predictions for lynx distribution during the current period (2014-2019) under the causal modeling framework were generally similar to the correlative model, but there were some notable differences ( Figure 4). For both models, predicted occupancy was highest in the northeastern and high-elevation areas of the region (Figure 4).
However, the correlative model also predicted high lynx occupancy in areas where they are absent (western and southeastern regions; Figure 4). In contrast, the causal model predicted the western and southeastern regions to have low occurrence due to high predicted occurrence of bobcats in these regions (Figure 4). Predictions for marten distribution, however, were similar between modeling approaches ( Figure 5). The primary difference was that the causal model predicted slightly lower occupancy in the region and sharper gradients than the occupancy model ( Figure 5). A notable pattern for both species was that predictions were more uncertain in areas with lower occupancy for the causal models, and the opposite pattern prevailed for correlative models (Figures S6 and S7).

| Model differences for the future (2080s)
Lynx distribution, as predicted using causal and correlative models, changed in the region under a high-emissions scenario for the 2080s (Figure 4). Both modeling frameworks predicted an increase in lynx occupancy by the 2080s in the northeastern part of the region (Figure 4), although the correlative model indicated higher uncertainty than other areas ( Figure S8). However, the causal lynx model predicted a lower increase in occupancy than the correlative model ( Figure 4). Lynx occupancy, according to the correlative model, was predicted to increase in areas they were unlikely to exist in the current period (e.g., western region; Figure 4). Although the causal model also predicted a slight increase in these regions, predicted lynx F I G U R E 4 Canada lynx (Lynx canadensis) and bobcat (Lynx rufus) distributions (occurrence probability) in northeastern US, as predicted using causal and correlative (corr) models for the current period (2014-2019) and the future (2080s) given projected changes in snow depth (mm) and forest biomass (g/m 2 ) under a RCP8.5 emission scenario. The projected snow depth for the 2080s was generated using the HadGEM2-ES global climate model. The third column indicates % difference in distributions between the current and future periods.
Positive differences indicate an increase in occupancy and negative values indicate a decrease occurrence was lower than that of the correlative model, low due to high predicted occurrence of bobcats for the 2080s (Figure 4). Both modeling approaches predicted an overall decline in lynx occupancy by the 2080s in the high-elevation regions, except for elevations >1200 m (Figure 4). However, the causal model predicted a slightly greater decrease due to increases in projected bobcat occupancy ( Figure 4).
Predicted marten distribution in the 2080s was similar for causal and correlative approaches; both predicted a decline in occupancy in the central and southern regions along its trailing edge and an increase in the northeastern region ( Figure 5). However, the correlative model indicated higher uncertainty in the northeastern region ( Figure S9). Similar to patterns in the current period, the causal models predicted more abrupt gradients in marten occupancy in the 2080s than the correlative model ( Figure 5).
We did not find evidence that predictions of future distributions under either modeling approach were sensitive to the choice of GCMs ( Figure S10). For lynx, there was little contrast between predicted occupancy using different GCMs; rather, the difference was attributed to the modeling framework ( Figure S10). A greater decline in marten occupancy was predicted for both modeling approaches using the HadGEM2-ES GCM that projects an extreme decrease in snow depth in the 2080s, relative to the GISS-E2-R GCM (Figures S2 and S10).
This result likely occurred because snow depth had a stronger positive effect on martens than lynx (Table 1). However, the same patterns emerged for each species, invariant of the GCM, indicating that differences in future distributions were due to the modeling frameworks.

| DISCUSS ION
There is overwhelming evidence that climate and land-use change alter wildlife species ranges even within short time frames (Sexton et al., 2009;Sirén & Morelli, 2017). However, the direct and indirect F I G U R E 5 American marten (Martes americana) and coyote (Canis latrans) distributions (occurrence probability) in the northeastern US, as predicted using causal and correlative (corr) models for the current period (2014-2019) and the future (2080s) given projected changes in snow depth (mm) and forest biomass (g/m 2 ) under a RCP8.5 emission scenario. The projected snow depth for the 2080s was generated using the HadGEM2-ES global climate model. The third column indicates % difference in distributions between the current and future periods.
Positive differences indicate an increase in occupancy and negative values indicate a decrease influence of these factors on biotic interactions has been largely unexplored due to limitations from correlative modeling frameworks.
Although other studies have used causal models to identify factors that influence species distributions (Mason et al., 2014;Mills et al., 2020;Wallingford & Sorte, 2019), including our previous study (Sirén et al., 2021), these studies were explanatory in nature. To our knowledge, our article is the first to generate spatially explicit predictions of species distributions using a causal modeling framework.
We show that prediction bias identified during the current period was related to the degree of complexity and strength of interactions of species embedded in larger systems. Importantly, our study indicates that ignoring indirect influences and biotic interactions risks erroneous forecasts for some species and provides continued support for taking a mechanistic approach when projecting ecological responses to climate change (Urban et al., 2016).
The causal lynx occupancy model that incorporated an indirect effect of snow through the phylogenetically similar bobcat had higher predictive power than the correlative one. These findings uphold previous work that suspected climate-mediated competition between these species (Hoving et al., 2005;Peers et al., 2013;Scully et al., 2018), yet could not rule out alternative hypotheses due to the use of correlative frameworks. Although the correlative model predicted lynx occurrence in similar regions, it also predicted occurrence in areas where lynxes do not occur (e.g., western Vermont, southeastern New Hampshire). These regions are highly developed (e.g., agriculture) with low forest biomass and shallow snow. Because forest biomass had a stronger effect than snow depth on lynx occurrence, according to the correlative model (see Figure 3), these regions were predicted to have high lynx occurrence. However, these regions were also predicted to have high bobcat occurrence because of shallow snow depth. Using a causal modeling framework, we were able to explicitly incorporate these interactions, and the model accurately predicted lower lynx occupancy in these regions. This finding has meaningful implications for lynx conservation, considering that this federally threatened species was recently recommended for delisting in the conterminous US (US Fish and Wildlife Service, 2017).
Contrastingly, predicted marten distribution was similar between modeling approaches during the current period, likely because the indirect effect of snow (through coyotes) was weak, and coyotes did not have a demonstrable direct effect on martens ( Figure S5).
However, the abrupt gradients in occupancy highlighted by the causal model that coincided with areas of deep snow (e.g., mountains) may indicate the marginal negative effect that coyotes have on martens and support previous work in the region that highlights the cascading effect that this top carnivore imparts on the species (e.g., coyotes → fisher → martens; Jensen & Humphries, 2019). The weaker negative interaction between these species contrasts with our previous work that included a causal model with several carnivore and prey species and different snow and forest biomass data (Sirén et al., 2021). However, we were unable to evaluate a complex model because it is currently unclear how to generate predictions and propagate errors for causal models with numerous paths (J. Lefcheck, personal communication). Regardless of this, the dynamics between martens and coyotes provides an important contrast with lynx and bobcats to illustrate how muted biotic interactions generate similar distributions using causal and correlative frameworks.
As hypothesized, the similarities and differences in lynx occupancy during the current period, as predicted using causal and correlative models, were propagated into the future. Both modeling approaches forecasted a decline in lynx occupancy at high elevation in the central region, and an increase in the northeast by 2080. However, as predicted, the changes highlighted by the causal model were less extreme than the correlative model. Furthermore, the correlative model forecasted a marked increase in the western region likely due to land-use change in southern New England that reduced forest biomass (Duveneck & Thompson, 2019). The causal lynx model, which predicted low occurrence in this region during the current period, only forecasted a slight increase in occupancy because bobcats were also predicted to increase in this region. A more likely scenario is the increase in lynx occupancy in the northeastern region which, despite a forecasted decline in snow depth, is projected to have more early successional habitat due to large-scale forest management (Duveneck & Thompson, 2019). These findings are consistent with those of previous work that highlight the influence of land management on lynx expansion in the region (Simons-Legaard et al., 2016).
The predicted future distribution of American martens, however, was similar between modeling approaches, likely due to the lack of a strong interaction with coyotes found during the current period.
Accordingly, our models predict that marten distribution will retract northward in the region due to the decline in snow depth projected by GCMs (see Figure 5, Figure S10). Previous research has indicated that martens have a strong positive association with snow due to its potential mediating effect on competitors (Jensen & Humphries, 2019;Sirén et al., 2021;Zielinski et al., 2017) and its role as habitat for foraging and thermoregulation (Buskirk et al., 1989). Thus, a decline in snow depth will likely have a detrimental impact and result in range contraction along its trailing edge (Carroll, 2007;Wasserman et al., 2012). Although this decline may be buffered by a projected increase in forest biomass (see Figure 2), we did not find a strong effect of forest biomass on predicted marten occupancy.
Our evaluation of the choice of GCMs for snow depth projections in the 2080s indicated that the primary differences in predicted distributions were still attributed to the modeling approach.
For example, the two GCMs provided different projections of future snow depth, yet the patterns generated by the casual and correlative modeling approach during the current period were propagated into the future for lynx and martens (see Figures 3 and 4; Figure S10). The only difference is that the GCM with higher warming (HadGEM2-ES) resulted in more of a decline in marten distribution than the GCM with lower warming (GISS-E2-R), for both modeling approaches. This finding can be explained given the strong positive effect that snow depth had on marten occupancy and the higher decline in snow depth projected by HadGEM2-ES. Interestingly, although there were differences in predicted lynx distribution that varied by causal and correlative models, that the choice of GCMs did not impart a strong effect as snow was less important for lynx, and forest biomass also imparted a strong effect compared to martens. Combined, these findings indicate that although there may be relative differences in the choice of GCMs, which varies based on the strength of that variable (e.g., snow depth) on a state parameter (e.g., distribution), the primary differences are still due to the choice of modeling approach and the degree of biotic interactions. These findings extend previous research that the choice of GCMs contributes toward differences in predicted distributions (Filazzola et al., 2020;Lyons & Kozak, 2020).
However, our study provides evidence that strong biotic interactions may nullify these differences.
We explored a new approach to predicting species distributions using a causal modeling approach that incorporated direct and in-  (Devarajan et al., 2020) will likely provide similar outputs to the causal model given that they all incorporate interactions explicitly in the model, although we note that given their relative infancy, have not yet proliferated into the standard workflow for projecting species distributions. It is also noteworthy, however, that these multi-species occupancy approaches generally require large sample sizes when compared to the causal models we apply here and are subject to multicollinearity issues when abiotic predictors are correlated with competitors (Godsoe et al., 2017), which limit the ability to separate out biotic and abiotic effects for the subordinate species. Finally, it is important to note that we cannot be certain that competitive interactions were present. Ideally, future research should couple experimental work and simulations with demographic and dietary studies to better understand the influence of climate and biotic interactions on species ranges.
Regardless of these uncertainties, our study represents the first, hopefully of many, spatially explicit species distribution model that uses a causal framework. We were able to show how climate potentially mediates interactions between similar species. These findings can be used to guide management and conservation decisions for Canada lynx, American martens, and other winter-adapted species along trailing edges that are predicted to decline due to climate change . Such actions could include land-use planning to protect or promote winter climate (Dickerson-Lange et al., 2017;Morelli et al., 2016) or predator management. However, we encourage others to use our approach to determine if indirect and direct factors differ in other regions of these species ranges (e.g., Hostetter et al., 2020;King et al., 2020;Zielinski et al., 2017) or to perform a range-wide examination using existing data from multiple camera networks. Most importantly, our study provides empirical support for recent studies that have advocated for using a mechanistic framework to forecast future distributions in light of climate and land-use change (Pacifici et al., 2015;Urban et al., 2016).

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interests.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13480.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data and R code used in this study are available at Dryad : https://doi.org/10.5061/dryad.k0p2n gf9j.