DOI: 10.5061/dryad.4j0zpc8f1
title: Data for: Predicting habitat suitability for Townsend’s big-eared bats across California in relation to climate change

description: Aim: Effective management decisions depend on knowledge of species
 distribution and habitat use. Maps generated from species distribution
 models are important in predicting previously unknown occurrences of
 protected species. However, if populations are seasonally dynamic or
 locally adapted, failing to consider population level differences could
 lead to erroneous determinations of occurrence probability and ineffective
 management. The study goal was to model the distribution of a species of
 special concern, Townsend’s big-eared bats (Corynorhinus townsendii), in
 California. We incorporate seasonal and spatial differences to estimate
 the distribution under current and future climate conditions. Methods: We
 built species distribution models using all records from statewide roost
 surveys and by subsetting data to seasonal colonies, representing
 different phenological stages, and to Environmental Protection Agency
 Level III Ecoregions to understand how environmental needs vary based on
 these factors. We projected species’ distribution for 2061-2080 in
 response to low and high emissions scenarios and calculated the expected
 range shifts. Results: The estimated distribution differed between the
 combined (full dataset) and phenologically-explicit models, while
 ecoregion-specific models were largely congruent with the combined model.
 Across the majority of models, precipitation was the most important
 variable predicting the presence of C. townsendii roosts. Under future
 climate scnearios, distribution of C. townsendii is expected to contract
 throughout the state, however suitable areas will expand within some
 ecoregions. Main conclusion: Comparison of phenologically-explicit models
 with combined models indicate the combined models better predict the
 extent of the known range of C. townsendii in California. However, life
 history-explicit models aid in understanding of different environmental
 needs and distribution of their major phenological stages. Differences
 between ecoregion-specific and statewide predictions of habitat
 contractions highlight the need to consider regional variation when
 forecasting species’ responses to climate change. These models can aid in
 directing seasonally explicit surveys and predicting regions most
 vulnerable under future climate conditions.
