Mapping fine-scale distribution of the northern pika Ochotona hyperborea considering duality in microhabitat thermal conditions

20 Species distributions are frequently modeled using predictors that exceed the spatial scale 21 experienced by the focal species. Incorporating fine-scale environmental conditions is 22 therefore expected to lead to more realistic model predictions. However, the importance of the 23 variety in the existing fine-scale conditions on species distribution remains poorly assessed 24 although species can effectively utilize multiple microhabitats for behavioral adaptation to 25 withstand climate change impacts. Here, we developed a fine-scale species distribution model 26 based on ambient air conditions for the northern pika (Ochotona hyperborea), a small 27 lagomorph found in rocky landforms, to first understand the improvement in the model 28 performance from the conventional coarse-scale model. We then assessed how model 29 predictions alter by incorporating the rock-interstice microclimates that are found in their 30 habitats for the baseline (1981–2010) and future periods (2041–2100). The fine-scale model 31 performed better and overall predicted lower habitat suitability across the study area than the 32 coarse-scale model. Incorporation of rock-interstice microclimate increased the habitat 33 suitability markedly relative to predictions based on ambient conditions, which resulted in 34 predicting more suitable areas in lower elevations and more areas remaining suitable into the 35 future. This result suggests that the northern pika may withstand the negative impacts from 36 rising ambient temperatures by effectively utilizing rock interstices via behavioral adaptation. 37 Our findings highlight the importance of analyzing species distribution at fine scales and 38 considering local environmental heterogeneity, which helps species mitigate the adverse 39 impacts of climate change, for conservation under climate change. 40 HIGHLIGHTS 41 1. Species use a wide variety of microhabitats and experience thermal conditions existing 42 locally. 43 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 3 2. In the northern pika, a cold-adapted lagomorph, fine-scale species distribution model 44 performed better than the coarse-scale model. 45 3. Complex topographical features that locally buffer ambient thermal conditions were 46 predicted to increase habitat suitability and elongate species persistence in the future. 47 4. Local environmental heterogeneity that helps species mitigate climate change impact will 48 be important for conservation. 49 50

behavioral thermoregulation, climate change, microclimate, northern pika, species 51 distribution model, topography 52 53 54 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 4

55 Climate change is having a profound influence on species distribution worldwide. 56 Evidence from historical resurveys shows that species distribution limits are generally 57 expanding at the leading edge and contracting at the trailing edge towards higher elevations and 58 latitudes (Lenoir et al., 2020; Moritz et al., 2008) . These range shifts are predicted not only to 59 affect species extinction risks (Urban, 2015), but also ecosystem functions and services (Pecl 60 et al., 2017). It is therefore important in contemporary biogeography to assess accurately how 61 species will change their spatial distribution in the future and identify conservation priority 62 areas to mitigate those i mpacts. Species distribution model s (SDMs) have become extremely 63 useful in this role by relating species occurrence with environmental data and enabling forecasts 64 of how suitable habitats and associated distribution ranges will change over time (Elith & 65 Leathwick, 2009) . Yet, SDMs have been criticiz ed for analyzing species occurrence and 66 environment relationships at spatial resolutions that are often too coarse relative to the actual 67 scales at which species utilize those habitats . A review by Potter et al. (2013) showed that the 68 grid lengths of the environmental data used in previous studies were approximately 10,000 69 times larger than the body length of the focal animal species and 1,000 times in plants, 70 suggesting a strong spatial mismatch between the environmental data and the actual conditions 71 experienced by species. While some researchers argue that species distribution models aim to 72 model population-level rather than individual-level responses (Bennie et al., 2014), the issue of 73 scale mismatch in distribution modelling is still recognized as important when considering the 74 population dynamics for a wide range of taxa, such as small mammals, herpetofauna, and plants 75 (Nadeau et al., 2017) . Therefore, endeavors to develop SDMs that incorporate environmental 76 conditions at biologically relevant scales are warranted (Lembrechts et al., 2019). 77 The spatial resolution of SDMs is of concern particularly in rugged terrains where coarse-78 scale climate data cannot capture appropriately the fine-scale environmental conditions created 79 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 5 by complex topographies (Austin & Van Niel, 2011). A growing body of literature has explored 80 the effect of resolution and spatial scale on SDMs (Franklin et al., 2013; Gillingham et al., 2012; 81 Guisan et al., 2007; Maclean & Early, 2023; Moudrý et al., 2023; Patiño et al., 2023; Randin et 82 al., 2009). Recent studies are further attempting to improve the performance and prediction 83 accuracy of fine -scale models by shifting from simply using downscaled climate data to 84 incorporating microclimates that are predicted based on in-situ thermal sensor measurements 85 (Haesen et al., 2023; Stark & Fridley, 2022) and remote sensing techniques (Stickley & 86 Fraterrigo, 2023). These studies show that incorporation of fine-scale conditions can improve 87 model performance and alter predictions of suitable areas, thus affecting the trajectory of 88 potential range shifts over time (e.g., Maclean & Early, 2023) . Therefore, relying solely on 89 coarse-scale models is problematic for conservation practice under climate change because it 90 can lead to poor decisions that overlook suitable areas when designing protected areas and 91 allocating resources for management. Further research is therefore required in this field because 92 a wide variety of environmental features create distinct climates at various spatial scales 93 (Morelli et al., 2017), and how species utilize and interact with them is system-specific (Shi et 94 al., 2016). 95 One of the perspectives still unexplored in this field , particularly for animals , is that a 96 given species can experience various fine -scale thermal conditions even within a single grid 97 cell as vertical structures also generate heterogenous conditions (e.g., Scheffers et al., 2014) . 98 Most previous studies have only considered the condition of a single vertical position, 99 overlooking the fine-grain heterogeneity in environmental condition s that likely influences 100 overall habitat suitability of a given location that an individual occupies . This is particularly 101 important in the context of species adaptation to climate change since species may mitigate the 102 negative impacts by altering their behavior effectively to exploit local heterogeneity in climatic 103 conditions (Beever, O’Leary, et al., 2016), which may cascade to have broader implications for 104 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 6 patterns in the overall distribution (Marske et al., 2023). Therefore, predicting habitat suitability 105 considering the variety of fine-scale conditions should improve the assessment of climate-106 induced species range shifts. 107 Pikas are small lagomorphs adapted to cold environments in montane ecosystems . 108 Various aspects of temperature, such as magnitude and variation, often act as strong drivers of 109 their distribution (Beever et al., 2010; Rodhouse et al., 2017) and abundance (Moyer-Horner et 110 al., 2016). Mounting evidence also suggest their vulnerability to recent climate change, with 111 field observations reporting local extinctions (Erb et al., 2011; Stewart et al., 2015) and range 112 contractions (Beever, Perrine, et al., 2016; Billman et al., 2021) . Other studies have predicted 113 future range contractions (Bhattacharyya et al., 2019; Schwalm et al., 2016). However, despite 114 occurring in complex t errain, previous studies modeling pika distribution have only used 115 coarse-scale climate data, possibly introducing bias into the habitat suitability predictions and 116 associated range shift estimates. 117 Some pika species belonging to the Pika subgenus dwell in rocky landforms such as talus 118 slopes, lava flows, and moraines. At the fine-scale, these rock-dwelling pikas experience duality 119 in thermal conditions as they utilize both the above-ground and sub-ground microhabitats. The 120 above-ground space is u sed for moving, vocalizing, foraging and collecting winter diet (i.e. 121 haypile), and pikas experience the local ambient air conditions during th ese activities. In 122 contrast, the sub-ground space in the rock interstices is used as the shelter (Gliwicz et al., 2005), 123 which provide pikas with cool and stable microclimates buffered from the ambient air 124 temperatures (Sakiyama et al., 2021; Varner & Dearing, 2014) . While both microhabitats are 125 essential for pikas, observations show that pikas respond flexibly to avoid the heat stress from 126 the local ambient air by utilizing rock interstices (Smith et al., 2016). Therefore, rock interstices 127 are considered a key f eature for their distribution by allowing the species to occur in 128 counterintuitive, low-elevation areas and for them to persist under climate change (Smith, 2020). 129 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 7 Predicting the habitat suitability at fine-scales for both the above- and sub-ground space should 130 thus facilitate a more accurate assessment of their vulnerability to climate change while also 131 taking their adaptive capacity into account. 132 Here, we incorporate fine-scale thermal conditions in rocky landforms into SDMs to 133 predict the distribution of the northe rn pika Ochotona hyperborea in Hokkaido , Japan. We 134 evaluate how this incorporation changes model performance, suitable area prediction, and range 135 shift forecast from the conventional coarse-scale approach. In particular, we use macroclimate 136 data (30 arcsecs) and downscaled macroclimate data (3 arcsecs) to build the coarse-scale model 137 (hereafter macroclimate model) and the fine-scale model ( local-climate model), respectively, 138 and compare their performances and habitat suitability predictions across the whole distribution. 139 While the local-climate model considers the effects of fine-scale, above-ground temperature, 140 we further assessed the habitat suitability for the sub -ground space by incorporating rock-141 interstice microclimates into the predictions, which we modeled based on local thermal 142 measurements (Sakiyama & García Molinos, 2024). We then projected future range shifts for 143 the species and compared the area of suitable habitat between local-climate and microclimate-144 based predictions. We hypothesized that (1) the local -climate model will have higher model 145 performance than the macroclimate model ; (2) the local -climate model will predict habitat 146 suitability distinctively from the macroclimate model reflecting the complex topography ; (3) 147 the rock-interstice microclimate enhances the habitat suitability from the local-climate model 148 prediction; and (4) the rock-interstices retain the suitable habitats in future predictions. 149

