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

Context: Sound taxonomy is the cornerstone of biodiversity conservation. Without a fundamental understanding of species delimitations, as well as their distributions and ecological requirements, our ability to conserve them is drastically impeded. Cryptic species – two or more distinct species currently classified as a single species – present a significant challenge to biodiversity conservation. How do we assess the conservation status and address potential drivers of extinction if we are unaware of a species’ existence? Here, we present a case where the reclassification of a species formerly considered widespread and secure – the sugar glider (Petaurus breviceps) – has dramatically increased our understanding of the potential impacts of the catastrophic 2019–20 Australian megafires to this species. Methods: We modelled and mapped the distribution of the former and reclassified sugar glider (Petaurus breviceps). We then compared the proportional overlap of fire severity classes between the former and reclassified distribution, and intersect habitat suitability and fire severity to help identify areas of important habitat following the 2019–20 fires. Key Results: Taxonomic revision means that the distribution of this iconic species appears to have been reduced to 8% of its formerly accepted range. Whereas the 2019–20 Australian megafires overlapped with 8% of the formerly accepted range, they overlapped with 33% of the proposed range of the redefined Petaurus breviceps. Conclusions: Our study serves as a sombre example of the substantial risk of underestimating impacts of mega-disturbance on cryptic species, and hence the urgent need for cataloguing Earth’s biodiversity in the age of megafire.,Methods Occurrence data Occurrence records of P. breviceps were collected from the Atlas of Living Australia (https://www.ala.org.au) and were subject to a filtering process. Because sugar gliders were introduced to Tasmania (Campbell et al. 2018), we excluded all Tasmanian records. We then removed dubious records by clipping all records to either the former P. breviceps range (based on IUCN maps; IUCN 2020) or the proposed reclassified P. breviceps range (Cremona et al. 2021) to create two sets of occurrence data (i.e., one each for the former and reclassified P. breviceps). It is worth noting, however, that the reclassified distribution of P. breviceps proposed by Cremona et al. (2021) is an estimate based on genetic and morphological data. Although evidence currently suggests that the Great Dividing Range acts as the western edge of the distribution of P. breviceps (Cremona et al. 2021), we cannot be certain of this. However, for the purposes of this study we have assumed it to be so. In both datasets, records were removed if: (i) they were missing date information or were collected before the year 2000; or (ii) they had high locational uncertainty (e.g., vague or inaccurate locations). Records within any 1 × 1 km grid cell were collapsed into a single record. The final filtered data base consisted of 7777 presence records within the formerly considered geographic range, and 5089 within the reclassified range (see Figure S1). Geographic range estimation We mapped the extent of occurrence (EOO) of the former and reclassified P. breviceps using the occurrence datasets. Extent of occurrence is defined as the area enclosed by the shortest possible boundary containing all sites in which a species is known to be present (IUCN 2021). We calculated EOO as α‐hulls (a generalisation of convex polygons that allow for breaks in species ranges), using the ‘alphahull’ package in R version 3.6.2 (R Core Team 2021), specifying a α value of two (IUCN 2021). We regarded EOO as preferable to area of occupancy (AOO) because maps of the latter showed clear spatial bias indicated by high densities of records surrounding major capital cities. Species distribution modelling Using the maxent algorithm, we developed species distribution models (SDMs) based on the two occurrence datasets outlined above (Phillips et al. 2006). We selected SDM environmental layers based on their likely importance to P. breviceps habitat suitability. All environmental layers were resampled to 1 × 1 km resolution prior to being included in models. A set of 10,000 background points were included within the SDM to compare densities in environmental values occupied by P. breviceps with those of the surrounding unoccupied environment. We addressed sample bias within the study area with a ‘target group’ background sampling approach (Phillips et al. 2009) (see Figure S1). We defined the target group as arboreal mammal species occurring within the study area, including P. breviceps. Sampling intensity for target group species was mapped by converting species presence records of the target group to a kernel density map using the kde2d function of the ‘MASS’ package (Venables and Ripley 2002) set with the default kernel bandwidth. Model performance was measured as area under the curve (AUC) of the receiver operating characteristic (ROC) plot, and the contribution of environmental variables to the response variable was measured as permutation importance (Phillips 2005). Fire overlap We overlapped the former and reclassified P. breviceps EOO with 2019–20 bushfire severity maps from the Google Earth Engine Burnt Area Mapping (GEEBAM; DIPE 2020). GEEBAM classifies the cells within the fire boundary as one of five fire severity classes: no data (cleared land, water etc.); unburnt (unburnt and lightly burnt); low and moderately burnt (some or moderate change post-fire); high severity (vegetation mostly scorched); and very high severity (vegetation clearly consumed). When calculating fire overlap, we considered only fires occurring within the Department of Agriculture, Water and Environment’s (2020) ‘preliminary area for environmental analysis’ (following Legge et al. 2020). This area encompasses bioregions that were deemed to have experienced anomalously substantial fire activity during the 2019–20 bushfire season. Overlap measures were calculated using QGIS version 3.14.1 (QGIS Development Team 2021). We created a fire severity × habitat quality matrix to help identify the spatial intersection between fire severity and habitat quality for the reclassified P. breviceps. First, we classified the continuous output of relative habitat quality derived from the SDM into four discrete classes: low quality (relative likelihood of occurrence 0–0.25); low–medium quality (relative likelihood of occurrence 0.25–0.50); medium–high quality (relative likelihood of occurrence 0.5–0.75); and high quality (relative likelihood of occurrence 0.75–1). We then combined the reclassified SDM with the GEEBAM fire severity layer to derive a layer with 16 unique combinations of all combinations of habitat quality and fire severity and mapped this across the range of P. breviceps.,

We explored the range shifts of alien and native breeding bird species of mainland Great Britain between the periods 1968–1972 and 2007–2011. We estimated β-diversity and its components (richness difference and replacement component) of alien and native communities per time period and explored the effect of abiotic factors (climate, land cover, human population, and elevation) on alien and native β-diversity and their components and of native diversity (species richness and functional diversity) on alien β-diversity by Generalized Dissimilarity Modelling. Breeding bird distribution data of Great Britain for the time periods 1968–1972 and 2007–2011 (data collected in hectads of 100 km2) were retrieved from The atlas of breeding birds in Britain and Ireland (Sharrock, 1976), Bird Atlas 2007–2011: the breeding and wintering birds of Britain and Ireland (Balmer et al., 2013), and Breeding and wintering bird distributions in Britain and Ireland from citizen science bird atlas (Gillings et al., 2019). Filtering of the species was applied based on their status that was retrieved from the British Trust for Ornithology's BirdFacts database (Robinson, 2005; https://www.bto.org/understanding-birds/birdfacts), Mcinerny et al. (2022), and Wayman et al. (2022). The trait dataset (body mass, clutch size, foraging location, habitat, activity time, and main diet) of breeding birds was retrieved from Tsianou et al. (2021). Specifically, the traits were compiled from datasets (Storchová et al., 2018) and electronic databases [Handbook of the Birds of the World Alive (https://www.hbw.com)]. Climatic factors (mean temperature, mean temperature range, mean temperature of the warmest month, mean temperature of the coldest month, and mean precipitation), percentage of different land cover types (forest, cropland, grassland, other, and, water), human population and elevation per hectad were calculated with using QGIS 3.18 using data available online. The detailed sources of abiotic factors are provided in the attached Readme.txt and the original paper, while the calculated variables are provided in attached .csv files (climaticA.csv, climaticC, HILDAA.csv, HILDAC.csv, human.popA.csv, human.popC.csv, elevation.csv). All analyses were performed using R version 4.1.0 (R Development Core Team 2021) and required packages are provided in the attached 1_Prepare environmental and diversity data.txt. Additional functions are provided in Nstar function.txt and Windrose function.txt. Species list (along with name changes), trait dataset, and hectads of mainland Great Britain are provided in the Species list.csv (name changes.csv), Traits.csv, and Mainland_GB.csv, respectively. Scripts for performing analyses are provided in the 2_Estimate range shifts metrics.txt and 3_Run GDMs.txt. Instructions for performing analyses are provided in the Readme.txt. References Balmer, D. E., Gillings, S., Caffrey, B., Swann, R., Downie, I., & Fuller, R. (2013). Bird Atlas 2007-11: the breeding and wintering birds of Britain and Ireland: BTO Thetford. Gillings, S., Balmer, D. E., Caffrey, B. J., Downie, I. S., Gibbons, D. W., Lack, P. C., ... & Fuller, R. J. (2019). Breeding and wintering bird distributions in Britain and Ireland from citizen science bird atlases. Global Ecology and Biogeography, 28(7), 866-874. https://doi.org/10.1111/geb.12906 McInerny, C. J., Musgrove, A. J., Stoddart, A., Harrop, A. H., & Dudley, S. P. (2022). The British List: a checklist of birds of Britain (10th edition). Ibis, 164, 860–910. https://doi.org/10.1111/ibi.13065 Robinson, R.A. (2005). BirdFacts: profiles of birds occurring in Britain & Ireland. BTO Research Report, 407, p.Thetford. (http://www.bto.org/birdfacts) Sharrock, J. T. R. (1976). The atlas of breeding birds in Britain and Ireland: A&C Black. Storchová, L., & Hořák, D. (2018). Life‐history characteristics of European birds. Global Ecology and Biogeography, 27(4), 400-406. https://doi.org/10.1111/geb.12709 Tsianou, M. A., Touloumis, K., & Kallimanis, A. S. (2021). Low spatial congruence between temporal functional β‐diversity and temporal taxonomic and phylogenetic β‐diversity in British avifauna. Ecological Research, 36(3), 491-505. https://doi.org/10.1111/1440-1703.12209 Wayman, J. P., Sadler, J. P., Pugh, T. A., Martin, T. E., Tobias, J. A., & Matthews, T. J. (2022). Assessing taxonomic and functional change in British breeding bird assemblages over time. Global Ecology and Biogeography, 31(5), 925-939. https://doi.org/10.1111/geb.13468 Aim: We explored the range shifts of alien and native birds, the responses of alien and native β-diversity to abiotic factors, and the effect of native diversity on alien β-diversity in two time periods. Location: Great Britain. Time period: 1968–1972, 2007–2011. Taxa studied: Breeding birds. Methods: We estimated range shifts of alien and native species between the periods 196...

We used data from 23 canyons in the Shoshone Mountains, Toiyabe Range, Toquima Range, and Monitor Range (central Great Basin, Nevada, USA). We derived riparian area and fragmentation at the canyon level from 2013 National Agricultural Inventory Program (NAIP) images. We did not attempt to distinguish the cause of fragmentation or the historical extent of riparian cover in these canyons. There were no substantial disturbances or changes in land use in these canyons from 2001-2018 that would drive changes in the extent of riparian area or its fragmentation. It is possible that changes occurred gradually, but we believe that the slow rate of growth of woody vegetation minimized these changes over the temporal extent of our work. To calculate riparian area, we delineated a buffer (500 meters from the center line of the canyon bottom) that was large enough to contain all riparian vegetation in most canyons. We mapped riparian cover within the buffer in QGIS version 3.0.2 (QGIS.org 2018...

Understanding species distribution and habitat preferences is crucial for effective conservation strategies. However, the lack of information about population responses to environmental change at different scales hinders effective conservation measures. In this study, we estimate the potential and realized distribution of Atractus lasallei, a semi-fossorial snake endemic to the northwestern region of Colombia. We modelled the potential distribution of A. lasallei based on ecological niche theory (using maxent), and habitat use was characterized while accounting for imperfect detection using a single-season occupancy model. Our results suggest that A. lasallei selects areas characterized by slopes below 10°, with high average annual precipitation (>2500mm/year) and herbaceous and shrubby vegetation. Its potential distribution encompasses the northern Central Cordillera and two smaller centers along the Western Cordillera, but its habitat is heavily fragmented within this potential distribution. When the two models are combined, the species’ realized distribution sums up to 935 km2, highlighting its vulnerability. We recommend approaches that focus on variability at different spatio-temporal scales to better comprehend the variables that affect species’ ranges and identify threats to vulnerable species. Prompt actions are needed to protect herbaceous and shrub vegetation in this region, highly demanded for agriculture and cattle grazing. Methods Ecological Niche Model and Potential Distribution Presence data were acquired from three sources: 1) specimens from biological collections, obtained from the Global Biodiversity Information Facility (accessed 22 March 2022) [35] and most of them revised in situ in the following collections: MHUA-Museo de Herpetología Universidad de Antioquia (Curator: J.M. Daza-Rojas), CSJ-h-Museo de Ciencias Naturales de La Salle (D.Z. Urrego), CBUCES-D-Colecciones Biológicas Universidad CES (J.C. Duque); 2) iNaturalist records obtained directly from them (accessed March-May 2022; we did not use iNaturalist records from GBIF) by searching Atractus records in the northwest of Colombia that included pictures that allowed verification through morphology (coloration patterns, and scale counts when possible), and 3) individuals encountered during the field phase for occupancy models. Identification of individuals was based on the original species description and taxonomic revisions of the genus [28, 33]. Further, a geographical filter was applied to presence records that were within 1 km of each other to reduce spatial autocorrelation [36, 37]. We used the final database to delimit the species accessible area or M [38] based on the intersection between the minimum convex polygon generated with 50 km buffers around each presence record using QGIS (v.3.10) [39], and the biogeographic regions of the northern Central Cordillera and the Western Cordillera of Colombia 40. The environmental variables for the niche model included topographic variables, atmospheric climate variables including temperature corrected to ground level [41] and soil variables 42. Climate variables represent long-term averages (1980-2010 in the case of atmospheric variables, and 2000-2020 in the case of ground-level temperature; see S1). These variables were selected based on previous research findings regarding the distribution of semi-fossorial reptiles [25, 43-45]. All variables were used in the models at a spatial resolution of 1 km. Variables with finer resolutions were resampled using the bilinear method, with the R-package “raster” v.3.5 [46] in R v.4.2.1 [47]. Subsequently, a Spearman correlation test (S2) was conducted to select non-correlated variables (< 0.8), using R-package “corrplot” v.0.92 [48]. Finally, two sets of variables were created, one that included two ground-level temperature variables estimated at five centimeters above the ground (S1) [41], and the second included the same two variables but measured at atmospheric level [49]. Models were calibrated with each data set independently, ensuring all variables used were not correlated. The ecological niche model was generated using the maximum entropy algorithm [15] through the R-package “kuenm” [50]. This methodology allows the evaluation of different sets of environmental variables (set 1 and set 2, S1) and various model parameterizations to ultimately identify the best model according to a set of criteria. We allowed the regularization parameter to vary from 0.1 (very strict in relation to observed values) to 5 (more flexible in relation to observed values), where 1 is the default value. We also evaluated across linear (l), quadratic (q), and linear-quadratic (lq) responses. The models were trained using a 20% random partition of the occurrence data for model evaluation. The evaluation criteria included omission rate (<5%), partial receiver operating characteristic (partial ROC), area under the curve (AUC ratio>1), and the Akaike Information Criterion corrected for sample size (AICc) [51]. In case several models achieved the evaluation criteria, we performed a consensus model using the median of the selected models. Finally, to obtain the geographic projection, a cutoff threshold was applied using the 10 percentile training presence criteria from the best model(s) to generate a presence/absence map. Occupancy Models To identify fine-scale factors influencing the occupancy of A. lasallei within its known distribution area, single-season occupancy models were employed [24]. The sampling design followed the recommendations of a previous study for semi-fossorial species [52], wherein 30 linear transects of 100 m x 2 m were established within the sampling area, spaced at a minimum of 200 m apart to ensure independence of detection histories across sites (Fig 1). Each transect was equipped with nine artificial cover objects (three roof tiles, three boards, and three plastic sheets), which were installed a minimum of two months prior to sampling for the organisms to habituate to their presence (S3). The transects were surveyed between October 2021 and January 2022 to ensure consistent occupancy status during the sampling period (closed-site assumption) between 8 AM and 4 PM. Surveys involved searching beneath leaf litter and under cover objects (both artificial and natural). Each transect was surveyed a minimum of four times, with visits spaced at least two weeks apart to satisfy the assumption of temporal independence. Animals were photographed and examined in the field to ensure correct identification (Approval Act No. 138, February 9, 2021, granted by the Committee on Ethics for Animal Experimentation, Universidad de Antioquia). Occupancy models were constructed using the R-package “unmarked” (v.1.2.5) [53] implemented in the R software. All covariates were standardized (mean=0, units in standard deviations) prior to modelling. To identify the best models, we first established the best detection model assuming constant occupancy, and then we used this detection model in all occupancy models [54]. To model detectability, we included as covariates, the number of cover objects, both natural and artificial (N_obj); vegetation height (Veg_H) [55]; soil moisture (Soil_moisture); and soil temperature (T_ground), both measured using a HOBO proV2 datalogger beneath a roof tile or under the object where an individual of the species was located at the time of each visit. As covariates for occupancy, we used vegetation height (Veg_H) [55]; terrain slope (Slope); topographic convergence (Con); compound topographic index (CTI) [56]; annual mean soil temperature (Tprom), maximum temperature of the warmest month (Tmax), and minimum temperature of the coldest month (Tmin)[41]; depth of leaf litter (Leaf_Dep) and depth of the 0 horizon (Hori0), both measured in the field using a soil auger; euclidean distance to the nearest house (D_house), nearest forest (D_forest), and nearest water body (D_water). These distances were estimated in QGIS [39], identifying the nearest houses and forests to the centroid of each transect using satellite imagery from GoogleEarth (https://www.google.com/intl/es/earth/). To calculate the distance to water bodies, it was necessary to construct a detailed hydrographic network for the area using a 12.5 m resolution DEM obtained from Alaska vertex (https://search.asf.alaska.edu/), utilizing the hydrology toolbox in ArcGIS Pro (v.2.7) [57]. A total of 87 biologically plausible and simple models were evaluated, each including one or two variables (S4), 20 of the models were for the detection component with constant occupancy, and the remaining models were for the occupancy component. Finally, to evaluate model fit to our data, we performed a parametric bootstrap test on the chosen model, using the parboot function of R package “unmarked” v.1.4.1 [53]. This test generates multiple sets of data iteratively from the best model and then compares these sets with the detection histories obtained in the field. A chi-squared test was employed to evaluate the null hypothesis that the observations are consistent with the proposed model. Integration of models To estimate the species’ realized distribution area [58], we used the binary (presence-absence) geographic projection from the consensus niche model to identify the areas where the macro conditions were suitable and applied the best occupancy model within those areas at a higher spatial resolution (0.00025° 27 meters). Finally, the resulting map was transformed into a binary outcome using a threshold of 0.78, based on the Q3 (third quartile) of the occupancy distribution values of that map; this threshold corresponds to 4 m of vegetation height according to the best occupancy model (Fig 2), which is biologically justified if we consider that all presence records obtained in the field phase were found in places with vegetation below 4 m. References 15.