Land snail microclimate niches identify suitable areas for climate refugia management on a montane landscape

Climate refugia management is an emerging natural resource sub-discipline but identifying which species would benefit, their climatic requirements


Introduction
The maintenance, management, or creation of climate change refugia are variations of an adaptation strategy to buffer ecological conditions from anthropogenic climate change (Morelli et al., 2016).Although the idea of climate refugia exists predominantly in the scientific literature, natural resource managers are beginning to implement strategies to adaptively manage landscapes for climate refugia (Morelli et al., 2020).However, identifying those areas and the species which would benefit remains challenging and is largely relegated to predictive models with limited on the ground validation (Barrows et al., 2020).Further complicating matters, data deficiency for most species (Bland et al., 2017) results in a lack of the basic ecological information and fine scale knowledge of where those species occur which would be needed to inform such models.Terrestrial gastropods, also known as land snails, are a case in point.
Land snails include both terrestrial snails and slugs and are sensitive to climatic conditions (Nicolai and Ansart, 2017).Major threats include temperature extremes, drought (Schweizer et al., 2019), and loss of snow cover in temperate regions (Nicolai and Ansart, 2017).Unfavorable climatic conditions can have severe consequences for land snails including thermal and desiccation deaths (McQuaid et al., 1979) which, even in areas of moderate climate, have been known to cause mass morality events (Nicolai et al., 2011).Closely related land snail species have evolved different physiological strategies adapted to microclimatic conditions (Schweizer et al., 2019).While there is even some evidence that within-species land snail populations may have different cold tolerance levels (Augspurger, 2013), adaptive or plastic capacity is difficult to prove (Nicolai and Ansart, 2017).Further complicating matters is that, although not without exception (Simonová et al., 2016), land snails have limited dispersal ability to modify range extent with climate change (Nicolai and Ansart, 2017).Therefore, instead of predicting future habitats for ectotherms of limited vagility, it is particularly important to conserve locations where these species currently occur.For instance, Wilson et al. (2019) identified cold sites of several meters square to conserve for persistence of Brazil's Araucaria tree (Araucaria angustifolia).These local microsites which are thought to provide locations for species to persist during climate change are called climate microrefugia (Keppel et al., 2015).
Climate refugia management, in the form of managing microrefugia or clusters of microrefugia, may be an effective path toward conservation of climate sensitive land snails but three immediate problems must be overcome.First, determining which species would benefit from climate microrefugia management; second, mapping where those species currently occur and; third, mapping where suitable microrefugia occur (Barrows et al., 2020).Gastropod species are associated with a M.K. Lucid et al. wide range of microclimate regimes (Nicolai and Ansart, 2017;Schweizer et al., 2019) but teasing out if and how different species occupy microclimate niches remains relatively untested.Successful climate refugia management would still rely on species occurrence and distribution data which is relatively limited for most gastropod species despite their being one of the most imperiled groups of animals (Lydeard et al., 2004).
We developed a study to address these obstacles on a large and diverse landscape covering portions of northern Idaho, northeastern Washington, and northwestern Montana, U.S.A. From 2010 to 14 we colocated air temperature data loggers with terrestrial gastropod surveys at 830 sites with the following objectives; 1) determine the distribution of terrestrial gastropod species across the study area, 2) determine if microclimate, macroclimate, or non-microclimate variables influence species occurrence, and 3) if microclimate is influential, identify microrefugia and clusters of microrefugia most suitable for climate refugia management.

Study area
The study area consists of a 22,975 km 2 area centered on the Idaho Panhandle (Fig. 1).It includes portions of five mountain ranges; the Selkirk, Purcell, West Cabinet, Coeur d'Alene, and Saint Joe Mountains.The topography is complex ranging from broad glacial valleys to mountainous areas with moderate to substantial relief.Elevation ranges from 702 to 2,326 m.Land uses include urban, rural developed, agricultural, and relatively pristine areas but is primarily managed temperate forest.The heavily forested area is dominated by a diverse mix of conifer species and is characterized as supporting inland temperate rainforest (DellaSala, 2011).The vast majority of survey sites were dominated by coniferous trees.
The study area climate is characterized by wet cool springs followed by warm and dry summers.Winters are wet and moderately cold with highly variable snowpack.Low elevation snowpack ranges from nonexistent to persisting for several months while higher elevations are characterized by deeper snowpacks that persist into early summer.Several climatic trends have been documented in the study area over past century.First, no discernable trend has been observed in maximum daily temperatures but the minimum daily temperature has increased 2.8 • F. Second, there has been no change in total annual precipitation but there has been a 33% increase in average annual stream run off due, in part, to a shift in timing due to earlier snowpack melt.Lastly, March 1st snowpack at lower elevations has decreased by 30% (Tinkham et al., 2015).

