Conservation implications of sex‐specific landscape suitability for a large generalist carnivore

Terrestrial mammal distribution models typically do not differentiate between sexes when making spatial predictions, which could have important conservation implications. As male carnivores are usually more risk tolerant and travel longer distances, male potential range should be larger and include more human‐modified areas than female range. To evaluate if differences between females and males could influence their conservation planning, we quantified sex‐specific suitable range for a recolonizing population of American black bears (Ursus americanus).

Differing body sizes and reproductive strategies can influence space use of mammalian females and males independently or interactively (Rode, Farley, & Robbins, 2006;Rubin & Bleich, 2006). For example, large carnivore males can be more tolerant of human-modified areas (Riley et al., 2003), are often involved in more vehicle collisions (Maehr, 1997) and are more likely to cause human-wildlife conflicts (Ditmer, Garshelis, Garshelis, Noyce, Haveles, & Fieberg, 2015;Odden, Smith, Aanes, & Swenson, 1999). When differences between sexes exist, lack of discrimination could result in underestimating important variables for one sex while overestimating them for the other (Da Conde et al., 2010), and loss or modification of a particular land cover may affect one sex more than the other, which in turn can influence population structure and growth (Bowyer, Pierce, Duffy, & Haggstrom, 2001;Kie & Bowyer, 1999). If land use differentially influences male and female mortality risk and sex ratios of a population (Reid & Peery, 2014;Rubin & Bleich, 2006), female decline could be masked by stable species occurrence. In ursids, an increased number of young males can be mistaken as population growth, even when female population is low and declining (Wielgus & Bunnell, 1994). Additionally, distribution models based on data from both sexes (or male-biased) could overestimate the female (i.e., reproductive) range (Fernández, Selva, Yuste, Okarma, & Jakubiec, 2012) and potentially misrepresent the species distribution range.
Species distribution models (SDMs), also referred to as habitat suitability models (Hirzel & Le Lay, 2008), combine species occurrence data with environmental variables to predict distributions usually across large spatial scales (Franklin, 2010;Shabani, Kumar, & Ahmadi, 2016). Species distribution models can be used to predict suitable conditions for species survival and to extrapolate areas where the species could occur based on known locations (Franklin 2010). A promising tool within SDMs is machine-learning methods, due to their nonparametric approaches and ability to overcome difficulties such as multicollinearity (Farrell et al., 2019).
Through the framework of machine learning and SDMs, we aim to quantify sex-specific distribution range for a low-density recolonizing population of American black bears (Ursus americanus; Family: Carnivora), evaluating if differences between females and males could influence their conservation planning. Black bears are large opportunistic omnivores that exhibit sexually dimorphic body size; males typically weight 50%-100% more than females (Hunter & Barrett, 2011). As sex ratio can strongly influence population growth (Caughley, 1977;Milner, Nilsen, & Andreassen, 2007), quantifying the occurrence of sex-specific distribution would be key in developing state-wide strategies to manage female black bear establishment, reproduction and recolonization. In addition, as male bears are more likely to occur in human-modified areas (Ditmer, Noyce, Fieberg, & Garshelis, 2018;Merkle, Robinson, Krausman, & Alaback, 2013) and involved in human-bear conflicts (Ditmer, Burk, & Garshelis, 2015;, by mapping male distribution we can better inform decision-makers concerning locations for outreach and conflict prevention.

| Study area
The study area was the state of Missouri, USA, where the recolonizing black bear population occurs primarily in the southern Ozark Highlands (Wilton, Belant, & Beringer, 2014). This ecoregion is characterized by karst topography, with elevation ranging from 76 to 274 m (Karstensen, 2010). Agriculture (e.g., corn, cattle) is more prevalent in the western region, while forest is more prevalent in the eastern region (Karstensen, 2010). The state has a humid temperate continental climate, with a mean daily minimum temperature (January) of −8 to −4°C, mean daily maximum temperature (July) of 29 to 32°C, and mean annual precipitation of 102 to 132 cm (National Oceanic and Atmospheric Administration, 2013).

| Data collection
Black bears were captured in the southern Ozark Highlands using modified Aldrich foot snares (Johnson and Pelton, 1980)

| Species distribution modelling
We used black bear telemetry data spanning from September 2010 to October 2018. To decrease spatial and temporal autocorrelation, we randomly subsampled telemetry locations to no more than one location per day per individual (Hiller, Belant, Beringer, Tyre 2015, Gantchoff, Wang, Beyer, & Belant, 2018. To minimize overrepresentation of individuals with greater number of telemetry locations (Edrén, Wisz, Teilmann, Dietz, & Söderkvist2010), we first calculated the median number of locations separately for females and males (371 and 92, respectively). For individuals with more than the median value, we randomly subsampled the median number of locations from the total for that individual and used the subsample for subsequent modelling.
We considered 11 variables: land cover (four categories), Shannon diversity index for land cover, vegetation productivity (mean and maximum), elevation, human density and distance to roads (minor and major). For land cover, we used data from the National Land Cover Database (NLCD; USGS, 2011). To reduce the number of variables for modelling, we re-grouped the 15 original land covers in four categories: forest cover, other natural cover (e.g., herbaceous, shrub), agriculture cover (crops and livestock) and developed cover.
We modified the 30-m cell size of NLCD data to fit the 250-m spatial resolution of other variables. We calculated the Shannon diversity index (Shannon & Weaver, 1949)  We used NDVI as a proxy for vegetation productivity (Bojarska & Selva, 2012;Pettorelli et al., 2006). We calculated the mean and maximum NDVI values for each cell in the study area during bear activity period (March to November). 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 (NASA Jet Propulsion Laboratory, 2009). 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 as 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 developed two distance to road layers: major roads (i.e., interstates and highways) and minor roads (i.e., local roads, rural roads and city streets).
By estimating distance to nearest road from the centroid of each cell, we created a 250-m raster for each layer.
We performed the following steps to develop sex-specific models of potential black bear distribution. We used the Spatial Analyst tool in ArcGIS Desktop (Environmental Systems Research Institute, Redlands, CA) to calculate a correlation matrix among all independent variables. No pair comparison had |r| > 0.70 (Dormann et al., 2013); therefore, we assumed multicollinearity did not compromise models results. As models require background data (e.g., pseudo-absence points), we generated a randomly drawn sample of background points for the extent of study area equal to the number of black bear locations. We developed three single distribution models (GBM, RF and Maxent) using the 11 spatial variables previously described (the specific modelling parameters and values used by biomod2 for GBM, RF and Maxent can be found in Appendix S1: Table S1.1). Models were calibrated using the locations of 80% of individuals as training data and 20% of individuals as test data. Models were evaluated using the area under the curve (AUC) of a receiver operating characteristic plot (ROC) and the true skill statistic (TSS) because of their independence from prevalence in the species data (Allouche, Tsoar, & Kadmon, 2006). The ROC is a curve of true positive rate (i.e., sensitivity) against false positive rate (i.e., 1-specificity), and values range from 0 to 1, with 0.5 being equivalent to random predictions (Hilden, 1991). The true skill statistic (TSS; traditionally used for assessing the accuracy of weather forecasts) ranges from 1 to −1, a value of 1 indicating perfect agreement and values of ≤0 equivalent to random predictions (Allouche et al., 2006). We created an ensemble model by weight-averaging all single models (GBM, RF, Maxent) proportionally to their evaluation metrics scores  which resulted in a final ensemble model representing a continuous likelihood of presence (because we used presence-only data, we were unable to calculate probability of presence; Franklin 2010). To compare distribution ranges for females and males, we transformed the ensemble model to a binary format (presence/absence) by defining a likelihood threshold to differentiate between cells predicted to be occupied and cells predicted to be unoccupied. We used an op-

| RE SULTS
We obtained 16,105 locations from 57 females and 2,976 locations from 43 males that were adequate for analysis (Appendix S2: Figure   S2.3). For females, 46 random individuals were used for developing models (13,562 locations) and 11 individuals for validation (2,543 locations). For males, 35 random individuals were used for modelling (2,335 locations) and 8 for validation (641 locations).
All female species distribution models had strong predictive performance; ROC and TSS values for testing data were 0.972 and 0.83, respectively, for GBM, 0.98 and 0.88 for RF, and 0.97 and 0.83 for Maxent (Appendix S2: Figure S2 and specificity of 87.8%. Variables most influential for males were forest cover (36%), max NDVI (21%), agriculture cover (20%) and elevation (12%; Appendix S3: Table S3.1). Similar to females, male locations were more likely to occur in forest cover, in areas with greater elevation, and in areas with greater mean and max NDVI values, and less likely in agriculture cover (Appendix S3: Figure S3.2). Variable permutation scores were greatest for elevation for males and females (Appendix S3: Table S3.2, Table S3.3).
The committee-averaged model for both sexes indicated greatest uncertainty in model prediction corresponded to the periphery of the potential distribution, whereas there was high agreement in the interior of the potential distribution (Appendix S2: Figure S2.2).
The threshold value to transform the continuous presence likelihood to a binary map was 0.53 for females (94% sensitivity) and 0.47 for males (89% sensitivity). The overlap of female and male binary models ( Figure 2) suggests that male distribution range (22,359.5 km 2 ) was 66.8% larger than female distribution range (14,939.6 km 2 ).
For both sexes, >90% of their distribution range consisted of forest cover (Table 1), and use of non-forested areas was relatively low.
However, the larger male range included 85% more developed areas and 54% more agriculture than female range (total area, Table 1).
Proportionally, males were 23% more likely to be present in developed areas, 11% more likely in cultivated crops and 56% more likely in non-forested natural areas (Table 1).

| D ISCUSS I ON
Our black bear sex-specific distribution models indicated that the distribution male range in Missouri was considerably larger than for females, and males were predicted to occur in more human-modified F I G U R E 2 Predicted sex-specific black bear (Ursus americanus) suitable range in Missouri, USA TA B L E 1 Land cover area (km 2 ) and proportion within female and male suitable range for American black bears (Ursus americanus) in Missouri, USA areas and more likely to be present in developed land and agriculture. In contrast, female range (i.e., reproductive range) was more constrained and included less human-modified areas. Since our study population is recolonizing and therefore below carrying capacity, it is possible that individuals, mostly females, are free of strong competition and only use the most suitable areas, a pattern that could change as these areas become saturated (Duprè, Genovesi, & Pedrotti, 2000). Selecting the most suitable areas could result in mapping a more restricted predicted range; however, similar to our results, the modelled distribution range for reproductive female brown bears (U. arctos) in Poland was smaller than the range based on occurrence of both sexes (Fernández et al., 2012), suggesting a general pattern may occur that is independent of density. Sex-specific behaviour has been observed in several ursid species, such as foraging behaviour in Asiatic black bears (U. thibetanus; Koike et al., 2012), road tolerance in black bears (Loosen, Morehouse, & Boyce, 2019) and brown bears (Mace et al., 1996), and spatial segregation in brown bears (Rode et al., 2006). Highlighting the importance of sex-specific modelling, a habitat-based conservation framework for brown bears in Canada was based on adult female habitat, as they are the most sensitive sex-age class for population trajectories and their habitat use differed from males (Nielsen, Stenhouse, & Boyce, 2006). A follow-up study, conducted in the same area but at a larger scale, focused on female-specific spatial management of the impact of roads on bear demographics (Boulanger & Stenhouse, 2014).
Modelling species distribution without considering potential differences between females and males, or the use of male-biased data, could result in the assumed species' reproductive range to be overestimated. This overestimation could lead to inaccurately assigning a species a better conservation status, overestimating a region's species richness or underestimating species range contractions and shifts.
Another potential problem from not differentiating between sexes is over-or underestimation of smaller scale presence-environment relationships. Previous species distribution models developed in Missouri from citizen reported black bear sightings found a strong positive effect of roads and negligible effect of forested areas (McFadden-Hiller & Belant, 2018), which contrasts with our results.
Missouri bear sightings are likely male-biased due to males moving larger distances (Gantchoff et al., 2018) and more likely to occur in human-modified areas (this study). Additionally, wildlife sightings are more likely to occur closer to human-modified areas due to higher human density. Even when multiple measures are taken to reduce sampling bias, the data available to develop species distribution models (e.g., sightings, tracks, museum specimens) can have constraints that cannot be completely eliminated; it is important to acknowledge those limitations as well as the potential for biased results. Not all data types or species distribution model approaches are valid for all applications, and using inadequate species distribution results may lead to misguided decisions and suboptimal use of resources (Guillera-Arroita et al., 2015).
Black bears are forest obligate species, and as expected, their distribution in this study was strongly associated with the distribution of forest cover (Sollmann, Gardner, Belant, Wilton, & Beringer, 2016;Wilton et al., 2014). The proportion of forest in the potential distribution range for both sexes (93%-95%) was similar to the individual proportion of forest within annual and seasonal home ranges for this same study area (93% F, 89% M; Gantchoff et al., 2018), which suggests non-forested land covers are a small but consistent and scale independent fraction of areas used by black bears in Missouri.
Second to land cover, elevation and vegetation productivity where most influential, suggesting areas of higher elevation and more productive are preferred likely because they offer more food resources and less human disturbance; human activity is often higher in areas of low elevation with less topographic relief (e.g., agriculture; Napton, Auch, Headley, & Taylor, 2010). Though the predicted range including human-modified areas was low for both sexes, the total area for males was considerably larger and males were more likely to use them, supporting with our predictions. Male black bears are more likely to use human-modified landscapes (Ditmer et al., 2018), use anthropogenic resources (Merkle et al., 2013) and are the primary crop depredators (Garshelis, 1989;Ditmer, Burk, et al., 2015).
Female black bears could occupy smaller areas and areas with adequate, but poorer, resources than males because of high dietary and physiological plasticity (Hilderbrand et al., 2018;McDonald & Fuller, 2005 Clapham, 1996) or when strong spatial segregation between sexes occurs for most of the year (e.g., grey seal Halichoerus grypus; Austin, Bowen, & McMillan, 2004).
Sex-specific space use and distribution range is possibly a widespread phenomenon that, while documented and acknowledged to exist, it is often not explicitly incorporated in applied ecology projects, particularly for terrestrial carnivores (Ahmadi et al., 2017;Klaassen & Broekhuis, 2018, but see Nielsen et al., 2006. By developing separate spatial models for each sex, we have identified this limitation for black bears and therefore suggest that traditional SDMs might be overestimating the biologically relevant distribution range of many species, particularly if both sexes have similar overall requirements but different tolerances to disturbance. In addition, models used to predict ranges, particularly for endangered species, are possibly not conservative, which in turn suggests that the conservation status of many species could be worse than currently believed or have different relationships with their environment than currently understood. Particularly for threatened and endangered species, not considering the distribution or needs of reproductive individuals could adversely affect long-term persistence of populations. As spatially explicit models of conservation priorities are becoming increasingly important (Brum et al., 2017;Rosauer et al., 2017), by acknowledging these limitations and developing improved models when data are available, we can improve our efforts to better focus resources in areas where population growth can occur, more accurately assess species' conservation status, and better identify areas vital for species persistence.

DATA ACCE SS I B I LIT Y
Data used in this manuscript will be available in Dryad: https ://doi. org/10.5061/dryad.j5b27q3.