Expertly validated models suggest responses to climate change are related to species traits: a phylogenetically-controlled analysis of the Order Lagomorpha

Climate change during the last five decades has impacted significantly on natural ecosystems and the rate of current climate change is of great concern among conservation biologists. Species Distribution Models (SDMs) have been used widely to project changes in species’ bioclimatic envelopes under future climate scenarios. Here, we aimed to advance this technique by assessing future changes in the bioclimatic envelopes of an entire mammalian Order, the Lagomorpha, using a novel framework for model validation based jointly on subjective expert evaluation and objective model evaluation statistics. SDMs were built using climatic, topographical and habitat variables for all 87 species under past and current climate scenarios. Expert evaluation and Kappa values were used to validate past and current distribution models and only those deemed ‘modellable’ through our framework were projected under future climate scenarios (58 species). We then used phylogenetically-controlled regressions to test whether species traits were correlated with predicted responses to climate change. Climate change will impact more than two-thirds of the Lagomorpha, with leporids (rabbits, hares and jackrabbits) likely to undertake poleward shifts with little overall change in range extent, whilst pikas are likely to show extreme shifts to higher altitudes associated with marked range declines, including the likely extinction of Kozlov’s Pika (Ochotona koslowi). Smaller-bodied species were more likely to exhibit range contractions and elevational increases, but showing little poleward movement, and fecund species were more likely to shift latitudinally and elevationally. Our results suggest that species traits may be important indicators of future climate change and we believe multi-species approaches, as demonstrated here, are likely to lead to more effective mitigation measures and conservation management.


Introduction
Changes in climate are predicted to have strong influences on the ecology and distribution of species [1], [2], with pronounced impacts on terrestrial biodiversity [3]. Although the climate changes naturally (and usually slowly), the rate of recent anthropogenically-induced change [4] is causing concern amongst conservation biologists [5]. Future climate change may have large effects on species niches i.e. the biotic and abiotic conditions in which a species can persist [6]. Species are predicted to adapt their bioclimatic niche, migrate to maintain their current niche, or become range restricted and undergo population decline, local or global extinction under future scenarios [7].
The Lagomorpha are an important mammalian Order economically and scientifically as a major human food resource, model laboratory animals, valued game species, pests of agricultural significance and key elements in food chains providing scientific insights into entire trophic systems. Lagomorphs are native on all continents except Antarctica, occurring from sea level to >5,000m and from the equator to 80°N spanning a huge range of environmental conditions, and also include some very successful invasive species [8].
The taxonomy of the Lagomorpha in recent decades has been in a state of flux but all species belong to two families: the Ochotonidae and the Leporidae. The Ochotonidae consists of one monotypic group in the genus Ochotona containing 25 species of small, social, high-altitude pikas. The Leporidae has 32 species of large, solitary, cursorial hares and jackrabbits in a single genus Lepus and 30 species of medium-sized, semi-social, fossorial rabbits comprising ten genera [8]. A quarter of lagomorphs are listed in the IUCN Red List of Threatened Species (www.iucnredlist.org) with a notable number of highly range-restricted species including fourteen listed under the IUCN Criteria B, with an extent of occurrence estimated to be less than 20,000 km 2 . In addition, pikas as high-altitude specialists with very high body temperatures of 39.3 -41.0°C [9] are extremely susceptible to changes in their environment, particularly ambient temperatures [10].
Species Distribution Models (SDMs) are widely used in ecology and relate species occurrences at known locations to environmental variables to produce models of environmental suitability, which can be spatially or temporally extrapolated to unsurveyed areas and into past or future conditions (e.g. [11]). Although SDMs have been highly influential in the field of ecology, their limitations have been widely reviewed (e.g. [12]). The impact of climate change on species distributions has been modelled in a wide range of studies and a number of responses have been described (e.g. [1], [5], [13]). Mammalian distributional changes have been well studied over the past decade and indicate that future climate change will have profound impacts (e.g. [5], [14], [15], [16]). Mammals in the Western hemisphere are unlikely to keep pace with climate change, with 87% expected to undergo range contractions [16], and mammals in Mediterranean regions, particularly endemic species, are predicted to be severely threatened by future climate change [15]; shrews are especially vulnerable to future changes [5], [15]. These distributional responses have also been noted in studies of past climatic changes, for example, Moritz et al. [14] found that from the early 20 th century to the present, small mammals in a North American national park substantially shifted their elevational range upwards corresponding to ~3°C increase in minimum temperatures.
The predicted impact of climate change on species distributions has only rarely been linked with species traits. Yet, species traits are widely accepted as potentially important indicators of responses to climate change and identifying such traits may be crucial for future conservation planning (e.g. [14], [17], [18]). Traits that directly impact climatic conditions experienced by a species, for example their activity cycle, are likely to be more important in mediating species responses to projected climate change than traits such as diet breadth. If species can broaden their occupied bioclimatic niche through trait plasticity, for example, altering their diel patterns of activity, then they may be less susceptible to future change [19].
Mammalian species active during certain times of the day will experience a limited range of climatic conditions, whereas more flexible species can select the conditions in which they are active [20], and may, therefore, be less susceptible to future change [19]. Small body size, nocturnal behaviour and burrowing may have allowed mammalian species to 'shelter' from climatic changes during the 66M year duration of the Cenozoic era [21]. Burrowing rabbits may be able to adapt to climate change by seeking underground refugia from extreme or fluctuating temperatures [22], whilst larger cursorial species, such as the hares and jackrabbits which, in the majority of cases, live above ground may have less variability in microclimate opportunities within which to shelter [23]. Thus, species with narrow environmental tolerances, poor dispersal ability and specialised habitats may be more susceptible to climate change [24]. Dispersal is also likely to be very important in future species distributions, with larger species being more mobile and, therefore, potentially better prepared to track climatic changes [16], Past studies have modelled the response of large numbers of species to predicted climate change [3] or dealt with a few key species at Order-level [16], lagomorphs due to their restricted diversity provide an opportunity to rigorously examine the response of every species yielding a detailed picture of change within an entire Order for the first time.
Crucially, the small number of lagomorph species, compared to other mammalian groups, means that datasets can be verified in detail, modelled individually and outputs expertly validated. Moreover, lagomorphs have a nearly global terrestrial distribution and occupy a wide range of biomes providing an opportunity to examine the response of similar species from tundra to desert and islands to mountain summits. 6 Here, we assess the projected change in the bioclimatic envelopes of all 'modellable' lagomorph species under future climate change using a framework for model validation based jointly on subjective expert evaluation cross-validated with objective model evaluation statistics. We predict lagomorph species distributions will increase in elevation and poleward movement under future climate change, but with significant differences between pikas, rabbits, hares and jackrabbits due to dissimilarities in species traits, for example body size.
Lagomorph morphological and life history traits will be correlated with the predicted responses to future climate change in order to test this hypothesis. We theorise that flexibility in activity cycle and larger body sizes, which may lead to greater mobility, will result in species being less vulnerable to future climatic changes and better able to track climate niches.