Sampling design
We stratified our study area into 5x5 km sampling cells and colocated gastropod surveys with microclimate data loggers at 1-2 sites per cell.We used ArcGIS 10.1 (Environmental Systems Research Institute, Redlands, California, USA) to generate a buffer around each road and trail from 50 to 150 m.We then generated a random point within this buffer for the survey location.This resulted in sites (n = 691) randomly located but biased to roads and trails to improve field efficiency.We also conducted surveys at randomly selected Forest Inventory and Analysis (FIA) plots (n = 139; Bechtold and Patterson, 2005) which were not biased to roads.In total, we co-located data loggers with gastropod surveys at 830 sites in 750 5x5km sampling cells.

Data collection and preparation
TRIX8, TRIX16, and HAXO8 LogTag® Transit data loggers were deployed within plastic radiation shields (Holden et al., 2013) designed to protect the logger from direct sunlight.The radiation shield was attached with nails on the north side of a conifer tree > 30 cm in diameter.Data loggers collected air temperature data every 90 min for 12-48 months.To determine microclimate variables for this study, we used only air temperature collected from each data logger for a continuous 12-month period when all sites with data loggers (n = 830) were being monitored (1 September 2013-31 August 2014).
Loggers could record erroneous data in some situations such as sunlight directly hitting the radiation shield or the logger being buried by deep snow.Therefore, we developed a set of rules (Table 1) which served as quality control checks to remove erroneous data from the data set.First, we identified observations that were likely to be erroneous because they were hotter or colder than temperatures likely to be observed in the region.Based on the temperature record from all NOAA weather stations within 50 km of any of the data loggers during the years of data collection, we removed any observations above 38.8• C or below 30 • C.
Next, we identified observations that were likely erroneous due to direct sun in or around the data logger creating a temperature spike.From visual inspection of the temperature records, we removed any observations that occurred while the temperature increased at a rate of >2.5 • C/hr.Finally, we identified observations that were likely erroneous because the data logger was buried under snow and thus not recording the air temperature.To do this, we removed any observations during winter and spring months corresponding to days where the daily temperature variation was <2 • C and the daily average temperature was near freezing (0C ± 1C).To minimize the impact of missing data on estimation of the microclimate variables, we removed daily records with >2 missing observations.For days with 1-2 missing observations, we imputed the missing values using a smoothing spline.We used R ( R Core Team, 2015) to perform all data manipulations and to derive eight microclimate variables (Table 2).
To determine macroclimate variables, we downloaded 30-year normal (1981-2010) mean annual temperature and precipitations at 800 m and 4 km resolutions from PRISM (PRISM Climate Group, 2021) for each of our survey locations.We performed a correlation analysis in the R environment (R Core Team, 2015) which showed these variables were highly correlated (r = 0.91) so we selected the 4 km resolution to provide a larger scale differential from the microclimate data.We used ArcGIS (Environmental Systems Research Institute, Redlands, California, USA) to determine elevation, aspect, distance from road, and percentage tree cover for each survey location.

Land snail surveys
We deployed land snail survey transects at each site consisting of cover board traps, leaf litter searches, timed searches, and pitfall traps (see Lucid et al., 2018a for full details).This method allows for detection of ground dwelling land snails but not arboreal or dendrophilic species (Lucid et al., 2018a).Specimens were stored in 95% ethanol and keyed to species according to characteristics outlined by Burke and Leonard (2013) and Lucid et al. (2018b).We detected 51 land snail species; however, we limited our statistical analysis to species we detected at a minimum of 10 sites.Therefore, our analysis includes 27 land snail species which we detected at a minimum of 10 and maximum of 264 sites.

Table 1
Microclimate air temperature algorithm rules.

