Habitat suitability and connectivity modeling reveal priority areas for Indiana bat (Myotis sodalis) conservation in a complex habitat mosaic

Conservation for the Indiana bat (Myotis sodalis), a federally endangered species in the United States of America, is typically focused on local maternity sites; however, the species is a regional migrant, interacting with the environment at multiple spatial scales. Hierarchical levels of management may be necessary, but we have limited knowledge of landscape-level ecology, distribution, and connectivity of suitable areas in complex landscapes. We sought to (1) identify factors influencing M. sodalis maternity colony distribution in a mosaic landscape, (2) map suitable maternity habitat, and (3) quantify connectivity importance of patches to direct conservation action. Using 3 decades of occurrence data, we tested a priori, hypothesis-driven habitat suitability models. We mapped suitable areas and quantified connectivity importance of habitat patches with probabilistic habitat availability metrics. Factors improving landscape-scale suitability included limited agriculture, more forest cover, forest edge, proximity to medium-sized water bodies, lower elevations, and limited urban development. Areas closer to hibernacula and rivers were suitable. Binary maps showed that 30% of the study area was suitable for M. sodalis and 29% was important for connectivity. Most suitable patches were important for intra-patch connectivity and far fewer contributed to inter-patch connectivity. While simple models may be effective for small, homogenous landscapes, complex models are needed to explain habitat suitability in large, mixed landscapes. Suitability modeling identified factors that made sites attractive as maternity areas. Connectivity analysis improved our understanding of important areas for bats and prioritized areas to target for restoration.