Species data
A total of 139,686 records including all 87 lagomorph species were either downloaded from the Global Biodiversity Information Facility (GBIF) Data Portal (data.gbif.org), collated from species experts or members of the IUCN Species Survival Commission (SSC) Lagomorph Specialist Group (LSG) and/or extracted from the literature for data deficient species as advised by experts. All past and current occurrence records, sorted by species, can be viewed on http://lagomorphclimatechange.wordpress.com/. Taxonomic accuracy was dealt with by checking all records against the latest IUCN taxonomy; if names did not match after crossreferencing with taxonomic synonyms and previous names they were rejected. Spatial data accuracy was dealt with by removing any obviously erroneous records for the target species if they fell outside the extent of the IUCN geographic range polygon. In addition, occurrences recorded with a spatial resolution of >2km were removed and duplicate records were eliminated unless they were recorded in different temporal periods (pre-1950 and post-1950).
This left 41,874 records of which 3,207 were pre-1950 and 38,667 were post-1950. In this study, spatial and temporal bias in sampling was eliminated by only selecting the background data (a random sample of the environmental layers describing pseudo-absences), 10,000 points, from sites at which any lagomorph species had been recorded and ensuring there were the same proportion of pre-1950 to post-1950 background points as there were pre-1950 to post-1950 species records.
combined to produce a single model reflecting the current classification of the Irish hare as an endemic sub-species.