Statistical analysis
We used Random Forest (RF) to identify variables most related to each species' occurrence.RF is a non-parametric machine learning algorithm that utilizes a Classification and Regression Trees (CART) and bagging approach (Breiman, 2001).Unlike many traditional parametric statistical models, RF does not make assumptions regarding data distribution, and thus performs well in assessing nonlinear relationships and is unsusceptible to overfitting and multicollinearity even when processing a large number of covariates.Briefly, in RF, a large number of decision trees are being generated.Each decision tree is randomly created by sampling two-thirds of the training data with replacement while the remaining third is used for performance evaluation (i.e., outof-bag [OOB] error estimate).Finally, the majority vote of all trees is used to create the final model.
We performed RF in the R environment using the RandomForest package (Liaw and Wiener, 2002).A separate RF model was developed for each species using the presence-absence occurrence data as the response variable.Predictor variables in the model included a suit of microclimate metrics, elevation, aspect, distance to road, and percent tree cover (Table 2).Variables were calculated study-area-wide for all species with the exception of seven species which were not detected in all mountain ranges (Fig. 2).For the seven species with restricted ranges we only used data from mountain ranges where the species is known to occur.For example, we only detected Polygyrella polygyrella in the Coeur d'Alene and Saint Joe Mountains.Therefore, we removed sites in the Purcell, Selkirk, and West Cabinet Mountains from the P. polygyrella models.Because we detected species at fewer sites than at which we collected microclimate data the number of presences was often low compared with absences.To deal with imbalanced presence-absence issue, we used the rf.classBalance function in the rfUtilities package when presences consisted of less than one-third of the presence-absence data set (Evans and Murphy, 2018).We used the rf.crossValidation in the rfUtilities package to calculate cross-validated metrics (i.e., OOB error estimate and percent correctly classified [PCC]) for evaluating model performances (Evans and Murphy, 2018).

Microclimate site scoring
To assign a microclimate score to each of the 830 survey sites we first calculated the median of each variable in Microsoft Excel®.For each variable and each site we assigned a value of − 1 or + 1 dependent on if the variable indicated site conditions which were cooler or warmer respectively than the median.For example, the median meanT for all sites was 6.4 • C. If a site had a meanT of 6.0 • C it received a score of − 1.If a site had a meanT of 7.0 • C it received a score of + 1.We did not calculate a score for TD because it was not a significant variable in any model.This left us with seven variables which were scored − 1 or + 1.We added those scores for each site which resulted in a microclimate score ranging from − 7 (coolest) to + 7 (warmest).

Species assignment to temperature groups
As described above, we used the median value of each microclimate variable to determine the 'cool' and 'warm' cutoffs.We excluded macroclimate and non-climate variables from the temperature assignment process.For each species and significant variable we categorized the variable as cool or warm if it fell entirely below or above the median value.We considered variables that spanned across the cutoff as 'overlap' values.For example, if a species temperature range was 6-8 • C, it would receive an 'overlap' for that variable because it includes values below and above the median meanT of 6.4 • C.

Gastropod surveys and microclimate data
We detected 51 terrestrial gastropod species (Table 3, see Lucid et al., 2018b for full list).We detected 27 species at ≥ 10 sites which was our minimum sample size for the Random Forest Analysis (Table 4).Twentyone gastropod species were fairly well distributed across the five mountain ranges in study area but six were not detected in at least one mountain range (Fig. 2).4).

Random forest analysis
With the exception of Zonitodes arboreus and Kootenai burkei, models for all species supported one (Anguispira kochi, Polygyrella, Punctum randolphi, and Zacoleus idahoensis) to five (M.magnipelta) microclimate variables.Models supported macroclimate variables for fifteen species with support for meanT4K for eleven species and ppt4K for six species.Elevation was supported in models for twelve species.With the exception of Tree Cover supporting the R. abietum model, non-climate variables other than elevation were not supported for any species.The most commonly supported microclimate variables were DDGT5 (n = 16) and meanT (n = 13) while the least commonly supported were minT (n = 3) maxT (n = 2), and TD (n = 0).Scaled variable importance support ranged from 0.01 to 0.08 (Table 5).

Microclimate site scoring
Forty-eight percent (n = 402) of sites had negative 'cool' scores and 52% (n = 428) had positive 'warm' scores.The Purcells had the highest ratio (64%) of 'cool' sites and the Saint Joe had the highest ratio (58%) of 'warm' sites (Fig. 1, Table 6).Microclimate scores were distributed across elevation groups with no discernable pattern.