150 Study area and species distribution data 151 We modeled the distribution of northern pika populations in central Hokkaido, Japan, 152 because this is the core distribution area of the species in Japan and detailed fine-scale habitat 153 thermal measurements have been conducted in this area (Sakiyama et al., 2021; Sakiyama & 154 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 8 García Molinos, 2024). The northern pika in Hokkaido is located at the southern margin of its 155 distribution and therefore its population is likely among the most vulnerable of the species to 156 climate warming. Although northerly neighboring populations are found across the Sea of Japan 157 in Sakhalin and Primorsky , the Hokkaido population is essentially a disjunct population . As 158 such, it is likely to respond to climate change by shifting upslope, but previous studies have not 159 predicted future distributions. 160 To model the distribution of northern pikas, we compiled occurrence data in central 161 Hokkaido from previous literature. Given almost all the compiled data comprised presence (but 162 not absence) locations, we decided to apply a presence-only algorithm for distribution modeling 163 (see Model Development section). We only used presence locations reported on high resolution 164 maps (i.e. 1/25,000 and 1/50,000 scale) for analysis since we were interested in modeling at 165 fine scales. Overall, a total of 305 presence points encompassing a wide elevational gradient 166 (50–2,210 m) were obtained from the literature (Table S1; Fig. S1a ). Survey years spanned 167 between 1963 and 2009 (Fig. S1a, b). Detection signs of northern pika presence were clearly 168 mentioned for almost all points (298 points, 97.7%), with the vocalization representing the main 169

