Spatiotemporal shifts in distribution of a recolonizing black bear population

. Following one of the most successful relocations of a large carnivore species, American black bears ( Ursus americanus ; hereafter bears) have expanded from Arkansas to the Ozark Highlands of southern Missouri, USA, where they were potentially extirpated in the early 1900s. Our ﬁ rst objective was to estimate spatial and temporal shifts in probability of bear presence to understand the population distribution in Missouri. Our second objective was to assess which factors might in ﬂ uence any detected shifts. We used public sightings of bears to develop statewide spatial distribution models for early (1991 – 2003) and late (2004 – 2015) bear presence probability in Missouri using maximum entropy modeling and determine the rate and direction of population expansion. We evaluated the effects of environmental variables (e.g., land cover, roads) on bear presence probability and predicted the total area of three levels of bear presence probabilities. The average annual density of bear sightings increased temporally by about 3.4/1000 km 2 in Missouri. The Ozark Highlands ecoregion contained the highest average annual density of sightings during both time periods. The bear presence probability increased spatially about 29 km in a north – northwest direction between the two time periods and the distribution became more dispersed over time. Our results corroborate prior ﬁ ndings suggesting bear presence is greatest in forests and rural areas. However, the increased importance of cultivated crops over time combined with the consistently positive relationship with roads indicates little negative effects from anthropogenic features on the bear presence probability in Missouri.