Model evaluation
A bespoke website (lagomorphclimatechange.wordpress.com) was created to allow each species expert to review the output of their allocated species. Expert evaluation, whereby an acknowledged expert on each species judges model predictions for current and past distributions, can be a useful tool prior to making future extrapolations [30]. A framework combining expert evaluation with reliable model evaluation metrics allows species distribution models to be assessed before they are used in future projections to infer likely future changes in distribution.
Forty-six lagomorph experts, including 20 members of the Lagomorph Specialist Group (LSG) of the IUCN Species Survival Commission (SSC) and 26 recognised lagomorph researchers selected based on recent publications, were paired to species (see Table S1 in Supporting Information) and asked to assess whether model projections accurately, roughly or did not capture the current and past range of each species i.e. good, medium or poor respectively according to the criteria in Anderson et al. [30]. Experts were asked to select the most accurate representation of the current and past range from the following models: i) pre- Independent model evaluation in SDM studies often uses the Area Under the Curve (AUC) value but this has been heavily criticised [31]. AUC is not advocated for model evaluation because Receiver Operating Characteristic (ROC) curves cannot be built for presence/absence or presence/pseudo-absence data [32] and AUC values can be influenced by the extent of model predictions [31]. There are also known limitations with using alternative metrics [32] such as sensitivity (proportion of presences which are correctly predicted), specificity (proportion of absences which are correctly predicted) or True Skill Statistic (a prevalence independent model metric calculated using sensitivity and specificity).
Misleading commission errors, which can arise from species not being at equilibrium with the environment, can affect such metrics. On the other hand, Kappa is an objective measure of prediction accuracy based on input species records and background points adjusted for the proportion of correct predictions expected by random chance [33]. Kappa utilises commission and omission errors [34], and although it does not take into account prevalence like the True Skill Statistic [32] and can sometimes produce misleading commission errors [35], it has widely accepted thresholds which can be useful in model evaluation [36], [37], [38]. There are documented flaws with all possible model evaluation statistics, however Kappa was chosen due to the reasons above and its common use in model evaluation (e.g. [39]), with the main focus of evaluation in this study being the expert validation approach.
The 'accuracy' function in the SDMTools package [40] in R (version 3.0.2) was used to calculate the Area under the Curve (AUC) value, omission rate, sensitivity, specificity, proportion correct, Kappa and True Skill Statistic for completeness, but only Kappa was taken forward for use in the model validation framework. The optimum threshold was taken as Kappa>0.4 as this value has been advocated in a range of previous studies [32], [36], [38].
Models that had a Kappa >0.4 and those that were ranked as either good or medium by expert evaluation were defined as "modellable" because they demonstrate good model fit and predictive ability; these species were carried forward for projection and further analysis.
Those models with a Kappa <0.4 or ranked as poor by expert evaluation were defined as "unmodellable", with poor model fit and predictive ability, and were rejected from further analysis (see Figure S1). Hereafter, species are referred to as "modellable" or "unmodellable" as explicitly defined above.

Future predictions
The model settings, for example, the input data used (pre-or post-1950) or restriction/ no restriction to occupied land classes, which provided the optimal outcome i.e. the highest Kappa value and best expert evaluation were used to project modellable species bioclimatic envelopes under future climate during the 2020s, 2050s and 2080s. Future predictions were clipped to the buffered minimum convex polygon of the target species further buffered by the dispersal distance (kilometres/year) of each species multiplied by the number of years elapsed from the present (1950-2000) taken as 1975. Predicted range size, mean latitude and minimum, mean and maximum elevation for each species and each time period were calculated. Model outputs for each species can be found in Appendix 2; unmodellable species projections are included only for reference.

Species traits
Species trait data were downloaded from the PanTHERIA database [41] and updated by searching the literature. Correlated traits were removed to reduce multicollinearity (i.e. tolerance <2 and VIF >5) and the final set of traits used to describe each modellable species were: activity cycle (nocturnal only, diurnal only, flexible), adult body mass (g), diet breadth (number of dietary categories eaten from: vertebrates, invertebrates, fruit, flowers/nectar/pollen, leaves/branches/bark, seeds, grass and roots/tubers), gestation length (days), habitat breadth (above ground-dwelling, aquatic, fossorial and ground-dwelling), home range (km 2 ), litters per year, litter size, population density (n/km 2 ) and age at sexual maturity (days).