enabling detection of 241 points (79.0%). This suggests that the detection methods were 170 relatively consistent among different literature consulted. 171 As calibration extent for this study, we considered the area corresponding to the known 172 geographical range of the species comprising the continuous mountainous areas in central 173 Hokkaido delimited by the flatlands (i.e., river valleys) as the extent borders. Specifically, we 174 considered the Abira and Ishikari River systems as the western, the Teshio River system as the 175 northwestern, the Shokotsu River system and the Sea of Okhotsk as the northeastern, the 176 Abashiri and the Tokachi River systems as the eastern, and the Pacific Ocean as the southern 177 borders of the extent (Fig. S1a). Next, we removed irrelevant land cover types from within the 178 study extent area such as urban land, agricultural fields, and aquatic area (e.g., lakes) since 179 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 9 northern pikas require pristine montane conditions for inhabitation. The resulting filtered study 180 extent area covered approximately 21,000 km2 in central Hokkaido (Fig. S1a). 181 Model settings 182 To understand how difference s in spatial resolution s of climate data affect model 183 performances and predictions of species distribution, we constructed macroclimate ( coarse-184 scale) and local-climate (fine-scale) models (Fig. 1; Table S2). In the macroclimate model, we 185 used environmental predictors at a resolution of 30 arcsecs, which converts to a grid cell of 691 186 m by 926 m (longitude x latitude) and size of 639,682 m 2. We consider this setting to be a 187 coarse resolution for the northern pika because the ir home ran ge is much smaller, typically 188 within the range 1,822 –11,530 m 2 (n = 4; Onoyama et al., 1991) . We also note that, the 189 mountainous terrains of our study area are significant ly smoothed at this resolution, with 190 median and maximum elevational differences within one 30-arcsec grid cell as large as 182 m 191 and 770 m, respectively, when computed based on 10-times finer (3 arcsecs) elevation data (Fig. 192 S2). Since species occurrence-environment relationships might be confounded by such a large 193 degree of simplification, we developed our local climate model using environmental data at 3 194 arcsecs, which converts to a grid cell of 69.1 m by 92.6 m and size of 6,397 m 2. We consider 195 this setting to reflect the local ambient air conditions within individual home ranges. Based on 196 the local-climate model, we then generated two habitat suitability predictions to account for the 197 variety of fine -scale thermal conditions . The fir st corresponds to the model raw predictions 198 based on local-climate data and represents the above-ground habitat suitability. The second 199 prediction incorporates information on rock-interstice microclimates to further predict sub-200 ground habitat suitability. 201 Topographical and geological variables 202 We used the terrain slope angle as topographical predictor in the SDMs because it has 203 been shown to have a positive relationship with the northern pika distribution (Sakiyama et al., 204 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 10 2021). We computed the slope angle raster in QGIS 3.34 (QGIS Development Team, 2023) 205 using the MERIT DEM (Yamazaki et al., 2017), a digital elevation model provided at a spatial 206 resolution of 3 arcsecs. Slope angles for the macroclimate setting were calculated using a DEM 207 aggregated to 30 arcsecs (Table S2). We also considered surface geology as a categorical binary 208 predictor because 99 % of the known distribution of the northern pika lies on igneous and 209 metamorphic rocks and moraines (Kawabe, 1992). To do so, we assigned whether any of these 210 geological types occurred per grid cell with in the study area based on the Seamless Digital 211 Geological Map of Japan (1:200,000) V2 (Geological Survey of Japan, National Institute of 212 Advanced Industrial Science and Technology, 2017) (Table S2). 213 Bioclimatic variables 214 To model and predict the distribution of northern pikas, w e computed the mean daily 215 temperature and mean diurnal range over July and August to account for the effects of thermal 216 magnitude and variation during the summer. We focused on summer conditions (i.e., July and 217 August) because the effects of climate change on cold-adapted species are likely to be 218 particularly strong during this season. Briefly (see Appendix S1 for a detailed explanation), we 219 first obtained climate data at 30 arcsecs using the chelsa-cmip6 R package, which provides 220 access to downscaled Global circulation model (GCM) projections (Karger et al., 2023). Four 221 GCMs that are known to reproduce well the climatic conditions of our study area were selected 222 (i.e. ACCESS-CM2, IPSL-CM6A-LR, MIROC6, and MRI-ESM2-0) (Shiogama et al., 2021). 223 We used climate data for 1981–2010 to model the baseline distribution. To create the 224 bioclimatic variables, we adapted the computation formulas used in ANUCLIM , a software 225 package used to build spatial climatic data (Xu & Hutchinson, 2011, 2013) , to the summer 226 season, and calculated mean summer temperature and mean summer diurnal range. The 227 resulting variables were then used for the macroclimate setting (Table S2). 228 For the local-climate setting, we downscaled the climate data obtained at 30 arcsecs to 3 229 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 11 arcsecs as follows. Temperature data was downscaled using kriging with the krigR R package 230 (Davy & Kusch, 2021; Kusch & Davy, 2022) . This approach first models the relationship 231 between temperature and elevation at the coarse scale (30 arcsecs), which is then used to 232 interpolate the fine-scale temperature using elevation at the same target scale (3 arcsecs). We 233 used the raw MERIT DEM (Yamazaki et al., 2017) dataset, provided at 3 arcsecs, as the 234 elevation data for the fine-scale interpolation and the aggregated version at 30 arcsecs to model 235 the coarse-scale relationship (Table S2). 236 Finally, for the rock-interstice microclimate, we considered that northern pikas can utilize 237 sub-ground thermal conditions in rocky landforms, which are known to be buffered from the 238 ambient air temperature (Sakiyama and García Molinos, 2024). However, given the absence of 239 sufficient spatial distribution information regarding rocky landforms throughout our study 240 extent, we selected three representative regions where we had sufficient data separated by 241 elevational ranges: Biei (1,100–2,000 m), Shikaribetsu (600–1,250 m), and Oketo (400–800 m) 242 (Fig. 4). For Biei, we were able to map the rocky landforms directly using aerial images. For 243 Shikaribetsu and Oketo, we referred to the mapped occurrence of rocky landforms in previous 244 literature (Kojima, 2004; Kurumada, 201 0). For these regions, we then used an empirical 245 relationship between ambient and interstice temperatures previously derived from northern pika 246 taluses in our study area (Sakiyama and Garcia Molinos, 2024) to convert downscaled ambient 247 temperature into interstice temperature for each of the 3 arc second grid cells containing rocky 248 landforms (Fig. S3; Appendix S1). Here, our assumption is that northern pikas experience local 249 ambient climatic conditions in cells devoid of rocks since these landforms are patchily 250 distributed within each region. We validated our downscaling and microclimate incorporation 251 approach by comparing the thermal conditions measured in northern pika habitats with the 252 values of each bioclimatic variable in the grid cells where those observations had been made 253 (see Appendix S2). Results confirmed that the local-climate data was moderately accurate in 254 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 12 reflecting observed ambient temperature and that the microclimate data was highly accurate in 255 reflecting observed rock interstice temperatures (Fig. S4, 5). 256 Model development 257 We performed the distribution modeling analysis of northern pikas in Hokkaido using 258 maximum entropy (MaxEnt; Phillips et al., 2006) . Models were developed separately for the 259 macroclimate and the local -climate settings considering the years 1981–2010 as the baseline 260 period, thus using the presence points during this period for model fitting. Prior to model fitting, 261 we reduced the effect of spatial sampling bias by thinning the presence data using the spThin R 262 package (Aiello-Lammens et al., 2015) at a distance of 500 m as presence locations could be 263 viewed as independent at this distance considering the complex mountainous terrain and the 264 limited mobility of the northern pika . Considering that the diagonal distance of macroclimate 265 grid cells is larger than this thinning distance and duplicate sampling thus could occur, we 266 randomly selected a single point for each cell after thinning. This pre-selection process resulted 267 in a total of 80 presence points remaining. We generated 10,000 background points throughout 268 the study area, which are used in MaxEnt for characterizing the extant environmental conditions 269 within the study area. To avoid multicollinearity issues (Dormann et al., 2013) , we computed 270 the pairwise Pearson’s correlation coefficients among the continuous predictors and confirmed 271 that all pairs did not have a strong correlation (|r| > 0.7) (Fig. S6). 272 Since MaxEnt models with default settings can lead to overfitting (Anderson & Gonzalez, 273 2011), we initially tuned the model hyperparameters to control the complexity while fitting well 274 to the data. To do this, we used various hyperparameter combinations of feature classes, which 275 determines the shape of the response curves, and regularization multipliers, which controls the 276 model complexity (Phillips et al., 2006; Phillips & Dudík, 2008), and selected the combination 277 that resulted in the highest predicative accuracy. For the feature class, we considered linear (L), 278 quadratic (Q), product (P), and/or hinge (H) characteristics, and f or the regularization 279 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 13 multipliers, we used values ranging from 0.2 to 2.0 with intervals of 0.2. The predictive 280 accuracy was assessed for each model by separating the presence data into training (80%) and 281 test (20%) data and computing the area under the receiver operating characteristic curve (AUC), 282 which indicates the accuracy condition on a range of 0 to 1 (Fielding & Bell, 1997). All tuning 283 processes were conducted using the optimizeModel function in the SDMtune R package 284 (Vignali et al., 2020) . After hyperparameter tuning, we evaluated the optimized model’s 285 predictive accuracy based on 5-fold cross validation using Continuous Boyce Index (Hirzel et 286 al., 2006), and compared these values among the model settings. 287 Comparing baseline predictions 288 For both the macroclimate and the local -climate m odels, w e predicted the habitat 289 suitability for the baseline period across the study area and compared the predicted suitability 290 by computing the difference. We also generated binarized maps of predicted species occurrence 291 (i.e. presence or absence) based on the habitat suitability prediction across the whole study area. 292 Since discriminant thresholds were likely to differ between the optimized models for the 293 macroclimate and local -climate and could be problematic for comparison, we instead used a 294 sensitivity analysis approach by applying a moving threshold from 0 to 1 at 0.01 increments to 295 discriminate presence or absence. Binary predictions were then compared between the 296 macroclimate and local-climate models by assessing the proportion of predicted presence areas 297 in the study area, calculated as the total amount of presence cells divided by the number of all 298 cells, for each moving threshold. We used p roportion instead of simply calculating the total 299 area of presence because the total study area was subtly different between the scales due to 300 processing procedures of environmental data. We then predicted the habitat suitability for the 301 sub-ground space for the three selected regions by replacing the local ambient temperatures by 302 the modeled rock interstice temperatures for the prediction of habitat suitability using the local-303 climate model . We consider this prediction to reflect the habitat suitability of the rock -304 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 14 interstices within individual home ranges defined based on the relationship between 305 occurrences and ambient air temperatures. The resultant prediction was compared with those 306 based on the macroclimate and local-climate. 307 Hindcasting 308 We use a hindcasting approach t o evaluate the transferability of the models to project 309 distributions over time (see Nogués‐Bravo, 2009) by projecting the distribution of the species 310 to the past and comparing it to known past presence points from the same past period for 311 validation (Hodel et al., 2022). Although this approach is still rarely used in SDM studies (Yates 312 et al., 2018) , it allowed us to assess how confident we can be using the models to predict 313 distributions for a (past or future) time period different to that the model was calibrated for. We 314 projected the distribution of northern pikas for years 1961–1980, a period for which early 315 presence records are available (Table S1; Fig S1). The same procedures as the baseline climate 316 data were used for processing past data. We evaluated the model’s predicative accuracy via the 317 Continuous Boyce Index (Hirzel et al., 2006). 318 Future projection 319 To project the future distribution of the northern pika, we obtained climate data for two 320 future time periods, 2041–2070 and 2071–2100, considering two emission scenarios, SSP1-2.6 321 and SSP5-8.5, representing a globally sustainable and a fossil-fuel based development scenario, 322 respectively. We predicted the future distribution for the study extent based on the macroclimate 323 and local-climate models using the future bioclimatic data, which were processed based on the 324 same procedures for the baseline climate da ta (i.e., we assume that future topographical and 325 geological variables will remain the same as present ). We binarized the suitability predictions 326 to presence or absence by applying the same approach used for the baseline period to calculate 327 the proportion of presence areas. 328 Moreover, we also projected the future distribution of the northern pika for the three 329 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 15 selected regions to assess how microclimate contributes to retaining suitable conditions in the 330 future. When processing the future microclimate data, we assumed that the buffering effect of 331 the rock interstices (i.e., the empirical relationship between local ambient and rock -interstice 332 temperatures) will remain invariable over time. We believe that this assumption is robust 333 because the field thermal measurements from which this relationship was derived contained a 334 large elevational gradient of 350–2200 m and the model residual was stable across the ambient 335 temperature gradient (Sakiyama and Garcia Molinos, 2024) , although we note that some 336 extrapolation did occur for some lower elevation areas where future conditions exceeded the 337 observed thermal gradient. We binarized the suitability predictions for both local-climate and 338 microclimate using the maximum true skill statistic computed for the local-climate model. The 339 consensus among the resulting binarized predictions from each of the four GCMs was then 340 computed by considering a cell as present whenever that cell was predicted as presen t by two 341 or more individual GCMs (i.e., ≥50% consensus) . The consensus predictions in the selected 342 regions were compared visually between the local-climate and microclimate-based predictions 343 for each period and emission scenarios. 344

