What drives spatially varying ecological relationships in a wide‐ranging species?

Decades of research on species distributions has revealed geographic variation in species‐environment relationships for a given species. That is, the way a species uses the local environment varies across geographic space. However, the drivers underlying this variation are contested and still largely unexplored. Niche traits that are conserved should reflect the evolutionary history of a species whereas more flexible ecological traits could vary at finer scales, reflecting local adaptation.


| INTRODUC TI ON
Understanding the ecological determinants of species distributions is a core area of wildlife research. Though studying species distributions initially provided basic ecological knowledge about a species' relationship with environmental factors, advancements in the field have made results useful for conservation planning, resource management, impact and risk assessments, and evaluating the implications of climate change on global biodiversity (Austin & Van Niel, 2011;Bled et al., 2011;Loiselle et al., 2003). The increasing number of applications of studying species distributions has in part been facilitated by larger, open-access species occurrence databases [e.g., iNaturalist ("iNaturalist.org", 2020), Global Biodiversity Information Facility ("GBIF.org", 2020)] coupled with the increasing availability of higher-resolution environmental predictors (e.g., remote sensing products) (Kays et al., 2020). This new era of studying species distributions has the potential to provide nuanced insight into how species use the environment across their entire geographic range (Maldonado et al., 2015).
Historically, studies of range-wide species distributions have assumed constant spatial relationships with the environment -that a species responds to a given environmental factor in the same way across their range (Koenig, 2001;Miller, 2012;Osborne et al., 2007).
However, a growing body of literature has challenged that assertion (e.g., Smith, Godsoe, et al., 2019;Sultaire et al., 2022;Ye et al., 2020), although the causes of within-species variation are not known and are often contested. For example, McNew et al. (2013) reported that greater prairie-chicken (Tympanuchus cupido) females varied spatially in their habitat selection, presumably due to individuals adapting to varying microsite conditions and quality across their three study regions; similar variable habitat selection has been documented in other taxa as well (Morellet et al., 2011).
The degree of ecological specialization within a species range reflects the extent to which its niche dimensions are evolutionary conserved vs locally adapted (Peterson, 2011). Most commonly, applications of species distribution models presume full niche conservatism by analysing all populations of a species together (i.e. 'global model'), and thus presume all animals will retain ecological characteristics of their specie's fundamental niche (Peterson, 2011;Peterson et al., 1999;Wiens et al., 2010). Some studies have considered within-species differences by modelling evolutionary lineages separately (e.g., from separate Pleistocene refugia, Myers et al., 2020).
However, most have not explicitly compared the performance of the lineage model to the global model (Chardon et al., 2020) and have not considered lineages of different ages. Understanding the extent to which the evolutionary history of the populations of a species predicts ecological relationships would be useful not only for creating more accurate SDMs but also for helping the populations persist through the rapidly changing climate and land use of the 21st century.
Ecological specialization could also occur at within or across genetic lineages through fine-scale local adaptation or differences in ecological context (Sexton et al., 2009, Smith, Beever, et al., 2019, Pease et al., 2022. In this case, the ecological context refers to environmental characteristics that might drive local adaptation including climate gradients, vegetation communities, human disturbance regimes, and competing species. Although evolutionary lineages of a species (i.e. subspecies) often map on major ecological zones, this is not necessarily the case for all aspects of ecological context (Waltari et al., 2007). Thus, evolutionary history and ecological context can offer independent predictions for patterns of ecological specialization of a species. For species where ecological context has a greater impact on intraspecific variation we might expect ecological boundaries such as ecoregion delineations or differing, competitive ecological communities may better describe how and where a given species is using the environment differently (Omernik, 1987).
Here we explore how American black bear (Ursus americanus; hereafter, black bear) varies in its ecological relationships across North America. Black bears are an ideal species for exploring spacevarying effects for several reasons. First, the species is a generalist with a large geographic range, suggesting a high potential for local adaptation and plasticity (Garshelis et al., 2016). Second, the genetic lineage boundaries of black bear are well characterized, allowing us to test two degrees of evolutionary history (Pleistocene refugia and subspecies) to explain patterns of non-constant ecological relationships (Pelton, 2003;Puckett et al., 2015). Finally, across their large range, black bears occupy many distinct ecoregions, coinciding with varying numbers of competing large carnivores and persist in many degrees of human disturbance, providing a means for exploring ecological context as a driver of local adaptation (Boitani et al., 2018;Omernik & Griffith, 2014). Although this species has been the subject of hundreds of local-level ecological studies (e.g., Gould, 2020), black bear ecology has not been synergized at a range-wide level.
Thus, the goal of our paper is to establish whether black bearenvironment relationships are constant across their range, and if not, which niche dimensions vary and rather the patterns of this variation are better explained by phylogenetic niche conservatism or ecological context.
To evaluate the variability of black bear ecology across their range, we first tested for evidence of spatial non-stationarity -the tendency for changes in the direction and/or magnitude of relationships of bear distribution to land cover, human densities, and climate. We evaluated non-stationarity by creating hundreds of local spatial models and comparing these findings to those from a range-wide model of black bear observations (Fink et al., 2010;Lloyd, 2010). While alternative approaches to K E Y W O R D S American black bear, niche conservatism, spatial non-stationarity, species-environment relationships, Ursus americanus assessing spatial non-stationarity exist (e.g., Spatially Varying Coefficients, Geographic Weighted Regression), we chose this approach as to ease the computational burden of the geographic extent and volume of data in this study (Fotheringham et al., 2003, Gelfand et al., 2003, Pease et al., 2022. We then used the results of the local spatial models to test four hypotheses about drivers of ecological non-stationary in black bear across their range by comparing the spatial pattern of non-stationarity to the boundaries of: (1) Pleistocene lineages, (2) subspecies, (3) ecoregions, and (4) patterns of competing species overlap (Figure 1). To test these questions and hypotheses, we built models using a publicly available database of opportunistic observations of black bear across North America from the iNaturalist platform and test model results with an independent national-scale camera trap dataset, Snapshot USA (Cove et al., 2021).

| Animal observations
We queried the iNaturalist database for observations of black bear across North America during 2015-2020. iNaturalist is an online social network for sharing biodiversity information to help others learn about nature ("iNaturalist.org", 2020). Observations are submitted to the database by iNaturalist users and the iNaturalist community identifies the observations to the highest taxonomic resolution possible (e.g., genus, species). iNaturalist users can submit observations across all taxa, anywhere on earth, with no required sampling protocol; there are, however, some coordinated efforts at local scales (e.g., "BioBlitz"). We used iNaturalist because no structured study on black bears matched the accessibility, quantity, and spatial extent of iNaturalist.
We queried the iNaturalist Application Programming Interface (API) for all research-grade observations (i.e., observations that include a photograph and have confirmed and agreed identification by more than 2 users) of the class Mammalia within North America (specific bounding box: −52 degrees northeast longitude, −170 degrees southwest longitude, 6 degrees southwest latitude, 75.5 degrees northeast latitude). We vetted 100 random black bear observations to check for species identification accuracy (100% accurate). We removed observations that were classified as scat, tracks, or signs, or associated with the North American Tracking Database Project on iNaturalist because they could be less distinguishable than photographs of the animals themselves (Lonsinger et al., 2015). We also filtered observations that had coordinates obscured (uncertainty >20 km) due to user preference or iNaturalist's protocol of obscuring the location of threatened or endangered F I G U R E 1 Visualization of four hypotheses driving non-stationarity in American black bear (Ursus americanus) ecological relationships across the species' current range, with an example of how a given spatial support window is assigned to a single subdivision of each hypothesis. The rectangle on the map shows one example spatial window, with the centroid of the spatial window shown as a black dot, while the grey points represent the centroid of all windows (n = 1373). The bear locations within the windows were used to build local models across the species' range. Each local model estimates ecological relationships based on observations within the window and was assigned one subdivision of each hypothesis based on its centroid. The example here shows a box centroid in the panhandle of Florida, USA, and is assigned to the: (1) "Eastern Temperate Forest" subdivision of the North America ecoregions; (2) "eastern" black bear lineage; (3) "floridanus" subspecies; and (4) no overlap with competing predators [brown bear (Ursus arctos, green), grey wolf (Canis lupus, brown), red wolf (Canis rufus), and cougar (Puma concolor)]. Colours for the competition subfigure represent the distribution maps for each of the top competitors species. From this North American mammal database, we then selected all observations of Ursus americanus that only occurred during our target time frame (01-Jan-2015 -31-Dec-2019). By using a fuzzy search function (Wickham, 2019), our filtering also included observations that were assigned a subspecies name (e.g., "Ursus americanus altifrontalis") but we did not differentiate among subspecies in any analysis. All remaining Mammalia observations were retained and used in addition to black bear observations as a proxy for iNaturalist search effort.

| Opportunistic reporting
Opportunistic observations of biodiversity, such as those reported in iNaturalist, may provide valuable insights into ecological relationships but the data typically prove challenging for unbiased statistical inference (Bird et al., 2014;Clare et al., 2019;Isaac et al., 2014;Yackulic et al., 2013). We recognize that these observations are often clustered and positively associated with human population density and the lack of standardized sampling protocols presents complications when evaluating survey search effort (Bird et al., 2014). Although human population density and accessible public lands such as National Parks are typically good proxies to account for the inherent spatial sampling biases in opportunistic observations (Geldmann et al., 2016), we chose to the use the iNaturalist database itself as a direct measure of iNaturalist user behaviour to evaluate the tendency for a species to be reported. We used the total number of iNaturalist mammal observations in each grid cell to reflect the relative search effort in that geographic space (Table 1; Phillips et al., 2009, Kéry et al., 2010. By using this measure of effort, we assume that the probability of a person reporting a black bear was related to the general probability of reporting any mammal species in the area. Some users in the iNaturalist database clearly violate this assumption if they tend to report observations of a single species only, so we removed their observations from our analysis. This type of behaviour may arise when groups (e.g., state agencies, non-profits) use iNaturalist as a platform for documenting observations associated with a given campaign such as the "Arkansas Bear Survey" promoted by the Arkansas Game and Fish Commission. We discovered these targeted efforts by calculating, for each iNaturalist user, the proportion of black bear observations out of their total number of mammal observations, and then removing data from the top 5th percentile.
This removed 2298 observations from 150 users in our dataset (~15 black bear observations/user).

| Study area
The spatial extent of black bear iNaturalist observations spanned the entire North American continent ( Figure S1). Rather than restricting to a specific country, we used data from across black bear's entire current range (IUCN; Garshelis et al., 2016); although only currently occupying about 65% of their historic range, this still spans Canada, Mexico, and the United States though sometimes in isolated pockets (Garshelis et al. (2016); Figure S1).

| Ecological covariates
We gathered environmental factors known to influence black bear distribution and abundance (Clark et al., 2020;Feldhamer et al., 2003;Pelton, 2003;Penteriani & Melletti, 2020), including land cover, climate, and human densities (Table 1). Ecological covariates were summarized within 50 sq km hexagonal grid cells over an area of 2.1 million sq km, resulting in 502, 961 grid cells ( Figure S1). We chose 50 sq km as a scale to reflect a typical home range of black bears (Pelton, 2003).
Black bear distribution is typically linked to areas of contiguous forest cover, but the forest types used are dependent on the geographic location (Pelton, 2003). To capture this variation in forest types used, we used the Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type (MCD12Q1) product from 2019 to differentiate forest types, where we included the proportion of evergreen needleleaf (i.e., conifer forests), deciduous broadleaf, and mixed forest types (Table 1). We also used the MODIS product to describe the proportion of cropland available as this is a dominant land cover type across much of the United States and black bear opportunistically forage in crop fields (Ditmer et al., 2018;Jones & Pelton, 2003).
We used the Gridded Population of the World (GPW 2015; Version 4) to calculate mean human population density per grid cell (Center for International Earth Science Information Network -CIESIN -Columbia University, 2018). Although other products exist (e.g., human footprint index), we chose the GPW product because we wanted to isolate the effect of humans and separately account for factors often incorporated into indices such as the human footprint index (e.g., cropland). Protected areas can be an important designation for a large, wide-ranging species, so we accounted for whether a given grid cell contained protected areas using information from the World Database of Protected Areas (UNEP-WCMC and IUCN 2020). Here, "protected areas" includes a range of designations including national forests, national parks, wildlife management areas, and wildlife sanctuaries with varying regulatory practices (e.g., hunting vs no hunting); these land uses typically fall into the category of "public lands" with varying levels of use and protection.
For climate variables, we used the Bioclimatic variables of annual average temperature and annual average precipitation from the WorldClim database (Fick & Hijmans, 2017). We selected these climatic variables as they reflected variation in climate conditions across the continent but also are correlated with productivity (e.g., gross primary productivity), which may be related to the wildlife carrying capacity of an area that may dictate species abundance and richness (Hobi et al., 2017). We considered including additional climatic variables to describe nuanced differences across space but found other WorldClim variables to be highly correlated with the existing set of predictors; this was also true for topographic characteristics such as elevation.

| Model description
There is a long history of fitting species distribution models using presence-only observations, which has become increasingly common with widespread availability of biodiversity occurrence records (Elith et al., 2011;Guillera-Arroita et al., 2015;Kays et al., 2020). Typically, presence-only observations are modelled by generating so-called "background points" -random locations within the area of interest where a presence-only observation did not occur -and then using logistic regression (1 for locations with observations, 0 otherwise) to make inference about ecological relationships (Renner et al., 2015;Warton et al., 2013). The described approach is typically used when 1) modelling is done in continuous space and 2) no relevant search effort information is available. This approach is conceptually equivalent to a resource selection function under a used-available framework, where black bear iNaturalist observations are considered "used" locations and the unused grid cells become the "available" background environmental conditions (Matthiopoulos et al., 2020).
In our case, we had discrete spatial units (i.e., hexagonal grids) and information on search effort (i.e., the number of iNaturalist observations in a grid cell). Because we had discrete spatial units, we were able to quantify the number of black bear observations within the a given grid cell -either 0 or more; therefore, there was no need to generate random background locations as we could simply assign a "0" to grid cells with no black bear observations. Additionally, the total number of iNaturalist observations recorded in each grid cell indicated whether the "0" assigned to a grid cell was a result of no black bear, no search effort, or both.
By using the iNaturalist database itself, we avoided the presencebackground modelling approach of opportunistic records and considered black bear observations as discrete count data (counts of black bear in a grid cell) and modelled using a generalized linear model with Poisson distributed errors. In this case, the predicted response (the Poisson intensity at grid cell s, λ s ) represented the expected number of black bear reports per-unit area (where the area is the size of the hexagon grid cells, 50 sq km). We fit our models with a log link using the glm function in R version 4.0.3 (R Core Team, 2020). The linear predictor of the model included all eight ecological covariates listed in Table 1. We included the search effort term (i.e., the number of iNaturalist mammal observations in a grid cell) in the linear predictor of the model instead of as an offset of the Poisson model so that a regression term is estimated instead of fixed at 1 (as is the case when used as an

| Non-stationarity in black bear-environment relationships
To explore non-stationarity in black bear ecological relationships we used an ensemble modelling procedure that used local spatial models within "spatial support windows" to parameterize local distribution models. A "spatial support window" is a user-defined geographically restricted sampling area (i.e., a square smaller than the total area of interest) where only observations and background points within each window are used to estimate regression coefficients. Estimates from overlapping spatial support windows were then averaged to effectively "scale up" and create an ensemble prediction across the entire landscape (Fink et al., 2010). We required a minimum number of 100 black bear observations in a given spatial support window for that window to be included in the ensemble model. This resulted in 1373 spatial windows that were approximately 10.5 degree longitude by 7.5 degree latitude in extent with at least 100 black bear observations (range = 100-1617; mean = 512) ( Figure 2) (Wisz et al., 2008). We determined the spatial window sizes to balance the geographic coverage of current black bear range with iNaturalist observation density; the spatial windows can be smaller when dense observations occur consistently across the study area, and window size increases with fewer observations. We then fit a model identical to the global model, with the same linear predictor that used the environmental covariates described in Table 1, but the black bear observations and background points were those only found in each window. That is, the ensemble and global models were identical in their specification (i.e., same use ecological predictors and response), but differed in their extent.
We compared absolute value of residual errors to quantify performance between the global and ensemble models (Rollinson et al., 2021). We

| Stationarity index
Following parameter estimation of the ensemble models, we calculated a stationary index (SI) for each ecological covariate to evaluate whether evidence of non-stationarity existed and, if so, the magnitude of the variation. In the ensemble modelling approach, each ecological covariate is estimated within each spatial window (n = 1373), so an interquartile range of the 1373 parameter estimates can be calculated; this represents the range or magnitude of variation in the parameter estimates across space (Osborne et al., 2007). The SI is calculated by dividing the interquartile range of a given ecological covariate by two times the standard error of the same ecological parameter estimated in the global model. Index values greater than one indicate evidence of non-stationarity, with higher values suggesting more spatial variation in a given ecological relationship (Osborne et al., 2007). Overall, the stationarity index is to indicate evidence of spatial non-stationarity and serve as a measure of strength of the non-stationarity.
We supplemented the stationarity index with a formal test of non-stationarity of model coefficients using a Moran test (Osborne et al., 2007). Given that we were using spatial support window centroids for our analysis, we defined neighbours for the Moran test based on distances among centroids (rather than defining neighbours on a grid or neighbouring polygons due to spatial support windows overlapping Finally, it is possible that a signal of spatial non-stationarity can arise from missing covariates or interactions in the linear predictor as well as the spatial pattern of observations. To explore the former, we included interactive effects between protected status and forest cover types and assessed how stationary indices changed with the inclusion of interactive effects. To test for the latter, we conducted a sensitivity test of the global model through bootstrapping.
Specifically, we subset the full dataset 1200 times to contain 100 black bear observations and a comparable number of grid cells as the local models, fit "global" models as described above, and then qualitatively compared the variation in the parameter estimates from this approach compared to those estimated in the local models.

| Evaluating drivers of non-stationarity
For any ecological covariate with evidence of non-stationarity (based on SI values >1) we separately tested four hypotheses (ecoregions, lineages, subspecies, or competing species) to explain the spatial pattern of non-stationarity (Figure 1). We did this for each covariate to allow that the drivers of non-stationarity may be different across land use, anthropogenic effects, and climatic responses. For the ecoregion hypothesis, we used the EPA Level 1 ecoregion delineations which are based on similar vegetation, climate, hydrology, land use, geology, landforms, and soils (Commission for Environmental Cooperation, 1997;Omernik, 1987;Omernik & Griffith, 2014). Pleistocene lineages were delineated following phylogeographic analyses by Puckett et al. (2015) and subspecies, representing evolutionary linages forming since the Pleistocene, described by Larivière (2001). The competition hypothesis describes whether a given spatial support centroid is within the current range of large North American predators [brown bears (Ursus arctos), grey wolf (Canis lupus), red wolf (Canis rufus), or cougar (Puma concolor)], which were determined using each species' IUCN range map (Boitani et al., 2018;McLellan et al., 2017); overlap with each species resulted in a nominal predictor of 1-4 (depending on the number of overlapped ranges). We included these four species because they are typically competitive with black bear, either dominantly (e.g., brown bear) or indirectly through kleptoparasitism (e.g., cougar) (Apps et al., 2006;Ballard et al., 2003;Latham & Boutin, 2011;McLellan, 2011).
We first generated centroids of each spatial support window and then assigned each centroid to the hypothesis subdivisions with each centroid, such that each centroid was assigned to one subdivision of each hypothesis (Figure 1). The hypotheses were then tested by fitting linear models with the response (i.e., dependent variable) as the ensemble estimated coefficient (n = 1373 coefficient estimates for each covariate, one from each spatial window) and the F I G U R E 2 Local ensemble models were developed using restricted spatial support windows across observations of American black bear (Ursus americanus) reported on iNaturalist during 2015-2019 (orange dots). Spatial support windows (black-outlined boxes) (total = 1373; n = 50 presented here for illustration) serve as bounds for observations and background points in the ensemble modelling approach and were included only if they included >100 bear observations. Results from overlapping boxes are averaged to effectively "scale up" to the range of black bear subdivision of each hypothesis as a categorical predictor variable.
Thus, for each covariate, we fit four linear models with different groupings of observations and determined hypothesis support using Akaike's Information Criterion (AIC) (Akaike, 1998).

| Assumptions
There are many assumptions inherent to wildlife sampling that are typically not met in opportunistic, presence-only databases (Bird et al., 2014;Hochachka et al., 2012;Yackulic et al., 2013). One key assumption is that sampling is either random or representative of the environmental conditions throughout the area of interest (Yackulic et al., 2013). We explored whether the spatial distribution of iNaturalist mammal observations was representative of land cover conditions across North America using the MODIS land cover product.
Using a used-available framework (akin to resource selection functions), we created a categorical classification indicating whether an iNaturalist observation occurred in a grid cell; cell values were 1 if an iNaturalist observation occurred and 0 otherwise. We then computed the most common land cover type within each hexagonal grid cell and assigned that land cover class to the hexagonal grid cell. Using this information, we compared use vs available across all 16 land cover classes described by the MODIS product using a chisquare test of homogeneity, testing the null hypothesis that different populations (used vs available) have the same proportions of land cover classifications (Thomas & Taylor, 1990).

| Verification
An important step in modelling animal occurrence records is determining how well a model describes where a species occurs. That is, predictive performance of such models is often an essential summary when evaluating model diagnostics (Elith et al., 2006). Typically, this is done by withholding a set of occurrence records, fitting a model without those observations, and then evaluating the model predictions to those locations ("testing points"). However, for the most reliable validation, the testing points should be independent of the observations used to train the model ("training points"); this standard is rarely met because an additional, independent data source is usually not available (Hijmans, 2012).
Here, however, we took advantage of a recent, large-scale animal monitoring effort -Snapshot USA -as an independent data source to verify parameter estimates from our modelling efforts (Cove et al., 2021 (Laboratory, 2015).

| iNaturalist observations
The iNaturalist query resulted in 676,007 terrestrial Mammalia observations during 01-January-2015 -31-October-2020, with 14,099 observations of black bear. We removed all iNaturalist observations that were flagged as animal tracks or scat, had an accuracy uncertainty >20,000 m, were missing coordinates, or were the result of a presumed targeted search effort. We additionally buffered the current IUCN black bear range map by 100 km and removed observations outside of this area (n = 5 observations). We chose be consistently found outside of their current IUCN range, based on home range, dispersal distances, and potential range shifts (Karelus et al., 2016;Koehler & Pierce, 2003;Schwartz & Franzmann, 1992 Overall, iNaturalist mammal observations were representative of the available land cover types across North America ( 2 14 = 0.101, p-value = 1; Figure S2).

| Range-wide iNaturalist model
Results from the global, range-wide model of black bear iNaturalist observations indicated that public land status had the strongest positive effect on the number of black bear observations relative to all other mammals across their range (Figure 3; Table S1). The global model indicated greater black bear occurrence with greater forest cover, regardless of forest cover type (conifer, deciduous, and mixed) but bear occurrence was lower where human population density and cropland cover were greater (Figure 3; Table S1). At the continental scale, black bear observations were more common in areas with greater precipitation, but there was no apparent effect of temperature ( Figure 3; Table S1).
Global regression coefficients from iNaturalist agreed directionally with Snapshot USA, which explicitly accounted for detection probability, in 88 percent of the covariates (Figure 3; Table S1). The greatest discrepancy between the iNaturalist and Snapshot USA model occurred in the estimated effect of mean annual temperature, where the iNaturalist model estimated no effect (i.e., point estimate = −0.01 with 95% CIs overlapping zero) but the Snapshot USA modelled estimated a negative effect, with non-overlapping 95% confidence intervals for the latter (Figure 3; Table S1). Despite this discrepancy, there was a prediction accuracy of 73% when using the iNaturalist global model to predict the Snapshot USA black bear occurrences (AUC = 0.73).

| Local ensemble model
We fit 1373 regression models geographically restricted to observations and environmental conditions within each spatial support window ( Figure 4). Through the local ensemble model, the iNaturalist database documented areas of black bear occurrence that range beyond the current IUCN range boundaries (Figure 4). Although all models converged, 200 had at least one standard error estimate greater than 10, which we deemed as erroneous and discarded from the ensemble predictions, bringing the total number of spatial F I G U R E 3 A comparison of regression coefficient estimated using a global Poisson regression model fit using American black bear (Ursus americanus) observations from the iNaturalist database during 2015-2020 (blue dots and lines) and a single-species, singleseason occupancy model fit using camera observations from the 2019 snapshot USA monitoring initiative (orange dots and lines). Dots represent point estimates and bars indicate 95% confidence intervals support windows to 1173. The percentage of Variance Inflation Factors exceeding our threshold of five had a maximum value of three percent across all covariates (temperature and precipitation), but six of the eight covariates had Variance Inflation Factors that exceeded five in less than 1 percent of the local models (Table S2). The local ensemble approach reduced residual sum of square errors by F I G U R E 4 Predicted intensity from an ensemble modelling approach estimated using 1173 geographically restricted spatial support windows and observations of American black bear (Ursus americanus) reported on iNaturalist during 2015-2019. Dotted outlines represent current known black bear species range. Areas of white colouring represent regions beyond the range of black bear or regions of the species range where windows had less than 100 black bear observations and were omitted from analysis F I G U R E 5 Regression coefficients estimates (violins; n = 1173) from the local Poisson regression models fit to iNaturalist black bear (Ursus americanus) observations during 2015-20 across North America. Violins represent the coefficient estimates across all spatial support windows. Also included for comparison is the 95% confidence intervals (dot and bars) from the global model. Solid horizontal black line is a marker for zero, and the dashed horizontal line is the median value of the coefficient estimates for the local models. All covariates were scaled and centered prior to analysis. iNaturalist records represents the regression coefficient for the search effort term in the model 11.3 percent compared to the global regression model. Additionally, the local ensemble modelling approach indicated strong support for spatially varying ecological relationships ( Figure 5). Each of the eight ecological covariates had stationary indices greater than one (range 9-64), with responses to climatic conditions (i.e., mean annual temperature and precipitation) being the most variable across black bear's range ( Figure 5). Additionally, we documented little change in the stationary indices when interactive effects were included in the linear predictor of the local models (Table S3). Finally, all but one (temperature) covariate failed to reject the null hypothesis of random dispersion ( Figure S3). The range at which the autocorrelation was no longer present varied across the covariates but ranged from 1800-3600 km ( Figure S3).
The ensemble model results showed notable variation in the spatial patterns of the eight coefficient estimates (Figure 6). The greatest variation in response to climatic conditions -the two ecological covariates with the greatest SI -tended to occur along the Midwest divide in the species' range, particularly in and around Arkansas, Oklahoma, Texas, and Louisiana ( Figure 6). The sensitivity test we conducted on the global model exhibited the same median directional effects as the local models, with the individual bootstrapped estimates spanning negative and positive effects -similar to our results from the local models -albeit less variation than the local models. The spatial pattern hypothesis comparison indicated the greatest support for the subspecies hypothesis in five of the eight ecological covariates; responses to climatic conditions (i.e., mean annual temperature and precipitation) were best explained by the ecoregion hypothesis and responses to human population density was best explained by the competition hypothesis (Table 2; Figure 7; Table S4).

| DISCUSS ION
Despite the mounting evidence that species-environment relationships can vary across space, the key drivers of the variation and the characteristics of the species exhibiting non-stationarity are relatively unknown (Miller, 2012, Smith, Beever, et al., 2019, Pease et al., 2022. Understanding these patterns is important because it may provide nuanced insight into a species' ecology, directly informing how regional management and conservation approaches could differ, and potentially informing how a species could respond to future environmental change. Our results show that black bear ecology varies across their range in their relationships with all types of covariates tested: land cover, human density, and climate. Black bear subspecies boundaries best explained differences in bear responses to five land cover variables, ecoregion boundaries best described variation in climatic conditions, and responses of bears to human population density was most related to the number of co-occurring competitors. We found no support for niche conservatism at the level of Pleistocene lineages. Our results demonstrate the complexity of spatial variation in ecological relationships with subspecieslevel local adaptation to landcover and finer-scale adaptation to climate. We additionally document the interplay of human population densities and the number of top competitors (e.g., wolves, grizzlies), and how that affects the way black bears use the landscape.
Future applications of large-scale SDMs (e.g., regional or continental) should work to account for varied relationships across space using biologically relevant delineations across the area of interest.
The ecological relationships of a local population reflect a balance between the population's potential to adapt to local conditions vs. the phylogenetic constraints associated with the population's genetic legacy. The ecology of black bears is variable across their range in terms of diet (Bull et al., 2001;Graber & White, 1983), hibernation (Gould, 2020;Johnson et al., 2018), and, as we show here, in their ecological relationships with landcover, humans, and climate. The fact that subspecies boundaries, and not Pleistocene lineages, were the best explanation for intraspecific variation in the response to land cover helps us understand the temporal scale of this local adaptation. The 16 subspecies represent evolutionary change after black bears expanded out of four Pleistocene refugia (2.6-0.012 Ma) into their current range. We propose that ecological adaptations during this time led to the differences we see today in the suitability of different forest types for different bear populations. The fact that the spatial pattern of these relationships did not better align with more fine-scale delineations of habitat type (ecoregion) suggests there is some genetic niche conservatism within these factors, and most are similar morphologically, with minor differences in cranial shape (Larivière, 2001;Puckett et al., 2015).
The two climatic variables we considered, mean annual temperature and precipitation, had the strongest signal of non-stationarity and the most localized spatial patterns, best explained by ecoregions. This is strong support that black bears are locally adapting to climatic conditions, regardless of their evolutionary history. The strongest relationships with these climate variables occurred on the boundaries of black bear's range (e.g., Missouri, southern Minnesota), and the magnitude of the effects suggested strong avoidance and attraction to climatic conditions, depending on the geographic region. For example, black bears were relatively more common in cooler conditions in Arkansas, Oklahoma, and east Texas, but at warmer conditions at the northern boundary of the species' range. Black bear's extensive geographical range and the inherent broad climatic niche associated with this range, coupled with evidence of local adaptation to climatic conditions, suggests a lack of vulnerability to forecasted climatic variation in much of their range.
The final variable we tested, human population density, was the only variable best explained by the competition hypothesis, where the varied effect of human population density was explained by the number of competing, sympatric predators in an area. This finding may be indicative of the situation black bears are forced into in modern day -one of simultaneous avoidance of competitors and, in many places, humans. In areas of high natural food availability, black bear tend to avoid humans and human-altered landscapes (Beckmann & Berger, 2003). In areas of range overlap with competing, sympatric grizzly bear, black bear's avoidance of grizzlies has been documented to result in increased temporal and spatial overlap | 1763 PEASE Et Al.
with human activity (Apps et al., 2006;Belant et al., 2010;Ladle et al., 2018). In contrast, in areas where black bear overlap with cougars, human avoidance may be heightened due to supplemental food availability from cached carcasses by cougars cached carcasses by cougars; black bear are known to exhibit kleptoparasitism and heavily scavenge on these caches (Elbroch & Kusler, 2018;Engebretsen et al., 2021). There are of course areas where black bear do not avoid increasing human population densities despite a lack of sympatric competitors (e.g., Urban areas in Appalachia; Gould, 2020) and additional factors contributing to differing black bear-environment relationships in these areas are likely at play (e.g., protected area status; hunting status of an area). We suggest that our finding that bears respond to human densities and that this variation is best explained by the existence of sympatric competitors is the result of a combination F I G U R E 6 Estimated regression coefficients from the local ensemble models (n = 1173) using American black bear (Ursus americanus) iNaturalist observations during 2015-2020 across North America. Each subfigure shows the spatial variation of each ecological relationship where the legend represents the range of the estimated coefficients for each covariate. Stationary index (SI) values greater than one indicate support for non-stationarity of niche conservatism (i.e., genetic predisposition) in how different populations respond to potential predators (such as humans or grizzlies), and local differences in bear management plans and human behaviour.
We also found an interesting spatial pattern in the relationship with public lands status. In most of their range, black bears are typically more common inside protected areas, particularly in the presence of high human disturbance, but in the Ouachita Mountain TA B L E 2 Model comparisons testing which non-stationarity hypothesis best explains variation in each ecological covariate's spatial pattern

| Testing assumptions of iNaturalist data
The use of iNaturalist observations to study wildlife ecology is in its infancy, and there is much debate over the validity of using such opportunistic, presence-only, citizen science observations to describe ecological relationships (Bird et al., 2014;Hochachka et al., 2012).
Given that we are one of the few to use these data to study North Using the presence-only observations of iNaturalist, we were unable to account for an important source of sampling bias -imperfect detection of black bear (Kéry & Schmidt, 2008). Imperfect detection of a species -that is, the tendency to not detect a species given it is present -can lead to biased parameter estimates and potentially conflate ecological effects with the observation process.
For black bear, we believe this source of bias is reduced (albeit not eliminated) because its large body size and diurnal habits make it conspicuous. We attempted to account for reporting probability through the use of sampling effort and occurrence probability using a suite of ecological covariates known to affect the species, thus we believe our ecological findings are reflective of the species' ecology.
Further, we were able to explicitly account for imperfect detection in the Snapshot USA through hierarchical modelling, and parameter estimates from those modelling efforts mirrored the ecological findings of our iNaturalist model, suggesting that the bias of imperfect detection may be minimal for iNaturalist black bear observations taken at the 50-km 2 scale. Additionally, we do not have any way of detecting non-independence of iNaturalist observations (e.g., pictures/observations of the same bear). We believe this may be accounted for by modelling of observation density -which is different from actual black bear density -but it is possible that this non-independence may be influencing covariate relationships.

| CON CLUS ION
Our research contributes to the increasing evidence that speciesenvironment relationships are not constant over space, and whether niche conservatism or local adaptation to ecological context explains the variation depends on the environmental conditions being evaluated (e.g., land use vs climatic) (Smith, Beever, et al., 2019, Pease et al., 2022. Although we documented support for subspecies delineations, local adaptation may be occurring at a finer scale than tested here and future research should look to fine-scale drivers such as human pressure or species-specific management regimes as drivers of local adaptation. Future work should also explore this question across taxa over large geographic space.
The strong ties we found to subspecies and ecoregions suggest that ecological relationships are dynamic over a relatively fine spatial scale. Our findings have implications for studying how species-environment relationships may change over time, and we encourage the acknowledgement of such variation in correlative assessments of species distribution, abundance, and richness; local spatial modelling and ensemble modelling are promising ways to explore this variation. We also believe that persistent research on the validity and use of opportunistic, presence-only databases such as iNaturalist will continue to lead to interesting, nuanced ecological findings not apparent in studies conducted over restricted geographic areas, and that the role of these databases in data fusion is likely to become increasingly important (Miller et al., 2019;Pacifici et al., 2017).
Looking forward, we believe ignoring spatial and/or temporal variation in large-scale modelling exercises is a missed opportunity to understand how species are changing their local behaviours with the changing world. Ignoring this variation may not only result in misleading inference, but could also result in overlooked species adaptations right before our eyes (Thompson, 1998). This is not to say that past range-wide modelling efforts are invalid, but rather to highlight the exciting path forward in understanding species distributions as we continue to accumulate billions of species occurrence vouchers with up-to-date, high-resolution environmental information.

ACK N OWLED G EM ENTS
We thank the iNaturalist community for submitting and reviewing thousands of American Black Bear observations. Thanks also to the Snapshot USA network for operating cameras across the United States. We are grateful for helpful comments on earlier drafts from Christopher Moorman. Three anonymous reviewers provided very helpful comments that substantially improved this work.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data used in this paper are downloaded from publicly accessible websites cited in the main text (https://emamm al.si.edu/; https:// inatu ralist.org).

PE E R 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.13594.