Combining past and contemporary species occurrences with ordinal species distribution modeling to investigate responses to climate change

Many organisms leave evidence of their former occurrence, such as scat, abandoned burrows, middens, ancient eDNA, or fossils, which indicate areas from which a species has since disappeared. However, combining this evidence with present-day occurrences within a single modeling framework remains challenging. Traditional binary species distribution modeling reduces occurrence to two temporally vague states (present/absent), so thus cannot leverage the information inherent in temporal sequences of evidence of past occurrence. In contrast, ordinal modeling can use the natural time-varying order of states (e.g., never occupied vs. occupied in the past vs. currently occupied) to provide greater insights into range shifts. We demonstrate the power of ordinal modeling for identifying the major influences of biogeographic and climatic, variables on current and past occupancy of the American pika (Ochotona princeps), a climate-sensitive mammal. Sampling over 5 years across the species’ southernmost, warm-edge range limit, we tested the effects of these variables at 570 habitat patches where occurrence was classified either as binary or ordinal. The two analyses produced different top models and predictors – ordinal modeling highlighted chronic cold as the most-important predictor of occurrence, whereas binary modeling indicated primacy of average summer-long temperatures. Colder wintertime temperatures were associated in ordinal models with higher likelihood of occurrence, which we hypothesize reflect longer retention of insulative and meltwater-provisioning snowpacks. Our binary results mirrored those of several other past pika investigations employing binary analysis. Because both ordinal- and binary-analysis top models included climatic and biogeographic factors, results constitute important considerations for climate-adaptation planning. Cross-time evidences of species occurrences are common, yet remain underutilized for assessing responses to climate change. Compared to multi-state occupancy modeling, which presumes all states occur in the same time period, ordinal models enable use of historical evidence of species’ occurrence to identify factors driving species’ distributions more finely across time.

Methods
Our dataset includes n = 570 habitat patches surveyed by >1 of the authors and/or field technicians and determined to be currently occupied, previously occupied, or having no evidence of occupancy by the study organism, the American pika (Ochotona princeps). The dataset spans 5 years, from 2016 to 2020 inclusive. Latitude and longitude were recorded using hand-held GPS devices accurate to 2-7 m (WAAS-enabled), and “survey Year” is the year in which the site was most recently surveyed. “Ordinal status” is parameterized as 0 (no evidence of pika occurrence), 1 (old evidence of pika occurrence), or 2 (evidence of current pika occurrence). Binary status is given as either 0 (not currently pika-occupied) or 1 (current pika occurrence detected). Our study domain is divided into northeast, southeast, northwest, and southwest “sub-regions” of mountain ranges that are separated from each other by ~45-100 km of lower-elevation topography. We recorded “elevation” (in meters) using handheld GPS units. The number of home ranges (“numHomeRanges”) is given for each site and was estimated in situ by observers, aided by laser rangefinders. If more than one value was given for the number of home ranges by different observers, then we used the mean value. The mean distance to the closest four patches is given in meters (m) and was calculated as the mean geographic distance to each of four nearest patches, which we used to determine patch-level isolation.