345 Model performance and variable importance 346 Both the macroclimate and local-climate models converged to the same hyperparameter 347 combinations for the best -fit model comprising an LQP feature class and a regularization 348 parameter value of 2. Using these hyperparameters, the Continuous Boyce Index was 0.73 and 349 0.80 for the macroclimate and local-climate models, respectively (Table 1). This indicates that 350 the continuous predictions reflected the species distribution well and that the local-climate 351 model was more accurate than the macroclimate model. 352 Mean summer temperature was the most important predictor explaining the northern pika 353 distribution in both models with a permutation importance of 74.5 % and 76.2 % for the 354 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 16 macroclimate and local-climate models, respectively (Fig. S 7). Mean summer temperature 355 decreased occurrence probability in both models, indicating that the northern pika is more likely 356 to be present in cooler areas. However, we detected a slightly steeper slope in the response 357 curve of the local-climate model, which indicates that the local -climate model fitted towards 358 cooler conditions than the macroclimate model (Fig. S8). In both models, surface geology was 359 the second most important predictor with 20% and 18 % of permutation importance , 360 respectively (Fig. S7), and a positive relationship with occurrence probability (Fig. S8). The 361 permutation contribution of daily thermal range and slope was less than 5% in both models (Fig. 362 S7). 363 Habitat suitability for the baseline period 364 The negative effect of mean summer temperature was reflected in the habitat suitability 365 predictions at the whole distribution scale in both models, with higher elevations predicted as 366 highly suitable areas and the suitability lowering gradually towards downslope (Fig. 2 a,b). 367 Specifically, large mountain ranges such as the Daisetsuzan and Hidaka Mountain Ranges , 368 which contain summits over 2 ,000 m , as well as mountains in the Eastern Daisetsu region 369 (elevation of the main peaks 1,800–2,000 m) were recognized as the most suitable areas. Mid- 370 to high -elevation mountain ranges, such as the Yubari-Ashibetsu (1,600–1,720 m) , Teshio 371 (1,400–1,600 m), Ashoro (1,200–1,300 m), and Shikaribestu Mountains (1,100–1,250 m), were 372 predicted as areas with moderate to high suitability . Finally, low-elevation mountains were 373 predicted to be low in habitat suitability for northern pikas in almost all regions although some 374 presence points were reported from such areas (e.g. central and northeastern region s). One 375 exception was the Samani-Erimo region, which was predicted to have high suitability despite 376 the region being predominated by low elevations (Fig. 2a,b). 377 As we compared the macroclimate-based and local-climate based predictions, the former 378 was generally higher in habitat suitability than the latter across the study area (Fig. 2c). The 379 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 17 difference in suitability was relatively large in the mountainous areas, particularly in the slopes 380 rather than the summits, with suitability differences up to 0.5. However, habitat suitability from 381 the local-climate based were also higher than the macroclimate -based predictions in some 382 mountainous areas as well as the Samani-Erimo region. The habitat suitability difference was 383 relatively small in low elevations. 384 The proportion of predicted presence areas for both models indicated similar curves along 385 the moving binarization threshold (Fig. 3a). The curve inverse shape means that, while the 386 whole study area is discriminated as present when the threshold is zero, a large proportion 387 becomes absent at relatively low thresholds, although a small proportion of the study area 388 remains present until the threshold becomes o ne. In other words, while most of the study area 389 is predicted to have low suitability, some areas have distinctively high values. Comparing the 390 curves between the models, the proportion of presence areas was higher in the macroclimate-391 based prediction than in the local -climate based prediction across the whole range of 392 binarization threshold (Fig. 3a). This suggests that the macroclimate model predicted a larger 393 area of suitable habitat than the local-climate model regardless of the binarization threshold. 394 In the three selected regions, the habitat suitability based on the macroclimate and local-395 climate showed a gradually changing pattern , likely reflecting the elevation difference. In 396 contrast, accounting for the rock interstice microclimate in local-climate predictions boosted 397 the habitat suitability distinctively within these regions (Fig. 4). For instance, while the habitat 398 suitability in Biei was relatively high given its high elevation in both the macroclimate - and 399 local-climate-based predictions, the rocky landforms further increased the suitability to mostly 400 reaching a suitability of 1 (Fig. 4 ). This trend was more evident in Shikaribetsu , the mid -401 elevation region, where the predicted suitability reached as high as 0.99 in the microclimate -402 based predictions despite the surrounding landscape having suitability less than 0.5 (Fig. 4). In 403 Oketo, the low-elevation region, the total area that increased in suitability was relatively small 404 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 18 given the limited area of rocky landforms, but the predicted suitability within rocky landforms 405 reached as high as 0.66 (Fig. 4). The maximum values of predicted suitability based on the 406 macroclimate and local-climate models were 0.55 and 0.41 in Shikaribetsu and 0.24 and 0.25 407 in Oketo, respectively. In all these three regional-scale predictions, the habitat suitability was 408 predicted in a more detailed and complex manner in the fine-scale (3 arc sec) predictions than 409 that in the coarse -scale (30 arc sec) predictions as evidenced by the variation in habitat 410 suitability of the local-climate and microclimate predict ions within the coarse-scale grid cells 411 (Fig. 4). 412 Hindcasting and future projection 413 In the hindcasting analysis, both macroclimate and local-climate models had Continuous 414 Boyce Index values well above zero (Table 1), suggesting that the models were a ccurate in 415 predicting past distributions. This result supports the transferability of these models to predict 416 distributional change over time . Accordingly, we were able to predict future distributions for 417 both model settings with confidence. 418 Both models indicated inverse curves in the proportion of presence areas in the future 419 (Fig. 3b). Comparing with the baseline proportions, future distribution was predicted to 420 decrease by more than 60% in in both time periods and emission scenarios considered (Fig S9). 421 Between the two future time periods, the difference in the habitat decrease rate was relatively 422 small under the SSP1-2.6 scenario, suggesting that presence areas will decrease until 2041–423 2070 but will then remain relatively stable until 2071–2100. However, the habitat decrease was 424 apparent under the SSP5-8.5 scenario, suggesting that suitable habitats will continue to decrease 425 until 2071–2100. (Fig. 3b; S9). 426 Patterns of predicted distributional change largely differed among the three selected 427 regions between the local-climate an d microclimate -based predictions . Based on the local -428 climate, predicted presences were found in Biei in both future time periods and scenarios and 429 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 19 in Shikaribetsu under the SSP1-2.6 scenario, while none were found in Oketo under all future 430 settings (Fig. 5). Conversely, predicted presences were found in almost all regions in all future 431 settings when accounting for rock interstice microclimates , except in Oketo for the far future 432 under the SSP5-8.5 scenario (Fig. 5). Discrepancies between local-climate and microclimate-433 based predictions were high towards the lower elevations and under the higher carbon emission 434 scenario (SSP5-8.5). 435