description: Study area and survey data The study area
 covers the U.S. state of California, which has steep environmental
 gradients that support an array of species (Dobrowski et al. 2011).
 Because California is ecologically diverse, with regions ranging from
 forested mountain ranges to deserts, we examined local environmental needs
 by modeling at both the state-wide and ecoregion scale, using U.S.
 Environmental Protection Agency (EPA) Level III ecoregion designations and
 there are thirteen Level III ecoregions in California (Table S1.1)
 (Griffith et al. 2016). Species occurrence data used in
 this study were from a statewide survey of <em>C.
 townsendii</em> in California conducted by Harris et al. (2019).
 Briefly, methods included field surveys from 2014-2017 following a
 modified bat survey protocol to create a stratified random sampling
 scheme. <em>Corynorhinus townsendii</em> presence at roost
 sites was based on visual bat sightings. From these survey efforts, we
 have visual occurrence data for 65 maternity roosts, 82 hibernation roosts
 (hibernacula), and 91 active-season non-maternity roosts (transition
 roosts) for a total of 238 occurrence records (Figure 1, Table
 S1.1). Ecogeographical factors We
 downloaded climatic variables from WorldClim 2.0 bioclimatic variables
 (Fick & Hijmans, 2017) at a resolution of 5 arcmin for broad-scale
 analysis and 30 arcsec for our ecoregion-specific analyses. To calculate
 elevation and slope, we used a digital elevation model (USGS 2022) in
 ArcGIS 10.8.1 (ESRI, 2006). The chosen set of environmental variables
 reflects knowledge on climatic conditions and habitat relevant to bat
 physiology, phenology, and life history (Rebelo et al. 2010, Razgour et
 al. 2011, Loeb and Winters 2013, Razgour 2015, Ancillotto et al. 2016). To
 trim the global environmental variables to the same extent (the state of
 California), we used the R package “raster” (Hijmans et al. 2022). We
 performed a correlation analysis on the raster layers using the
 “layerStats” function and removed variables with a Pearson’s coefficient
 &gt; 0.7  (see Table 1 for final model variables). For future climate
 conditions, we selected three general circulation models (GCMs) based on
 previous species distribution models of temperate bat species (Razgour et
 al. 2019) [Hadley Centre Global Environment Model version 2 Earth Systems
 model (HadGEM3-GC31_LL; Webb, 2019), Institut Pierre-Simon Laplace Coupled
 Model 6th Assessment Low Resolution (IPSL-CM6A-LR; Boucher et al., 2018),
 and Max Planck Institute for Meteorology Earth System Model Low Resolution
 (MPI-ESM1-2-LR; Brovkin et al., 2019)] and two contrasting greenhouse
 concentration trajectories (Shared Socio-economic Pathways (SSPs): a
 steady decline pathway with CO2 concentrations of 360 ppmv (SSP1-2.6) and
 an increasing pathway with CO2 reaching around 2,000 ppmv (SSP5-8.5)
 (IPCC6). We modeled distribution for present conditions future (2061-2080)
 time periods. Because one aim of our study was to determine the
 consequences of changing climate, we changed only the climatic data when
 projecting future distributions, while keeping the other variables
 constant over time (elevation, slope). Species
 distribution modeling We generated distribution maps
 for total occurrences (maternity + hibernacula + transition, hereafter
 defined as “combined models”), maternity colonies , hibernacula, and
 transition roosts. To estimate the present and future habitat suitability
 for <em>C. townsendii</em> in California, we used the maximum
 entropy (MaxEnt) algorithm in the “dismo” R package (Hijmans et al. 2021)
 through the advanced computing resources provided by Texas A&M
 High Performance Research Computing. We chose MaxEnt to aid in the
 comparisons of state-wide and ecoregion-specific models as MaxEnt
 outperforms other approaches when using small datasets (as is the case in
 our ecoregion-specific models). We created 1,000 background points from
 random points in the environmental layers and performed a 5-fold cross
 validation approach, which divided the occurrence records into training
 (80%) and testing (20%) datasets. We assessed the performance of our
 models by measuring the area under the receiver operating characteristic
 curve (AUC; Hanley & McNeil, 1982), where values &gt;0.5
 indicate that the model is performing better than random, values 0.5-0.7
 indicating poor performance, 0.7-0.9 moderate performance and values of
 0.9-1 excellent performance (BCCVL, Hallgren et al., 2016). We also
 measured the maximum true skill statistic (TSS; Allouche, Tsoar, &
 Kadmon, 2006) to assess model performance. The maxTSS ranges from -1 to
 +1:values &lt;0.4 indicate a model that performs no better than
 random, 0.4-0.55 indicates poor performance, (0.55-0.7) moderate
 performance, (0.7-0.85) good performance, and values &gt;0.80 indicate
 excellent performance (Samadi et al. 2022). Final distribution maps  were
 generated using all occurrence records for each region (rather than the
 training/testing subset), and the models were projected onto present and
 future climate conditions. Additionally, because the climatic conditions
 of the different ecoregions of California vary widely, we generated
 separate models for each ecoregion in an attempt to capture potential
 local effects of climate change. A general rule in species distribution
 modeling is that the occurrence points should be 10 times the number of
 predictors included in the model, meaning that we would need 50
 occurrences in each ecoregion. One common way to overcome this limitation
 is through the ensemble of small models (ESMs) (Breiner et al. 2015.,
 2018; Virtanen et al. 2018; Scherrer et al. 2019; Song et al. 2019)
 included in ecospat R package (references). For our ESMs we implemented
 MaxEnt modeling, and the final ensemble model was created by averaging
 individual bivariate models by weighted performance (AUC &gt; 0.5). 
 We also used null model significance testing with to evaluate the
 performance of our ESMs (Raes and Ter Steege 2007). To perform null model
 testing we compared AUC scores from 100 null models using randomly
 generated presence locations equal to the number used in the developed
 distribution model. All ecoregion models outperformed the null expectation
 (p&lt;0.002). Estimating range shifts
 For each of the three GCMs and each RCP scenario, we converted
 the probability distribution map into a binary map (0=unsuitable,
 1=suitable) using the threshold that maximizes sensitivity and specificity
 (Liu et al. 2016). To create the final maps for each SSP scenario, we
 summed the three binary GCM layers and took a consensus approach, meaning
 climatically suitable areas were pixels where at least two of the three
 models predicted species presence (Araújo and New 2007, Piccioli Cappelli
 et al. 2021). We combined the future binary maps
 (<em>fmap</em>) and the present binary maps
 (<em>pmap</em>) following the formula <em>fmap
 </em>x 2 + <em>pmap</em> (from Huang et al., 2017) to
 produce maps with values of 0 (areas not suitable), 1 (areas that are
 suitable in the present but not the future), 2 (areas that are not
 suitable in the present but suitable in the future), and 3 (areas
 currently suitable that will remain suitable) using the raster calculator
 function in QGIS. We then calculated the total area of suitability, area
 of maintenance, area of expansion, and area of contraction for each binary
 model using the “BIOMOD_RangeSize” function in R package “biomod2”
 (Thuiller et al. 2021).
description: We cannot provide original bat occurrence points used in analyses
 in the interest of protecting sensitive colony locations from disturbance
 or vandalism.