INTRODUCTION
Many large carnivore species are expanding and recolonizing historical range distributions globally in response to conservation actions and improved human attitudes (Gompper et al. 2015). Though variable public perceptions of certain carnivore species, such as grizzly bears (Ursus arctos) and cougars (Puma concolor) in North America, persist in some areas, cultural acceptance of positive human-carnivore coexistence is increasing, particularly in areas with favorable wildlife management policies (Kellert et al. 1996, Linnell et al. 2001, LaRue et al. 2012. Conservation measures such as implementation of appropriate harvest regulations or protections afforded by threatened or endangered species designations reduced black bear mortality rates and supported the recovery and increase of black bear (Ursus americanus) populations by the late 20th century in many portions of their historical range within the United States (Hristienko andMcDonald 2007, Beckmann andLackey 2008).
Range expansion and recovery rates of carnivore species can be spatially and temporally inhibited by environmental variables (Gaston 2003). Very large rivers, mountain ranges, and major roads comprise ecological and anthropogenic barriers that result in habitat fragmentation that can inhibit range expansion (White et al. 2000, Riley et al. 2006, Long et al. 2013, Simek et al. 2015. Increased land cover heterogeneity also can result in increased land use area for carnivores by affecting the distribution of available resources (Hiller et al. 2015a). Specifically, black bears often display greater selection for forested areas relative to non-forested areas and may become limited in areas with high nonforest-to-forest ratio (Hiller et al. 2015b, Sollmann et al. 2016. Topographic features, such as high elevation, are often positively associated with black bear presence because high-elevation terrain offers unique seasonal food availability Pelton 1981, Frary et al. 2011). Additionally, lower elevations are more likely fragmented by transportation networks and agriculture than higher elevations (Nelleman et al. 2007). Thus, when high-elevation areas are limited and black bear population densities are low, population expansion can be restricted (Frary et al. 2011). Physical barriers, availability of land covers, and topographic features present myriad landscapes for carnivores to navigate, and in fragmented landscapes, habitat connectivity is vital for large carnivore species to persist and disperse (Noss et al. 1996).
The ideal free distribution model postulates that individual fitness and habitat selection are correlated, potentially resulting in greater densities of animals in higher-quality habitats (Fretwell and Calver 1969, Fretwell 1972, Loegering and Fraser 1995. Through the connection of such higher-quality habitat patches, dispersing individuals can contribute to population expansion. In the Midwest United States, cougar populations appear to be expanding through steppingstone dispersal where the largest habitat patches closest to source populations are likely colonized first (MacAthur and Wilson 1967, LaRue et al. 2012). Based on the affinity for forest cover by black bears, we would expect selection for larger forest patches by black bears and, thus, higher population densities in forested patches (Herrero 1972, Hiller et al. 2015b). Further, for coarse-scale assessments of colonizing species, we would expect presence probabilities to temporally increase in selected land covers.
By 1931, black bears were considered extirpated from much of the midwestern United States (Smith and Clark 1994). Conservation efforts (e.g., successful reintroductions, improved habitat connectivity, regulated hunting) beginning in 1959 and continuing into the 1960s facilitated the recolonization of black bears in the Interior Highlands of Arkansas and Oklahoma Clark 1994, Bales et al. 2005). The population expanded northward into Missouri as evidenced by increased reports of bear sightings (Wilton et al. 2014).
Our goal was to estimate spatial and temporal shifts in probability of black bear presence in Missouri, and to assess what factors might influence these shifts. We hypothesized that black bear distribution would be correlated with forest cover and expand northward over time due to spatial availability of resources.
Because of the relatively low density of the black bear population in southern Missouri, we expected an increased probability of black bear presence over time in forested areas relative to non-forested areas. However, an increase of presence in non-forested areas may indicate the black bear population is increasing beyond the carrying capacity of forested areas, requiring individuals to utilize areas with potentially less suitable habitat. Lastly, due to the opportunistic nature of the sampling method, we expected black bear presence to be positively biased by local roads that provide opportunities for humans to observe wildlife when traveling.

STUDY AREA
Our study area comprised the state of Missouri (181,188 km 2 ) with about 17% (10,667 km 2 ) in public ownership and a human population of 5,988,927 (Raeker et al. 2010, U.S (Table 1; USGS 2011). Of the four ecoregions in Missouri, the Ozark Highlands contains 80% (53,011 km 2 ) of the state's forest cover and consists primarily of upland oak (Quercus spp.) and oak-pine (Pinus spp.) woodlands ( Fig. 1; Nigh and Schroeder 2002, Raeker et al. 2010, USGS 2011. Elevations ranged from 40 m in the southeastern Mississippi River Alluvial Basin ecoregion to 540 m in the Ozark Highlands ecoregion (Nigh andSchroeder 2002, NASA Jet Propulsion Laboratory 2009).

METHODS
To understand potential spatial and temporal shifts in the probability of black bear presence (hereafter, bear presence probability) in Missouri, we established five mechanical objectives: (1) develop and compare the statewide spatial distribution models for bear presence probabilities across two equal time periods (early, 1991-2003; late, 2004-2015); (2) evaluate the effects of various environmental factors on bear presence probabilities for both time periods; (3) predict the total area of various levels of bear presence probabilities for both time periods; (4) determine the direction of and quantify the change in the distribution of bear presence probabilities; and (5) assess the relative dispersion patterns for the distribution of bear sightings for each time period.

Data collection
We used public reports of black bear sightings in Missouri from 1991 to 2015 obtained from the Missouri Department of Conservation's (MDC) Report a Bear Sighting program (Fig. 1;MDC 2008MDC , 2016. Sightings included observation type (direct observations, e.g., sighting, remote camera image; signs, e.g., track, scat; vehicle collisions, illegal harvests, and bear incidents, e.g., bear-related property damage), date, location (i.e., geographic coordinates), and number of observed adults and dependent young. As McKelvey et al. (2008) found, analyses based on insufficiently vetted anecdotal occurrence data could have substantial conservation implications due to the prevalence of incorrectly identified species observations provided by non-experts. Additionally, for consistency, we used methods described by Wilton et al. (2014) to identify and remove unreliable reports from the dataset. Hence, we excluded reports that lacked information regarding observation type and those with an observation type classified as sign (e.g., tracks, scat, hair) from analyses because black bears have widely variable scat signatures and black bear tracks and hair can easily be mistaken for other species (e.g., domestic dogs, deer) by non-experts. We considered other observation types to be reliable since there are no species similar in physical size and appearance to black bears in Missouri. Lastly, reports with incomplete location information were excluded from analyses.
We considered 11 independent variables to estimate bear presence probability, including land cover, Shannon diversity index of land cover,  vegetation productivity, elevation, human population, and roads (Table 2; Shannon and Weaver 1949). We used spatial land cover data from the 2006 National Land Cover Database (Fry et al. 2011, USGS 2011. We modified the cell size of the 15 land covers to fit the 250 m spatial resolution of other variables (Table 3). We also estimated a relative index of land cover diversity (e.g., patch richness, patch evenness) by calculating the Shannon diversity index at the 250 m spatial resolution using the 15 land covers at the original 30 m spatial resolution (USGS 2011).
We obtained Normalized Difference Vegetation Index (NDVI) data from the Earth Resources Observation and Science (EROS) Center's Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Aqua satellite launched in 2003 and with a spatial resolution of 250 m and a 16-bit radiometric resolution (i.e., À2000 to 10,001 scale; USGS 2016). We used NDVI as a proxy for natural sources of vegetative food available during bear activity (non-hibernation) period and the statewide growing season (March-October 2003Pettorelli et al. 2006, Bojarska andSelva 2012). Though some vegetation productivity changes may have occurred between the two time periods (possibly due to succession or variation in precipitation), we acknowledge that larger scale changes caused by climate change are not detectable at the timescale of our study (Easterling et al. 2000a, b). We converted the NDVI data to an 8-bit radiometric resolution (i.e., 0-255 scale) more commonly used (D'Elia et al. 2015, Mohammadi et al. 2015 and estimated the seasonal minimum, mean, and maximum NDVI values for each pixel in the study area. We obtained digital elevation models acquired by the Advanced Spaceborne Thermal Emission and Reflection Radiometer to extract elevation and modified the cell size to meet the 250 m spatial resolution requirement for our analysis (NASA Jet Propulsion Laboratory 2009). Notes: NASA JPL, National Aeronautics and Space Administration Jet Propulsion Laboratory; NDVI, Normalized Difference Vegetation Index; SD, standard deviation; USGS, U.S. Geological Survey; USCB, U.S. Census Bureau. En dash indicates that the given metric was not able to be calculated or provided.
† Each pixel is at a resolution of 250 m ‡ Land cover was the only categorical variable, and the rest were continuous.
We used 2010 human population and road data from the U.S. Census Bureau 2010, 2012, respectively) as additional anthropogenic variables. As human population data are provided at a spatial resolution of the number of people per block and blocks varied in area, we calculated the density of people by block and modified the cell size by density to a raster with 250 m spatial resolution. We extracted roads according to their designated Feature Classification Codes defined by the Master Address File/Topologically Integrated Geographic Encoding and Referencing system. We selected four road classes including (1) primary roads (i.e., generally divided, limited-access highways within the interstate highway system or under state management distinguished by the presence of interchanges), (2) secondary roads (i.e., main arteries in the U.S. Highway, State Highway, or County Highway systems with ≥1 traffic lane), (3) local roads, rural roads, or city streets (i.e., paved nonarterial single lane roads), and (4) vehicular trails (i.e., unpaved dirt trails in rural areas requiring vehicles equipped with four-wheel drive). By estimating distance to nearest road from the centroid of each cell, we created a 250-m raster for each road class.

Statistical analyses
We modeled the bear presence probability in response to 11 independent variables, all of which were continuous except the categorical land cover variable. We estimated the bear presence probability over the early time period of 1991-2003 and the late time period of 2004-2015 using maximum entropy modeling (MaxEnt, version 3.3.3k; Phillips et al. 2006, Merow et al. 2013. We used MaxEnt to compare presence locations (i.e., bear sightings) with background locations (i.e., pseudoabsences) as our species presence data lacked absence locations (Phillips and Dud ık 2008). Though sighting location accuracy may contain error, MaxEnt is apparently resilient to spatial error up to 5 km (Graham et al. 2008, Baldwin 2009). Because the public Land cover (ID no.) Definition Open water (11) Water areas (e.g., rivers, ponds, lakes) with <25% cover (e.g., vegetation, soil) Open space development (21) Vegetated areas with <20% impervious surfaces (e.g., large-lot single-family housing units, parks, golf courses) Low-intensity development (22) Mixed areas of constructed materials and vegetation with 20-49% impervious surfaces (e.g., single-family housing units) Medium-intensity development (23) Mixed areas of constructed materials and vegetation with 50-79% impervious surfaces (e.g., single-family housing units) High-intensity development (24) Areas where people reside or work in high numbers with 80-100% impervious surfaces (e.g., apartment complexes, row houses, commercial or industrial areas) Barren land (31) Bare earthen areas (i.e., rock, sand, clay) with <15% vegetation cover (e.g., bedrock, desert pavement, gravel pits) Deciduous forest (41) Areas with >20% of vegetation cover dominated by trees >5 m tall; >75% of tree species are deciduous Evergreen forest (42) Areas with >20% of vegetation cover dominated by trees >5 m tall; >75% of tree species maintain their leaves year-round resulting in a continuously green canopy Mixed forest (43) Areas with >20% of vegetation cover dominated by trees >5 m tall; neither deciduous nor evergreen species are >75% of total tree cover Shrub/scrub (52) Areas dominated by shrubs (e.g., true shrubs, young trees) <5 m tall with canopy composed of >20% shrubs Grassland/herbaceous (71) Areas with >80% gramminoid or herbaceous vegetation; not subject to intensive management but can be grazed Pasture/hay (81) Areas with >20% of grasses, legumes, or grass-legume mixtures planted for seed or hay crop production or livestock grazing Cultivated crops (82) Areas with >20% crop vegetation cover (e.g., soybeans, corn, orchards); includes all actively tilled land Woody wetlands (90) Areas with >20% forest or shrub land vegetation cover and substrate is periodically saturated or covered with water Emergent herbaceous wetlands (95) Areas with >80% perennial herbaceous vegetation cover and substrate is periodically saturated or covered with water reports wildlife sightings opportunistically, we assumed some degree of spatial sampling bias was associated with the bear sightings data and compared eight combinations of methods to understand the potential impacts of the spatial sampling bias (Table 4). For four of the eight combinations, we included a systematic sampling correction method (i.e., spatial filtering) advised by Fourcade et al. (2014). When compared to other methods (e.g., background manipulation, splitting biased datasets) by Kramer-Schadt et al. (2013) and Fourcade et al. (2014), systematic sampling correction consistently performed better with minimizing omission and commission errors across a wide array of sampling bias types and intensities (Reddy andD avalos 2003, Merow et al. 2013). Therefore, we placed a grid with 6.49 km 2 cells (the average 50% utilization distribution for adult male and female black bear space use in southern Missouri and northern Arkansas, USA; Hiller et al. 2015a) over our sightings data for each time period and randomly reduced the number of sightings to one per cell. The second bias correction measure we employed was a bias file that attempts to emulate the sampling bias evident in the presence records with the background locations in MaxEnt by limiting the selection area of background locations (Phillips et al. 2009, Kramer-Schadt et al. 2013, Fourcade et al. 2014). Our bias file restricted background location to areas within 500 km 2 of a sighting location (the average 95% utilization distribution of adult male black bears in southern Missouri and northern Arkansas, USA; Hiller et al. 2015a). Lastly, because black bears avoid roads in some regions and not others, we included roads as independent variables in four of the eight combinations (Reynolds-Hogland and Mitchell 2007, Hiller et al. 2015b, Stillfried et al. 2015. We applied the eight bias correction combinations to the early and late time periods, resulting in a total of 16 models. Using ENMTools (version 1.4.4; Warren et al. 2010), we tested all possible pairwise combinations of independent variables for multicollinearity. We assumed multicollinearity did not compromise model results if |r| < 0.70 for any pair of variables, and for any pairs with | r| > 0.70, we retained the variable we considered most ecologically relevant (Dormann et al. 2013).
For each time period and all eight bias correction combinations, we ran MaxEnt with the following customized settings: auto features, create response curves, perform jackknife tests, logistic output format, regularization multiplier = 1, maximum number of background points = 10,000, replicates (i.e., model runs) = 10, and maximum iterations = 5000. We selected 10 replicates as opposed to 1 to assess average model behavior and performance by producing a single probability model (Phillips et al. 2006). We assessed model generality, or the ability of the model to identify attributes of the species distribution rather than artifacts of noisy sampling procedures, by conducting 10-fold cross-validation (Elith et al. 2011, Merow et al. 2013, Radosavljevic and Anderson 2014. We designated 90% of the data for model training and 10% for model validation to produce a single probability model, and then estimated average model behavior and performance (Phillips et al. 2006). For the four bias correction combinations that contained the bias file, we allowed for extrapolation to the remained of the study area.
We assessed model performance via thresholdindependent measures using receiver operating characteristic plots. Receiver operating characteristic plots compare sensitivity and 1-specificity, where sensitivity represents the ability of the model to correctly predict black bear presence and specificity is a measure of correctly predicted black bear background locations (Fielding and Bell 1997). This relationship is quantified by the ❖ www.esajournals.org area under the curve (AUC) statistic with values ranging from 0.5 to 1.0, where 0.5 indicates no greater fit than randomly expected and 1.0 indicates a perfect model fit (Baldwin 2009). The AUC statistic assesses model fit based on the probability that a black bear presence location will be ranked higher than a background location (Phillips et al. 2006). MaxEnt models are considered to have good discriminatory power with AUC values >0.7 (Hosmer and Lemeshow 1989). We quantified the spatiotemporal shift in black bear presence locations across three regions between the two time periods. We parsed the presence locations into three bins of equal sample size by latitude and identified the mean geographic centers for both time periods in each of the three regions. By measuring the directional distance between the early and late geographic centers in each region, we determined the average distance and direction in which the distribution of black bear sightings shifted in each region during our study. Since the distribution of the data was heavily skewed left with more sightings occurring in the late time period than the early time period, we used nonparametric bootstrapping in the boot package in program R to compute the bootstrapped mean (Davison andHinkley 1997, Canty andRipley 2015). To more accurately reflect the non-normal distribution of the data, we computed the bias-corrected and accelerated (BC a ) bootstrapped 95% confidence intervals (DiCiccio and Efron 1996).
After spatial filtering, our final dataset contained 971 black bear sightings (Data S1). Because the NDVI minimum and NDVI mean variables were correlated (|r| = 0.77), we excluded NDVI minimum from further analyses. Of the eight 10fold cross-validation MaxEnt species distribution models (SDMs) of bear presence probability, five had good discrimination ability for across both time periods (Fig. 3A). Of the combinations that included a method to correct for spatial sampling bias, the SC_rd combination resulted in the highest mean AUC value (early, 0.83 AE 0.05 [95% CL]; late, 0.81 AE 0.04). We selected this method because the AUC values and associated 95% CL for both time periods were greater than or did not overlap with 0.7. We also sought to avoid any potential over-prediction of areas of high and medium bear presence probability while maintaining high discriminatory power. Based on the mean percent contributions, distance to nearest local road was the most important variable for both time periods for predicting bear presence probability, followed by the categorical land cover and distance to nearest secondary road (Table 5). Remaining variables in combination contributed ≤30.6% to either predictive model.
Holding all other variables at their respective mean sample value, the bear presence probability was greatest (>0.3) proximate to local (early, <77 m; late, <160 m) and secondary (early, <152 m; late, <1225 m) roads (Appendix S1:   Fig. 4). Notable differences between the two time periods include increased bear presence probability from early to late for open water, deciduous forest, and cultivated crop covers.
For the SC_rd bias combination, we predicted areas of low (0.0-0.5), medium (0.5-0.75), and high (>0.75) bear presence probability to be 164,645, 12,615, and 1983 km 2 , respectively, for the early time period and 160,657, 16,465, and 1795 km 2 , respectively, for the late time period 5). Most increases in medium and high levels of bear presence probabilities from the early to late time periods occurred in deciduous forest (2615 km 2 ; Fig. 6) followed by herbaceous (324 km 2 ), open water (252 km 2 ), and low-intensity development (234 km 2 ) with a total increase of 3661 km 2 across all land covers. These changes resulted in spatiotemporal shifts between the geographic centers for the two northern regions of 25 km to the north-northwest for the middle region and 60 km to the northwest for the northern region with an overall statewide shift of 28.7 km to the north-northwest (332°) between the two time periods (Fig. 6).  Among the three regions, the frequency of sightings was greater during the early timer period for distances <80 km from the geographic centers but was greater for distances >80 km from the geographic centers for the late time period (Fig. 7A-C). Overall, sighting locations were 15.9 km further away from the geographic center during the late time than during the early time period (Fig. 8).

DISCUSSION
To our knowledge, our study presents the first attempt to quantify expansion rates of a recolonizing carnivore population through assessing changes in presence probabilities. We predicted an approximately 2% increase in the area of Missouri considered medium or high bear presence probability between time periods, generally occurring in a north-northwest direction with the majority of sightings occurring in the forested Missouri Ozark Highlands (MOH) ecoregion of the state. In 1991, the MOH ecoregion was estimated to have a black bear population density of <1.6 bears/100 km 2 (Smith and Clark 1994). The slight increase in density by 2012 of 1-10 bears/ 100 km 2 demonstrates slow expansion and indicates the population is likely below ecological carrying capacity in the MOH ecoregion (Sinclair et al. 1992, Sollmann et al. 2016. Bear densities in Missouri are low compared to other locations in the Interior Highlands considered to have stable black bear populations (e.g., Oklahoma: 27 individuals/100 km 2 , Arkansas: 14-25 individuals/100 km 2 ; Bales et al. 2005, Kristensen 2013. Low densities combined with few human-bear conflicts and large areas of suitable habitat may likely support higher densities of bears in the future (Sollmann et al. 2016). As black bears are habitat generalists with relative high reproductive potential, we expect black bears are still re-populating suitable habitat in the MOH ecoregion before expanding further north into sparsely forested landscapes fragmented and dominated by anthropogenic features (e.g., roads, agriculture).
Black bear presence appeared to increase with decreasing distance to nearest local and secondary roads for both time periods, with greater importance during the late time period. While bears avoid paved highways in western Maryland and unpaved roads in western North Carolina, we found a positive relationship potentially indicating habituation to human-influenced habitats (Fescke et al. 2002, Beckmann and Berger 2003, Reynolds-Hogland and Mitchell 2007. Alternatively, unhunted bear populations may use lowtraffic-volume roads as travel corridors without the seasonal mortality risk associated with hunting that can alter movements (Reynolds-Hogland andMitchell 2007, Stillfried et al. 2015). Bears Fig. 6. Change in bear presence probability within deciduous forest cover. Predicted areas of deciduous forest cover where bear presence probability increased (red), did not change (white), or decreased (blue) between early (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) and late (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) time periods in Missouri, USA. Additionally, the regional geographic centers for early (green dots) and late (red dots) time periods suggest an average shift in distribution of 28.7 km in a north-northwest direction (332°). may be more likely to use roads as travel corridors through unsuitable cover, particularly during the breeding season and for dispersal (Hiller et al. 2015b). Though not all carnivore species benefit from human-altered landscapes, some carnivores, such as brown bears (Ursus arctos), gray wolves (Canis lupus), and leopards (Panthera pardus), have recolonized or recovered portions of their historical ranges now considered urban or rural (Bateman andFleming 2012, Athreya et al. 2013). Individual behavioral modifications may balance anthropogenic mortality risk, such as from vehicular collisions, and anthropogenic resource acquisition (e.g., forage from agricultural crops and human refuse, cover from culverts and other human structures, and travel corridors such as roads; Brody and Pelton 1989, Martin et al. 2010, Duquette et al. 2017. The relationship between wildlife and anthropogenic landscapes presents a complex management scenario likely to increase with human population expansion and corresponding development. The observed positive relationship between bear presence and roads may appear stronger than expected because observation data based on public reports are likely to contain higher spatial and sampling biases in landscapes with high road densities than systematically collected data (McShea et al. 2011, Garc ıa-Rangel andPettorelli 2013). Akin to other SDMs based on public sightings datasets, a more rigorous assessment of the complex relationship between wildlife and landscape features, particularly roads, will require the development of systematic survey datasets. However, advancements in statistical methods combined with increasing public interest in wildlife conservation have allowed for greater utility of presence-only datasets (Rest et al. 2015). By using machine learning modeling techniques, opportunistic public sightings datasets can have great utility for assessing widespread species distributions at large spatial scales and can even outperform regression statistical approaches based on intensive field survey data (Hochachka et al. 2012, Santika et al. 2014, Jackson et al. 2015. Though we took multiple measures (e.g., conducting quality control, correcting for sampling bias, using machine learning modeling techniques) to reduce effects from sampling bias and because it is impossible to determine if sampling bias is completely eliminated, we assume some bias effects remain and potentially inflate the relationship between bears and roads. However, we suggest that the general relationship is valid resulting from bears opportunistically utilizing roads as travel corridors.
Land cover also influenced bear presence probability during both time periods, though slightly less during the late time period. While we expected forest cover to be correlated with black Fig. 7. Distance from geographic centers. Greater proportions of black bear sightings were located further (>80 km) from the regional geographic centers in the northern (A), middle (B), and southern (C) regions during the late time period (2004-2015; light gray bars) than in the early time period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003); dark gray bars) in Missouri, USA. Statewide, greater proportions of sightings were further (>200 km) from the geographic center in the late time period than in the early time period. bear distribution based on previous research suggesting deciduous forests are an important habitat feature for black bears, our results suggest that black bears use multiple land covers as opportunistic omnivores (Herrero 1972, Beeman and Pelton 1980, Hiller et al. 2015b). Benson and Chamberlain (2007) found female black bears selected diverse habitat types in Louisiana that varied across scales, seasons, and reproductive status. Though all three forest covers (i.e., deciduous forest, evergreen forest, and mixed forest) are important for black bears in Missouri, grasslands and emergent herbaceous wetlands were of similar importance with a significant increase in the importance of cultivated crops from early to late time periods (Wilton et al. 2014). Edible crops, such as corn, sunflower, and oats, present abundant clumped food sources where bears can conserve search energy by limiting travel time (Ditmer et al. 2015). Expanding black bear populations in northern Minnesota and northern Lower Peninsula of Michigan also use land covers with a wide range of human impacts and forage options (Carter et al. 2010, Ditmer et al. 2015. Opportunistic strategies allow black bears to alter behaviors to meet life requisites in fragmented landscapes, resulting in black bear populations expanding into and persisting in areas previously considered uninhabitable (Carter et al. 2010). Thus, we can expect colonization of similar areas in Missouri, though the rate of expansion would likely be low.
Our results suggest the greatest bear presence probability (i.e., probability of a black bear sighting) occurs in predominantly anthropogenic landscapes (e.g., high road densities, rural development, agriculture). Supporting evidence from previous studies suggests that forest cover fragmented with anthropogenic activities, such as urban development, presents an environment that likely facilitates human-bear interactions, including sightings (McFadden-Hiller et al. 2016). Concatenated with the findings from Wilton et al. (2014), our results suggest that the black bear population range is expanding in multiple directions (e.g., north, east, west) and the distribution has become more dispersed in Missouri. We observed a continuance of the bimodal Fig. 8. Mean distance from regional geographic centers. The regional mean distances from the geographic center of bear sighting locations in the northern, middle, and southern regions increased from 65.2 km (bias-corrected and accelerated [BC a ] bootstrapped 95% confidence limit [CL] = 57.6-74.9), 73.9 km (BC a 95% CL = 65.9-82.7), and 76.8 km (BCa 95% CL = 67.6-88.7), respectively, in the early time period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003); dark gray bars) to 80.7 km (BCa 95% CL = 74.5-87.2), 92.5 km (BCa 95% CL = 85.7-99.0), and 104.0 km (BCa 95% CL = 96.4-112.1), respectively, in the late time period (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015); light gray bars) in Missouri, USA. Statewide, sighting locations were 15.9 km further from the geographic center during the late time period than the early time period. distribution in the number of black bear sightings reported annually as Wilton et al. (2014), which may be more reflective of public outreach efforts by MDC than true fluctuations in bear population abundance , Danielsen et al. 2009). Missouri Department of Conservation initiated campaigns for public information requests in 1990 and 1991 likely resulting in the first bimodal peak from 1991 to 1993 (MDC 1993). The second bump in public sightings starting in 2009 may likely have resulted from MDC's Report a Bear Sighting System going online in 2008 (MDC 2008). While we diligently corrected for various biases often associated with public sightings datasets (e.g., species misidentification, location error, sampling error), though unquantifiable, we acknowledge the existence of the novelty bias where the likelihood of reporting wildlife presence decreases with increasing commonality of the species , Danielsen et al. 2009, Lee et al. 2010). Specifically, it is possible that the decline in probability of presence in southeastern Missouri observed between the two time periods was a consequence of a reduction in reporting bear observations, and further a consequence of increasing familiarity of black bears by people (Bowman et al. 2001, Howe et al. 2010). However, there is no direct evidence available to support or refute this statement. While we estimated potential black bear range expansion at a rate of 1.1 km/yr, suggesting colonization of this species across Missouri is very low, we suggest that any deviation from the potential rate of expansion as a consequence of reporting bias would be ecologically unimportant. Additionally, given that public outreach programs occurred during both time periods and that our results are based on a combination of years within a time period rather than on individual years, we assert that our results remain accurately reflective of relative black bear population trends in Missouri. With that, we expect public sightings to continue increasing as bears expand into areas of northern and eastern Missouri with increased agriculture and human development.
As large carnivore populations are low density and difficult to observe, monitoring across large spatial extents poses many challenges for wildlife agencies. Though numerous survey methods are available for monitoring carnivores, such as radiotelemetry, camera surveys, and genetic sampling, such techniques are often impractical across large spatial scales and limited by personnel, time, and budget restrictions (Gompper et al. 2006). Public sightings of wildlife can serve as an alternative for monitoring carnivore populations across large spatial extents when appropriately assessed for data quality and biases (McKelvey et al. 2008, Kramer-Schadt et al. 2013).

ACKNOWLEDGMENTS
We thank W. Cooke, Q. Meng, T. Hiller, and A. Mercer for geospatial and analytical discussions that improved our study and W. Cooke for providing early reviews of our manuscript. We thank the numerous biologists from the Missouri Department of Conservation, particularly J. Beringer, for collecting statewide black bear sightings. This study was supported by the Forest and Wildlife Research Center at Mississippi State University and the Missouri Department of Conservation.