436 We developed distribution model s for the northern pika that incorporate d fine-scale 437 temperatures experienced by the species to assess how model performances and predictions 438 change from a conventional model based on coarse-scale conditions . Our results show that 439 predictive accuracy was higher for the fine-scale model. We also found that the predicted 440 habitat suitability varied among regions and at smaller spatial scales in the fine-scale prediction. 441 These results confirm that accounting for fine-scale environments can lead to distinct and more 442 complex predictions than those obtained from a coarse-scale model and that solely depending 443 on coarse-scale predictions can result in biased predictions that overlook potential suitable 444 habitats. 445 There are two possible reasons that likely contributed to the higher predictive accuracy 446 in the fine-scale model. First, model response curves showed some important differences in the 447 modeled relationship between occurrence and the environmental predictors . While mean 448 temperature had a dominant negative effect on the distribution in both models (Fig. S7; S8), a 449 response also observed in a previous study (Sakiyama et al., 2021) , the response curve was 450 steeper in the local-climate model. This indicates that the model fitted towards cooler conditions, 451 which may be reflecting the actual response of the northern pika more accurate ly than th e 452 macroclimate model. Interestingly, our input data for the distribution model indicated that many 453 of the occurrence points in higher elevations were located in areas higher than its surroundings 454 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 20 and thus had cooler temperatures at the fine-scale (Fig. S10). Although we referred to presence 455 points recorded in previous literature, it remains unclear whether this trend is due to a potential 456 sampling bias in convex topographies such as ridgelines and mountain tops. We believe that 457 assessing species distribution at fine-scales could introduce its scale-dependent issues, which 458 requires attention as more fine-scale SDMs are developed (Mitchell et al., 2017). Second, the 459 quality of the predictors may have also had an effect on model accuracy . Even when coarse- 460 and fine-scale models exhibit the same relationship between occurrence and environment, the 461 validated accuracy is likely to decrease if predictors are generalized at the coarse-scale. Instead, 462 fine-scale predictors are more likely to accurately represent true environmental conditions, 463 which could have resulted in higher predictive accuracy by resolving better local complex 464 topographies. Thus, the overall increase in accuracy of the local-climate model can be attributed 465 to these small refinements in the modeled relationship and predictors. In general, our results are 466 in agreement with recent studies finding higher performance in SDMs incorporating fine-scale 467 climates (Haesen et al., 2023; Stark & Fridley, 2022). 468 Our results showed that incorporating fine-scale environments into SDMs resulted in a 469 decrease of suitable areas compared to the coarse-scale prediction. This is likely due to 470 differences in the modeled response of mean summer temperature between the models. Such 471 differences likely emerged from the local -climate model being fitted towards cooler 472 temperatures relative to the coarser macroclimate model resulting in lower predicted occurrence 473 probabilities at certain temperatures or elevations. Consequently, this leads to a greater area of 474 predicted presence in the macroclimate model along the moving threshold when binarizing the 475 habitat suitability predictions . In previous literature, fine-scale models have resulted in both 476 greater (Stickley & Fraterrigo, 2023) and smaller suitable areas (Franklin et al., 2013; Haesen 477 et al., 2023) . Hence, the difference in the predicted habitat area resulting from incorporating 478 fine-scale environments into S DMs may be s pecific to each study system and the scale -479 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 21 dependent relationship between the occurrence records and environmental predictors used in 480 the model . Furthermore, in our study, both models projected drastic future declines in the 481 northern pika distribution, with over half of the suitable habitats predicted to disappear by mid-482 century. This suggests that, similar to other pika species, northern pikas in Hokkaido are 483 vulnerable to climate change (Beever et al., 2011; Bhattacharyya et al., 2019). Importantly, the 484 predicted change in species distribution depended strongly on the future period and emission 485 scenario. This result stresses the need to enforce effective policies for greenhouse gas emission 486 mitigation and to guide effective conservation of the species by protecting areas that will 487 account for the effect of climate change and maintain its distribution in the future. 488 At the regional scale, the incorporation of the rock-interstice microclimate provided 489 distinct predictions for both the baseline and future. Overall, the cooler conditions provided by 490 rock-interstices resulted in the microclimate predictions markedly enhancing habitat suitability 491 compared to t he local-climate based prediction in all three regions. This demonstrates that a 492 single location (grid cell) can exhibit dual habitat suitability due to the existence of diverse 493 thermal conditions in complex topograph ies. However , what are the implications for the 494 northern pika distribution? In rocky landforms, the habitat suitability can simply be categorized 495 into three combinations based on the local -climate and microclimate conditions: high -high, 496 low-high, and low-low. The combination of high-low suitability was not observed because rock 497 interstices were always cooler than local ambient conditions , hence more suitable. Northern 498 pikas are expected to be present in areas with high-high suitability, and absent from areas with 499 low-low suitability . However, whether areas with low-high suitability serve as functional 500 habitats for northern pikas will likely depend on the extent to which the rock interstice condition 501 mitigates the impact of poor thermal ambient conditions, which should be strongly linked to the 502 capacity of northern pika s to adapt their thermoregulation behavior. Previous studies have 503 indicated that northern pika s utilize above-ground space during the cooler times of the day 504 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 22 (Onoyama, 1 991); a behavior also observed in American pika stud ies (Smith et al., 2016) . 505 Additionally, research on the American pika suggests that it regulates body temperature by 506 reducing above-ground activities and retreating into rock interstices (MacArthur & Wang, 507 1974). Thus, the northern pika may exhibit high behavioral flexibility to respond to unfavorable 508 ambient conditions, enabling them to inhabit areas with low-high suitability. Consequently, the 509 actual overall suitability of a given rocky landform should be influenced by a more complex 510 combination of these two suitability conditions (above and below ground). 511 Despite its simplicity, our approach to predict microclimate-based habitat suitability is 512 innovative because it enables detecting potential suitable habitats that are likely to be 513 overlooked in predictions based solely on ambient conditions, especially in the lower elevations. 514 Therefore, this approach can be used as a first -cut approach to identify potential priority areas 515 for conservation. Assessing occupancy patterns in areas with low-high suitability will be 516 therefore important for resource allocation , although the large elevational range of observed 517 presence of the northern pika, despite their cold adaptation, already suggests they are inhabiting 518 such counterintuitive areas perhaps through behavioral adaptation. Moreover, this approach has 519 implications for identifying climate change refugia for the northern pika as all three regions 520 were predicted to contain areas supported by rock interstice microclimate in the future despite 521 the deterioration of local-climate conditions. Areas with favorable microclimates that support 522 species persistence under a changing climate (i.e., microrefugia) are often cryptic (Provan & 523 Bennett, 2008), and therefore, extensive identification of such areas based on spatial data should 524 facilitate effective protection of vulnerable populations. In our case, however, the current lack 525 of an extensive dataset reporting presence of rocky landforms in our study area limits such 526 endeavors. Therefore, further refinements in geospatial data are warranted for predicting 527 potential microrefugia. Deeper understanding of the behavioral responses of the northern pika 528 to ambient and rock interstice thermal conditions will also be important to assess the reliability 529 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 23 of behavioral adaptation in driving species distribution patterns (Marske et al., 2023). 530 The importance of f ine-scale SDMs are being increasingly recognized as they often 531 provide distinct predictions from conventional scales and information in fine-scale predictions 532 are essential for improving downstream objectives, such as designing effective conservation 533 plans under climate change (Nadeau et al., 2022) . In this study, the local-climate model 534 prediction proved to be valuable for achieving such objectives as it provided distinct habitat 535 suitability predictions. Additionally, an important implication of fine-scale SDMs is their role 536 in filling the gap between predictions and observations of ongoing species range shifts (Chen 537 et al., 2011). This is because conventional approaches have been criticized for overestimating 538 range shifts due to the lack of consideration for the fine-scale climates that may facilitate species 539 persistence in warmer regions (Pearson, 2006). Therefore, it is important to validate whether 540 fine-scale models effectively reduce d overestimation by conducting ground-truthing surveys. 541 A study by Maclean and Early (2023) exemplifies this, demonstrating that fine-scale predictions 542 reduced the number of species at extinction risk and were in agreement with the actual observed 543 patterns, while coarse-scale predictions estimated several species to become extinct. 544 One of t he limitations of our study is that, despite using the downscaled macroclimate 545 data, we still rely on a simplified representation of time and space to develop the fine-scale 546 models. Rocky landforms exhibit high heterogeneity, with ambient air and rock interstice 547 temperatures varying depending on height and depth, respectively, and changing across seasons 548 (Benedict et al., 2020; Millar et al., 2016; Smith, 2020) . In this regard, our study may still be 549 considered as a coarse-scale representation of reality in time and space. Like other SDM studies, 550 our model relies on the correlative relationships between presences and climatic and 551 topographical predictors. Studies have demonstrated that predicted distributions are altered by 552 incorporating biotic interactions (Bateman et al., 2012), dispersal settings (Zurell et al., 2016), 553 and lineage information (Zhang et al., 2021) into the models. Therefore , further research is 554 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 24 needed to incorporate these aspects into the prediction of the distribution of the northern pika. 555 In conclusion, the incorporation of fine-scale climate data into SDMs yielded contrasting 556 predictions to the coarse -scale prediction . Future distribution s were predicted to decrease 557 increasingly with both further away time horizons and higher emission scenarios considered , 558 stressing the impo rtance of enforcing effective mitigation measures to curve greenhouse gas 559 emissions. Together, this information is essential for the conservation of northern pika s in 560 Hokkaido, as previous studies have not assessed such impacts. Incorporation of rock-interstice 561 microclimate data significantly increased the habitat suitability and predicted more suitable 562 habitats for both the baseline and future. Behavioral adaptation is considered crucial for the 563 inhabitance of the northern pikas at lower elevations, which is important for the conservation 564 of the northern pika because its distribution is often considered to be restricted to the alpine 565 zone. Based on these results, we advocate for the identification of rocky landforms and their 566 integration into conservation planning to effectively leverage microclimates in conservation 567 efforts. 568 569 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 25

570 This study was supported by the JSPS KAKENHI Grant Number JP22J10191, grant-in-aid of 571 The Zoshinkai Fund for Protection of Endangered Animals, and grant-in-aid of The Inui 572 Memorial Trust for Research on Animal Science, Japan. 573 DECLARATION OF COMPETING INTEREST 574 The authors have no competing interests to declare. 575 CRediT AUTHORSHIP CONTRIBUTION STATEMENT 576 TS contributed to conceptualization and design o f methodology, data analysis, funding 577 acquisition, and writing the original draft; TS and JGM contributed to discussion of results, 578 writing the final manuscript and gave final approval for publication. 579 DATA AVAILABILITY 580 Data will be made available on request. 581 SUPPLEMENTAL MATERIAL 582 The following materials are available as part of the online article at 583 https://escholarship.org/uc/fb. 584 Appendix S1. Derivation of the bioclimatic variables. 585 Appendix S2. Validation of the bioclimatic variables. 586 Table S1. List of previous literature reporting the northern pika distribution in Hokkaido. 587 Table S2. Procedures for generating the environmental variables. 588 Fig. S1. Geographical map of the presence points used in this study. 589 Fig. S2. A boxplot indicating the elevational difference within macroclimate grid cells. 590 Fig. S3. The modeled relationship between mean ambient and rock interstice temperatures. 591 Fig. S4. Scatterplots used for validating the bioclimatic variables. 592 Fig. S5. Validation plots of the bioclimatic variables used in the study. 593 Fig. S6. Result of the correlation analysis. 594 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 26 Fig. S7. Permutation importance of the model predictors. 595 Fig. S8. Response plots of the model predictors. 596 Fig. S9. Percentage of presence areas predicted to change in the future from the baseline period. 597 Fig. S10. Elevational difference between the macroclimate and local-scale cells at the presence 598 points. 599 600 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 27