Introduction
Bats are critical components of functioning ecosystems and, in the United States of America (USA), are at risk from many factors. Hibernating populations are impacted by human disturbance and adversely affected by White-Nose Syndrome (WNS), a fungal disease caused by a non-native pathogen, Pseudogymnoascus destructans (USFWS 2007). Following hibernation, bats migrating to summer maternity areas are likely to encounter fragmented habitat (Carter et al. 2002;USFWS 2007), increased wind energy development (Arnett and Baerwald 2013;Roscioni et al. 2013), decreased prey availability via agricultural pesticide application (IDNR 2017 andSanchez-Bayo andWyckhuys 2019), and exposure to environmental contaminants (O'Shea and Clark 2002), all of which may impact survival.
The Indiana bat (Myotis sodalis) is an insectivorous species listed as ''Endangered'' under the USA Endangered Species Act (ESA) and ''Near Threatened'' on the International Union for Conservation of Nature's Red List of Threatened Species. The species communally hibernates in caves and mines (Whitaker and Brack 2002) and migrates regionally in spring (Krauel et al. 2018), traveling away from hibernacula to summer grounds, using stopover habitat to roost and forage along the way (Roby et al. 2019). At summer grounds, females form large maternity colonies, usually under the bark of dead trees, and remain in the area, rearing young until migrating back to hibernacula in the fall (Guthrie 1933). Myotis sodalis are philopatric and return to maternity areas (Gumbert et al. 2002), suggesting that maintaining suitable maternity habitat enhances fitness (USFWS 2007).
Maternity habitat conservation focuses on the protection of roosts and surrounding trees by assigning a seasonal conservation buffer around identified sites (USFWS 2018). However, ample data suggest additional strategies are essential for long-term conservation of M. sodalis. Protection beyond identified roosts would be effective in the long-term, as trees naturally deteriorate after 1-2 years, rendering them unusable (Gumbert et al. 2002;Kurta et al. 2002;O'Keefe and Loeb 2017), or within seasons when destroyed by weather (Gardner et al. 1991). Detecting bat presence and locating new roosts requires substantial time and resources. Moreover, M. sodalis use a network of trees every year (Humphrey et al. 1977;Gumbert et al. 2002), making locating all of them difficult. Numbers of hibernating M. sodalis (* 537,297 bats; USFWS 2017) suggest that many maternity colonies have yet to be located on the landscape and not all suitable areas are protected.
Myotis sodalis has been studied intensively due to its endangered status; thus, summer roost and home range-level needs of the species are well-described (e.g., Humphrey et al. 1977;Gumbert et al. 2002;Kurta et al. 2002;Carter and Feldhamer 2005;Menzel et al. 2005;Sparks et al. 2005;Britzke et al. 2006). However, M. sodalis interacts with the environment at multiple scales-fine scale for roost habitat, stand scale for foraging, and larger landscape contexts for migratory decisions. Hierarchical management may be necessary and landscape-level conservation is likely important. Migration is energetically demanding, and functionally connected habitat may allow bats to allocate less energy to searching for roosts and foraging areas. In addition, policy, regulation, and management decisions regarding endangered species in the USA are often made at state or regional government scales covering large mosaic areas. For these reasons, understanding the landscape-level habitat needs of M. sodalis is crucial to promote the recovery of the species.
Myotis sodalis habitat suitability has been modeled for relatively small study areas with little landscape variation (Pauli et al. 2015;De La Cruz and Ward 2016;Hammond et al. 2016), though Loeb and Winters (2013) assessed suitability range-wide using only climate and elevation variables. Habitat suitability models (HSMs) may be informative only for landscapes and scales in which they were created and may not extrapolate well to larger areas or mosaic landscapes (Pauli et al. 2015;Hammond et al. 2016). We sought to identify factors that are important for M. sodalis maternity distribution across a large, complex landscape and to identify areas to target conservation efforts.

Overview
Using Maxent ) and a long-term occurrence dataset to identify factors influencing the distribution of M. sodalis maternity colonies, we tested an array of simple and complex hypotheses, including some gleaned from the literature. As M. sodalis show habitat affinities at smaller scales, we also proposed the hypothesis that maternity habitat needs are complex and require ''just right'' (i.e., ''Goldilocks'') conditions at landscape scales, with more forest and water available, less agriculture and urban disturbance, lower elevations, farther from roads, and closer to hibernacula.
Using an information theoretic approach and plausible models, we mapped maternity habitat and identified suitable patches. We quantified connectivity importance for suitable patches using graph theory and used probabilistic habitat availability metrics to target areas for conservation action.

Study area
We studied M. sodalis in Illinois, USA, a large landscape (14,594,235 ha) with mixed land cover types. Historically dominated by prairie and some forest, contemporary row-crop agriculture and urbanization has transformed the landscape into a mosaic of agricultural, urban, sub/ex-urban, and fragmented natural habitat (Iverson 1988 Missouri: 195,157 bats, 36%;USFWS 2017). The Mississippi River Floodplain provides important roosting habitat for bats in this study area, likely due to snag creation via flood events (Carter and Feldhamer 2005;Bergeson et al. 2015). The southern quarter of the state encompasses the Shawnee National Forest, * 107,485 ha managed by the United States Forest Service.

Occurrence data
We acquired 294 M. sodalis summer occurrence records via United States Fish and Wildlife Service and Illinois Department of Natural Resources (data compiled from IDNR, Illinois Nature Preserves Commission, Illinois Endangered Species Board, and the Natural Heritage Database). We removed records labeled ''low positional accuracy'' and kept only maternity-associated records (captured reproductively active females, juveniles, and maternity roost locations). This left records from years 1989-2017.
To supplement these records and sample lessstudied areas, we completed 39 nights of mist-netting in 10 areas in 2017 and 31 nights in 12 areas in 2018, focusing on areas with distribution data gaps. A sampling area consisted of at least one mist-net site surveyed for at least two nights. Reproductively active females and juveniles were fitted with radio transmitters (0.27 g, Holohil Systems Ltd., Ontario, Canada) and tracked with a receiver and antenna for C 7 days (or until the transmitter failed) to locate maternity roosts.
Field-collected data were combined with the USFWS and IDNR datasets; the resulting data set spanned nearly 3 decades . In ArcMap 10.6.1 (Environmental Systems Research Institute, Redlands, California, USA), we randomly removed records within 1 km of others to address sampling bias and spatial autocorrelation, consistent with Pauli et al. (2015). To determine if suitable sites from the late 1980s and the 1990s were still informative for modern suitability, we evaluated the percent of land use change for cover types (barren, urban, water, grassland, agriculture, and forest) from 1991 to 2018 at occurrence locations and a surrounding 500-m buffer (Online Resource 1). We used a paired t-test to compare distance from major roads in 1996 versus 2018. The amount of forest around records increased by 7.2% during the time span, and amount of agriculture decreased by 4.8%. Distance to major roads significantly decreased during the span (5,798 m ± 343 SE in 1996; 6,400 m ± 359 SE in 2017; p = 0.00002). These changes around occurrence locations are potentially attributed to the ESA protections around sites and development of newer roads rendering old major roads less traveled. Given that forest increased during our occurrence time span, we retained all points and believe our models represent a conservative estimate of available suitable habitat.

Environmental variables
We initially considered 36 environmental variables predicted to be important to the ecology of M. sodalis. Environmental variables were related to climate, elevation, surrounding land cover, total forest edge, and distance to features such as hibernacula, roads, urban areas, and water sources. We resampled the resolution of land cover and all other raster variables (except climate data) from 30-m to 100-m (100 9 100-m grid cell dimensions) to reduce computational power needed for the analysis while still accounting for landscape complexity.
We sourced land cover variables using the Land Cover of Illinois (LCOI) 1999-2000 dataset (IDNR et al. 1999). We simplified land classes from 29 to 6 major classes, which were forest, agriculture, grassland, urban, water, and other. To examine the influence of specific forest types, we created a second layer of minor classes: bottomland forest, deciduous forest, coniferous forest, agriculture, grassland, urban, water, and other.
We used moving window analysis in FRAGSTATS (version 4.2;McGarigal et al. 2012) to calculate landscape metrics from major (water, agriculture, forest, urban) and minor (bottomland, coniferous, and deciduous forest) land cover classes. Per Bellamy et al. (2013), we used a multi-scale approach. A 100-m, 500-m, and 1-km radius window was centered on each raster cell to calculate total edge, total area, and number of patches for each land cover class within the window, so that each cell in the resulting rasters contained a continuous value for the surrounding landscape (* 7 ha for 100-m radius; * 95 ha for 500-m; * 345 ha for 1-km). The 500-m radius moving window roughly corresponds to M. sodalis dispersal distances from roosts to foraging sites (* 1 km; Timpone et al. 2010;Pauli et al. 2015).
We did not consider spatial scales with radii [ 1 km because this was the spatial autocorrelation distance we used to filter the occurrence data.
We used elevation data (ISGS 2003) in our models because elevation is important to M. sodalis distribution at large (Loeb and Winters 2013) and smaller spatial scales (Hammond et al. 2016). Temperature is important to distribution and maternity colony formation (Loeb and Winters 2013;Pettit and O'Keefe 2017); therefore, we assessed the effect of average maximum temperature. Temperature data were sourced from the 30-year normal (1981-2010) PRISM dataset (800-m resampled to 100-m resolution, PRISM Climate Group 2004).We used May temperatures because M. sodalis begin forming maternity colonies * 8 April-14 May in Indiana, USA (Pettit and O'Keefe 2017) and USFWS recognizes the maternity season to begin mid-May (USFWS 2018).
For all ''distance-to'' variables, we used ArcMap to calculate distance to a feature from the center of each raster cell. We created a distance-to-water variable (USGS 2018), eliminating water bodies with area \ 1 ha. We derived a distance-to-major roads variable (Illinois Department of Transportation 2018), selecting roads with annual average daily traffic rate of 2 cars/minute as major roads (per Pauli et al. 2015). We calculated distance to closest hibernation site using locations for sites surveyed within the last 10 years where C 1 M. sodalis was observed (n = 47, RA King, USFWS, Bloomington, IN Field Office, unpublished data, 29 November 2017). We predicted hibernacula would influence distributions in a study area of this size. Hibernacula are spatially clustered in Illinois' limited karst regions; therefore, knowledge of undocumented hibernation sites would be unlikely to change our results substantially.

Habitat modeling
We developed models via Maxent version 3.4.1 (Phillips et al. 2017), quantifying relative habitat suitability by finding the most uniform distribution given a set of constraints (environmental conditions; Elith et al. 2006;Phillips et al. 2006). Maxent performs well with presence-only data compared to other modeling approaches (Elith et al. 2006), improving our understanding of habitat associations for rare and elusive species (Pearson et al. 2007) for which true absences are difficult to confirm. To reduce the number of variables and identify appropriate spatial scales, we used a univariate vetting process to identify variables to populate the candidate model set, similar to methods used by Bellamy et al. (2013). We tested 36 univariate Maxent models (default settings in Maxent, but only linear and quadratic features selected). We withheld 10% of random records for testing, partitioning 142 records for training and 16 for testing; thus, AUC Test scores from Maxent output were comparable among models. This allowed us to determine discriminant ability of variables at differing scales (100-m, 500-m, and 1-km) and decide which of 2 correlated variables to retain in subsequent candidate models (Online Resource 2). For variables created at multiple scales we retained the scale with the highest discriminant ability (AUC Test score). We removed variables with AUC Test B 0.5, a score that suggested the variable performed worse than random at discriminating between occurrence locations and background points; we kept variables with AUC Test [0.5 to use in the candidate model set.
Habitat suitability models and map Using an information theoretic approach (Burnham and Anderson 2002), we tested a series of 26 competing models, including both null and global models (Table 1). Each model contained 1 to 7 variables and each variable was used at least twice in the candidate model set. These candidate models were multi-scale and included variables at the cell resolution (e.g., elevation), at the best performing scale for land-cover derived variables (100-m, 500-m, or 1-km), or that were non-scalar (e.g., ''distance-to'' variables). We calculated a Pearson's correlation matrix using SDMtoolbox (Brown et al. 2017). When 2 variables were highly correlated (r [ 0.7), we used only the variable with the higher univariate AUC Test score within a candidate model.
Of the 26 models, 7 were based on important variables from previous M. sodalis studies (Carter et al. 2002;Britzke et al. 2006;Carter 2006;Loeb and Winters 2013;Womack et al. 2013;De La Cruz and Ward 2016;Hammond et al. 2016). Sixteen models represented novel hypotheses, including 2 models accounting for potential sampling bias (e.g., more surveys may have occurred closer to hibernacula, water, and roads) and land use history bias (e.g., bats may occur in areas that are difficult to access). The global model contained 1 variable from every general environmental category, minus correlated variables. The ''Goldilocks'' model also contained variables from every general category but contained 1 less parameter than the global. We created an effective null model (AUC Test = 0.49) by using an ASCII grid generated in ArcMap, in which integer values from 1 to 100 were randomly distributed across the study area extent.
As bats have complex relationships with environmental conditions, we selected linear and quadratic features in Maxent (Elith et al. 2011). We increased the regularization multiplier from 1 (default) to 3 to reduce overparameterization (Warren and Seifert 2011). All other Maxent settings were default, including the number of background points.
We implemented ENMTools (Warren et al. 2010) in Perl to calculate AIC c values (Burnham and Anderson 2002) from raw Maxent output ASCII and LAMBDAS files. Models were ranked on deviation from the model with the lowest AIC c score (DAIC c ). Models with DAIC c B 2 are typically considered to have substantial and equivalent support given the candidate set (Burnham and Anderson 2002); we also considered a model plausible if its AIC c weight was within 1% of top model's weight. To evaluate model predictive performance, we followed Hovick et al. (2015) to create binary models using the 10% error threshold of suitability values of the training data and calculate the omission rate (proportion of test occurrences that were predicted absent). A 10% tolerance accounts for potential spatial and species identification errors in the dataset by omitting suitability scores in the lowest 10th percentile of the training data (Hovick et al. 2015). Thresholding based on occurrence data produces maps that are biologically relevant. We determined the proportion of the study area predicted to be suitable to check the efficacy of models with low omission rates (i.e., a model with a low omission rate that predicted most of the study area suitable was not informative for our objectives). Despite recent concerns of assessing presence-background models with discriminant metrics (Leroy et al. 2018;Warren et al. 2020), AUC is reported frequently in the HSM literature; thus, we provide AUC Test scores for candidate models as a threshold-independent evaluation metric. We calculated parameter importance values by summing weights across all plausible models in which each parameter appeared, with values closest to 1 considered most important. We also calculated modelaveraged parameter estimates and direction of influence for each of the variables within plausible models, using values from all plausible models.
To define suitable areas, we summed binary maps from the plausible model set using the Raster Calculator in ArcMap. In the final raster, cells with summed values [ 1 were considered suitable habitat. For example, if two models were plausible, then 2 binary maps would be summed together. The possible values in the resulting raster would range from 0 to 2, and cells with values [1 would be considered suitable.

Habitat connectivity
We quantified patch importance to functional connectivity of the habitat network (Saura and Rubio 2010) in Conefor Sensinode 2.6 (CS26; Saura and Torne 2009). We converted suitable raster cells into simplified polygons using the Raster to Polygon tool in ArcMap to represent maternity habitat patches. CS26 requires substantial computational power (Saura and Torne 2009), so we removed 6201 patches \ 40 ha in area to reduce the number of nodes. The average home range area for females in Illinois is 161.1 ha (Menzel et al. 2005) and, therefore, 40 ha is likely a conservative value. We retained 1573 patches ranging from 40 to 3,328,276 ha (mean = 2,765 ha). Using Conefor Inputs for ArcGIS Extension (Jenness 2016), we converted the landscape to a graph of nodes (suitable patches) and all possible links (Euclidean distances between patches).
In CS26, we calculated the Probability of Connectivity (PC) of the overall patch network and the connectivity importance of individual patches by the deviation in the PC when a patch was removed from the landscape (dPC k ; see Saura and Rubio (2010) and Saura and Pascual-Hortal (2007) for extensive details on PC and dPC k metrics). Three fractions sum to determine the patch importance (dPC k ): dPCintra k , dPCflux k , and dPCconnector k . The dPCintra k considers the movement within the patch k itself and is influenced by the attribute (area in this study). The dPCflux k is determined by the area of the patch and the topological position in the network. A patch only contributes to the dPC k with a dPCconnector k value if it is the most likely path between patches. The dPCconnector k fraction can identify patches that are likely important stepping-stones for reaching other patches.
Bats are capable of traversing a matrix of suitable and unsuitable patches; however, they do require protective cover, foraging habitat, and roosting stops along migration routes. We used the median distance traveled by M. sodalis between foraging bouts during spring migration (Roby et al. 2019) to define dispersal distance in CS26; if a patch was 10.1 km from another patch (based on Euclidean distance of the links), it was defined as 50% likely to be connected to another patch. Distances [ 10.1 km were assigned a lower probability of dispersal. Similar to Saura and Pascual-Hortal (2007), this process allowed for prioritization of suitable areas (patches with higher relative dPC k values) important for connectivity and areas with opportunity to improve dPC k values. To assess if available habitat patches were isolated from the network and identify those important for inter-patch connectivity, we calculated the number of patches important for each fraction of dPC, and the relative importance of each fraction by dividing the sum of dPCintra k , dPCflux k , and dPCconnector k values by the sum of dPC k .

Occurrence data
Our 2017 and 2018 mist-netting efforts yielded 28 new occurrence records: 8 capture sites and 20 roost sites. We combined these with the USFWS and IDNR records and then reduced the total dataset to 158 records after culling points with low positional accuracy, within 1 km of other records, or not associated with a maternity colony (Fig. 1). The resulting dataset spanned 45 of Illinois' 102 counties. Eighty-seven of the records were roost occurrences and 71 were capture records.

Univariate models
Most land cover variables had more discriminant ability at the 500-m scale (* 95 ha area). The exceptions were area of coniferous and bottomland forest (best at 1-km) and number of forest patches (best at 100-m). Three variables had little discriminant ability-number of urban patches in 1-km radius, (AUC Test = 0.31), area of urban in 1-km radius (AUC Test = 0.38), and area of coniferous forest in 500-m radius (AUC Test = 0.48).

Habitat suitability models
The ''Goldilocks'' model, hypothesizing that M. sodalis require complex conditions, had the lowest AIC c score and held 78% of AIC c weight ( Table 2). The global model ranked second (DAIC c = 2.54) and was also plausible based on AIC c weight (22%). Both plausible models had moderate omission rates and the b Fig. 1 Maternity-associated records in Illinois (data compiled from USFWS, IDNR, Illinois Nature Preserves Commission, Illinois Endangered Species Board, and the Natural Heritage Database, and the Illinois Bat Conservation Program; years 1989-2018; hexagons) and hibernacula in Illinois and neighboring states provided by USFWS and IDNR (diamonds). Map depicts the 158 maternity records retained after filtering out spatially autocorrelated (\ 1 km away from other records) and low positional accuracy records. Areas predicted to be suitable by plausible Maxent models are depicted in gray

Parameter importance and estimates
Five variables appeared in both plausible models and had importance values of 1-area of agriculture in 500-m radius, area of water in 500-m radius, elevation, distance to closest hibernacula, and area of urban in 500-m radius (Table 3). Response curves for these variables and area of forest in 500-m radius, which was also in the top-ranked ''Goldilocks'' model, are plotted Importance values were calculated by summing the AICc weight across the plausible models in which the parameter appeared. Parameters with ''/_q'' refer to quadratic effects. Variables are defined in Table 1 and Table S2 Fig in Fig. 2. When we considered area in terms of proportion of the window, habitat suitability was highest when the landscape was less than * 20% agriculture in a 500-m radius, and above this value suitability steadily declined (Fig. 2c). Suitability was also highest when the proportion of water was * 30-50% of the landscape (Fig. 2e), and proportion of forest was * 60-80% (Fig. 2f). Suitability decreased farther from hibernacula, reaching an asymptote at distances [ 125 km (Fig. 2a). Elevation and area of urban cover had negative linear effects (Fig. 2b, d).

Suitability map
In the final summed raster, 30% of the study area was suitable habitat (Fig. 1). There was a high concentration of suitable areas in southern Illinois, which included the Shawnee National Forest (Fig. 1). Further north, suitable areas were generally adjacent to rivers (e.g., Illinois River Valley). Fourteen of 16 test occurrences were in areas predicted to be suitable. Two test records were in unsuitable areas: a roost from 1989 and a capture location from 2011. Of the more recent test records (years 2011-2017), 9 of 10 were successfully predicted.

Habitat connectivity
We analyzed 1573 habitat patches in CS26 (Fig. 3). The most critical patch for connectivity was a single large patch covering the southern Illinois landscape and extending across the Illinois River Valley. When this patch was removed, the overall functional connectivity of the habitat network decreased by 99.7% (dPC k = 99.67; dPCintra k = 98.96; dPCflux k = 0.71; dPCconnector k = 0.001). This critical patch comprised 22.8% of the total state area and 76% of the suitable habitat available in the state. Most patches contributed to importance (dPC k ) via the dPCintra k fraction (n = 1132; Fig. 4), and fewer patches contributed through dPCflux k (n = 130) and dPCconnector k (n = 8). The relative importance of each fraction to the dPC were as follows: dPCintra (98%), dPCflux (1%), and dPCconnector (1%). There were 1,443 patches that only contributed to importance through the dPCintra fraction and were not important for inter-patch connectivity.

Discussion
Two complex models successfully identified landscape-scale factors explaining M. sodalis maternity habitat site suitability in a large mosaic landscape (88% of test occurrences were predicted present). Models derived for smaller, more homogenous study areas were implausible at the larger spatial scale of our study area and had either higher omission errors or lower discriminant ability than plausible models, supporting the idea that M. sodalis habitat needs are multi-scale. When applied to Illinois, the range-wide model created by Loeb and Winters (2013) for a climatic suitability study had lower discriminant ability compared to our plausible models and likely over-predicted suitability (67% of the study area predicted suitable), which suggests it did not adequately account for M. sodalis ecology. Therefore, effective habitat suitability models should consider hierarchical levels of bat ecology (i.e., roost habitat, foraging areas, and migration corridors connecting summer habitat and hibernacula). At the regional landscape scale considered herein, M. sodalis inhabited complex maternity habitat near medium-to large-sized water bodies, with limited agriculture and more forest cover, at lower elevations, and closer to hibernacula. For land cover variables tested at 3 scales, conditions available in the surrounding * 95 ha area (500-m radius), rather than smaller (* 7 ha; 100-m radius) or larger (* 346 ha; 1-km radius) scales, were generally more influential. The 500-m radius moving window corresponds to typical M. sodalis dispersal distances from roosts to foraging sites (* 1 km; Timpone et al. 2010;Pauli et al. 2015). However, consistent with findings for multiple bat species in Europe (Bellamy and Altringham 2015), different scales were important for different variables (e.g., bottomland forest area was best at b Fig. 3 Connectivity importance of patches (dPC k ) of modeled maternity habitat for Myotis sodalis in Illinois, USA. Larger dPC k values are more important than lower values to the overall habitat connectivity (PC) 1-km radius and number of forest patches was best at 100-m radius). This suggests that, at landscape scales, M. sodalis distribution is a function of foraging needs, not just local roosting needs and that habitat needs are likely multi-scale. The concentration of suitable areas around rivers was not surprising; low-elevation, forested riparian areas may flood regularly, which is important for snag creation (Carter 2006). Elevation is an important predictor of habitat suitability for M. sodalis in the heavily forested Southern Appalachian Mountains; however, there bats roost at higher elevations than in other parts of the species' range (Hammond et al. 2016). In forested regions of the Midwest where remote sensing data were available, Pauli et al. (2015) determined that proximity to hibernacula was not an important predictor of M. sodalis roosting habitat. However, ours is the first study to assess the relationship between M. sodalis maternity habitat suitability and distance to hibernacula across an entire state, and we found that suitability steadily decreased as distance increased. The average distance from occurrence records to the nearest hibernaculum was * 55 km and areas within 125 km of hibernacula were more suitable. While agricultural lands may provide some foraging opportunities for M. sodalis (Sparks et al. 2005;Kaiser and O'Keefe 2015), we found that bats used areas with less agriculture than is typically available in this large mosaic landscape. Our models identified a threshold where suitability decreased as proportion of agriculture within 95 ha increased above 20%. Clearing large tracts of land for agriculture may contribute to forest fragmentation and isolation of suitable maternity roost habitat. Large agricultural areas may also indirectly affect habitat suitability by decreasing biodiversity and abundance of insect prey (Sanchez-Bayo and Wyckhuys 2019). This result of our study has implications for bats globally, as agriculture is estimated to be a threat to half of bat species (Frick et al. 2019).
To maximize return on conservation investment, improvements in both habitat suitability and functional connectivity are essential. Our study was the first to examine connectivity of the distribution of M. sodalis habitat across a large landscape. Similar modeling methods were used in a study of horseshoe bats (Rhinolophus spp.) to identify key corridors and areas for conservation in a 176,000-ha area in Europe (Le Roux et al. 2017). Metrics such as dPC and dPC fractions enhanced our knowledge of the habitat network and patches important for maintaining or restoring connectivity. Based on quality of habitat patches alone, 1573 patches were suitable in our * 14,594,000-ha study area. Of these suitable patches, 72% contributed to connectivity importance, owing mostly to connectivity within the area of a patch (dPCintra); however, intra-patch connectivity is likely unimportant to migratory movements. Only 130 patches (8%) were important for inter-patch connectivity (dPCflux k and dPCconnector k ), indicating there is relatively low connectivity among suitable patches across the larger landscape we studied.
Our results have landscape-scale implications for M. sodalis management. For example, if a patch only contributes to importance via the dPCintra fraction (1,443 patches in our study; red patches in Fig. 5), this is indicative of isolation from other patches. Accordingly, in Illinois, many patches of M. sodalis habitat are isolated. The most influential landscape variables were typically at the 95-ha scale in our models, while 10.1 km is the median distance migrating M. sodalis travel before they stop in habitat to forage or roost (Roby et al. 2019). Our results suggest that allocating resources close to hibernacula may yield the greatest return on investment; thus distance to hibernacula should be considered in studies on other hibernating bat species. We recommend first preserving and restoring native prairies and forests in areas that are B 10.1 km away from existing habitat patches and \ 125 km away from known hibernacula (Fig. 5). Improving conditions in the matrix at the 95-ha scale would improve landscape-scale maternity habitat suitability, and likely increase inter-patch connectivity importance values (dPCflux and dPCconnector).
More specifically, we propose two ways to improve conservation outcomes for M. sodalis in Illinois. First, we suggest restoring habitat around existing patches near a ''Priority 2'' hibernaculum in the Illinois River Valley (''Priority'' designation by USFWS where 1000 to 10,000 individual M. sodalis have been observed hibernating; USFWS 2007; IDNR 2017; b Fig. 4 Myotis sodalis maternity habitat patches important for intra-patch connectivity (dPCintra k ) and inter-patch connectivity (dPCflux k and dPCconnector k ). These three fractions sum to the patch connectivity importance (dPC k ) Figs. 1 and 5). This is an area of conservation interest, as it is near the northernmost edge of the species' range and is predicted to be a vulnerable area under future climate change (Loeb and Winters 2013) and wind energy development scenarios (Erickson et al. 2016). Second, as wind farms threaten migratory bats (Arnett and Baerwald 2013), we suggest siting new wind farms farther from hibernacula and away from the maternity habitat network (i.e., outside of the orange suitable areas and the gray priority restoration areas in Fig. 5). Areas within the range of M. sodalis are under increasing pressure for wind development; placing farms outside of the main habitat network may reduce mortality risk during migration. Wind farm risk maps have been created for two bat species in Europe and are likely a valuable tool for conservation planning in other areas (Roscioni et al. 2013). While these models have numerous applications, important caveats exist. Sampling bias is a common problem with HSMs, potentially skewing variable importance (Elith et al. 2011). Although we incorporated a research bias hypothesis into our candidate model set and reduced bias by spatially filtering records, our occurrence dataset was acquired from multiple sources; therefore, possible errors and some sampling bias may still exist. However, because we combined data from more than one source and conducted field surveys to address distribution gaps, we likely alleviated some biases. Further, the research bias model was not plausible ( Table 2). As is the case with most HSM studies, we did not model abundance data for summer or winter populations; incorporating abundance data might alter the relative importance of some patches. Our models were created with currently available data, but as climate and land use patterns change, M. sodalis may alter its behavior and distribution (Loeb and Winters 2013); therefore, these models should be retested in the future and in different landscapes. Although we found that land use has changed in the span of our occurrence data, there is evidence that M. sodalis are loyal to maternity areas in the Midwest for decades, even in cases of extensive urbanization; e.g., the Indianapolis Airport colony was discovered in 1994 and is still present in the area 26 years later (Whitaker and Sparks 2008;Bergeson et al. 2020). Though areas around maternity colonies in Illinois may have become more suitable (increased forest) during the span of our occurrence data, potentially via habitat protections under the ESA, this may not be the case for unprotected areas in the state.
Models created at larger, range-wide scales might provide more insight and could help managers understand if range-wide management is needed or if statewide scales are sufficient to promote recovery. Myotis sodalis cross geopolitical boundaries, so larger-scale models may be necessary. Additional work addressing migration-related questions may lead to better landscape-scale models as we gain more insight into the behavior, dispersal capabilities, and habitat needs of migrating bats. We recognize that M. sodalis is but one species of concern and that resources are limited; therefore, we recommend future studies that overlay HSMs for multiple sensitive species (Cooper-Bohannon et al. 2016) and identify important patches for a variety of organisms to maximize return on conservation investment. Though we model habitat and connectivity for one bat species, our methods could be applied to numerous species in myriad ecological contexts to inform conservation action for imperiled bats worldwide. We conclude that including data collected across several decades identifies factors influencing suitability and incorporating multi-scale variables in complex models is crucial when modeling bat habitat at larger spatial scales.
Ethical approval Work was conducted under IACUC Protocol 16074 approved by the University of Illinois Committee and all appropriate state and federal research permits (USFWS Recovery Permit TE11170C). We strictly followed all current White-Nose Syndrome decontamination protocols.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.