Discussion
Our results demonstrate a clear pattern that the climate variables we tested are more important than the non-climate habitat variables we tested for land snails and that species can be grouped into cool, warm, and generalist temperature groups.Elevation is often used as a surrogate for microclimate variables (e.g.Mumladze et al., 2017), however, our data do not show a clear pattern linking elevation with microsite air temperatures (Fig. 3) and elevation was only marginally important across species.The variability in temperature within elevation groups is likely due to cold air pooling and midslope thermals, which are what cause microclimate to diverge from elevation lapse rates (e.g.Holden et al., 2011), particularly in cirque and deep valley locations which are common throughout our study system.This suggests the need to use microclimate in addition to just elevation as a surrogate to modelled climate products that use elevation lapse rates.This is further evidenced by the differing performance of the microclimate meanT and macroclimate meanT4K variables.At least one of these variables supported 18 (67%) of the 27 study species models, however, both variables were important for only six (22%) of those 27 species.This is likely an artifact of scale in which either 1) some study sites were not represented well by the meantT4K macroclimate data and/or 2) some species were particularly sensitive to this variable on the microsite scale.Examining multiple scales is a critical, but often overlooked, component of habitat modelling (McGarigal et al., 2016) and our results indicate air temperature on both macro and micro scale improve modelling efforts.

Cool air associate species
Five (C.sanburni, H. camelus, H. skadei, M. mycophaga, and P. wascoense) of the twelve species we identify as cool air associates are currently considered species of conservation concern in our study area based both on distribution and unpublished data which suggested climate sensitivity (IDFG, 2017).Our study confirms this assumption and provides evidence D. whitneyi, E. fulvus, M. ingersolli, P. idahoense, P. humile, R. abietum, and U. lyrata may also be climate sensitive.These species should be evaluated as potential species of conservation concern within the same criteria as the preceding five species along with assessments of their climate adaptive capacity (i.e.Thurman et al., 2020).
With the exception of P. humile and R. abietum these species were typically not detected in valley locations and were restricted to mountainous portions of the study area.M. mycophaga is a large bodied slug (Burke and Leonard, 2013) which, prior to our study, had not been officially documented in the state of Idaho for over half a century (Lucid et al., 2016;Pilsbry, 1953).MeanT and DDLT0 are the most important variables for this species (0.06, Fig. 3).Although meanT has not shown discernable change over the past century, the steady increase in minT (Tinkham et al., 2015) suggests DDLT0 will also increase and likely reduce favorable microclimate conditions for M. mycophaga.D. whitneyi is a small snail with a ribbed shell (Burke and Leonard, 2013).Elevation was not an important variable and we detected this species predominantly at elevations lower than the mean study area elevation.However, it has a low minT requirement (− 26--23 • C) and narrow range of suitable DDGT18 (0-60).C. sanburni is a medium, range restricted (Fig. 2), snail which also has a low minT requirement (-29 -− 27 • C).This argues that the sites suitable for these species are more dependent on physical attributes creating cool microclimate rather than elevation.
E. Fulvus, M. ingersolli, and U. lyrata are relatively common snails and slugs.R. abietum is exceptionally common and, second only to A. kochi, was the most commonly detected species in our survey.Although these species are currently common, they should still be evaluated as potentially climate sensitive species in the context of adaptive capacity (Thurman et al., 2020).
H. camelus and H. skadei are large 'jumping' slugs that have mostly allopatric distributions with well-defined contact zones where limited sympatry occurs (Lucid et al., 2018b;Rankin et al., 2019).That they share the same four modelled variables of importance suggests they occupy similar ecological and microclimate niches in different geographic areas.Although interbreeding between these two species is likely prevented by distinctly different penis morphology, the reasons behind geographic disparity of the two species remains unclear (Lucid et al., 2018b).
P. idahoense and P. wascoense are small snails with overlapping distributions.P. idahoense occurs primarily at higher elevations (also see Hendricks, 2016) while P. wascoense is less elevationally dependent.This is supported by the models which indicate elevation is the most important variable for P. idahoense while unimportant for P. wascoense for which meanT is the most important variable (Fig. 3).
Both of the Hemphillia and Pristiloma species pairs demonstrate members of the same genus using similar microclimate niches.We observed a different pattern in the Prophysaon, large tail-dropper slugs, species pair where P. humile was a cool air associate and P. andersoni is associated with warm air indicating microclimate niche separation.

Warm air associate species
P. andersoni is one of seven warm air associates along with D. reticulatum, H. vancouverense, L. maximus, P. minutissimum, S. pugetensis and V. pellucida.L. maximus is a non-native large bodied slug which is largely restricted to anthropogenically dominated valley locations (Fig. 2).The models for this species were the most strongly supported in our study indicating this importance of warm local air temperatures for this species and suggests it may have the capacity to expand its distribution into more natural areas as the climate warms.