601 Aiello-Lammens, M. E., Boria, R. A., Radosavljevic, A., Vilela, B., & Anderson, R. P. (2015). 602 spThin: An R package for spatial thinning of species occurrence records for use in 603 ecological niche models. Ecography, 38(5), 541 –545. 604 https://doi.org/10.1111/ecog.01132 605 Anderson, R. P., & Gonzalez, I. (2011). Species -specific tuning increases robustness to 606 sampling bias in models of species distributions: An implementation with Maxent. 607 Ecological Modelling , 222(15), 2796 –2811. 608 https://doi.org/10.1016/j.ecolmodel.2011.04.011 609 Austin, M. P., & Van Niel, K. P. (2011). Improving species distribution models for climate 610 change studies: Variable selection and scale. Journal of Biogeography , 38(1), 1–8. 611 https://doi.org/10.1111/j.1365-2699.2010.02416.x 612 Bateman, B. L., VanDerWal, J., Williams, S. E., & Johnson, C. N. (2012). Biotic interactions 613 influence the projected distribution of a specialist mammal under climate change. 614 Diversity and Distributions , 18(9), 861 –872. https://doi.org/10.1111/j.1472 -615 4642.2012.00922.x 616 Beever, E. A., O’Leary, J., Mengelt, C., West, J. M., Julius, S., Green, N., Magness, D., Petes, 617 L., Stein, B., Nicotra, A. B., Hellmann, J. J., Robertson, A. L., Staudinger, M. D., 618 Rosenberg, A. A., Babij, E., Brennan, J., Schuurman, G. W., & Hofmann, G. E. (2016). 619 Improving Conservation Outcomes with a New Paradigm for Understanding Species’ 620 Fundamental and Realized Adaptive Capacity: A new paradigm for defining adaptive 621 capacity. Conservation Letters, 9(2), 131–137. https://doi.org/10.1111/conl.12190 622 Beever, E. A., Perrine, J. D., Rickman, T., Flores, M., Clark, J. P., Waters, C., Weber, S. S., 623 Yardley, B., Thoma, D., Chesley -Preston, T., Goehring, K. E., Magnuson, M., 624 Nordensten, N., Nelson, M., & Collins, G. H. (2016). Pika ( Ochotona princeps ) losses 625 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 28 from two isolated regions reflect temperature and water balance, but reflect habitat area 626 in a mainland region. Journal of Mammalogy , 97(6), 1495 –1511. 627 https://doi.org/10.1093/jmammal/gyw128 628 Beever, E. A., Ray, C., Mote, P. W., & Wilkening, J. L. (2010). Testing alternative models of 629 climate-mediated extirpations. Ecological Applications , 20(1), 164 –178. 630 https://doi.org/10.1890/08-1011.1 631 Beever, E. A., Ray, C., Wilkening, J. L., Brussard, P. F., & Mote, P. W. (2011). Contemporary 632 climate change alters the pace and drivers of extinction. Global Change Biology, 17(6), 633 2054–2070. https://doi.org/10.1111/j.1365-2486.2010.02389.x 634 Benedict, L. M., Wiebe, M., Plichta, M., Batts, H., Johnson, J., Mon k, E., & Ray, C. (2020). 635 Microclimate and Summer Surface Activity in the American Pika (Ochotona princeps). 636 Western North American Naturalist, 80(3). https://doi.org/10.3398/064.080.0303 637 Bennie, J., Wilson, R. J., Maclean, I. M. D., & Suggitt, A. J. (2014). Seeing the woods for the 638 trees – when is microclimate important in species distribution models? Global Change 639 Biology, 20(9), 2699–2700. https://doi.org/10.1111/gcb.12525 640 Bhattacharyya, S., Mungi, N. A., Kawamichi, T., Rawat, G. S., Adhikari, B. S., & Wilkening, 641 J. L. (2019). Insights from present distribution of an alpine mammal Royle’s pika 642 (Ochotona roylei) to predict future climate change impacts in the Himalaya. Regional 643 Environmental Change, 19(8), 2423–2435. https://doi.org/10.1007/s10113-019-01556-644 x 645 Billman, P. D., Beever, E. A., McWethy, D. B., Thurman, L. L., & Wilson, K. C. (2021). 646 Factors influencing distributional shifts and abundance at the range core of a climate‐647 sensitive mammal. Global Change Biology , 27(19), 4498 –4515. 648 https://doi.org/10.1111/gcb.15793 649 Chen, I.-C., Hill, J. K., Ohlemüller, R., Roy, D. B., & Thomas, C. D. (2011). Rapid range shifts 650 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 29 of species associated with high levels of climate warming. Science, 333(6045), 1024–651 1026. https://doi.org/10.1126/science.1206432 652 Davy, R., & Kusch, E. (2021). Reconciling high resolution climate datasets using KrigR. 653 Environmental Research Letters , 16(12), 124040. https://doi.org/10.1088/1748-654 9326/ac39bf 655 Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., Marquéz, J. R. G., 656 Gruber, B., Lafourcade, B., Leitão, P. J., Münkemüller, T., McClean, C., Osborne, P. 657 E., Reineking, B., Schröder, B., Skidmore, A. K., Zurell, D., & Lautenbach, S. (2013). 658 Collinearity: A review of methods to deal with it and a simulation study evaluating their 659 performance. Ecography, 36(1), 27 –46. https://doi.org/10.1111/j.1600 -660 0587.2012.07348.x 661 Elith, J., & Leathwick, J. R. (2009). Species distribution models: Ecological explanation and 662 prediction across space and time. Annual Review of Ecology, Evolution, and Systematics, 663 40(1), 677–697. https://doi.org/10.1146/annurev.ecolsys.110308.120159 664 Erb, L. P., Ray, C., & Guralnick, R. (2011). On the generality of a climate-mediated shift in the 665 distribution of the American pika ( Ochotona princeps ). Ecology, 92(9), 1730 –1735. 666 https://doi.org/10.1890/11-0175.1 667 Fielding, A. H., & Bell, J. F. (1997). A review of methods for the assessme nt of prediction 668 errors in conservation presence/absence models. Environmental Conservation, 24(1), 669 38–49. https://doi.org/10.1017/S0376892997000088 670 Franklin, J., Davis, F. W., Ikegami, M., Syphard, A. D., Flint, L. E., Flint, A. L., & Hannah, L. 671 (2013). Modeling plant species distributions under future climates: How fine scale do 672 climate projections need to be? Global Change Biology , 19(2), 473 –483. 673 https://doi.org/10.1111/gcb.12051 674 Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology. 675 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 30 (2017). Seamless Digital Geological Map of Japan V2 1: 200,000 [dataset]. 676 https://gbank.gsj.jp/seamless/ 677 Gillingham, P. K., Huntley, B., Kunin, W. E., & Thomas, C. D. (2012). The effect of spatial 678 resolution on projected responses to cl imate warming. Diversity and Distributions , 679 18(10), 990–1000. https://doi.org/10.1111/j.1472-4642.2012.00933.x 680 Gliwicz, J., Witczuk, J., & Pagacz, S. (2005). Spatial behaviour of the rock‐dwelling pika 681 (Ochotona hyperborea ). Journal of Zoology , 267(2), 113 –120. 682 https://doi.org/10.1017/S0952836905007211 683 Guisan, A., Graham, C. H., Elith, J., Huettmann, F., & the NCEAS Species Distribution 684 Modelling Group. (2007). Sensitivity of predictive species distribution models to 685 change in grain size. Diversity and Dist ributions, 13(3), 332 –340. 686 https://doi.org/10.1111/j.1472-4642.2007.00342.x 687 Haesen, S., Lenoir, J., Gril, E., De Frenne, P., Lembrechts, J. J., Kopecký, M., Macek, M., Man, 688 M., Wild, J., & Van Meerbeek, K. (2023). Microclimate reveals the true thermal nich e 689 of forest plant species. Ecology Letters , 26(12), 2043 –2055. 690 https://doi.org/10.1111/ele.14312 691 Hirzel, A. H., Le Lay, G., Helfer, V., Randin, C., & Guisan, A. (2006). Evaluating the ability 692 of habitat suitability models to predict species presences. Ecological Modelling, 199(2), 693 142–152. https://doi.org/10.1016/j.ecolmodel.2006.05.017 694 Hodel, R. G. J., Soltis, D. E., & Soltis, P. S. (2022). Hindcast‐validated species distribution 695 models reveal future vulnerabilities of mangroves and salt marsh species. Ecology and 696 Evolution, 12(9), e9252. https://doi.org/10.1002/ece3.9252 697 Karger, D. N., Chauvier, Y., & Zimmermann, N. E. (2023). chelsa‐cmip6 1.0: A python 698 package to create high resolution bioclimatic variables based on CHELSA ver. 2.1 and 699 CMIP6 data. Ecography, e06535. https://doi.org/10.1111/ecog.06535 700 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 31 Kawabe, M. (1992). Geological factor of distribution of Northern Pika Ochotona hyperborea 701 in Hokkaido. Bulletin of the Higashi Taisetsu Museum of Natural History, 14, 103–106. 702 Kojima, N. (2004). Conservation biology of the Japanese pika Ochotona hyperborea yesoensis. 703 Iwate University. 704 Kurumada, T. (2010). The habitat patches distribution and utilization of the northern pika in 705 Mt. Nakayama, Hokkaido. Report of Hokkaido Institute of Environmental Sciences, 36, 706 65–69. 707 Kusch, E., & Davy, R. (2022). KrigR —a tool for downloading and statistically downscaling 708 climate reanalysis data. Environmental Research Letters , 17(2), 024005. 709 https://doi.org/10.1088/1748-9326/ac48b3 710 Lembrechts, J. J., Nijs, I., & Lenoir, J. (2019) . Incorporating microclimate into species 711 distribution models. Ecography, 42(7), 1267–1279. https://doi.org/10.1111/ecog.03947 712 Lenoir, J., Bertrand, R., Comte, L., Bourgeaud, L., Hattab, T., Murienne, J., & Grenouillet, G. 713 (2020). Species better track climate warming in the oceans than on land. Nature Ecology 714 & Evolution, 4(8), 1044–1059. https://doi.org/10.1038/s41559-020-1198-2 715 MacArthur, R. A., & Wang, L. C. H. (1974). Behavioral thermoregulation in the pika Ochotona 716 princeps: A field study using radiotelemetry. Canadian Journal of Zoology, 52(3), 353–717 358. https://doi.org/10.1139/z74-042 718 Maclean, I. M. D., & Early, R. (2023). Macroclimate data overestimate range shifts of plants in 719 response to climate change. Nature Climate Change , 13(5), 484 –490. 720 https://doi.org/10.1038/s41558-023-01650-3 721 Marske, K. A., Lanier, H. C., Siler, C. D., Rowe, A. H., & Stein, L. R. (2023). Integrating 722 biogeography and behavioral ecology to rapidly address biodiversity loss. Proceedings 723 of the National Academy of Sciences , 120(15), e2110866120. 724 https://doi.org/10.1073/pnas.2110866120 725 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 32 Millar, C. I., Westfall, R. D., & Delany, D. L. (2016). Thermal components of American pika 726 habitat—How does a small lagomorph encounter climate? Arctic, Antarctic, and Alpine 727 Research, 48(2), 327–343. https://doi.org/10.1657/AAAR0015-046 728 Mitchell, P. J., Monk, J., & Laurenson, L. (2017). Sensitivity of fine‐scale species distribution 729 models to locational uncertainty in occurrence data across multiple sample sizes. 730

