Aim: Geographic variation in metabolic resources necessary for functional trait expression can set limits on species distributions. For species that need to produce and maintain biomineralized traits for survival, spatial variation in mineral macronutrients may constrain species’ distributions by limiting the expression of biomineralized traits. Here, we examine whether threatened, heavily biomineralized Oreohelix land snails are restricted to CaCO₃ rock regions, if they incorporate greater amounts of CaCO₃ rock carbon in their shell than less biomineralized smooth forms, and if ornamentation increases shell strength.

Location: Western United States

Methods: We used random forest (RF) classification models at multiple spatial resolutions to evaluate the contribution of topographic, vegetation, climate, and geologic variables in predicting the presence of heavily biomineralized shell ornaments. We then measured and compared shell biometric variables, ¹⁴C/¹²C ratios, and peak force for fracture for ornamented and smooth forms from calcareous and non-calcareous substrates.

Results: Distance to CaCO₃ rock was the most important variable in all trait distribution models and was highly associated with local ornamentation classification and forecasted distribution. Pairwise comparisons of ¹⁴C/¹²C ratios in closely occurring ornamented vs. smooth population pairs revealed ornamented forms incorporate greater CaCO₃ rock carbon than smooth forms. Ornamented types measured in this study were generally heavier and required greater peak force for fracture than smooth snails, except when compared to smooth forms sampled from CaCO₃ rock.

Main conclusions: Biomineralization expression, species distribution, and trait function appear to be constrained by mineral supply in a highly diverse group of land snails. This trait-environment relationship suggests similar CaCO₃ macronutrient constraints may modulate biomineralization expression and restrict species distribution in other terrestrial molluscs and has a direct impact on the management of Oreohelix species.

RF predictors:

Predictor Dataset Creation

The predictors used in this study came from a variety of sources (Supplementary Table A1). In this section, we will detail how they were made to facilitate replication of our results. All predictors were reprojected in ArcGIS Pro v.2.6.0 to WGS1984 and clipped to the same raster resolution. Predictor names used in the R code are shown in parentheses. See Supplemental Table A1 for references.

• Elevation (elevation): This layer was sourced from the publicly available ASTER Global Digital Elevation data reprojected to 90m resolution using the Project tool and clipped to the desired extent using the Clip Raster tool.
• Slope (slope): This layer was created using the Slope tool in ArcGIS on the 90m elevation data using a z-factor of 0.00001171 appropriate for 40 degrees latitude (https://pro.arcgis.com/en/pro-app/latest/tool-reference/3d-analyst/applying-a-z-factor.htm) which is close to the mean latitude of our study area.
• Compound topographic index (CTI): This layer was created using the GradientMetrics ArcGIS toolbox (https://github.com/jeffreyevans/GradientMetrics) compound topographic index tool on the 90m resolution elevation data. 
• Heat load index (HLI): This layer was created using the GradientMetrics ArcGIS toolbox (https://github.com/jeffreyevans/GradientMetrics) heat load index tool on the 90m resolution elevation data.
• Height distance above nearest drainage (HAND): This layer was created by first creating a stream layer using the Con tool in ArcGIS pro on the perennial and intermittent streams present in the USGS National Hydrography Dataset (https://www.usgs.gov/core-science-systems/ngp/national-hydrography). Then, we used the Flow Distance tool using the stream layer and 90m elevation data as inputs and set the function to measure only vertical flow distance.
• Horizontal distance to nearest drainage e(DTND): This layer was created by first creating a stream layer from the perennial and intermittent streams present in the National Hydrography Dataset (https://www.usgs.gov/core-science-systems/ngp/national-hydrography) using the Con tool in ArcGIS pro. Then, we used the Flow Distance tool with the stream layer and 90m elevation data as inputs and set the function to measure only horizontal flow distance.
• Global horizontal irradiance (GHI): This layer was sourced from the Global Solar Atlas (Solar Atlas 2020). No modifications were made other than reprojecting and clipping was necessary.
• Soil ph (soilph): We downloaded tiled POLARIS probabilistic mean soil ph data from across the western USA from http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/ph/mean/0_5/.
• Tiled data were combined using the Mosaic to New Raster tool in ArcGIS Pro. Missing data were filled in based on a sliding window average of 20 x 20 pixels using the Cell Statistics tool and the Mosaic to New Raster tool.
• Soil clay content (soilclay): We downloaded tiled POLARIS probabilistic mean soil clay content data across the western USA from http://hydrology.cee.duke.edu/POLARIS/PROPERTIES/v1.0/clay/mean/0_5/.

Tiled data were combined using the Mosaic to New Raster tool in ArcGIS Pro. Missing data were filled in based on a sliding window average of 20 x 20 pixels using the Cell Statistics tool and the Mosaic to New Raster tool.

• Normalized Difference Vegetation Index (NDVI): We used modified SEBALIGEE Google Earth Engine code (Mhawej and Faour 2020; original code used in Mhawej and Faour here: https://code.earthengine.google.com/48200ed2b76ff4acc530c618bb047635; code used in this paper is provided on Dryad) to take the mean NDVI across the month of July for the entirety of LANDSAT8’s available data (2013 – 2020). We chose the month of July to gather NDVI data as cloud cover artifacts are usually less prevalent in summer months, and this time of the year represents some of the greatest extremes in temperature/desiccation land snails will experience. Regional segments were combined using the Mosaic to New Raster tool in ArcGIS Pro.
• LANDSAT8 Surface Temperature (LST): We used modified SEBALIGEE Google Earth Engine code (Mhawej and Faour 2020; original code used in Mhawej and Faour here: https://code.earthengine.google.com/48200ed2b76ff4acc530c618bb047635; code used in this paper is provided on Dryad) to download regional mean NDVI across the month of July for the entirety of LANDSAT8’s available data (2013 – 2020). We chose the month of July to gather NDVI data as cloud cover artifacts are usually less prevalent in summer months, and this time of the year represents some of the greatest extremes in temperature/desiccation land snails will experience. Regional segments were combined using the Mosaic to New Raster tool in ArcGIS Pro.
• Distance to nearest developed area (developed): We used the Con tool in ArcGIS Pro on the National Land Cover Database (https://pubs.er.usgs.gov/publication/fs20123020; Homer et al. 2012) pasture, crops, moderately developed, and highly developed areas to create a layer of developed areas. We then input the developed area as the feature of interest in the Path Distance tool and the cost raster of equal values (in this case, 1) for the entire Western USA. This value was then multiplied by the distance of the original projection (90m) to generate real distance away from developed areas.
• Distance to CaCO₃ rock (LS): We used the Con tool in ArcGIS pro on the National Karst Map’s (https://pubs.usgs.gov/of/2014/1156/; Weary and Doctor 2014) near-surface carbonate layers from wet and dry environments to create our initial CaCO₃ rock layer. This layer was then improved locally in the Salmon and Snake river area (Idaho, USA) with carbonate layers from Kauffman et al. (2014) as this information was already available from a preliminary study. The two CaCO₃ rock layers were combined using the Mosaic to New Raster tool to create the final CaCO₃ rock layer. We then input the CaCO₃ rock areas as the feature of interest in the Path Distance tool and the cost raster of equal values (in this case, 1) for the entire Western USA. This value was then multiplied by the distance of the original projection (90m) to generate distance away from CaCO₃ rock.

Dataset Composition

Number of records for classification and regression models were taken from a combination of government and private surveys. All records are less than 500m spatial uncertainty. All records on Dryad are obfuscated to greater than 1km in accordance with state and federal government data use agreements.

See publication for crushing/radiosiotope data.

ArcGIS Pro/QGIS to modify layers

R for scripts