Statistical analysis
To illustrate projected changes in the distribution of lagomorph species, the difference in predicted species richness per cell was calculated between the 1930s (1900-1949) and 2080s . The difference in model output metrics (range size, mean latitude and minimum, mean and maximum elevation) was calculated between the 1930s and 2080s.
Change in range size was expressed in percentage change but change in latitude was represented as degree movement towards the poles and change in elevation in metres.
Generalised Least Square (GLS) models, with an autoregressive moving average (ARMA) correlation structure, in the 'nlme' package in R version 3.0.2 [42] were used to test the differences in temporal trends for range change, poleward movement and elevation between lagomorph groups: 1) pikas, 2) rabbits and 3) hares and jackrabbits. ARMA was used to explicitly account for the non-independent nature of the time-series periods. Phylogenetic Generalized Least-Squares (PGLS) regressions were performed to test whether changes in model predictions varied with morphological or life history traits. PGLS analysis was run using the 'caper' package in R [43]. A lagomorph phylogeny was extracted from the mammalian supertree provided by Fritz et al. [44]. Likely clade membership for five species not included in this phylogeny was determined from Ge et al. [45], and then missing tips were grafted on using an expanded tree approach [46]. Outliers were removed prior to analysis, because species with large residuals may overly influence the results of the regression, and they were identified as those with a studentized residual >3 units following Cooper et al. [47]. All subsequent models exhibited normally distributed residuals tested using Shapiro-Wilk. The scaling parameter lambda varies from 0, where traits are independent of phylogeny, to 1 where species' traits covary in direct proportion to their shared evolutionary history [48]. We estimated lambda for each model and tested whether it was significantly different from 0 or 1 during the PGLS analysis. All subset regressions were 14 run using the 'dredge' function in the 'MuMIn' [49] package in R, using AICc as the rank estimate, and then model averaging was used to describe the effect of each variable.

Results
There was a high degree of agreement between expert model classification and Kappa values, however experts were often more critical because, for example, they classified 25 species as 'poor' but these species had a mean Kappa of ~0.6 ( Fig. 1). Fifty-eight species (67%) were deemed modellable with expert validation classed as medium or good and Kappa >0.4 and 29 species (33%) rejected as unmodellable with expert validation classed as poor and/or Kappa <0.4. Unmodellable species were 4 times more likely to be listed by the IUCN as Data Deficient than modellable species, with 8 unmodellable species (28%) listed as threatened.
The majority of species with small sample sizes were classed as unmodellable; the median number of occurrence points for modellable species was 36 and for unmodellable 13.
Hereafter, all results are for modellable species only and, therefore, are a highly conservative estimate of the impact of climate change on the Order.
Global changes in predicted lagomorph species richness suggest that almost a third of the Earth's terrestrial surface (31.5 million km 2 ) is predicted to experience loss of lagomorph species by the 2080s (Fig. 2). The desert regions of North-eastern China and hills of Sichuan become increasingly unsuitable, potentially losing up to 10 species by 2080, including the woolly hare (Lepus oiostolus) and Glover's pika (Ochotona gloveri) which are predicted to undergo dramatic movements to higher elevations. In contrast, predicted species gains were notable across: (a) northern Eurasia, due to poleward movement of the mountain hare (Lepus timidus); and, (b) North America, where some regions e.g. the Upper Missouri catchment of Montana and North Dakota were predicted to gain up to 5 species. This includes the desert and eastern cottontails (Sylvilagus audubonii and S. floridanus respectively) which are predicted to exhibit strong poleward movements. The majority of African lagomorph species were classed as unmodellable and as such Fig. 2 is largely data deficient for the continent.
Pikas exhibited the most substantial mean increase in elevation becoming increasingly isolated on mountain summits (e.g. the Rockies in North America and the Tibetan Plateau and high Himalayas) resulting in a significant 31% range contraction (F df=1,234 =3.7, p=0.03).
In addition, lagomorph species occupying islands (n=6) will, on average, lose 8km 2 of their ranges compared to 0.2km 2 gain for continental species (n=52), whilst mountain dwelling species (n=24) will lose 37km 2 of their ranges compared to 25km 2 gain for lowland species (n=34).
PGLS models indicate that members of each group were capable of showing a variety of responses i.e. species of pika, rabbits, hares or jackrabbits exhibited both increases and decreases in each response variable (Fig. 4). All traits used in the PGLS models are described in the methods and full results can be found in Table S3. Here we present only the significant results. There was a significant positive relationship between predicted range change and adult body mass (β=0.258, F df=4,52 =2.308, p=0.021) with hares and jackrabbits generally increasing their range by 2080 and pikas exhibiting range contraction ( Fig. 4; also see Table   S3). Adult body mass was positively associated with predicted poleward movement were found between activity cycle and predicted changes in species distribution.