in Ecology and Evolution , 8(1), 12 –21. https://doi.org/10.1111/2041 -731 210X.12645 732 Morelli, T. L., Maher, S. P., Lim, M. C. W., Kastely, C., Eastman, L. M., Flint, L. E., Flint, A. 733 L., Beissinger, S. R., & Moritz, C. (2017). Climate change refugia and habitat 734 connectivity promote species persistence. Climate Change Responses , 4(1), 8. 735 https://doi.org/10.1186/s40665-017-0036-5 736 Moritz, C., Patton, J. L., Conroy, C. J., Parra, J. L., White, G. C., & Beissinger, S. R. (2008). 737 Impact of a Century of Climate Change on Small -Mammal Communities in Yosemite 738 National Park, USA. Science, 322(5899), 261–264. 739 https://doi.org/10.1126/science.1163428 740 Moudrý, V., Keil, P., Gábor, L., Lecours, V., Zarzo -Arias, A., Barták, V., Malavasi, M., 741 Rocchini, D., Torresani, M., Gdulová, K., Grattarola, F., Leroy, F., Marchetto, E., 742 Thouverai, E., Prošek, J., Wild, J ., & Šímová, P. (2023). Scale mismatches between 743 predictor and response variables in species distribution modelling: A review of practices 744 for appropriate grain selection. Progress in Physical Geography: Earth and 745 Environment, 47(3), 467–482. https://doi.org/10.1177/03091333231156362 746 Moyer-Horner, L., Beever, E. A., Johnson, D. H., Biel, M., & Belt, J. (2016). Predictors of 747 Current and Longer -Term Patterns of Abundance of American Pikas ( Ochotona 748 princeps) across a Leading -Edge Protected Area. PLOS ONE , 11(11), e0167051. 749 https://doi.org/10.1371/journal.pone.0167051 750 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 33 Nadeau, C. P., Giacomazzo, A., & Urban, M. C. (2022). Cool microrefugia accumulate and 751 conserve biodiversity under climate change. Global Change Biology , 28(10), 3222 –752 3235. https://doi.org/10.1111/gcb.16143 753 Nadeau, C. P., Urban, M. C., & Bridle, J. R. (2017). Coarse climate change projections for 754 species living in a fine‐scaled world. Global Change Biology , 23(1), 12 –24. 755 https://doi.org/10.1111/gcb.13475 756 Nogués‐Bravo, D. (2009). Predicting the past distribution of species climatic niches. Global 757 Ecology and Biogeography , 18(5), 521 –531. https://doi.org/10.1111/j.1466 -758 8238.2009.00476.x 759 Onoyama, K. (1991). Diurnal behaviors. In Hokkaido Government (Ed.), Survey Report of 760 Wildlife Distribution: Pika Ecology (pp. 56–65). 761 Onoyama, K., Kurumada, T., & Ohmija, S. (1991). Home range. In Hokkaido Government 762 (Ed.), Survey Report of Wildlife Distribution: Pika Ecology (pp. 66–94). 763 Patiño, J., Collart, F., Vanderpoorten, A., Martin‐Esquivel, J. L., Naranjo‐Cigal a, A., Mirolo, 764 S., & Karger, D. N. (2023). Spatial resolution impacts projected plant responses to 765 climate change on topographically complex islands. Diversity and Distributions, 29(10), 766 1245–1262. https://doi.org/10.1111/ddi.13757 767 Pearson, R. G. (2006). Climate change and the migration capacity of species. Trends in Ecology 768 & Evolution, 21(3), 111–113. https://doi.org/10.1016/j.tree.2005.11.022 769 Pecl, G. T., Araújo, M. B., Bell, J. D., Blanchard, J., Bonebrake, T. C., Chen, I. -C., Clark, T. 770 D., Colwell, R. K., Danielsen, F., Evengård, B., Falconi, L., Ferrier, S., Frusher, S., 771 Garcia, R. A., Griffis, R. B., Hobday, A. J., Janion -Scheepers, C., Jarzyna, M. A., 772 Jennings, S., … Williams, S. E. (2017). Biodiversity redistribution under climate 773 change: Impacts on ecosystems and human well -being. Science, 355(6332), eaai9214. 774 https://doi.org/10.1126/science.aai9214 775 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 34 Phillips, S. J., Anderson, R. P., & Schapire, R. E. (2006). Maximum entropy modeling of 776 species geographic distributions. Ecological Modelling , 190(3–4), 231 –259. 777 https://doi.org/10.1016/j.ecolmodel.2005.03.026 778 Phillips, S. J., & Dudík, M. (2008). Modeling of species distributions with Maxent: New 779 extensions and a comprehensive evaluation. Ecography, 31(2), 161 –175. 780 https://doi.org/10.1111/j.0906-7590.2008.5203.x 781 Potter, K. A., Arthur Woods, H., & Pincebourde, S. (2013). Microclimatic challenges in global 782 change biology. Global Change Biology , 19(10), 2932 –2939. 783 https://doi.org/10.1111/gcb.12257 784 QGIS Development Team. (2023). QGIS Geographic Information Sy stem (3.34) [Computer 785 software]. http://www.qgis.org 786 Randin, C. F., Engler, R., Normand, S., Zappa, M., Zimmermann, N. E., Pearman, P. B., Vittoz, 787 P., Thuiller, W., & Guisan, A. (2009). Climate change and plant distribution: Local 788 models predict high-elevation persistence. Global Change Biology, 15(6), 1557–1569. 789 https://doi.org/10.1111/j.1365-2486.2008.01766.x 790 Rodhouse, T. J., Hovland, M., & Jeffress, M. R. (2017). Variation in subsurface thermal 791 characteristics of microrefuges used by range core and perip heral populations of the 792 American pika ( Ochotona princeps ). Ecology and Evolution , 7(5), 1514 –1526. 793 https://doi.org/10.1002/ece3.2763 794 Sakiyama, T., & García Molinos, J. (2024). Northern pikas experience reduced occupancy due 795 to surrounding human land use d espite the occurrence of suitable microclimates. 796 Journal of Biogeography, jbi.14815. https://doi.org/10.1111/jbi.14815 797 Sakiyama, T., Morimoto, J., Watanabe, O., Watanabe, N., & Nakamura, F. (2021). Occurrence 798 of favorable local habitat conditions in an atypical landscape: Evidence of Japanese pika 799 microrefugia. Global Ecology and Conservation , 27, e01509. 800 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 35 https://doi.org/10.1016/j.gecco.2021.e01509 801 Scheffers, B. R., Evans, T. A., Williams, S. E., & Edwards, D. P. (2014). Microhabitats in the 802 tropics buffer temperature in a globally coherent manner. Biology Letters , 10(12), 803 20140819. https://doi.org/10.1098/rsbl.2014.0819 804 Schwalm, D., Epps, C. W., Rodhouse, T. J., Monahan, W. B., Castillo, J. A., Ray, C., & Jeffress, 805 M. R. (2016). Habitat availability and gen e flow influence diverging local population 806 trajectories under scenarios of climate change: A place-based approach. Global Change 807 Biology, 22(4), 1572–1584. https://doi.org/10.1111/gcb.13189 808 Shi, H., Wen, Z., Paull, D., & Guo, M. (2016). A framework for qu antifying the thermal 809 buffering effect of microhabitats. Biological Conservation , 204, 175 –180. 810 https://doi.org/10.1016/j.biocon.2016.11.006 811 Shiogama, H., Ishizaki, N. N., Hanasaki, N., Takahashi, K., Emori, S., Ito, R., Nakaegawa, T., 812 Takayabu, I., Hijiok a, Y., Takayabu, Y. N., & Shibuya, R. (2021). Selecting CMIP6 -813 Based Future Climate Scenarios for Impact and Adaptation Studies. SOLA, 17(0), 57–814 62. https://doi.org/10.2151/sola.2021-009 815 Smith, A. T. (2020). Conservation status of American pikas ( Ochotona princeps). Journal of 816 Mammalogy, 101(6), 1466–1488. https://doi.org/10.1093/jmammal/gyaa110 817 Smith, A. T., Nagy, J. D., & Millar, C. I. (2016). Behavioral Ecology of American Pikas 818 (Ochotona princeps) at Mono Craters, California: Living on the Edge. Western North 819 American Naturalist, 76(4), 459. https://doi.org/10.3398/064.076.0408 820 Stark, J. R., & Fridley, J. D. (2022). Microclimate‐based species distribution models in complex 821 forested terrain indicate widespread cryptic refugia under climate change. Global 822 Ecology and Biogeography, 31(3), 562–575. https://doi.org/10.1111/geb.13447 823 Stewart, J. A. E ., Perrine, J. D., Nichols, L. B., Thorne, J. H., Millar, C. I., Goehring, K. E., 824 Massing, C. P., & Wright, D. H. (2015). Revisiting the past to foretell the future: 825 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 36 Summer temperature and habitat area predict pika extirpations in California. Journal of 826 Biogeography, 42(5), 880–890. https://doi.org/10.1111/jbi.12466 827 Stickley, S. F., & Fraterrigo, J. M. (2023). Microclimate species distribution models estimate 828 lower levels of climate -related habitat loss for salamanders. Journal for Nature 829 Conservation, 72, 126333. https://doi.org/10.1016/j.jnc.2023.126333 830 Urban, M. C. (2015). Accelerating extinction risk from climate change. Science, 348(6234), 831 571–573. https://doi.org/10.1126/science.aaa4984 832 Varner, J., & Dearing, M. D. (2014). The importance of biologically relevant microclimates in 833 habitat suitability assessments. PLoS ONE , 9(8), e104648. 834 https://doi.org/10.1371/journal.pone.0104648 835 Vignali, S., Barras, A. G., Arlettaz, R., & Braunisch, V. ( 2020). SDMtune: An R package to 836 tune and evaluate species distribution models. Ecology and Evolution, 10(20), 11488–837 11506. https://doi.org/10.1002/ece3.6786 838 Xu, T., & Hutchinson, M. F. (2011). ANUCLIM Version 6.1 User Guide. 839 Xu, T., & Hutchinson, M. F. (20 13). New developments and applications in the ANUCLIM 840 spatial climatic and bioclimatic modelling package. Environmental Modelling & 841 Software, 40, 267–279. https://doi.org/10.1016/j.envsoft.2012.10.003 842 Yamazaki, D., Ikeshima, D., Tawatari, R., Yamaguchi, T., O’Loughlin, F., Neal, J. C., Sampson, 843 C. C., Kanae, S., & Bates, P. D. (2017). A high‐accuracy map of global terrain 844 elevations. Geophysical Research Letters , 44(11), 5844 –5853. 845 https://doi.org/10.1002/2017GL072874 846 Yates, K. L., Bouchet, P. J., Caley, M. J., Mengersen, K., Randin, C. F., Parnell, S., Fielding, 847 A. H., Bamford, A. J., Ban, S., Barbosa, A. M., Dormann, C. F., Elith, J., Embling, C. 848 B., Ervin, G. N., Fisher, R., Gould, S., Graf, R. F., Gregr, E. J., Halpin, P. N., … 849 Sequeira, A. M. M. (2018). Outstanding Challenges in the Transferability of Ecological 850 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 37 Models. Trends in Ecology & Evolution , 33(10), 790 –802. 851 https://doi.org/10.1016/j.tree.2018.08.001 852 Zhang, Z., Kass, J. M., Mammola, S., Koizumi, I., Li, X., Tanaka, K., Ikeda, K., Suzuki, T., 853 Yokota, M., & Usio, N. (2021). Lineage‐level distribution models lead to more realistic 854 climate change predictions for a threatened crayfish. Diversity and Distributions, 27(4), 855 684–695. https://doi.org/10.1111/ddi.13225 856 Zurell, D., Thuiller, W., Pagel, J., C abral, J. S., Münkemüller, T., Gravel, D., Dullinger, S., 857 Normand, S., Schiffers, K. H., Moore, K. A., & Zimmermann, N. E. (2016). 858 Benchmarking novel approaches for modelling species range dynamics. Global Change 859 Biology, 22(8), 2651–2664. https://doi.org/10.1111/gcb.13251 860 861 862 863 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 38 Tables 864 Table 1. Optimized hyperparameters and Continuous Boyce Index values of the best macroclimate and local-climate models. 865 Model Feature class Regularization multiplier Continuous Boyce index Hindcast Continuous Boyce index Macroclimate LQP 2.0 0.73 0.82 Local-climate LQP 2.0 0.80 0.96 866 867 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 39 Figures 868 869 870 Figure 1. Schematic diagram describing the different model and prediction settings of the 871 study. First, distribution modeling was conducted at two spatial resolutions: coarse (30 872 arcsecs) and fine (3 arcsecs). We used the macroclimate data in the coarse-scale setting to 873 reflect the conditions often used in traditional SDM studies. In the fine-scale setting, we used 874 the downscaled macroclimate data or the local-climate data, which we considered to represent 875 the local ambient conditions. Second, we predicted the baseline habitat suitability based on 876 these models. At the fine-scale, we further predicted the habitat suitability using the rock-877 interstice microclimate data, which we predicted based on empirical thermal measurements. 878 Third, we conducted a hindcasting analysis to evaluate the original models’ capability of 879 predicting species distribution over time. Fourth, we projected the future distribution based on 880 macroclimate, local-climate, and microclimate, considering four Global circulation models 881 and two emission scenarios. 882 883 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 40 884 Figure 2. Predicted habitat suitability of the whole distribution of the northern pika in 885 Hokkaido, Japan based on (a) macroclimate and (b) local-climate models. In (c), the habitat 886 suitability difference between these two predictions (Local – Macro) is provided. Black points 887 represent presence locations. In (d), the elevation of central Hokkaido is shown in the 888 background. The black border line denotes the model’s calibration extent with blue areas 889 within representing areas that were excluded from the analysis (see Methods for details). 890 891 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 41 892 Figure 3. Proportion of the study extent area predicted as present area along the binarization 893 threshold applied as a moving window at 0.01 increments (see Methods for details) for the (a) 894 baseline and (b) future periods. The lines and shade colors represent model settings (orange = 895 macroclimate, green = local-climate). The future panels are divided by time period (2041–896 2070 and 2071–2100) and emission scenario (SSP1-2.6 and SSP5-8.5), while the filled 897 envelopes represent the range of the future proportions of presence area predicted for the four 898 Global circulation models. 899 900 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 42 901 Figure 4. Predicted habitat suitability of the northern pika in (starting from top row) the Biei, 902 Shikaribetsu, and Oketo regions based on (starting from left column) the macroclimate, local-903 climate, and microclimate. Areas with rocky landforms are depicted by black lines in the 904 microclimate panels. The elevation for each region is shown in the last column. 905 906 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956 43 907 Figure 5 Predicted future distribution of the northern pika in the Biei, Shikaribetsu, and 908 Oketo regions based on the consensus approach. The 1st and 2nd columns represent predictions 909 for the 2041–2070 and 2071–2100 under the SSP1-2.6 scenario, and the 3rd and 4th those 910 under the SSP5-8.5 scenario. Dark gray areas represent areas suitable for the baseline period 911 that are predicted to become unsuitable in the future, while light gray areas represent areas 912 unsuitable for the baseline period. Areas with rocky landforms are depicted by the black lines. 913 The light blue areas represent suitable habitats supported by the local climate outside rocky 914 landforms. Within rocky landforms, the dark blue areas represent suitable habitats supported 915 by both the local-climate and microclimate, whereas the pink areas represent those supported 916 only by microclimate. 917 918 Author-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956