{"paper_id":"17cb73c9-62dc-4fa6-b45e-4851435e350a","body_text":"PREPRINT\nAuthor-formatted, not peer-reviewed document posted on 15/07/2024\nDOI: https://doi.org/10.3897/arphapreprints.e131956\nMapping fine-scale distribution of the northern pika \nOchotona hyperborea considering duality in\nmicrohabitat thermal conditions\n Tomoki Sakiyama, Jorge García Molinos\n\n1 \n \nArticle type: Research article 1 \n 2 \nMapping fine-scale distribution  of the northern pika Ochotona hyperborea considering 3 \nduality in microhabitat thermal conditions 4 \n 5 \nRunning head: Fine-scale SDM in northern pikas 6 \n 7 \nTomoki Sakiyama* 1 and Jorge García Molinos 2 8 \n 9 \n1 Graduate School of Environmental Science, Hokkaido University, Kita-10, Nishi-5, Kita-ku, 10 \nSapporo, Hokkaido, 060-0810, Japan 11 \n2 Arctic Research Center, Hokkaido University, Kita-21, Nishi-11, Kita-ku, Sapporo, Hokkaido, 12 \n001-0021, Japan 13 \n 14 \n* To whom correspondence should be addressed: wisdom.tree.1994@gmail.com 15 \n 16 \nTomoki Sakiyama: https://orcid.org/0000-0003-3655-0832 17 \nJorge García Molinos: https://orcid.org/0000-0001-7516-1835 18 \n  19 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n2 \n \nABSTRACT 20 \nSpecies distributions are frequently modeled using predictors that exceed the spatial scale 21 \nexperienced by the focal species. Incorporating fine-scale environmental conditions is 22 \ntherefore expected to lead to more realistic model predictions. However, the importance of the 23 \nvariety in the existing fine-scale conditions on species distribution remains poorly assessed 24 \nalthough species can effectively utilize multiple microhabitats for behavioral adaptation to 25 \nwithstand climate change impacts. Here, we developed a fine-scale species distribution model 26 \nbased on ambient air conditions for the northern pika (Ochotona hyperborea), a small 27 \nlagomorph found in rocky landforms, to first understand the improvement in the model 28 \nperformance from the conventional coarse-scale model. We then assessed how model 29 \npredictions alter by incorporating the rock-interstice microclimates that are found in their 30 \nhabitats for the baseline (1981–2010) and future periods (2041–2100). The fine-scale model 31 \nperformed better and overall predicted lower habitat suitability across the study area than the 32 \ncoarse-scale model. Incorporation of rock-interstice microclimate increased the habitat 33 \nsuitability markedly relative to predictions based on ambient conditions, which resulted in 34 \npredicting more suitable areas in lower elevations and more areas remaining suitable into the 35 \nfuture. This result suggests that the northern pika may withstand the negative impacts from 36 \nrising ambient temperatures by effectively utilizing rock interstices via behavioral adaptation. 37 \nOur findings highlight the importance of analyzing species distribution at fine scales and 38 \nconsidering local environmental heterogeneity, which helps species mitigate the adverse 39 \nimpacts of climate change, for conservation under climate change. 40 \nHIGHLIGHTS 41 \n1. Species use a wide variety of microhabitats and experience thermal conditions existing 42 \nlocally. 43 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n3 \n \n2. In the northern pika, a cold-adapted lagomorph, fine-scale species distribution model 44 \nperformed better than the coarse-scale model. 45 \n3. Complex topographical features that locally buffer ambient thermal conditions were 46 \npredicted to increase habitat suitability and elongate species persistence in the future. 47 \n4. Local environmental heterogeneity that helps species mitigate climate change impact will 48 \nbe important for conservation. 49 \n 50 \nKeywords: behavioral thermoregulation, climate change, microclimate, northern pika, species 51 \ndistribution model, topography 52 \n 53 \n  54 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n4 \n \nINTRODUCTION 55 \nClimate change is having a profound influence on species distribution worldwide. 56 \nEvidence from historical resurveys shows that species distribution limits are generally 57 \nexpanding at the leading edge and contracting at the trailing edge towards higher elevations and 58 \nlatitudes (Lenoir et al., 2020; Moritz et al., 2008) . These range shifts are predicted not only to 59 \naffect species extinction risks  (Urban, 2015), but also ecosystem functions and services (Pecl 60 \net al., 2017). It is therefore important in contemporary biogeography to assess accurately how 61 \nspecies will change their spatial distribution in the future  and identify conservation priority 62 \nareas to mitigate those i mpacts. Species distribution model s (SDMs) have become extremely 63 \nuseful in this role by relating species occurrence with environmental data and enabling forecasts 64 \nof how suitable habitats and associated distribution ranges will change  over time (Elith & 65 \nLeathwick, 2009) . Yet, SDMs have been criticiz ed for analyzing species occurrence and 66 \nenvironment relationships at spatial resolutions that are often too coarse relative to the actual 67 \nscales at which species utilize those habitats . A review by Potter et al. (2013) showed that the 68 \ngrid lengths of the environmental data used in previous studies were approximately 10,000 69 \ntimes larger than the body length of the focal animal species and 1,000 times in plants, 70 \nsuggesting a strong spatial mismatch between the environmental data and the actual conditions 71 \nexperienced by species. While some researchers argue that species distribution models aim to 72 \nmodel population-level rather than individual-level responses (Bennie et al., 2014), the issue of 73 \nscale mismatch in distribution modelling is still recognized as important when considering the 74 \npopulation dynamics for a wide range of taxa, such as small mammals, herpetofauna, and plants 75 \n(Nadeau et al., 2017) . Therefore, endeavors to develop SDMs that incorporate environmental 76 \nconditions at biologically relevant scales are warranted (Lembrechts et al., 2019). 77 \nThe spatial resolution of SDMs is of concern particularly in rugged terrains where coarse-78 \nscale climate data cannot capture appropriately the fine-scale environmental conditions created 79 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n5 \n \nby complex topographies (Austin & Van Niel, 2011). A growing body of literature has explored 80 \nthe effect of resolution and spatial scale on SDMs (Franklin et al., 2013; Gillingham et al., 2012; 81 \nGuisan et al., 2007; Maclean & Early, 2023; Moudrý et al., 2023; Patiño et al., 2023; Randin et 82 \nal., 2009). Recent studies are further attempting to improve the performance and prediction 83 \naccuracy of fine -scale models by shifting from simply using downscaled climate data to 84 \nincorporating microclimates that are predicted based on in-situ thermal sensor measurements 85 \n(Haesen et al., 2023; Stark & Fridley, 2022)  and remote sensing techniques (Stickley & 86 \nFraterrigo, 2023). These studies show that incorporation of fine-scale conditions can improve 87 \nmodel performance and alter predictions of suitable areas, thus affecting the trajectory of 88 \npotential range shifts over time (e.g., Maclean & Early, 2023) . Therefore, relying solely on 89 \ncoarse-scale models is problematic for conservation practice under climate change because it 90 \ncan lead to poor decisions that overlook suitable areas when designing protected areas and 91 \nallocating resources for management. Further research is therefore required in this field because 92 \na wide variety of environmental features create distinct climates at various spatial scales 93 \n(Morelli et al., 2017), and how species utilize and interact with them is system-specific (Shi et 94 \nal., 2016). 95 \nOne of the perspectives  still unexplored in this field , particularly for animals , is that a  96 \ngiven species can experience various fine -scale thermal conditions even within a single grid 97 \ncell as vertical structures also generate heterogenous conditions (e.g., Scheffers et al., 2014) . 98 \nMost previous studies have only considered the condition of a single vertical position, 99 \noverlooking the fine-grain heterogeneity  in environmental condition s that likely influences 100 \noverall habitat suitability of a given location that an individual occupies . This is particularly 101 \nimportant in the context of species adaptation to climate change since species may mitigate the 102 \nnegative impacts by altering their behavior effectively to exploit local heterogeneity in climatic 103 \nconditions (Beever, O’Leary, et al., 2016), which may cascade to have broader implications for 104 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n6 \n \npatterns in the overall distribution (Marske et al., 2023). Therefore, predicting habitat suitability 105 \nconsidering the variety of  fine-scale conditions should improve the  assessment of climate-106 \ninduced species range shifts. 107 \nPikas are small lagomorphs adapted to cold environments  in montane ecosystems . 108 \nVarious aspects of temperature, such as magnitude and variation, often act as strong drivers of 109 \ntheir distribution (Beever et al., 2010; Rodhouse et al., 2017) and abundance (Moyer-Horner et 110 \nal., 2016). Mounting evidence also suggest their vulnerability to recent climate change, with 111 \nfield observations reporting local extinctions (Erb et al., 2011; Stewart et al., 2015)  and range 112 \ncontractions (Beever, Perrine, et al., 2016; Billman et al., 2021) . Other studies have predicted 113 \nfuture range contractions (Bhattacharyya et al., 2019; Schwalm et al., 2016). However, despite 114 \noccurring in complex t errain, previous studies modeling pika distribution have only used 115 \ncoarse-scale climate data, possibly introducing bias into the habitat suitability predictions and 116 \nassociated range shift estimates.  117 \nSome pika species belonging to the Pika subgenus dwell in rocky landforms such as talus 118 \nslopes, lava flows, and moraines. At the fine-scale, these rock-dwelling pikas experience duality 119 \nin thermal conditions as they utilize both the above-ground and sub-ground microhabitats. The 120 \nabove-ground space is u sed for moving, vocalizing, foraging and collecting winter diet (i.e. 121 \nhaypile), and pikas experience the local  ambient air conditions during th ese activities. In 122 \ncontrast, the sub-ground space in the rock interstices is used as the shelter (Gliwicz et al., 2005), 123 \nwhich provide pikas with cool and stable microclimates buffered from the ambient air 124 \ntemperatures (Sakiyama et al., 2021; Varner & Dearing, 2014) . While both microhabitats are 125 \nessential for pikas, observations show that pikas respond flexibly to avoid the heat stress from 126 \nthe local ambient air by utilizing rock interstices (Smith et al., 2016). Therefore, rock interstices 127 \nare considered a key f eature for their distribution  by allowing the species to occur in 128 \ncounterintuitive, low-elevation areas and for them to persist under climate change (Smith, 2020). 129 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n7 \n \nPredicting the habitat suitability at fine-scales for both the above- and sub-ground space should 130 \nthus facilitate a more accurate assessment of  their vulnerability to climate change while also 131 \ntaking their adaptive capacity into account. 132 \nHere, we incorporate fine-scale thermal conditions  in rocky landforms into SDMs to 133 \npredict the distribution of the northe rn pika Ochotona hyperborea  in Hokkaido , Japan. We 134 \nevaluate how this incorporation changes model performance, suitable area prediction, and range 135 \nshift forecast from the conventional coarse-scale approach. In particular, we use macroclimate 136 \ndata (30 arcsecs) and downscaled macroclimate data (3 arcsecs) to build the coarse-scale model 137 \n(hereafter macroclimate model) and the fine-scale model ( local-climate model), respectively, 138 \nand compare their performances and habitat suitability predictions across the whole distribution. 139 \nWhile the local-climate model considers the effects of fine-scale, above-ground temperature, 140 \nwe further assessed the habitat suitability for the sub -ground space by incorporating rock-141 \ninterstice microclimates  into the predictions, which we modeled based on local thermal 142 \nmeasurements (Sakiyama & García Molinos, 2024). We then projected future range shifts for 143 \nthe species and compared the area of suitable habitat between local-climate and microclimate-144 \nbased predictions. We hypothesized that (1) the local -climate model will have higher model 145 \nperformance than the macroclimate model ; (2) the local -climate model  will predict habitat 146 \nsuitability distinctively from the macroclimate model reflecting the complex topography ; (3) 147 \nthe rock-interstice microclimate enhances the habitat suitability  from the local-climate model 148 \nprediction; and (4) the rock-interstices retain the suitable habitats in future predictions. 149 \nMATERIALS AND METHODS 150 \nStudy area and species distribution data 151 \nWe modeled the distribution of  northern pika populations in central Hokkaido, Japan, 152 \nbecause this is the core distribution area of the species in Japan and detailed fine-scale habitat 153 \nthermal measurements have been conducted in this area (Sakiyama et al., 2021; Sakiyama & 154 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n8 \n \nGarcía Molinos, 2024). The northern pika in Hokkaido is located at the southern margin of its 155 \ndistribution and therefore its population is likely among the most vulnerable of the species to 156 \nclimate warming. Although northerly neighboring populations are found across the Sea of Japan 157 \nin Sakhalin and Primorsky , the Hokkaido population is essentially a disjunct population . As 158 \nsuch, it is likely to respond to climate change by shifting upslope, but previous studies have not 159 \npredicted future distributions. 160 \nTo model the distribution of northern pikas, we compiled occurrence data  in central 161 \nHokkaido from previous literature. Given almost all the compiled data comprised presence (but 162 \nnot absence) locations, we decided to apply a presence-only algorithm for distribution modeling 163 \n(see Model Development section). We only used presence locations reported on high resolution 164 \nmaps (i.e. 1/25,000 and 1/50,000 scale)  for analysis since we were interested in modeling at 165 \nfine scales. Overall, a total of 305 presence points encompassing a wide elevational gradient 166 \n(50–2,210 m) were obtained from the literature (Table S1; Fig. S1a ). Survey years spanned 167 \nbetween 1963 and 2009 (Fig. S1a, b). Detection signs of northern pika presence were clearly 168 \nmentioned for almost all points (298 points, 97.7%), with the vocalization representing the main 169 \nmethod enabling detection of 241 points (79.0%). This suggests that the detection methods were 170 \nrelatively consistent among different literature consulted. 171 \nAs calibration extent for this study, we considered the area corresponding to the known 172 \ngeographical range of the species comprising the continuous mountainous areas in central 173 \nHokkaido delimited by the flatlands (i.e., river valleys) as the extent borders. Specifically, we 174 \nconsidered the Abira and Ishikari River systems as the western, the Teshio River system as the 175 \nnorthwestern, the Shokotsu River system and the Sea of Okhotsk as the northeastern, the 176 \nAbashiri and the Tokachi River systems as the eastern, and the Pacific Ocean as the southern 177 \nborders of the extent (Fig. S1a). Next, we removed irrelevant land cover types from within the 178 \nstudy extent area such as urban land, agricultural fields, and aquatic area (e.g., lakes) since 179 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n9 \n \nnorthern pikas require pristine montane conditions for inhabitation. The resulting filtered study 180 \nextent area covered approximately 21,000 km2 in central Hokkaido (Fig. S1a). 181 \nModel settings 182 \nTo understand how difference s in spatial resolution s of climate data  affect model 183 \nperformances and predictions of species distribution, we constructed macroclimate ( coarse-184 \nscale) and local-climate (fine-scale) models (Fig. 1; Table S2). In the macroclimate model, we 185 \nused environmental predictors at a resolution of 30 arcsecs, which converts to a grid cell of 691 186 \nm by 926 m (longitude x latitude) and size of 639,682 m 2. We consider this setting to be a 187 \ncoarse resolution for the northern pika because the ir home ran ge is much smaller, typically 188 \nwithin the range 1,822 –11,530 m 2 (n = 4; Onoyama et al., 1991) . We also note that, the 189 \nmountainous terrains of our study area are significant ly smoothed at this resolution, with 190 \nmedian and maximum elevational differences within one 30-arcsec grid cell as large as 182 m 191 \nand 770 m, respectively, when computed based on 10-times finer (3 arcsecs) elevation data (Fig. 192 \nS2). Since species occurrence-environment relationships might be confounded by such a large 193 \ndegree of simplification, we developed our local climate model using environmental data at 3 194 \narcsecs, which converts to a grid cell of 69.1 m by 92.6 m and size of 6,397 m 2. We consider 195 \nthis setting to reflect the local ambient air conditions within individual home ranges. Based on 196 \nthe local-climate model, we then generated two habitat suitability predictions to account for the 197 \nvariety of fine -scale thermal conditions . The fir st corresponds to the model raw predictions 198 \nbased on local-climate data and represents the above-ground habitat suitability. The second 199 \nprediction incorporates information on rock-interstice microclimates to further predict  sub-200 \nground habitat suitability. 201 \nTopographical and geological variables 202 \nWe used the terrain slope angle as topographical predictor in the SDMs because  it has 203 \nbeen shown to have a positive relationship with the northern pika distribution (Sakiyama et al., 204 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n10 \n \n2021). We computed the slope angle raster in QGIS 3.34 (QGIS Development Team, 2023)  205 \nusing the MERIT DEM (Yamazaki et al., 2017), a digital elevation model provided at a spatial 206 \nresolution of 3 arcsecs. Slope angles for the macroclimate setting were calculated using a DEM 207 \naggregated to 30 arcsecs (Table S2). We also considered surface geology as a categorical binary 208 \npredictor because 99 % of the known distribution of the northern pika lies on igneous and 209 \nmetamorphic rocks and moraines (Kawabe, 1992). To do so, we assigned whether any of these 210 \ngeological types occurred per grid cell with in the study area  based on the Seamless Digital 211 \nGeological Map of Japan (1:200,000) V2 (Geological Survey of Japan, National Institute of 212 \nAdvanced Industrial Science and Technology, 2017) (Table S2). 213 \nBioclimatic variables 214 \nTo model and predict the distribution of northern pikas, w e computed the mean daily 215 \ntemperature and mean diurnal range over July and August to account for the effects of thermal 216 \nmagnitude and variation during the summer. We focused on summer conditions (i.e., July and 217 \nAugust) because the effects of climate change on cold-adapted species are likely to be 218 \nparticularly strong during this season. Briefly (see Appendix S1 for a detailed explanation), we 219 \nfirst obtained climate data at 30 arcsecs using the chelsa-cmip6 R package, which provides 220 \naccess to downscaled Global circulation model (GCM) projections (Karger et al., 2023). Four 221 \nGCMs that are known to reproduce well the climatic conditions of our study area were selected 222 \n(i.e. ACCESS-CM2, IPSL-CM6A-LR, MIROC6, and MRI-ESM2-0) (Shiogama et al., 2021). 223 \nWe used climate data for 1981–2010 to model the baseline distribution. To create the 224 \nbioclimatic variables, we adapted the computation formulas used in ANUCLIM , a software 225 \npackage used to build spatial climatic data  (Xu & Hutchinson, 2011, 2013) , to the summer 226 \nseason, and calculated mean summer temperature  and mean summer diurnal range.  The 227 \nresulting variables were then used for the macroclimate setting (Table S2). 228 \nFor the local-climate setting, we downscaled the climate data obtained at 30 arcsecs to 3 229 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n11 \n \narcsecs as follows. Temperature data was downscaled using kriging with the krigR R package 230 \n(Davy & Kusch, 2021; Kusch & Davy, 2022) . This approach first  models the relationship 231 \nbetween temperature and elevation at the coarse scale (30 arcsecs), which is  then used to 232 \ninterpolate the fine-scale temperature using elevation at the same target scale (3 arcsecs). We 233 \nused the raw MERIT DEM (Yamazaki et al., 2017)  dataset, provided at 3 arcsecs, as the 234 \nelevation data for the fine-scale interpolation and the aggregated version at 30 arcsecs to model 235 \nthe coarse-scale relationship (Table S2).  236 \nFinally, for the rock-interstice microclimate, we considered that northern pikas can utilize 237 \nsub-ground thermal conditions in rocky landforms, which are known to be buffered from the 238 \nambient air temperature (Sakiyama and García Molinos, 2024). However, given the absence of 239 \nsufficient spatial distribution information regarding rocky landforms throughout our study 240 \nextent, we selected three representative regions  where we had sufficient data separated by 241 \nelevational ranges: Biei (1,100–2,000 m), Shikaribetsu (600–1,250 m), and Oketo (400–800 m) 242 \n(Fig. 4). For Biei, we were able to map the rocky landforms directly using aerial images. For 243 \nShikaribetsu and Oketo, we referred to the mapped occurrence of rocky landforms in previous 244 \nliterature (Kojima, 2004; Kurumada, 201 0). For these regions, we then used an empirical 245 \nrelationship between ambient and interstice temperatures previously derived from northern pika 246 \ntaluses in our study area (Sakiyama and Garcia Molinos, 2024) to convert downscaled ambient 247 \ntemperature into interstice temperature for each of the 3 arc second grid cells containing rocky 248 \nlandforms (Fig. S3; Appendix S1). Here, our assumption is that northern pikas experience local 249 \nambient climatic conditions  in cells devoid of  rocks since these landforms are patchily 250 \ndistributed within each region. We validated our downscaling and microclimate incorporation 251 \napproach by comparing the thermal conditions measured in northern pika habitats with the 252 \nvalues of each bioclimatic variable in the grid cells where those observations had been made 253 \n(see Appendix S2). Results confirmed that  the local-climate data was moderately accurate in 254 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n12 \n \nreflecting observed ambient temperature and that the microclimate data was highly accurate in 255 \nreflecting observed rock interstice temperatures (Fig. S4, 5). 256 \nModel development 257 \nWe performed the distribution modeling analysis of northern pikas in Hokkaido using 258 \nmaximum entropy (MaxEnt; Phillips et al., 2006) . Models were developed separately for the 259 \nmacroclimate and the local -climate settings considering the years 1981–2010 as the baseline 260 \nperiod, thus using the presence points during this period for model fitting. Prior to model fitting, 261 \nwe reduced the effect of spatial sampling bias by thinning the presence data using the spThin R 262 \npackage (Aiello-Lammens et al., 2015)  at a distance of 500 m as presence locations could be 263 \nviewed as independent at this distance considering the complex mountainous terrain and the  264 \nlimited mobility of the northern pika . Considering that the diagonal distance of macroclimate 265 \ngrid cells  is larger than this thinning distance  and duplicate sampling thus could occur, we 266 \nrandomly selected a single point for each cell after thinning. This pre-selection process resulted 267 \nin a total of 80 presence points remaining. We generated 10,000 background points throughout 268 \nthe study area, which are used in MaxEnt for characterizing the extant environmental conditions 269 \nwithin the study area. To avoid multicollinearity issues (Dormann et al., 2013) , we computed 270 \nthe pairwise Pearson’s correlation coefficients among the continuous predictors and confirmed 271 \nthat all pairs did not have a strong correlation (|r| > 0.7) (Fig. S6). 272 \nSince MaxEnt models with default settings can lead to overfitting (Anderson & Gonzalez, 273 \n2011), we initially tuned the model hyperparameters to control the complexity while fitting well 274 \nto the data. To do this, we used various hyperparameter combinations of feature classes, which 275 \ndetermines the shape of the response curves, and regularization multipliers, which controls the 276 \nmodel complexity (Phillips et al., 2006; Phillips & Dudík, 2008), and selected the combination 277 \nthat resulted in the highest predicative accuracy. For the feature class, we considered linear (L), 278 \nquadratic (Q), product (P), and/or hinge (H) characteristics, and f or the regularization 279 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n13 \n \nmultipliers, we used values ranging from 0.2 to 2.0 with intervals of 0.2.  The predictive 280 \naccuracy was assessed for each model by separating the presence data into training (80%) and 281 \ntest (20%) data and computing the area under the receiver operating characteristic curve (AUC), 282 \nwhich indicates the accuracy condition on a range of 0 to 1 (Fielding & Bell, 1997). All tuning 283 \nprocesses were  conducted using the optimizeModel function in  the SDMtune R package 284 \n(Vignali et al., 2020) . After hyperparameter tuning, we evaluated the optimized model’s 285 \npredictive accuracy based on 5-fold cross validation using Continuous Boyce Index (Hirzel et 286 \nal., 2006), and compared these values among the model settings. 287 \nComparing baseline predictions 288 \nFor both the macroclimate and the local -climate m odels, w e predicted the habitat 289 \nsuitability for the baseline period across the study area and compared the predicted suitability 290 \nby computing the difference. We also generated binarized maps of predicted species occurrence 291 \n(i.e. presence or absence) based on the habitat suitability prediction across the whole study area. 292 \nSince discriminant thresholds were likely to differ between the optimized models  for the 293 \nmacroclimate and local -climate and could be problematic for comparison, we instead used a 294 \nsensitivity analysis approach by applying a moving threshold from 0 to 1 at 0.01 increments to 295 \ndiscriminate presence or absence. Binary predictions were then  compared between the 296 \nmacroclimate and local-climate models by assessing the proportion of predicted presence areas 297 \nin the study area, calculated as the total amount of presence cells divided by the number of all 298 \ncells, for each moving threshold. We used p roportion instead of simply calculating the total 299 \narea of presence because the total study area was subtly different between the scales due to  300 \nprocessing procedures of environmental data.  We then predicted the habitat suitability for the 301 \nsub-ground space for the three selected regions by replacing the local ambient temperatures by 302 \nthe modeled rock interstice temperatures for the prediction of habitat suitability using the local-303 \nclimate model . We consider this prediction to reflect the habitat suitability of the rock -304 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n14 \n \ninterstices within individual home ranges defined based on the relationship between 305 \noccurrences and ambient air  temperatures. The resultant prediction was compared with those 306 \nbased on the macroclimate and local-climate. 307 \nHindcasting 308 \nWe use a hindcasting approach t o evaluate the transferability of the models to project 309 \ndistributions over time (see Nogués‐Bravo, 2009) by projecting the distribution of the species 310 \nto the past and  comparing it to  known past presence points from the same past period for 311 \nvalidation (Hodel et al., 2022). Although this approach is still rarely used in SDM studies (Yates 312 \net al., 2018) , it allowed us to assess how confident we can be using the models to predict 313 \ndistributions for a (past or future) time period different to that the model was calibrated for. We 314 \nprojected the distribution of northern pikas for years 1961–1980, a period for which early 315 \npresence records are available (Table S1; Fig S1). The same procedures as the baseline climate 316 \ndata were used for processing past data. We evaluated the model’s predicative accuracy via the 317 \nContinuous Boyce Index (Hirzel et al., 2006). 318 \nFuture projection 319 \nTo project the future distribution of the northern pika, we obtained climate data for two 320 \nfuture time periods, 2041–2070 and 2071–2100, considering two emission scenarios, SSP1-2.6 321 \nand SSP5-8.5, representing a globally sustainable and a fossil-fuel based development scenario, 322 \nrespectively. We predicted the future distribution for the study extent based on the macroclimate 323 \nand local-climate models using the future bioclimatic data, which were processed based on the 324 \nsame procedures for the baseline climate da ta (i.e., we assume that future topographical and 325 \ngeological variables will remain the same as present ). We binarized the suitability predictions 326 \nto presence or absence by applying the same approach used for the baseline period to calculate 327 \nthe proportion of presence areas.  328 \nMoreover, we also projected the future distribution of the northern pika for the three 329 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n15 \n \nselected regions to assess how microclimate contributes to retaining suitable conditions in the 330 \nfuture. When processing the future microclimate data, we assumed that the buffering effect of 331 \nthe rock interstices  (i.e., the empirical relationship between local ambient and rock -interstice 332 \ntemperatures) will remain invariable over time. We believe that this assumption  is robust 333 \nbecause the field thermal measurements from which this relationship was derived contained a 334 \nlarge elevational gradient of 350–2200 m and the model residual was stable across the ambient 335 \ntemperature gradient (Sakiyama and Garcia Molinos, 2024) , although we note that some 336 \nextrapolation did occur for some lower elevation areas where future conditions exceeded  the 337 \nobserved thermal gradient. We binarized the suitability predictions for both local-climate and 338 \nmicroclimate using the maximum true skill statistic computed for the local-climate model. The 339 \nconsensus among the resulting  binarized predictions from each of the four GCMs was then 340 \ncomputed by considering a cell as present whenever that cell was predicted as presen t by two 341 \nor more individual GCMs (i.e., ≥50% consensus) . The consensus predictions in the selected 342 \nregions were compared visually between the local-climate and microclimate-based predictions 343 \nfor each period and emission scenarios. 344 \nRESULTS 345 \nModel performance and variable importance 346 \nBoth the macroclimate and local-climate models converged to the same hyperparameter 347 \ncombinations for the best -fit model comprising an LQP feature class and a regularization 348 \nparameter value of 2. Using these hyperparameters, the Continuous Boyce Index was 0.73 and 349 \n0.80 for the macroclimate and local-climate models, respectively (Table 1). This indicates that 350 \nthe continuous  predictions reflected the species distribution well  and that the local-climate 351 \nmodel was more accurate than the macroclimate model. 352 \nMean summer temperature was the most important predictor explaining the northern pika 353 \ndistribution in both models  with a permutation importance of 74.5 % and 76.2 % for the 354 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n16 \n \nmacroclimate and local-climate models, respectively (Fig. S 7). Mean summer temperature 355 \ndecreased occurrence probability in both models, indicating that the northern pika is more likely 356 \nto be present in cooler areas. However, we detected a slightly steeper slope in the response 357 \ncurve of the local-climate model, which indicates that the local -climate model fitted  towards 358 \ncooler conditions than the macroclimate model (Fig. S8). In both models, surface geology was 359 \nthe second most important predictor  with 20% and 18 % of permutation importance , 360 \nrespectively (Fig. S7), and a positive relationship  with occurrence probability (Fig. S8). The 361 \npermutation contribution of daily thermal range and slope was less than 5% in both models (Fig. 362 \nS7). 363 \nHabitat suitability for the baseline period 364 \nThe negative effect of mean summer temperature  was reflected in the habitat suitability 365 \npredictions at the whole distribution scale in both models, with higher elevations predicted as 366 \nhighly suitable areas and the  suitability lowering gradually towards downslope (Fig. 2 a,b). 367 \nSpecifically, large mountain ranges such as the Daisetsuzan and Hidaka Mountain Ranges , 368 \nwhich contain summits over 2 ,000 m , as well as mountains in the Eastern Daisetsu region  369 \n(elevation of the main peaks 1,800–2,000 m) were recognized as the most suitable areas. Mid- 370 \nto high -elevation mountain ranges, such as the Yubari-Ashibetsu (1,600–1,720 m) , Teshio 371 \n(1,400–1,600 m), Ashoro (1,200–1,300 m), and Shikaribestu Mountains (1,100–1,250 m), were 372 \npredicted as areas  with moderate to high suitability . Finally, low-elevation mountains  were 373 \npredicted to be low in habitat suitability for northern pikas in almost all regions although some 374 \npresence points were reported from such areas (e.g. central and northeastern region s). One 375 \nexception was the Samani-Erimo region, which was predicted to have high suitability despite 376 \nthe region being predominated by low elevations (Fig. 2a,b).  377 \nAs we compared the macroclimate-based and local-climate based predictions, the former 378 \nwas generally higher in habitat suitability than the latter across the study area  (Fig. 2c). The 379 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n17 \n \ndifference in suitability was relatively large in the mountainous areas, particularly in the slopes 380 \nrather than the summits, with suitability differences up to 0.5. However, habitat suitability from 381 \nthe local-climate based were also higher than the macroclimate -based predictions in some 382 \nmountainous areas as well as the Samani-Erimo region. The habitat suitability difference was 383 \nrelatively small in low elevations. 384 \nThe proportion of predicted presence areas for both models indicated similar curves along 385 \nthe moving binarization threshold  (Fig. 3a). The curve inverse shape means that, while the 386 \nwhole study area is discriminated as present when the threshold is zero, a large proportion 387 \nbecomes absent at relatively low thresholds, although a small proportion of the study area 388 \nremains present until the threshold becomes o ne. In other words, while most of the study area 389 \nis predicted to have low suitability, some areas have distinctively high values. Comparing the 390 \ncurves between the models, the proportion of presence areas was higher in the macroclimate-391 \nbased prediction than  in the local -climate based prediction  across the whole range of 392 \nbinarization threshold (Fig. 3a). This suggests that the macroclimate model predicted a larger 393 \narea of suitable habitat than the local-climate model regardless of the binarization threshold. 394 \nIn the three selected regions, the habitat suitability based on the macroclimate and local-395 \nclimate showed a gradually changing pattern , likely  reflecting the elevation difference. In 396 \ncontrast, accounting for the rock interstice microclimate in local-climate predictions boosted 397 \nthe habitat suitability distinctively within these regions (Fig. 4). For instance, while the habitat 398 \nsuitability in Biei was relatively high given its high elevation in both the macroclimate - and 399 \nlocal-climate-based predictions, the rocky landforms further increased the suitability to mostly 400 \nreaching a suitability of 1  (Fig. 4 ). This trend was  more evident in Shikaribetsu , the mid -401 \nelevation region, where the predicted suitability reached as high as 0.99  in the microclimate -402 \nbased predictions despite the surrounding landscape having suitability less than 0.5 (Fig. 4). In 403 \nOketo, the low-elevation region, the total area that increased in suitability was relatively small 404 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n18 \n \ngiven the limited area of rocky landforms, but the predicted suitability within rocky landforms 405 \nreached as high as 0.66 (Fig. 4). The maximum values of predicted suitability based on the 406 \nmacroclimate and local-climate models were 0.55 and 0.41 in Shikaribetsu and 0.24 and 0.25 407 \nin Oketo, respectively. In all these three regional-scale predictions, the habitat suitability was 408 \npredicted in a more detailed and complex manner in the fine-scale (3 arc sec) predictions than 409 \nthat in the coarse -scale (30 arc sec) predictions as evidenced by the variation in habitat 410 \nsuitability of the local-climate and microclimate predict ions within the coarse-scale grid cells 411 \n(Fig. 4). 412 \nHindcasting and future projection 413 \nIn the hindcasting analysis, both macroclimate and local-climate models had Continuous 414 \nBoyce Index values well above zero  (Table 1), suggesting that the models were a ccurate in 415 \npredicting past distributions. This result supports the transferability of these models to predict 416 \ndistributional change over time . Accordingly, we were able to predict future distributions for 417 \nboth model settings with confidence.  418 \nBoth models indicated inverse curves  in the proportion of presence areas  in the future 419 \n(Fig. 3b). Comparing with the baseline proportions, future distribution was predicted to 420 \ndecrease by more than 60% in in both time periods and emission scenarios considered (Fig S9). 421 \nBetween the two future time periods, the difference in the habitat decrease rate was relatively 422 \nsmall under the SSP1-2.6 scenario, suggesting that presence areas will decrease until 2041–423 \n2070 but will then remain relatively stable until 2071–2100. However, the habitat decrease was 424 \napparent under the SSP5-8.5 scenario, suggesting that suitable habitats will continue to decrease 425 \nuntil 2071–2100. (Fig. 3b; S9). 426 \nPatterns of predicted distributional change largely differed among the three selected 427 \nregions between the local-climate an d microclimate -based predictions . Based on the local -428 \nclimate, predicted presences were found in Biei in both future time periods and scenarios and 429 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n19 \n \nin Shikaribetsu under the SSP1-2.6 scenario, while none were found in Oketo under all future 430 \nsettings (Fig. 5). Conversely, predicted presences were found in almost all regions in all future 431 \nsettings when accounting for rock interstice microclimates , except in Oketo for the far future 432 \nunder the SSP5-8.5 scenario (Fig. 5). Discrepancies between local-climate and microclimate-433 \nbased predictions were high towards the lower elevations and under the higher carbon emission 434 \nscenario (SSP5-8.5). 435 \nDISCUSSION 436 \nWe developed distribution model s for the northern pika that incorporate d fine-scale 437 \ntemperatures experienced by the species to assess how model performances and predictions 438 \nchange from a conventional model  based on coarse-scale conditions . Our results show that 439 \npredictive accuracy was higher for the fine-scale model. We also found that  the predicted 440 \nhabitat suitability varied among regions and at smaller spatial scales in the fine-scale prediction. 441 \nThese results confirm that accounting for fine-scale environments can lead to distinct and more 442 \ncomplex predictions than those obtained from a coarse-scale model and that solely depending 443 \non coarse-scale predictions can  result in biased predictions that overlook potential suitable 444 \nhabitats. 445 \nThere are two possible reasons that likely contributed to the higher predictive accuracy 446 \nin the fine-scale model. First, model response curves showed some important differences in the 447 \nmodeled relationship between occurrence and the environmental predictors . While mean 448 \ntemperature had a dominant negative effect on the distribution in both models (Fig. S7; S8), a 449 \nresponse also observed in a previous study (Sakiyama et al., 2021) , the response curve was 450 \nsteeper in the local-climate model. This indicates that the model fitted towards cooler conditions, 451 \nwhich may be reflecting the actual response of the northern pika more accurate ly than th e 452 \nmacroclimate model. Interestingly, our input data for the distribution model indicated that many 453 \nof the occurrence points in higher elevations were located in areas higher than its surroundings 454 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n20 \n \nand thus had cooler temperatures at the fine-scale (Fig. S10). Although we referred to presence 455 \npoints recorded in previous literature, it remains unclear whether this trend is due to a potential 456 \nsampling bias in convex topographies such as ridgelines and mountain tops. We believe that 457 \nassessing species distribution at fine-scales could introduce its scale-dependent issues, which 458 \nrequires attention as more fine-scale SDMs are developed  (Mitchell et al., 2017). Second, the 459 \nquality of the predictors may have also had an effect on model accuracy . Even when coarse- 460 \nand fine-scale models exhibit the same relationship between occurrence and environment, the 461 \nvalidated accuracy is likely to decrease if predictors are generalized at the coarse-scale. Instead, 462 \nfine-scale predictors are more likely to accurately represent true environmental conditions, 463 \nwhich could have resulted in higher predictive accuracy  by resolving better  local complex 464 \ntopographies. Thus, the overall increase in accuracy of the local-climate model can be attributed 465 \nto these small refinements in the modeled relationship and predictors. In general, our results are 466 \nin agreement with recent studies finding higher performance in SDMs incorporating fine-scale 467 \nclimates (Haesen et al., 2023; Stark & Fridley, 2022). 468 \nOur results showed that incorporating fine-scale environments into SDMs resulted in a 469 \ndecrease of suitable areas compared to the coarse-scale prediction. This is likely due  to 470 \ndifferences in the modeled response of mean summer temperature between the models. Such 471 \ndifferences likely emerged from the local -climate model being fitted towards cooler 472 \ntemperatures relative to the coarser macroclimate model resulting in lower predicted occurrence 473 \nprobabilities at certain temperatures or elevations. Consequently, this leads to a greater area of 474 \npredicted presence in the macroclimate model along the moving threshold when binarizing the 475 \nhabitat suitability predictions . In previous literature, fine-scale models have resulted in both 476 \ngreater (Stickley & Fraterrigo, 2023) and smaller suitable areas (Franklin et al., 2013; Haesen 477 \net al., 2023) . Hence, the difference in the predicted habitat  area resulting from incorporating 478 \nfine-scale environments into S DMs may be s pecific to each study system and the scale -479 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n21 \n \ndependent relationship between the occurrence records and environmental predictors used in 480 \nthe model . Furthermore, in our study, both models projected drastic future declines in the 481 \nnorthern pika distribution, with over half of the suitable habitats predicted to disappear by mid-482 \ncentury. This suggests that, similar to other pika species, northern pikas in Hokkaido are 483 \nvulnerable to climate change (Beever et al., 2011; Bhattacharyya et al., 2019). Importantly, the 484 \npredicted change in species  distribution depended strongly on the future period and emission 485 \nscenario. This result stresses the need to enforce effective policies for greenhouse gas emission 486 \nmitigation and  to guide effective conservation of the species by protecting areas that will 487 \naccount for the effect of climate change and maintain its distribution in the future. 488 \nAt the regional scale, the incorporation of the rock-interstice microclimate provided 489 \ndistinct predictions for both the baseline and future. Overall, the cooler conditions provided by 490 \nrock-interstices resulted in the microclimate predictions markedly enhancing habitat suitability 491 \ncompared to t he local-climate based prediction in all three regions. This demonstrates that a 492 \nsingle location (grid cell) can exhibit dual habitat suitability due to the existence of diverse 493 \nthermal conditions in complex topograph ies. However , what are the implications for the 494 \nnorthern pika distribution? In rocky landforms, the habitat suitability can simply be categorized 495 \ninto three combinations  based on the local -climate and microclimate  conditions: high -high, 496 \nlow-high, and low-low. The combination of high-low suitability was not observed because rock 497 \ninterstices were always cooler  than local ambient conditions , hence more suitable. Northern 498 \npikas are expected to be present in areas with high-high suitability, and absent from areas with 499 \nlow-low suitability . However, whether areas with low-high suitability  serve as functional 500 \nhabitats for northern pikas will likely depend on the extent to which the rock interstice condition 501 \nmitigates the impact of poor thermal ambient conditions, which should be strongly linked to the 502 \ncapacity of northern pika s to adapt their thermoregulation behavior. Previous studies have 503 \nindicated that northern pika s utilize above-ground space during the cooler times of the day  504 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n22 \n \n(Onoyama, 1 991); a behavior also observed in American pika stud ies (Smith et al., 2016) . 505 \nAdditionally, research on the American pika suggests that it regulates body temperature by 506 \nreducing above-ground activities  and retreating into rock interstices (MacArthur & Wang, 507 \n1974). Thus, the northern pika may exhibit high behavioral flexibility to respond to unfavorable 508 \nambient conditions, enabling them to inhabit areas with low-high suitability. Consequently, the 509 \nactual overall suitability of a given rocky landform should be influenced by a more complex 510 \ncombination of these two suitability conditions (above and below ground).  511 \nDespite its simplicity, our  approach to predict microclimate-based habitat suitability is 512 \ninnovative because it enables detecting potential suitable habitats that are likely to be  513 \noverlooked in predictions based solely on ambient conditions, especially in the lower elevations. 514 \nTherefore, this approach can be used as a first -cut approach to identify potential priority areas 515 \nfor conservation. Assessing occupancy patterns in areas with low-high suitability will be 516 \ntherefore important for resource allocation , although the large elevational range of observed 517 \npresence of the northern pika, despite their cold adaptation, already suggests they are inhabiting 518 \nsuch counterintuitive areas perhaps through behavioral adaptation. Moreover, this approach has 519 \nimplications for identifying climate change refugia for the northern pika  as all three regions 520 \nwere predicted to contain areas supported by rock interstice microclimate in the future despite 521 \nthe deterioration of local-climate conditions. Areas with favorable microclimates that support 522 \nspecies persistence under a changing climate (i.e., microrefugia) are often cryptic (Provan & 523 \nBennett, 2008), and therefore, extensive identification of such areas based on spatial data should 524 \nfacilitate effective protection of vulnerable populations. In our case, however, the current lack 525 \nof an extensive dataset reporting presence of rocky landforms  in our study area  limits such 526 \nendeavors. Therefore, further refinements in geospatial data are warranted for predicting 527 \npotential microrefugia. Deeper understanding of the behavioral responses of the northern pika 528 \nto ambient and rock interstice thermal conditions will also be important to assess the reliability 529 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n23 \n \nof behavioral adaptation in driving species distribution patterns (Marske et al., 2023).  530 \nThe importance of f ine-scale SDMs are being increasingly recognized  as they often 531 \nprovide distinct predictions from conventional scales and information in fine-scale predictions 532 \nare essential for improving downstream objectives, such as designing effective conservation 533 \nplans under climate change  (Nadeau et al., 2022) . In this study, the local-climate model 534 \nprediction proved to be valuable for achieving such objectives as it provided distinct habitat 535 \nsuitability predictions. Additionally, an important implication of fine-scale SDMs is their role 536 \nin filling the gap between predictions and observations of ongoing species range shifts (Chen 537 \net al., 2011). This is because conventional approaches have been criticized for overestimating 538 \nrange shifts due to the lack of consideration for the fine-scale climates that may facilitate species 539 \npersistence in warmer regions (Pearson, 2006). Therefore, it is important to validate whether 540 \nfine-scale models effectively reduce d overestimation by conducting ground-truthing surveys. 541 \nA study by Maclean and Early (2023) exemplifies this, demonstrating that fine-scale predictions 542 \nreduced the number of species at extinction risk and were in agreement with the actual observed 543 \npatterns, while coarse-scale predictions estimated several species to become extinct. 544 \nOne of t he limitations of our study is that, despite using the downscaled macroclimate 545 \ndata, we still rely on a simplified representation of time and space to develop the fine-scale 546 \nmodels. Rocky landforms exhibit high heterogeneity, with ambient air and rock interstice 547 \ntemperatures varying depending on height and depth, respectively, and changing across seasons 548 \n(Benedict et al., 2020; Millar et al., 2016; Smith, 2020) . In this regard, our study may still be 549 \nconsidered as a coarse-scale representation of reality in time and space. Like other SDM studies, 550 \nour model relies on  the correlative relationships between presences and  climatic and 551 \ntopographical predictors. Studies have demonstrated that predicted distributions are altered by 552 \nincorporating biotic interactions (Bateman et al., 2012), dispersal settings (Zurell et al., 2016), 553 \nand lineage information (Zhang et al., 2021)  into the models. Therefore , further research is 554 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n24 \n \nneeded to incorporate these aspects into the prediction of the distribution of the northern pika. 555 \nIn conclusion, the incorporation of fine-scale climate data into SDMs yielded contrasting 556 \npredictions to the coarse -scale prediction . Future distribution s were  predicted to decrease  557 \nincreasingly with both further away time horizons and higher emission scenarios considered , 558 \nstressing the impo rtance of enforcing effective mitigation measures to curve greenhouse gas 559 \nemissions. Together, this information is essential for the conservation of northern pika s in 560 \nHokkaido, as previous studies have not assessed such impacts. Incorporation of rock-interstice 561 \nmicroclimate data significantly increased the habitat suitability and predicted more suitable 562 \nhabitats for both the baseline and future. Behavioral adaptation is considered crucial for the 563 \ninhabitance of the northern pikas at lower elevations, which is important for the conservation 564 \nof the northern pika because its distribution is often considered to be restricted to the alpine  565 \nzone. Based on these results, we advocate for the  identification of rocky landforms and their 566 \nintegration into conservation planning  to effectively leverage microclimates in conservation 567 \nefforts.  568 \n  569 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n25 \n \nACKNOWLEDGEMENTS 570 \nThis study was supported by the JSPS KAKENHI Grant Number JP22J10191, grant-in-aid of 571 \nThe Zoshinkai Fund for Protection of Endangered Animals, and grant-in-aid of The Inui 572 \nMemorial Trust for Research on Animal Science, Japan. 573 \nDECLARATION OF COMPETING INTEREST 574 \nThe authors have no competing interests to declare. 575 \nCRediT AUTHORSHIP CONTRIBUTION STATEMENT 576 \nTS contributed to conceptualization and design o f methodology, data analysis, funding 577 \nacquisition, and writing the original draft; TS and JGM contributed to discussion of results, 578 \nwriting the final manuscript and gave final approval for publication. 579 \nDATA AVAILABILITY 580 \nData will be made available on request. 581 \nSUPPLEMENTAL MATERIAL 582 \nThe following materials are available as part of the online article at  583 \nhttps://escholarship.org/uc/fb. 584 \nAppendix S1. Derivation of the bioclimatic variables. 585 \nAppendix S2. Validation of the bioclimatic variables. 586 \nTable S1. List of previous literature reporting the northern pika distribution in Hokkaido. 587 \nTable S2. Procedures for generating the environmental variables. 588 \nFig. S1. Geographical map of the presence points used in this study. 589 \nFig. S2. A boxplot indicating the elevational difference within macroclimate grid cells. 590 \nFig. S3. The modeled relationship between mean ambient and rock interstice temperatures. 591 \nFig. S4. Scatterplots used for validating the bioclimatic variables. 592 \nFig. S5. Validation plots of the bioclimatic variables used in the study. 593 \nFig. S6. Result of the correlation analysis.  594 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n26 \n \nFig. S7. Permutation importance of the model predictors. 595 \nFig. S8. Response plots of the model predictors. 596 \nFig. S9. Percentage of presence areas predicted to change in the future from the baseline period. 597 \nFig. S10. Elevational difference between the macroclimate and local-scale cells at the presence 598 \npoints. 599 \n  600 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n27 \n \nREFERENCES 601 \nAiello-Lammens, M. E., Boria, R. A., Radosavljevic, A., Vilela, B., & Anderson, R. P. (2015). 602 \nspThin: An R package for spatial thinning of species occurrence records for use in 603 \necological niche models. Ecography, 38(5), 541 –545. 604 \nhttps://doi.org/10.1111/ecog.01132 605 \nAnderson, R. P., & Gonzalez, I. (2011). Species -specific tuning increases robustness to 606 \nsampling bias in models of species distributions: An implementation with Maxent. 607 \nEcological Modelling , 222(15), 2796 –2811. 608 \nhttps://doi.org/10.1016/j.ecolmodel.2011.04.011 609 \nAustin, M. P., & Van Niel, K. P. (2011). Improving species distribution models for climate 610 \nchange studies: Variable selection and scale. Journal of Biogeography , 38(1), 1–8. 611 \nhttps://doi.org/10.1111/j.1365-2699.2010.02416.x 612 \nBateman, B. L., VanDerWal, J., Williams, S. E., & Johnson, C. N. (2012). Biotic interactions 613 \ninfluence the projected distribution of a specialist mammal under climate change. 614 \nDiversity and Distributions , 18(9), 861 –872. https://doi.org/10.1111/j.1472 -615 \n4642.2012.00922.x 616 \nBeever, E. A., O’Leary, J., Mengelt, C., West, J. M., Julius, S., Green, N., Magness, D., Petes, 617 \nL., Stein, B., Nicotra, A. B., Hellmann, J. J., Robertson, A. L., Staudinger, M. D., 618 \nRosenberg, A. A., Babij, E., Brennan, J., Schuurman, G. W., & Hofmann, G. E. (2016). 619 \nImproving Conservation Outcomes with a New Paradigm for Understanding Species’ 620 \nFundamental and Realized Adaptive Capacity: A new paradigm for defining adaptive 621 \ncapacity. Conservation Letters, 9(2), 131–137. https://doi.org/10.1111/conl.12190 622 \nBeever, E. A., Perrine, J. D., Rickman, T., Flores, M., Clark, J. P., Waters, C., Weber, S. S., 623 \nYardley, B., Thoma, D., Chesley -Preston, T., Goehring, K. E., Magnuson, M., 624 \nNordensten, N., Nelson, M., & Collins, G. H. (2016). Pika ( Ochotona princeps ) losses 625 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n28 \n \nfrom two isolated regions reflect temperature and water balance, but reflect habitat area 626 \nin a mainland region. Journal of Mammalogy , 97(6), 1495 –1511. 627 \nhttps://doi.org/10.1093/jmammal/gyw128 628 \nBeever, E. A., Ray, C., Mote, P. W., & Wilkening, J. L. (2010). Testing alternative models of 629 \nclimate-mediated extirpations. Ecological Applications , 20(1), 164 –178. 630 \nhttps://doi.org/10.1890/08-1011.1 631 \nBeever, E. A., Ray, C., Wilkening, J. L., Brussard, P. F., & Mote, P. W. (2011). Contemporary 632 \nclimate change alters the pace and drivers of extinction. Global Change Biology, 17(6), 633 \n2054–2070. https://doi.org/10.1111/j.1365-2486.2010.02389.x 634 \nBenedict, L. M., Wiebe, M., Plichta, M., Batts, H., Johnson, J., Mon k, E., & Ray, C. (2020). 635 \nMicroclimate and Summer Surface Activity in the American Pika (Ochotona princeps). 636 \nWestern North American Naturalist, 80(3). https://doi.org/10.3398/064.080.0303 637 \nBennie, J., Wilson, R. J., Maclean, I. M. D., & Suggitt, A. J. (2014). Seeing the woods for the 638 \ntrees – when is microclimate important in species distribution models? Global Change 639 \nBiology, 20(9), 2699–2700. https://doi.org/10.1111/gcb.12525 640 \nBhattacharyya, S., Mungi, N. A., Kawamichi, T., Rawat, G. S., Adhikari, B. S., & Wilkening, 641 \nJ. L. (2019). Insights from present distribution of an alpine mammal Royle’s pika 642 \n(Ochotona roylei) to predict future climate change impacts in the Himalaya. Regional 643 \nEnvironmental Change, 19(8), 2423–2435. https://doi.org/10.1007/s10113-019-01556-644 \nx 645 \nBillman, P. D., Beever, E. A., McWethy, D. B., Thurman, L. L., & Wilson, K. C. (2021). 646 \nFactors influencing distributional shifts and abundance at the range core of a climate‐647 \nsensitive mammal. Global Change Biology , 27(19), 4498 –4515. 648 \nhttps://doi.org/10.1111/gcb.15793 649 \nChen, I.-C., Hill, J. K., Ohlemüller, R., Roy, D. B., & Thomas, C. D. (2011). Rapid range shifts 650 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n29 \n \nof species associated with high levels of climate warming. Science, 333(6045), 1024–651 \n1026. https://doi.org/10.1126/science.1206432 652 \nDavy, R., & Kusch, E. (2021). Reconciling high resolution climate datasets using KrigR. 653 \nEnvironmental Research Letters , 16(12), 124040. https://doi.org/10.1088/1748-654 \n9326/ac39bf 655 \nDormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., Marquéz, J. R. G., 656 \nGruber, B., Lafourcade, B., Leitão, P. J., Münkemüller, T., McClean, C., Osborne, P. 657 \nE., Reineking, B., Schröder, B., Skidmore, A. K., Zurell, D., & Lautenbach, S. (2013). 658 \nCollinearity: A review of methods to deal with it and a simulation study evaluating their 659 \nperformance. Ecography, 36(1), 27 –46. https://doi.org/10.1111/j.1600 -660 \n0587.2012.07348.x 661 \nElith, J., & Leathwick, J. R.  (2009). Species distribution models: Ecological explanation and 662 \nprediction across space and time. Annual Review of Ecology, Evolution, and Systematics, 663 \n40(1), 677–697. https://doi.org/10.1146/annurev.ecolsys.110308.120159 664 \nErb, L. P., Ray, C., & Guralnick, R. (2011). On the generality of a climate-mediated shift in the 665 \ndistribution of the American pika ( Ochotona princeps ). Ecology, 92(9), 1730 –1735. 666 \nhttps://doi.org/10.1890/11-0175.1 667 \nFielding, A. H., & Bell, J. F. (1997). A review of methods for the assessme nt of prediction 668 \nerrors in conservation presence/absence models. Environmental Conservation, 24(1), 669 \n38–49. https://doi.org/10.1017/S0376892997000088 670 \nFranklin, J., Davis, F. W., Ikegami, M., Syphard, A. D., Flint, L. E., Flint, A. L., & Hannah, L. 671 \n(2013). Modeling plant species distributions under future climates: How fine scale do 672 \nclimate projections need to be? Global Change Biology , 19(2), 473 –483. 673 \nhttps://doi.org/10.1111/gcb.12051 674 \nGeological Survey of Japan, National Institute of Advanced Industrial Science and Technology. 675 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n30 \n \n(2017). Seamless Digital Geological Map of Japan V2 1: 200,000  [dataset]. 676 \nhttps://gbank.gsj.jp/seamless/ 677 \nGillingham, P. K., Huntley, B., Kunin, W. E., & Thomas, C. D. (2012). The effect of spatial 678 \nresolution on projected responses to cl imate warming. Diversity and Distributions , 679 \n18(10), 990–1000. https://doi.org/10.1111/j.1472-4642.2012.00933.x 680 \nGliwicz, J., Witczuk, J., & Pagacz, S. (2005). Spatial behaviour of the rock‐dwelling pika 681 \n(Ochotona hyperborea ). Journal of Zoology , 267(2), 113 –120. 682 \nhttps://doi.org/10.1017/S0952836905007211 683 \nGuisan, A., Graham, C. H., Elith, J., Huettmann, F., & the NCEAS Species Distribution 684 \nModelling Group. (2007). Sensitivity of predictive species distribution models to 685 \nchange in grain size. Diversity and Dist ributions, 13(3), 332 –340. 686 \nhttps://doi.org/10.1111/j.1472-4642.2007.00342.x 687 \nHaesen, S., Lenoir, J., Gril, E., De Frenne, P., Lembrechts, J. J., Kopecký, M., Macek, M., Man, 688 \nM., Wild, J., & Van Meerbeek, K. (2023). Microclimate reveals the true thermal nich e 689 \nof forest plant species. Ecology Letters , 26(12), 2043 –2055. 690 \nhttps://doi.org/10.1111/ele.14312 691 \nHirzel, A. H., Le Lay, G., Helfer, V., Randin, C., & Guisan, A. (2006). Evaluating the ability 692 \nof habitat suitability models to predict species presences. Ecological Modelling, 199(2), 693 \n142–152. https://doi.org/10.1016/j.ecolmodel.2006.05.017 694 \nHodel, R. G. J., Soltis, D. E., & Soltis, P. S. (2022). Hindcast‐validated species distribution 695 \nmodels reveal future vulnerabilities of mangroves and salt marsh species. Ecology and 696 \nEvolution, 12(9), e9252. https://doi.org/10.1002/ece3.9252 697 \nKarger, D. N., Chauvier, Y., & Zimmermann, N. E. (2023). chelsa‐cmip6 1.0: A python 698 \npackage to create high resolution bioclimatic variables based on CHELSA ver. 2.1 and 699 \nCMIP6 data. Ecography, e06535. https://doi.org/10.1111/ecog.06535 700 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n31 \n \nKawabe, M. (1992). Geological factor of distribution of Northern Pika Ochotona hyperborea 701 \nin Hokkaido. Bulletin of the Higashi Taisetsu Museum of Natural History, 14, 103–106. 702 \nKojima, N. (2004). Conservation biology of the Japanese pika Ochotona hyperborea yesoensis. 703 \nIwate University. 704 \nKurumada, T. (2010). The habitat patches distribution and utilization of the northern pika in 705 \nMt. Nakayama, Hokkaido. Report of Hokkaido Institute of Environmental Sciences, 36, 706 \n65–69. 707 \nKusch, E., & Davy, R. (2022). KrigR —a tool for downloading and statistically downscaling 708 \nclimate reanalysis data. Environmental Research Letters , 17(2), 024005. 709 \nhttps://doi.org/10.1088/1748-9326/ac48b3 710 \nLembrechts, J. J., Nijs, I., & Lenoir, J. (2019) . Incorporating microclimate into species 711 \ndistribution models. Ecography, 42(7), 1267–1279. https://doi.org/10.1111/ecog.03947 712 \nLenoir, J., Bertrand, R., Comte, L., Bourgeaud, L., Hattab, T., Murienne, J., & Grenouillet, G. 713 \n(2020). Species better track climate warming in the oceans than on land. Nature Ecology 714 \n& Evolution, 4(8), 1044–1059. https://doi.org/10.1038/s41559-020-1198-2 715 \nMacArthur, R. A., & Wang, L. C. H. (1974). Behavioral thermoregulation in the pika Ochotona 716 \nprinceps: A field study using radiotelemetry. Canadian Journal of Zoology, 52(3), 353–717 \n358. https://doi.org/10.1139/z74-042 718 \nMaclean, I. M. D., & Early, R. (2023). Macroclimate data overestimate range shifts of plants in 719 \nresponse to climate change. Nature Climate Change , 13(5), 484 –490. 720 \nhttps://doi.org/10.1038/s41558-023-01650-3 721 \nMarske, K. A., Lanier, H. C., Siler, C. D., Rowe, A. H., & Stein, L. R. (2023). Integrating 722 \nbiogeography and behavioral ecology to rapidly address biodiversity loss. Proceedings 723 \nof the National Academy of Sciences , 120(15), e2110866120. 724 \nhttps://doi.org/10.1073/pnas.2110866120 725 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n32 \n \nMillar, C. I., Westfall, R. D., & Delany, D. L. (2016). Thermal components of American pika 726 \nhabitat—How does a small lagomorph encounter climate? Arctic, Antarctic, and Alpine 727 \nResearch, 48(2), 327–343. https://doi.org/10.1657/AAAR0015-046 728 \nMitchell, P. J., Monk, J., & Laurenson, L. (2017). Sensitivity of fine‐scale species distribution 729 \nmodels to locational uncertainty in occurrence data across multiple sample sizes. 730 \nMethods in Ecology and Evolution , 8(1), 12 –21. https://doi.org/10.1111/2041 -731 \n210X.12645 732 \nMorelli, T. L., Maher, S. P., Lim, M. C. W., Kastely, C., Eastman, L. M., Flint, L. E., Flint, A. 733 \nL., Beissinger, S. R., & Moritz, C. (2017). Climate change refugia and habitat 734 \nconnectivity promote species persistence. Climate Change Responses , 4(1), 8. 735 \nhttps://doi.org/10.1186/s40665-017-0036-5 736 \nMoritz, C., Patton, J. L., Conroy, C. J., Parra, J. L., White, G. C., & Beissinger, S. R. (2008). 737 \nImpact of a Century of Climate Change on Small -Mammal Communities in Yosemite 738 \nNational Park, USA. Science, 322(5899), 261–264. 739 \nhttps://doi.org/10.1126/science.1163428 740 \nMoudrý, V., Keil, P., Gábor, L., Lecours, V., Zarzo -Arias, A., Barták, V., Malavasi, M., 741 \nRocchini, D., Torresani, M., Gdulová, K., Grattarola, F., Leroy, F., Marchetto, E., 742 \nThouverai, E., Prošek, J., Wild, J ., & Šímová, P. (2023). Scale mismatches between 743 \npredictor and response variables in species distribution modelling: A review of practices 744 \nfor appropriate grain selection. Progress in Physical Geography: Earth and 745 \nEnvironment, 47(3), 467–482. https://doi.org/10.1177/03091333231156362 746 \nMoyer-Horner, L., Beever, E. A., Johnson, D. H., Biel, M., & Belt, J. (2016). Predictors of 747 \nCurrent and Longer -Term Patterns of Abundance of American Pikas ( Ochotona 748 \nprinceps) across a Leading -Edge Protected Area. PLOS ONE , 11(11), e0167051. 749 \nhttps://doi.org/10.1371/journal.pone.0167051 750 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n33 \n \nNadeau, C. P., Giacomazzo, A., & Urban, M. C. (2022). Cool microrefugia accumulate and 751 \nconserve biodiversity under climate change. Global Change Biology , 28(10), 3222 –752 \n3235. https://doi.org/10.1111/gcb.16143 753 \nNadeau, C. P., Urban, M. C., & Bridle, J. R. (2017). Coarse climate change projections for 754 \nspecies living in a fine‐scaled world. Global Change Biology , 23(1), 12 –24. 755 \nhttps://doi.org/10.1111/gcb.13475 756 \nNogués‐Bravo, D. (2009). Predicting the past  distribution of species climatic niches. Global 757 \nEcology and Biogeography , 18(5), 521 –531. https://doi.org/10.1111/j.1466 -758 \n8238.2009.00476.x 759 \nOnoyama, K. (1991). Diurnal behaviors. In Hokkaido Government (Ed.), Survey Report of 760 \nWildlife Distribution: Pika Ecology (pp. 56–65). 761 \nOnoyama, K., Kurumada, T., & Ohmija, S. (1991). Home range. In Hokkaido Government 762 \n(Ed.), Survey Report of Wildlife Distribution: Pika Ecology (pp. 66–94). 763 \nPatiño, J., Collart, F., Vanderpoorten, A., Martin‐Esquivel, J. L., Naranjo‐Cigal a, A., Mirolo, 764 \nS., & Karger, D. N. (2023). Spatial resolution impacts projected plant responses to 765 \nclimate change on topographically complex islands. Diversity and Distributions, 29(10), 766 \n1245–1262. https://doi.org/10.1111/ddi.13757 767 \nPearson, R. G. (2006). Climate change and the migration capacity of species. Trends in Ecology 768 \n& Evolution, 21(3), 111–113. https://doi.org/10.1016/j.tree.2005.11.022 769 \nPecl, G. T., Araújo, M. B., Bell, J. D., Blanchard, J., Bonebrake, T. C., Chen, I. -C., Clark, T. 770 \nD., Colwell, R. K., Danielsen, F., Evengård, B., Falconi, L., Ferrier, S., Frusher, S., 771 \nGarcia, R. A., Griffis, R. B., Hobday, A. J., Janion -Scheepers, C., Jarzyna, M. A., 772 \nJennings, S., … Williams, S. E. (2017). Biodiversity redistribution under climate 773 \nchange: Impacts on ecosystems and human well -being. Science, 355(6332), eaai9214. 774 \nhttps://doi.org/10.1126/science.aai9214 775 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n34 \n \nPhillips, S. J., Anderson, R. P., & Schapire, R. E. (2006). Maximum entropy modeling of 776 \nspecies geographic distributions. Ecological Modelling , 190(3–4), 231 –259. 777 \nhttps://doi.org/10.1016/j.ecolmodel.2005.03.026 778 \nPhillips, S. J., & Dudík, M. (2008). Modeling of species distributions with Maxent: New 779 \nextensions and a comprehensive evaluation. Ecography, 31(2), 161 –175. 780 \nhttps://doi.org/10.1111/j.0906-7590.2008.5203.x 781 \nPotter, K. A., Arthur Woods, H., & Pincebourde, S. (2013). Microclimatic challenges in global 782 \nchange biology. Global Change Biology , 19(10), 2932 –2939. 783 \nhttps://doi.org/10.1111/gcb.12257 784 \nQGIS Development Team. (2023). QGIS Geographic Information Sy stem (3.34) [Computer 785 \nsoftware]. http://www.qgis.org 786 \nRandin, C. F., Engler, R., Normand, S., Zappa, M., Zimmermann, N. E., Pearman, P. B., Vittoz, 787 \nP., Thuiller, W., & Guisan, A. (2009). Climate change and plant distribution: Local 788 \nmodels predict high-elevation persistence. Global Change Biology, 15(6), 1557–1569. 789 \nhttps://doi.org/10.1111/j.1365-2486.2008.01766.x 790 \nRodhouse, T. J., Hovland, M., & Jeffress, M. R. (2017). Variation in subsurface thermal 791 \ncharacteristics of microrefuges used by range core and perip heral populations of the 792 \nAmerican pika ( Ochotona princeps ). Ecology and Evolution , 7(5), 1514 –1526. 793 \nhttps://doi.org/10.1002/ece3.2763 794 \nSakiyama, T., & García Molinos, J. (2024). Northern pikas experience reduced occupancy due 795 \nto surrounding human land use d espite the occurrence of suitable microclimates. 796 \nJournal of Biogeography, jbi.14815. https://doi.org/10.1111/jbi.14815 797 \nSakiyama, T., Morimoto, J., Watanabe, O., Watanabe, N., & Nakamura, F. (2021). Occurrence 798 \nof favorable local habitat conditions in an atypical landscape: Evidence of Japanese pika 799 \nmicrorefugia. Global Ecology and Conservation , 27, e01509. 800 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n35 \n \nhttps://doi.org/10.1016/j.gecco.2021.e01509 801 \nScheffers, B. R., Evans, T. A., Williams, S. E., & Edwards, D. P. (2014). Microhabitats in the 802 \ntropics buffer temperature in a globally coherent manner. Biology Letters , 10(12), 803 \n20140819. https://doi.org/10.1098/rsbl.2014.0819 804 \nSchwalm, D., Epps, C. W., Rodhouse, T. J., Monahan, W. B., Castillo, J. A., Ray, C., & Jeffress, 805 \nM. R. (2016). Habitat availability and gen e flow influence diverging local population 806 \ntrajectories under scenarios of climate change: A place-based approach. Global Change 807 \nBiology, 22(4), 1572–1584. https://doi.org/10.1111/gcb.13189 808 \nShi, H., Wen, Z., Paull, D., & Guo, M. (2016). A framework for qu antifying the thermal 809 \nbuffering effect of microhabitats. Biological Conservation , 204, 175 –180. 810 \nhttps://doi.org/10.1016/j.biocon.2016.11.006 811 \nShiogama, H., Ishizaki, N. N., Hanasaki, N., Takahashi, K., Emori, S., Ito, R., Nakaegawa, T., 812 \nTakayabu, I., Hijiok a, Y., Takayabu, Y. N., & Shibuya, R. (2021). Selecting CMIP6 -813 \nBased Future Climate Scenarios for Impact and Adaptation Studies. SOLA, 17(0), 57–814 \n62. https://doi.org/10.2151/sola.2021-009 815 \nSmith, A. T. (2020). Conservation status of American pikas ( Ochotona princeps). Journal of 816 \nMammalogy, 101(6), 1466–1488. https://doi.org/10.1093/jmammal/gyaa110 817 \nSmith, A. T., Nagy, J. D., & Millar, C. I. (2016). Behavioral Ecology of American Pikas 818 \n(Ochotona princeps) at Mono Craters, California: Living on the Edge. Western North 819 \nAmerican Naturalist, 76(4), 459. https://doi.org/10.3398/064.076.0408 820 \nStark, J. R., & Fridley, J. D. (2022). Microclimate‐based species distribution models in complex 821 \nforested terrain indicate widespread cryptic refugia under climate change. Global 822 \nEcology and Biogeography, 31(3), 562–575. https://doi.org/10.1111/geb.13447 823 \nStewart, J. A. E ., Perrine, J. D., Nichols, L. B., Thorne, J. H., Millar, C. I., Goehring, K. E., 824 \nMassing, C. P., & Wright, D. H. (2015). Revisiting the past to foretell the future: 825 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n36 \n \nSummer temperature and habitat area predict pika extirpations in California. Journal of 826 \nBiogeography, 42(5), 880–890. https://doi.org/10.1111/jbi.12466 827 \nStickley, S. F., & Fraterrigo, J. M. (2023). Microclimate species distribution models estimate 828 \nlower levels of climate -related habitat loss for salamanders. Journal for Nature 829 \nConservation, 72, 126333. https://doi.org/10.1016/j.jnc.2023.126333 830 \nUrban, M. C. (2015). Accelerating extinction risk from climate change. Science, 348(6234), 831 \n571–573. https://doi.org/10.1126/science.aaa4984 832 \nVarner, J., & Dearing, M. D. (2014). The importance of biologically relevant microclimates in 833 \nhabitat suitability assessments. PLoS ONE , 9(8), e104648. 834 \nhttps://doi.org/10.1371/journal.pone.0104648 835 \nVignali, S., Barras, A. G., Arlettaz, R., & Braunisch, V. ( 2020). SDMtune: An R package to 836 \ntune and evaluate species distribution models. Ecology and Evolution, 10(20), 11488–837 \n11506. https://doi.org/10.1002/ece3.6786 838 \nXu, T., & Hutchinson, M. F. (2011). ANUCLIM Version 6.1 User Guide. 839 \nXu, T., & Hutchinson, M. F. (20 13). New developments and applications in the ANUCLIM 840 \nspatial climatic and bioclimatic modelling package. Environmental Modelling & 841 \nSoftware, 40, 267–279. https://doi.org/10.1016/j.envsoft.2012.10.003 842 \nYamazaki, D., Ikeshima, D., Tawatari, R., Yamaguchi, T., O’Loughlin, F., Neal, J. C., Sampson, 843 \nC. C., Kanae, S., & Bates, P. D. (2017). A high‐accuracy map of global terrain 844 \nelevations. Geophysical Research Letters , 44(11), 5844 –5853. 845 \nhttps://doi.org/10.1002/2017GL072874 846 \nYates, K. L., Bouchet, P. J., Caley, M.  J., Mengersen, K., Randin, C. F., Parnell, S., Fielding, 847 \nA. H., Bamford, A. J., Ban, S., Barbosa, A. M., Dormann, C. F., Elith, J., Embling, C. 848 \nB., Ervin, G. N., Fisher, R., Gould, S., Graf, R. F., Gregr, E. J., Halpin, P. N., … 849 \nSequeira, A. M. M. (2018). Outstanding Challenges in the Transferability of Ecological 850 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n37 \n \nModels. Trends in Ecology & Evolution , 33(10), 790 –802. 851 \nhttps://doi.org/10.1016/j.tree.2018.08.001 852 \nZhang, Z., Kass, J. M., Mammola, S., Koizumi, I., Li, X., Tanaka, K., Ikeda, K., Suzuki, T., 853 \nYokota, M., & Usio, N. (2021). Lineage‐level distribution models lead to more realistic 854 \nclimate change predictions for a threatened crayfish. Diversity and Distributions, 27(4), 855 \n684–695. https://doi.org/10.1111/ddi.13225 856 \nZurell, D., Thuiller, W., Pagel, J., C abral, J. S., Münkemüller, T., Gravel, D., Dullinger, S., 857 \nNormand, S., Schiffers, K. H., Moore, K. A., & Zimmermann, N. E. (2016). 858 \nBenchmarking novel approaches for modelling species range dynamics. Global Change 859 \nBiology, 22(8), 2651–2664. https://doi.org/10.1111/gcb.13251 860 \n 861 \n 862 \n  863 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n38 \n \nTables 864 \nTable 1. Optimized hyperparameters and Continuous Boyce Index values of the best macroclimate and local-climate models.  865 \nModel Feature class Regularization multiplier Continuous Boyce index Hindcast Continuous Boyce index \nMacroclimate LQP 2.0 0.73 0.82 \nLocal-climate LQP 2.0 0.80 0.96 \n 866 \n  867 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n39 \n \nFigures 868 \n 869 \n 870 \nFigure 1. Schematic diagram describing the different model and prediction settings of the 871 \nstudy. First, distribution modeling was conducted at two spatial resolutions: coarse (30 872 \narcsecs) and fine (3 arcsecs). We used the macroclimate data in the coarse-scale setting to 873 \nreflect the conditions often used in traditional SDM studies. In the fine-scale setting, we used 874 \nthe downscaled macroclimate data or the local-climate data, which we considered to represent 875 \nthe local ambient conditions. Second, we predicted the baseline habitat suitability based on 876 \nthese models. At the fine-scale, we further predicted the habitat suitability using the rock-877 \ninterstice microclimate data, which we predicted based on empirical thermal measurements. 878 \nThird, we conducted a hindcasting analysis to evaluate the original models’ capability of 879 \npredicting species distribution over time. Fourth, we projected the future distribution based on 880 \nmacroclimate, local-climate, and microclimate, considering four Global circulation models 881 \nand two emission scenarios. 882 \n  883 \n      \n          \n                       \n            \n                        \n                                          \n                           \n                    \n                         \n                         \n            \n                               \n                                 \n            \n             \n            \n                  \n              \n         \n               \n                                       \n                               \n                         \n                           \n                 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n40 \n \n 884 \nFigure 2. Predicted habitat suitability of the whole distribution of the northern pika in 885 \nHokkaido, Japan based on (a) macroclimate and (b) local-climate models. In (c), the habitat 886 \nsuitability difference between these two predictions (Local – Macro) is provided. Black points 887 \nrepresent presence locations. In (d), the elevation of central Hokkaido is shown in the 888 \nbackground. The black border line denotes the model’s calibration extent with blue areas 889 \nwithin representing areas that were excluded from the analysis (see Methods for details). 890 \n  891 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n41 \n \n 892 \nFigure 3. Proportion of the study extent area predicted as present area along the binarization 893 \nthreshold applied as a moving window at 0.01 increments (see Methods for details) for the (a) 894 \nbaseline and (b) future periods. The lines and shade colors represent model settings (orange = 895 \nmacroclimate, green = local-climate). The future panels are divided by time period (2041–896 \n2070 and 2071–2100) and emission scenario (SSP1-2.6 and SSP5-8.5), while the filled 897 \nenvelopes represent the range of the future proportions of presence area predicted for the four 898 \nGlobal circulation models. 899 \n  900 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n42 \n \n 901 \nFigure 4. Predicted habitat suitability of the northern pika in (starting from top row) the Biei, 902 \nShikaribetsu, and Oketo regions based on (starting from left column) the macroclimate, local-903 \nclimate, and microclimate. Areas with rocky landforms are depicted by black lines in the 904 \nmicroclimate panels. The elevation for each region is shown in the last column. 905 \n  906 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956\n\n43 \n \n 907 \nFigure 5 Predicted future distribution of the northern pika in the Biei, Shikaribetsu, and 908 \nOketo regions based on the consensus approach. The 1st and 2nd columns represent predictions 909 \nfor the 2041–2070 and 2071–2100 under the SSP1-2.6 scenario, and the 3rd and 4th those 910 \nunder the SSP5-8.5 scenario. Dark gray areas represent areas suitable for the baseline period 911 \nthat are predicted to become unsuitable in the future, while light gray areas represent areas 912 \nunsuitable for the baseline period. Areas with rocky landforms are depicted by the black lines. 913 \nThe light blue areas represent suitable habitats supported by the local climate outside rocky 914 \nlandforms. Within rocky landforms, the dark blue areas represent suitable habitats supported 915 \nby both the local-climate and microclimate, whereas the pink areas represent those supported 916 \nonly by microclimate. 917 \n 918 \nAuthor-formatted, not peer-reviewed document posted on 15/07/2024. DOI:  https://doi.org/10.3897/arphapreprints.e131956","source_license":"CC-BY-4.0","license_restricted":false}