Discussion
The Lagomorpha as a whole are predicted to exhibit much greater poleward (mean ± 95%CI = 1.1° ± 0.5°) and elevational shifts (165m ± 73m) by 2100 than calculated in a recent metaanalysis collating information on a wide variety of taxonomic groups [1]. On average birds, butterflies and alpine herbs were predicted to shift ~0.8° poleward and increase in elevation by ~90m by 2100. Thuiller et al. [5] found that mammals were less vulnerable to change than other groups, but more than 50% of shrews could lose more than 30% of their suitable climatic space. In comparison, our study shows that only 34% of modellable lagomorphs will lose more than 30% of their suitable climatic space, but lagomorphs appear to show more notable changes in elevation and poleward movement. A study by Schloss et al. [16] found that, in some scenarios, 50% of mammal species in North and South America will be unable to keep pace with future climate change. Our results, therefore, indicate that lagomorphs follow similar trends to other mammalian climate change studies, but with more substantial poleward and elevational shifts. Furthermore, our results are conservative estimates due to the exclusion of unmodellable species, most notably African species. Parts of Africa are expected to become drier and warmer under future climate, with substantial increases in arid land [4] which will likely lead to negative consequences for lagomorphs not assessed here.
In contrast, a recent study of empirical climate studies by McCain & King [19] has shown that only about 50% of, mostly North American, mammal species respond to climate change by shifting in latitude and elevation. Results of four lagomorph studies are included which show that the pygmy rabbit (Brachylagus idahoensis) undergoes extirpation and contraction due to climatic changes [50], the range of the snowshoe hare (Lepus americanus) and collared pika (Ochotona collaris) does not change [51], [52] and the American pika (Ochotona princeps) mostly undergoes extirpation and upslope contraction, but some sites within the range show no change [14], [53]. Although we find that more than 50% of lagomorphs respond to climate change, our results are largely congruent with these empirical studies, indicating range contraction in the American pika and pygmy rabbit, and little change in the ranges of the snowshoe hare and collared pika.
Pikas are predicted to show elevational rather than latitudinal shifts because these highaltitude specialists are known to be extremely susceptible to increases in temperature or unpredictable seasonality, which could potentially lead to heat stress [9]. On the other hand leporids, typically being lowland species, exhibited less substantial increases in elevation but greater poleward shifts, for example, the mountain hare (Lepus timidus) which is predicted to shift poleward. This is probably due to the high sensitivity of its boreo-alpine niche to changes in temperature [54], [55]. Indeed, even the Irish hare (L. t. hibernicus) which inhabits temperate grasslands, unlike other mountain hares, is predicted to experience a contraction in the south-east of Ireland. As global temperatures increase, northern latitudes will become more climatically suitable for southern leporids and, therefore, species bioclimatic envelopes will track poleward to match. Rabbits, hares and jackrabbits were thus predicted to exhibit little overall change in the total extent of their ranges. The majority of the modellable rabbit species were from the genus Sylvilagus inhabiting south and eastern USA and Mexico; a region projected to become generally drier [56] inducing latitudinal shifts in species to track suitable habitats or vegetation. Thus, by shifting their ranges poleward, leporids are predicted to be able to maintain or increase their range extent suggesting they are considerably less sensitive to projected climate change than pikas.
PGLS models estimated lambda to be close to zero for changes in range, poleward movement and elevation indicating that observed trends were independent of phylogeny [48].
Smaller species, principally pikas and some rabbits, were typically more likely to make elevational shifts in range whilst larger species, principally hares and jackrabbits, had a greater tendency for poleward shifts. This relationship was also discovered in empirical 20 climate studies (e.g. [19]) and is probably due to the large number of small-bodied species that dwell in mountainous regions (n=24) which have more opportunity to shift altitudinally.
A 1m elevational shift is equivalent to a 6km latitudinal shift [1], so it is easier for smaller species to shift in elevation rather than latitudinally. This could be explained due to a relationship between adult body mass and dispersal distance, which was used to clip model projections, but a Spearman's Rank correlation suggested no correlation between these two variables in our dataset (see Table S4).
Trait analyses also showed a significant positive relationship between fecundity and extreme elevational or poleward movement which is comparable to a study by Moritz et al. [14], who found that fecund small mammals in Yosemite National Park, USA, were more likely to expand their ranges upward than less fecund counterparts. We found no association between activity cycle and vulnerability to climate change in the Lagomorpha as hypothesized and reported in McCain & King [19], but this may be due to the varied nature of lagomorph activity. The activity of the European hare is known to be less consistent and partially diurnal in the summer [57] and the American pika also shows seasonal changes in activity patterns [58]. However, we do acknowledge the potential drawbacks of linking traits to modelled distributional changes rather than empirical-based studies, i.e. non-independence of trait results; nevertheless the analysis presented here enables understanding of the traits which could potentially lead to vulnerability to future climate change.
The uncertainties in SDM using projected future climate scenarios are well described [12].
Models are vulnerable to sampling bias [59], spatial scale issues [60], lack of data for rare species [26], uncertainty regarding future climate conditions [61] and insufficient independent evaluation [60]. We have tackled these explicitly by accounting for sampling bias by restricting the set of background points used, using data with the highest spatial resolution available (30 arc-second) and selecting species records to match, bootstrapping models for rare species with few records, averaging climatic data from five global circulation models and using a framework of joint expert and metric-driven model validation to segregate current distribution model outputs into those that were modellable and unmodellable before subsequent projection and analysis. Nevertheless, our predictions could still be potentially confounded by species-area relationships [62] and biological interactions [63]. Regardless of individual model outcomes (see Appendix 2), the overall trends observed across the Order Lagomorpha as a whole are compelling. This study did not take account of shifts in habitats, vegetation or human impact in response to climate change, but we have shown that adaptation to future climate conditions may be possible as some species were predicted to exhibit poleward movements, with only modest shifts in range or elevation e.g.