Generalist microclimate species
The generalist species with the highest support was A. kochi, a large and genetically diverse (Rankin et al., 2021) native snail, which was the most commonly detected species in our survey.P. polygyrella, a range restricted snail (Fig. 2), would have met our cool air association cutoff had we allowed meanT4K as a variable in the selection process.The remaining generalist species had a variety of cool, warm, and overlapping variables none of which showed a strong indication of precise temperature requirements.

Microclimate mapping
Scoring sites for microclimate enables the identification of larger clusters of areas which might serve as a buffer to climate change and be the most useful to manage as climate refugia.The largest cluster of − 7 sites is in the northern Selkirk Mountains.This is unsurprising as this area is known to have complex air temperature variability and significant amounts of cold air drainage that occurs nightly from high to low elevations (Holden et al., 2011).The average elevation of the Selkirks changes 600 m every 1 km (Holden et al., 2011) and this extreme topographic variation is likely responsible for the largest cluster of the coolest sites.This extreme relief is not as pronounced in the remainder of the study area but other areas of clustering tend to be in the areas of highest topographic relief for individual mountain ranges.Additional clusters of − 7 sites occur in the northeastern Purcells, central West Cabinets, and southeastern Saint Joe.The Coeur d'Alenes do not show as obvious clustering patterns however there is a weak cluster of cool sites in the northeastern Coeur d'Alenes.This area almost entirely overlaps the range of cool air associate C. sanburni and P. polygyrella (Fig. 2).P. polygyrella did not meet our cut-off for 'cool air associate' becuase only one microclimate variable (DDLT0) was important in the RF models.However, the one important microclimate variable (DDLT0) and one important macroclimate variable (meanT4K) did meet the defined 'cool' cutoffs.Therefore, P. polygyrella does have an association with cool air but not strong enough to meet the cutoffs we established for this study.Exploration of this area as a cool air refugia in the Coeur d'Alenes would be worthwhile.

Climate refugia management
Refugia management is an emerging sub-discipline but guidance, theory, and solid examples exist that demonstrate how refugia might be managed on the ground (Barrows et al., 2020).A wide range of tools and decision making frameworks are available for refugia management (Morelli et al., 2020) ranging from preservation of late successional forests (Krawchuk et al., 2020), selective thinning of younger forests to expand snowpack melt period (Ellis et al., 2013), to topographic alteration to reduce local air temperatures (Greenwood et al., 2016).
The ecological structure of our study area is heavily influenced by past glacial activity and ice-free refugia during the last glacial maximum (Shafer et al., 2010).For instance, the West Cabinet and Coeur d'Alene mountain ranges host relatively low populations American (Martes americana) and Pacific Coast (Martes caurina) Marten which abut at a strong contact zone in the valley where these mountain ranges meet (Lucid et al., 2020a).Western Toads (Anaxyrus boreus) primarily occur in the northern part of the study area and are largely restricted areas predominated by cool air refugia sites in the Selkirk Mountains (Lucid et al., 2020b).We recommend natural resource managers use the coldest microrefugia sites identified in this study as land snail climate refugia.Additionally, we recommend using the clusters of microrefugia we identify, in concert with knowledge of other species ecological requirements, to prioritize larger areas to be managed as wildlife climate refugium.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Study area.Sampling cells were stratified by a 5x5 km grid.Air temperature data loggers were co-located with land snail surveys at 830 sites.Colors represent air temperature microclimate score ranging from − 7 (coolest) to + 7 (warmest).

Fig. 2 .
Fig. 2. Detection maps of each of the 27 land snail species used in this study.See Table3for full genus names.*Species with restricted ranges.For these species we used only data from ranges where they were detected in the RF models.

Table 2
Variables used in Random Forest analysis and median, minimum, and maximum values for microclimate and macroclimate variables.

Table 3
Land snail species, common name, sample size (n = # detection sites used in Random Forest analysis), small (<5mm) or large (>5mm) snail (external shell) or slug (no external shell), native to study area, and microclimate ecological niche assignment (see Table4).

Table 5
Random Forest scaled variables, as defined in Table4, of importance and total number of species for which variable was important.

Table 6
Microclimate scores by mountain range including total number of sites (above) and percentage of sites within each mountain range (below).