eastern (Sylvilagus floridanus) and Appalachian cottontail (Sylvilagus obscurus).
The predicted changes in climate conditions are likely to have greater impacts on isolated lagomorph populations, i.e. those on islands and in high-altitude systems. If changes continue at the rate currently predicted until the 2080s, then there may be no climatically suitable range available for some montane or isolated species e.g. the Tres Marias cottontail (Sylvilagus graysoni) or black jackrabbit (Lepus insularis). Conservation strategies, such as assisted migration, could be the only option for these highly range-restricted species. Furthermore, conservation management will need to focus on small-bodied mammals as these are predicted to show more dramatic responses to changing climate. Small mammals are key in food webs sustaining predator populations, impacting plant communities by grazing and soil biology and hydrology by burrowing. Thus, fundamental shifts in lagomorphs globally may cause trophic cascade effects, especially in northern latitudes such as the cyclic systems of the Arctic.
The advancing knowledge of past extinction rates for the Lagomorpha [45], along with significantly better bioclimatic envelope modelling, could aid the prediction and prevention of future extinctions. Our models suggest that Kozlov's pika (Ochotona koslowi) may become extinct by the 2080s as the elevational increases required to maintain its current bioclimatic envelope disappear as it reaches the maximum elevation available. We have shown here how expert validation can be effectively integrated into the model evaluation process in order to improve model predictions and advocate use of this framework in future SDM studies. Assessment of vulnerability to climate change is needed urgently to identify how and to what extent species, taxa, communities and ecosystems are susceptible to future changes, taking into account likely changes to vegetation, human pressures and interspecific interactions, and to direct conservation management in an efficient and effective manner.
Multi-species approaches are likely to lead to more effective mitigation measures and contribute to our understanding of the general principles underpinning the biogeographical and ecological consequences of climate change impacts.            Annual evapotranspiration was taken as the sum of all monthly values. Annual water balance was calculated by subtracting annual evapotranspiration from mean annual precipitation. The number of months with a positive water balance was calculated by subtracting each monthly evapotranspiration from its corresponding monthly precipitation, and then converting these into a binary format, where a value greater than zero was given a value of one and a value less than zero was kept at zero [2]. The twelve binary files were then summed to calculate the number of months with a positive water balance.
Mean annual Normalised Difference Vegetation Index (NDVI) was calculated from monthly values which were downloaded from the EDIT Geoplatform (http://edit.csic.es/Soil-Vegetation-LandCover.html). NDVI is commonly used to measure the density of plant growth and is obtained from NOAA AVHRR satellite images. A negative value indicates snow or ice, a value around 0 indicates barren areas, values between 0.2 and 0.3 indicate grassland, and values near 1 indicate rainforests [3]. Human influence index was downloaded from the SEDAC website (http://sedac.ciesin.columbia.edu/data/set/wildareas-v2-humaninfluence-index-geographic). This was a composite of human population density, railways, roads, navigable rivers, coastlines, night-time lights, built-up areas and agricultural and urban land cover. Values within the index range from 0 to 64, where zero equalled no human influence and 64 represented maximum human influence [4]. Solar radiation was calculated using the Spatial Analyst function in ArcGIS 10.1 (ESRI, California, USA). Solar radiation is defined as the total amount of incoming solar insolation (direct and diffuse), or global radiation, and was measured in watt hours per square meter or WH/m 2 [5]. An index of surface roughness was also calculated by finding the difference between maximum and minimum gradient values, based on a global Digital Elevation Model at 30 arc-second resolution downloaded from WorldClim [6]. Altitude was not included as a variable independently because organisms perceive climatic and habitat variables as proxies for altitude [7].

Future climate data
Pierce et al. [8] report that using data averaged across five global circulation models (GCMs) is substantially better than any one individual model and significantly reduces model error. We averaged all variables for the following five GCMs: CCCma-CGCM3.1/T47, CNRM-CM3, CSIRO-Mk3.0, HadCM3 and NASA-GISS-ER. Future projections adopted the time periods: 2020s (2010-2039), 2050s (2040-2069) and 2080s (2070-2099). Data from the A2 future climate change scenario was used because although it was originally described as "extreme climate change" it now appears to best represent the trend in observed climate. The Pygmy rabbit's bioclimatic envelope is predicted to decline by 87% with a 1 o mean latitudinal poleward shift and mean increase in elevation of ~300m driven predominately by an increase in mean minimum elevation (>600m) with little change in mean maximum elevation (~50m). 95% of the permutation importance of the model was contributed to by mean annual temperature (64.5%), maximum temperature (25.2%) and annual water balance (5.9%).  The Riverine rabbit's bioclimatic envelope is predicted to decline by 85% with a ~1 o mean latitudinal poleward shift and mean increase in elevation of ~200m driven by similar increases in both minimum and maximum elevation. 95% of the permutation importance of the model was contributed to by minimum temperature (33.1%) and precipitation (22.5%), temperature seasonality (22.3%), mean annual temperature (6.7%) and precipitation (5.0%), annual evapotranspiration (4.3%) and precipitation seasonality (2.6%).

Summary:
The Hispid hare's bioclimatic envelope is predicted to increase by 21% with a ~1.5 o mean latitudinal poleward shift and mean increase in elevation of ~70m driven by increases in maximum elevation. 95% of the permutation importance of the model was contributed to by mean annual temperature (52.7%), precipitation seasonality (29.0%), annual evapotranspiration (6.6%), number of months with a positive water balance (2.9%), maximum precipitation (2.2%) and minimum temperature (1.7%).

Summary:
The Japanese hare's bioclimatic envelope is predicted to increase by 9% with no latitudinal poleward shift and a mean increase in elevation of ~10m driven by a decrease in minimum elevation. 95% of the permutation importance of the model was contributed to by temperature seasonality (31.8%), annual water balance (20.6%), human influence index (13.5%), mean annual precipitation (8.9%), maximum precipitation (6.2%), precipitation seasonality (5.4%) and minimum precipitation (4.4%).

Summary:
The Iberian hare's bioclimatic envelope is predicted to increase by 40% with a ~1 o mean latitudinal poleward shift and a mean increase in elevation of ~10m driven by an increase in maximum elevation. 95% of the permutation importance of the model was contributed to by maximum precipitation (39.0%), annual evapotranspiration (38.0%), minimum temperature (15.0%) and maximum temperature (3.0%).

Summary:
The Kozlov's pika's bioclimatic envelope is predicted to decrease by 100% (total extinction). 95% of the permutation importance of the model was contributed to by mean annual precipitation (66.8%), minimum precipitation (24.2%) and human influence index (5.8%).

Summary:
The Large-eared pika's bioclimatic envelope is predicted to decrease by 40% with a ~1 o mean latitudinal shift towards the Equator and a mean increase in elevation of ~880m driven by an increase in minimum elevation. 95% of the permutation importance of the model was contributed to by minimum precipitation (87.3%), minimum temperature (5.8%) and mean annual temperature (2.1%).

Summary:
The Nubra's pika's bioclimatic envelope is predicted to increase by 1% with no latitudinal polewards shift, but a mean increase in elevation of ~200m driven by an increase in minimum and maximum elevation. 95% of the permutation importance of the model was contributed to by minimum precipitation (76.1%) and mean annual temperature (20.9%).

Summary:
The Pallas's pika's bioclimatic envelope is predicted to decrease by 60% with a ~2 o mean latitudinal polewards shift and a mean increase in elevation of ~40m driven by an increase in minimum elevation. 95% of the permutation importance of the model was contributed to by minimum precipitation (46.1%), mean annual precipitation (35.5%), mean annual temperature (9.4%) and minimum temperature (6.7%).

Summary:
The American pika's bioclimatic envelope is predicted to decrease by 25% with a ~1 o mean latitudinal polewards shift and a mean decrease in elevation of ~10m driven by an decrease in minimum elevation. 95% of the permutation importance of the model was contributed to by mean annual temperature (88.1%), annual evapotranspiration (5.0%) and maximum temperature (2.8%).

Summary:
The Afghan pika's bioclimatic envelope is predicted to increase by 5% with a ~1 o mean latitudinal polewards shift and a mean increase in elevation of ~380m driven by an increase in maximum elevation. 95% of the permutation importance of the model was contributed to by minimum precipitation (89.2%), minimum temperature (5.5%) and normalised difference vegetation index (2.3%).

Summary:
The Turkestan red pika's bioclimatic envelope is predicted to decrease by 10% with a ~1 o mean latitudinal shift towards the Equator and a mean increase in elevation of ~630m driven by an increase in maximum elevation. 95% of the permutation importance of the model was contributed to by minimum precipitation (82.5%), minimum temperature (5.9%), human influence index (3.0%), precipitation seasonality (2.5%) and surface roughness index (1.7%).        Status: MODELLABLE; Included in final analysis: √ Summary: The San Jose brush rabbit's bioclimatic envelope is predicted to decrease by 25% with a no latitudinal polewards shift and a mean increase in elevation of ~30m driven by an increase in minimum elevation. 95% of the permutation importance of the model was contributed to by temperature seasonality (36.7%), human influence index (34.0%), minimum precipitation (21.8%) and precipitation seasonality (3.2%).

Summary:
The Marsh rabbit's bioclimatic envelope is predicted to increase by 90% with a ~1° mean latitudinal polewards shift and a mean increase in elevation of ~2m driven by an increase in maximum elevation. 95% of the permutation importance of the model was contributed to by surface roughness index (73.1%), minimum precipitation (18.3%) and precipitation seasonality (4.9%).

Summary:
The Robust cottontail's bioclimatic envelope is predicted to decrease by 90% with a ~4° mean latitudinal polewards shift and a mean increase in elevation of ~480m driven by an increase in minimum elevation. 95% of the permutation importance of the model was contributed to by precipitation seasonality (68.5%), minimum precipitation (8.2%), temperature seasonality (7.1%), minimum temperature (4.5%), mean annual precipitation (3.6%), annual water balance (2.6%) and number of months with a positive water balance (1.6%).

Summary:
The Venezuelan lowland rabbit's bioclimatic envelope is predicted to decrease by 100% with a ~1.5° mean latitudinal polewards shift and a mean increase in elevation of ~275m driven by an increase in minimum elevation. 95% of the permutation importance of the model was contributed to by temperature seasonality (97.7%).   Figure S1. Framework for assessing whether species were "modellable" or 5 "unmodellable" based on Kappa values and expert evaluation classification. The 6 optimum threshold for Kappa was taken as 0.4 [9], [10], [11]. Expert evaluations were 7 classified according to Anderson et al. [12].   Table S3. Results of phylogenetically-controlled generalised least square regressions. 16 Significant p values for model-averaged coefficients are in bold. F and p values for the top 17 model are listed under each response variables; asterisks (*) indicate traits in the top model. 18 Lambda (λ) confidence intervals and significance from 0 and 1 are also shown.