Arctic biodiversity: increasing richness ... - Wiley Online Library

Arctic biodiversity: increasing richness accompanies shrinking refugia for a cold-associated tundra fauna A. G. HOPE,1,6, E. WALTARI,2 J. L. MALANEY,3 D. C. PAYER,4 J. A. COOK,5

AND

S. L. TALBOT1

1

U.S. Geological Survey, Alaska Science Center, 4210 University Drive, Anchorage, Alaska 99508 USA 2 City College of New York, Department of Biology, Marshak Science Building 814, 160 Convent Avenue, New York, New York 10031 USA 3 Department of Natural Resources and Environmental Science, University of Nevada–Reno, Reno, Nevada 89557 USA 4 Arctic Landscape Conservation Cooperative, 101 12th Avenue, Room 216, Fairbanks, Alaska 99701 USA 5 Museum of Southwestern Biology and Department of Biology, University of New Mexico, Albuquerque, New Mexico 87131 USA Citation: Hope, A. G., E. Waltari, J. L. Malaney, D. C. Payer, J. A. Cook, and S. L. Talbot. 2015. Arctic biodiversity: increasing richness accompanies shrinking refugia for a cold-associated tundra fauna. Ecosphere 6(9):159. http://dx.doi. org/10.1890/ES15-00104.1

Abstract. As ancestral biodiversity responded dynamically to late-Quaternary climate changes, so are extant organisms responding to the warming trajectory of the Anthropocene. Ecological predictive modeling, statistical hypothesis tests, and genetic signatures of demographic change can provide a powerful integrated toolset for investigating these biodiversity responses to climate change, and relative resiliency across different communities. Within the biotic province of Beringia, we analyzed specimen localities and DNA sequences from 28 mammal species associated with boreal forest and Arctic tundra biomes to assess both historical distributional and evolutionary responses and then forecasted future changes based on statistical assessments of past and present trajectories, and quantified distributional and demographic changes in relation to major management regions within the study area. We addressed three sets of hypotheses associated with aspects of methodological, biological, and socio-political importance by asking (1) what is the consistency among implications of predicted changes based on the results of both ecological and evolutionary analyses; (2) what are the ecological and evolutionary implications of climate change considering either total regional diversity or distinct communities associated with major biomes; and (3) are there differences in management implications across regions? Our results indicate increasing Arctic richness through time that highlights a potential state shift across the Arctic landscape. However, within distinct ecological communities, we found a predicted decline in the range and effective population size of tundra species into several discrete refugial areas. Consistency in results based on a combination of both ecological and evolutionary approaches demonstrates increased statistical confidence by applying cross-discipline comparative analyses to conservation of biodiversity, particularly considering variable management regimes that seek to balance sustainable ecosystems with other anthropogenic values. Refugial areas for cold-adapted taxa appear to be persistent across both warm and cold climate phases and although fragmented, constitute vital regions for persistence of Arctic mammals. Key words: boreal-tundra ecotone; climate change; community turnover; conservation; ecological niche prediction; small mammals; species distribution model; statistical phylogeography; wildlife management. Received 20 February 2015; revised 6 March 2015; accepted 23 March 2015; final version received 3 May 2015; published 28 September 2015. Corresponding Editor: D. P. C. Peters. Copyright: Ó 2015 Hope et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. http://creativecommons.org/licenses/by/3.0/ 6

Present address: Division of Biology, 116 Ackert Hall, Kansas State University, Manhattan, Kansas 66506 USA.

E-mail: [email protected]

v www.esajournals.org

1

September 2015 v Volume 6(9) v Article 159

HOPE ET AL.

et al. 2007, Eidesen et al. 2013). Communities now occurring within Beringia reflect constituent species, or distinct lineages within species, that have evolutionary origins elsewhere in either Eurasia or North America, as well as taxa that diverged and adapted in situ. Community turnover in Beringia through time has thus been pervasive (Post et al. 2009, Eidesen et al. 2013) and this dynamic region connecting processes throughout the Holarctic represents a valuable arena for investigating how species with diverse ecological requirements and evolutionary histories respond to dramatic climate fluctuations (Callaghan et al. 2004). Here we test hypotheses related to regional changes in biodiversity, a critical step for understanding both global atmospheric influences and impacts across distinct biomes that result in novel communities (Barnosky et al. 2012). We apply genetic and niche-based approaches in a comparative context, reflecting 15 years of phylogeographic studies across 28 species of high-latitude mammals to statistically address hypotheses concerning (1) consistency of inferences across both ecological and evolutionary analyses; (2) differences in total diversity versus diversity within distinct biomes; and (3) different management implications reflecting heterogeneity in landscape characteristics and land uses. The questions that frame these hypotheses, as outlined below, are associated with assessing demographic change of multiple species (Hope et al. 2013b, Wood et al. 2013) in response to lateQuaternary climate trends as a means of building confidence in forecasts of contemporary trends. Biogeographic responses to climate change may be estimated by implementing climatebased distribution models (Bellard et al. 2012, Hof et al. 2012) or through interpretations based on evolutionary signatures (Avise 2000, Hewitt 2004). Increasingly, coupling species distribution models (SDMs) with testing alternative phylogeographic scenarios allows for a more rigorous exploration and synthesis of both ecological and evolutionary approaches (Carnaval et al. 2009, Knowles and Alvarado-Serrano 2010, Chan et al. 2011, Lorenzen et al. 2011, Prost et al. 2013). Development of comparative phylogeographic methods that incorporate coalescent theory also facilitates increased analytical power through simultaneous analyses of multiple co-distributed

INTRODUCTION Prediction is difficult, especially if it’s about the future. —Neils Bohr Understanding how biodiversity at high latitudes, including species richness and associated genetic diversity, will respond to rapidly changing environments is critical considering that Arctic climate predictions suggest continued warming into the next century (Diffenbaugh and Giorgi 2012, Overland et al. 2013). Warmer and wetter conditions are shifting organismal distributions and changing demography (Prowse et al. 2009, Pauls et al. 2013). For instance, poleward latitudinal climate shifts will decrease the extent of the tundra biome as the Arctic Ocean delimits a northern barrier (Hof et al. 2012). Known as the Arctic squeeze (Meltofte 2013), decreasing area may translate to fewer individuals, decreasing b diversity, loss of genetic variants, and potentially decreased adaptive potential for tundra specialists (Hofreiter and Barnes 2010, Pauls et al. 2013, Yannic et al. 2014). As a consequence, the Arctic tundra biome is of conservation concern because non-migratory species associated with this relatively young and depauperate biome are adapted to cold environments that have prevailed at high latitudes through much of the Quaternary (2.6 Myr; Abbott and Brochmann 2003, Callaghan et al. 2004). Historical biogeographic studies revealed that Arctic species persisted in refugial areas through previous climatic warm phases (e.g., Fedorov and Stenseth 2002, Abbott and Brochmann 2003, Brunhoff et al. 2003, Hewitt 2004, Hope et al. 2013a), but continued persistence under current climate trends may depend on our ability to locate and accommodate refugial regions for Arctic biodiversity in an increasingly anthropic landscape (Chapin et al. 2006, Kuemmerle et al. 2014, Raynolds et al. 2014). In the northern refugium known as Beringia, large areas spanning northeastern Asia and northwestern North America remained ice-free throughout the Quaternary (Hulteń 1937). During glacial climate phases, the Beringian refugium represented both a center of endemic genetic diversity and a corridor for movement between the northern continents (Cook et al. 2005, Waltari v www.esajournals.org

2


HOPE ET AL.

taxa (Hickerson et al. 2010, Palsbøll et al. 2013). Our first question therefore assesses the effectiveness of integrating multiple disciplines for comparative community analyses. Specifically, is there consistency among predicted changes in biodiversity using alternate statistical perspectives from ecological and evolutionary data analyses? Making fundamental links between ecological responses of species and evolutionary processes associated with environmental change should provide a powerful approach for investigating changes in biodiversity. We incorporate data from both SDMs and genetic demographic analyses across multiple communities to reciprocally verify and refine inferences of distributional and demographic changes in biodiversity. We show that integrating methods can significantly increase confidence in forecasting biodiversity responses to climate, producing robust distributional hypotheses (predictions) that may be empirically tested and refined into the next century—a valuable management need. Since onset of the current interglacial phase (;11 ka), tundra has receded and the temperateadapted boreal biome has shifted northward (Lessa et al. 2003, Hewitt 2004). Secondary contact between boreal and tundra biomes in Beringia during Quaternary interglacial phases resulted in novel species interactions across ecotonal areas (Stewart and Lister 2001, Hof et al. 2012). External climate forcing is moving the current boreal-tundra ecotone northward (Frost and Epstein 2014) and creating challenges for conservationists and land managers who seek to maintain ecological function and promote overall biodiversity while protecting endemic tundra components (Usher et al. 2005). Therefore, our second question is what can we expect in the near-future concerning changes in total biodiversity at high latitudes as related to differential changes among distinct tundra (Arctic) and boreal (sub-Arctic) communities? This question reflects the possibility of fundamental differences in responses among communities with different ecological affinities, particularly with respect to their relative velocities of change (Schloss et al. 2012). We therefore test whether boreal and tundra mammals are responding to contemporary climate with H0: the same, or H1: different velocities of change. Marked differences between boreal and tundra mammalian responses emv www.esajournals.org

phasize different life histories among communities that are experiencing increasing overlap through time, with potentially severe implications for future resilience of tundra species as communities reform. Extended persistence of tundra species through multiple glacial phases indicates resiliency within this biome to climate cycling through the Quaternary and suggests that at least some elements of tundra communities will continue to persist in the face of modern changes. However, the most recent Wisconsinan glacial cycle (;110–11 ka) was distinctive because humans became an integral part of the biosphere, significantly influencing wildlife abundance and distribution through hunting and land-use changes (Duerden 2004, Cook et al. 2006, Bellard et al. 2012). Atmospheric changes and associated climate warming of the Anthropocene have surpassed levels of previous Quaternary interglacial maxima, and in the future the northern high latitudes will be characterized by novel environments (Overland et al. 2013, Reu et al. 2014). In addition to global-scale climate changes, Arctic amplification is prevalent at northern high-latitudes and is also accompanied by increased localized activities related to industrial development, harvest of biological resources, science, and tourism (Duerden 2004, Hinzman et al. 2005, Raynolds et al. 2014). The footprint of these activities includes areas recognized as warm-phase refugia for tundra species and essential for persistence through their evolutionary history (Hope et al. 2013b). Within Alaska (constituting a large portion of eastern Beringia), 89% of the land area (;134M ha) is administered by multiple state or federal agencies with significant differences in mission reflected by differing management regimes (http://www. alaskacenters.gov/partners.cfm). Understanding the ecological and spatiotemporal evolutionary responses of species to past and contemporary climate changes will help inform future planning. Balancing potentially competing objectives, for instance of wildlife conservation and resource extraction (Ford and Smit 2004, Kuemmerle et al. 2014), will be complicated by significant uncertainty concerning the long-term stability of ecosystems (Hinzman et al. 2005, Barnosky et al. 2012). Our final question assesses how knowledge of shifting communities may influ3


HOPE ET AL.

ence land-use decisions across a landscape characterized by multiple units with differing objectives. Specifically, what are the differences in major regional management implications considering a broad-scale community approach? We tested whether distributional changes among major regional land-use areas across the Alaskan landscape are expected to be H0: homogeneous, reflecting a broad cross-site trend in biodiversity change that can be extrapolated through the Holarctic, or H1: heterogeneous, indicating a need for careful consideration of regional processes, and to assess related communities distributed across the Holarctic independently. This study provides a comprehensive mammalian comparative phylogeographic assessment that integrates SDMs with genetic demographics over a suite of 28 species to provide both individual and whole-community assessments. We utilize a molecular dataset and analyze all taxa simultaneously to calibrate evolutionary rates for more accurate inference of relative timing and both SDMs and genetic demographics have been refined by incorporating reciprocal knowledge from ecological and evolutionary standpoints.

MATERIAL

AND

also often translate to rapid demographic and evolutionary change (Bromham 2011).

Species distribution modeling For SDM development, we used current, past, and future monthly climate data at 2.5 0 (4 km) spatial resolution, using bioclimatic variables from the WorldClim data set (Hijmans et al. 2005), version 1.4. Both future and Last Glacial Maximum (LGM) climate data are from CMIP5 and were based on three General Circulation Model simulations: NCAR-CCSM, version 4.0, MIROC, version ESM 2010, and MPI, version ESM-LR, while Last Interglacial (LIG) data were based only on the CCSM model (Otto-Bliesner et al. 2006). Future GCM data were used based on the rcp85 scenario, a business-as-usual atmospheric projection, at two time intervals: 2050s, and 2070s. Compiled geospatial occurrences per species ranged from 21-221 (mean 102). To correct for potential sampling bias which can lead to over-fit models (Reddy and Da´valos 2003, Veloz 2009, Kramer-Schadt et al. 2013), we refined occurrence data by removing redundant localities within 20km, and increased thinning up to 120km for well-sampled species (Appendix: Table A1). To test for consistency of model predictions among taxa given disparity in number of locality records across taxa, we chose one taxon from boreal (Sorex cinereus) and tundra (Lemmus trimucronatus) biomes that were represented by numerous records, and then randomly selected 20, 30, 40, and 50 locality records for each, producing SDM predictions for each random dataset for all timeframes as described for other taxa. We used correlative modelling to explore niche-based variation and assess temporal trends in northern biodiversity. We used Maxent, version 3.3.3k, to construct present-day SDMs and for temporal projections. Recently the inherent difficulties associated with spatiotemporal projection of species distributions (Anderson and Gonzalez 2011, Peterson et al. 2011) have resulted in criticisms of SDMs, particularly when SDMs over-fit to present day climates are projected across timescales (Davis et al. 2014, Worth et al. 2014). In light of these developments, we used a series of techniques based on operational criteria from the literature to ensure models were not overly fit to the present

METHODS

Study system and species Our investigation centers on the mammals of northern Alaska, a fauna that extends broadly throughout the Holarctic. Because the range of some taxa shifted into or out of this region through time or potentially could in the future, specimens used for both niche and genetic analyses in this study included distinct lineages with distributions that extend beyond Alaska’s borders in western Beringia (i.e., Far East Siberia) and far eastern Beringia (i.e., Yukon Territory and northern British Columbia; Appendix: Table A1). The Alexander Archipelago of Southeast Alaska was excluded from sample representation as many specimens from that region represent endemic lineages that reflect complex and unresolved post-glacial phylogeographic histories (Cook et al. 2006). Small and medium-sized mammals are ideal focal taxa as they are yearround residents, relatively abundant, and extensive sampling is feasible within this logistically challenging region. Their short generation times v www.esajournals.org

4


HOPE ET AL.

conditions, but rather constituted more objective estimates of fundamental niches (Peterson et al. 2011). In addition to the thinning of occurrences, our present-day SDMs were developed in masks bounding each lineage in question (Appendix: Table A1) to avoid inclusion of regions where species are absent for non-climatic reasons (Anderson and Raza 2010, Barve et al. 2011), for example the presence of congeners. Correlated variables can also result in overly complex models, so we reduced the 19 bioclimatic variables by removing redundancy (variables with a Pearson’s correlation coefficient R 2 . 0.8). Species-specific parameter tuning often enhances model performance (Anderson and Gonzalez 2011) and can greatly reduce overfitting, therefore we conducted an analysis of model selection (Warren and Seifert 2011) for each species to assess the importance of the bregularization parameter. This parameter is a model penalty proportional to variance of features observed at presence localities but is also inversely proportional to the square root of sample sizes (Merow et al. 2013). We optimized the values (1, 3, 5, 7, 9, 11, 13, and 15) using the model selection process in ENMTools v. 1.3 (Warren et al. 2010) but held features constant. Following model optimization, we projected the initial present-day SDM with limited extent first onto the entire Holarctic region for present climate, then to both past and future periods for each species. Each model prediction is an average of five replicates under the ‘‘crossvalidate’’ option. We used thresholds to represent a binary distribution of our results that summarize species-specific models across multiple future climate predictions. Whenever possible we used LIG or LGM fossil records as calibrations for models, using a minimum fossil presence threshold for the ‘hindcast’ models. When fossil records were unavailable, we used present-day SDMs to identify optimal thresholds that included minimum training presence for shrew species and equal training sensitivity and specificity thresholds for all others. Importantly, SDMs were created for all taxa considering only the distribution of distinct genetic units of analysis (lineages). As a consequence, uncertainty surrounding distributional predictions is controlled for by incorporating knowledge of the evolutionary history of study taxa. In general, genetic v www.esajournals.org

lineages included are well defined by geographic extent, making inclusion of only representative occurrence data for each lineage possible. These genetically and geographically refined datasets, including many individual specimens used in both genetic and SDM analyses, also helped to ensure that inclusion of mis-identified specimen records (particularly common for cryptic taxa such as Sorex, but rarely assessed when refining SDM datasets) was minimized. Three taxa (the continental lineage of Tamiasciurus hudsonicus and two distinct lineages of Mustela erminea) were excluded from SDM projections as definitive lineage associations of specimen occurrences could not be substantiated.

Genetic demographics Genetic analyses utilized sequences of variable lengths (541–1143bp) of the mitochondrial cytochrome b gene that we retrieved from GenBank, for a total of 1005 specimens representing 24 of 28 taxa (Table 1). Sample sizes for independent taxa, each representing a distinct lineage for analysis, ranged from 6 to 176 and sequence length was held constant within each taxon. GenBank accession numbers and associated references to original phylogeographic studies are listed in Hope et al. (2014). Considering that taxa with disparate phylogeographic histories likely responded differentially to widespread environmental trends, we used DnaSP v5 (Librado and Rozas 2009) to summarize and test population demography with Tajima’s D (Tajima 1989) and Ramos-Onsins and Rozas’ R2 (Ramos-Onsins and Rozas 2002). We assessed the power of these summary statistics with 10,000 coalescent simulations. These statistics consider the number of segregating sites and are sensitive to both large and small sample sizes, respectively. We assumed panmixia given that individuals of each taxon constituted well-supported clades without substantial intraclade substructure (Sta¨dler et al. 2009) and no evidence of inter-clade gene flow, based on independently published phylogeographic assessments. Estimation of relative timing of demographic changes in response to past climates is critical for interpreting the direction of responses (increase vs. decline) to climatic trends and provides a basis for predicting future population trajecto5


HOPE ET AL. Table 1. Summary statistics (n ¼ sample size; l ¼ length in base pairs of the mtDNA Cyt b gene), demographic expansion statistics (D ¼ Tajima’s D; R2 ¼ Ramos-Onsins and Rozas R2), and predicted changes in size of range through time. Species distribution model predictions indicate positive or negative change in predicted percentage occupancy within the study area for each taxon since the Last Glacial Maximum, and up to the 2070s timeframe. Percentages were used in order to accommodate changes in sea level and associated changes in total size of the study area between timeframes. Groups include T ¼ tundra, B ¼ boreal, and W ¼ widespread. SDM predictions and percent change (%) Species

Group

n

l (bp)

D (P)

R2 (P)

LGM–present

Present–2070s

Dicrostonyx groenlandicus Lemmus trimucronatus Microtus longicaudus Microtus miurus Microtus oeconomus Microtus pennsylvanicus Myodes rutilus Synaptomys borealis Zapus hudsonius Zapus princeps Marmota broweri Marmota caligata Urocitellus parryii Tamiasciurus hudsonicus (N) Tamiasciurus hudsonicus (S) Ochotona collaris Lepus americanus Lepus othus Gulo gulo Mustela erminea (H) Mustela erminea (C) Sorex cinereus Sorex hoyi Sorex minutissimus Sorex monticolus Sorex navigator Sorex tundrensis Sorex ugyunak

T T B T W B W B B B T B W B

68 18 16 81 78 27 26 n/a 28 28 49 6 n/a 28

1075 789 1143 1072 1143 672 1140 n/a 1140 1140 1140 1140 n/a 771

1.89 (0.007) 1.27 (0.089) 1.95 (0.012) 1.30 (0.087) 1.69 (0.018) 1.71 (0.023) 2.19 (0.003) n/a 2.06 (0.003) 1.91 (0.010) 0.09 (0.559) n/a n/a 2.33 (,0.001)

0.039 (0.007) 0.078 (0.011) 0.061 (,0.001) 0.060 (0.099) 0.051 (0.039) 0.054 (,0.001) 0.038 (,0.001) n/a 0.055 (0.003) 0.051 (,0.001) 0.106 (0.508) n/a n/a 0.039 (,0.001)

26.4 (82.9–56.5) þ20.3 (13.2–33.4) þ22.8 (0.9–23.7) 44.7 (85.9–41.2) þ24.4(42.4–66.8) þ44.3 (4.5–48.7) þ39.1 (24.4–63.5) þ28.3 (2.2–30.5) þ22.6 (11.0–33.6) þ0.2 (0.0–0.2) 45.7 (68.9–23.2) 0.1 (10.3–10.2) 12.3 (48.6–60.9) þ30.2 (3.2–33.3)

37.4 (56.5–19.1) þ11.1 (33.4–44.5) þ42.1 (23.7–65.8) 30.7 (41.2–10.5) þ2.9 (66.8–69.8) þ36.5 (48.7–85.2) þ28.5 (63.5–92.0) þ47.2 (30.5–77.7) þ52.7 (33.6–86.3) þ3.7 (0.2–4.0) 22.9 (23.2–0.3) 6.9 (10.2–3.3) 31.5 (60.9–29.4) þ60.1 (33.3–93.4)

B

25

771

2.46 (,0.001)

0.049 (,0.001)

n/a

n/a

W B T T W B B B W B B W T

176 n/a n/a 18 29 25 114 20 27 35 44 27 68

1119 n/a n/a 1140 541 541 1075 1011 1056 801 1140 1050 1075

1.86 (0.006) n/a n/a 1.16 (0.135) 1.98 (0.007) 1.81 (0.011) 2.68 (,0.001) 1.68 (0.032) 2.17 (0.002) 2.28 (0.001) 1.85 (0.016) 2.17 (0.002) 1.89 (0.007)

0.036 (0.007) n/a n/a 0.116 (0.160) 0.051 (,0.001) 0.086 (0.065) 0.015 (,0.001) 0.061 (,0.001) 0.046 (,0.001) 0.038 (,0.001) 0.049 (0.005) 0.044 (,0.001) 0.039 (0.007)

0.9 (64.5–63.7) þ38.1 (13.3–51.4) þ6.5 (0.2–6.7) þ14.2 (40.2–54.4) n/a n/a þ40.7 (14.5–55.2) þ48.5 (24.1–72.6) þ55.3 (24.7–80.1) þ33.8 (62.4–96.2) þ30.6 (22.9–53.5) þ56.9 (33.6–90.5) þ5.9 (36.6–42.5)

35.4 (63.7–28.3) þ4.2 (51.4–55.6) þ8.1 (6.7–14.8) 33.7 (54.4–20.7) n/a n/a þ43.9 (55.2–99.0) þ17.4 (72.6–90.1) þ1.1 (80.1–81.2) þ3.8 (96.2–100.0) þ36.8 (53.5–90.3) þ0.8 (90.5–91.3) 4.4 (42.5–38.1)

ries. The efficacy of Bayesian skyline plots (BSPs) for estimating demographic histories has been demonstrated, particularly for detecting prolonged changes in effective population size (Ne; e.g., Ho and Shapiro 2011, Palsbøll et al. 2013). However, accuracy of timing relies on precise taxon-specific substitution rates (Ho et al. 2008). Relative substitution rates among multiple taxa, and subsequent taxon-specific BSPs, were estimated previously (Hope et al. 2014). BSPs were used to provide estimates of Ne through time and develop alternative scenarios for hypothesis testing using coalescent simulations. Species with limited genetic sampling were excluded from tests of demographic expansion (Table 1), but were included in BSP estimates considering the added rigor of coalescent simulations for testing relative significance of BSP trends. v www.esajournals.org

A continuing limitation of comparative phylogeographic methods is a lack of comparative multi-locus data across multiple taxa in a region. Therefore to statistically assess the strength of observed demographic changes considering broad credible intervals from single locus BSPs, we utilized approximate Bayesian computation (ABC; Beaumont et al. 2002) to test competing demographic hypotheses (scenarios; Appendix: Fig. A1). All competing scenarios considered single non-bifurcating evolutionary lineages, for which one million simulated datasets were generated in the program DIYABC v2.0 (Cornuet et al. 2014). All parameters (Appendix: Table A3) were assigned a uniform prior (interval of 60.5 3 parameter value). Mutation rates for all scenarios were based on empirical values co-estimated for each taxon (Appendix: Table A3). Transition/ 6


HOPE ET AL.

transversion ratios (K) and proportion of invariant sites (PINV) represent mean values from empirical model-test results or were set to default values where no empirical values were provided. Statistical evaluation of simulated data from each set of models included pre-evaluation of the suitability of model-prior combinations, and computation of posterior probabilities of models using single point logistic regressions with default values. Confidence in scenarios was also tested using single point logistic regression methods by simulating 30,000 pseudo-observed datasets for the top two competing scenarios based on results of posterior probability model tests and drawing parameter values from distributions. The first set of three scenarios considered the shape of BSPs (i.e., how population sizes change through time; ‘‘BSP scenarios’’) for each taxon where values of Ne and coalescent-based time estimates from empirical median BSP curves were used as mean prior values for a Null scenario of demographic change (Appendix: Table A3). For competing scenarios, we reversed the demographic BSP curve along the y-axis from the Null scenario, and finally tested a scenario of no change in average population size through time (Appendix: Fig. A1). The second set of three scenarios tested relative statistical support for observed recent population size changes from BSPs post-LGM, considering the curves often have broader credible intervals near the present (‘‘LGM scenarios;’’ Drummond et al. 2005). Mostrecent historic population size trends were therefore used to loosely infer potential future trajectories. Although future trends are difficult to accurately predict, contemporary genetic signatures compared over multiple species translates to multiple lines of evidence for concerted responses to the most recent episode of warming climate, thereby increasing confidence of continued future trajectories among common community components. Importantly, the LGM scenarios statistically incorporated information from SDM predictions for each taxon. The percentage change in predicted distributional extent between present and 2070s timeframes was used as the percentage contemporary increase or decrease in effective population size for ABC scenario tests in order to assess if such population changes are feasible given the most v www.esajournals.org

recent genetic demographic trends. The Null scenario therefore utilized parameter values from empirical BSPs up to the maximum estimated Ne, and then increased Ne following the LGM (coincident with the most recent interglacial warming trend; Appendix: Fig. A1) by a percentage equal to predicted range expansion (or contraction) from SDMs (Table 1, values in present-to-2070s column). If these values were less than 20% predicted change, or if SDM predictions were not available for a taxon, then a standard 20% increase/decrease was used for ABC scenarios. The two competing scenarios included a corresponding decrease in Ne from the maximum population size prediction, and a scenario of no change. For taxa where SDM predictions indicate only minimal change in geographic extent, we would expect a scenario of no recent change in population size as most likely based on coalescent simulations. Interpretation of ABC simulations is therefore statistically dependent on both empirical sequence data and distributional predictions based on georeferenced specimen records, reflecting strongly integrated eco-evolutionary methods.

Community assessment among management regions We investigated community-level response to late-Quaternary climate and potential future climate change considering all taxa simultaneously, and also by splitting the fauna into ecological groups that may reflect common spatial and temporal phylogeographic responses associated with the major biomes (boreal forest and alpine/tundra). As an example of both regional differences in shifting patterns of species richness and to inform multiple land-management regimes with differing priorities, we focused on three U.S. federal management regions associated with northern portions of the study area (Fig. 1). The Arctic Network of National Parks (ARCN), administered by the U.S. National Park Service, extends from the Seward Peninsula in the west to the central Brooks Range. The mission of the National Park Service is to preserve unimpaired the natural and cultural resources and values of the National Park System. The ARCN comprises the Bering Land Bridge National Preserve (1.1M ha), Cape Krusenstern 7


HOPE ET AL.

Fig. 1. Land use, ecoregions, and species distribution model (SDM) predictions for all taxa. (A) Distributions of the four largest federal land-use administrations within Alaska. The three quantified management regions here


8


HOPE ET AL. (continuation of Fig. 1 legend) include the Arctic Network of National Parks (ARCN; administered by the NPS; square), Arctic National Wildlife Refuge (Arctic NWR; administered by the USFWS; circle) and National Petroleum Reserve–Alaska (NPR-A; administered by the BLM; triangle). Management region borders and assigned shapes (circle, square, triangle) are included on all plots for comparison of regional distribution of diversity through time (Fig. 3). (B) Distribution of major ecoregions within the study area (Level I; Environmental Protection Agency). SDMs were based on (C) Last Interglacial (LIG), (D) Last Glacial Maximum (LGM), (E) current (Now), (F) 2050s, and (G) 2070s climate projections. The color gradient reflects areas of low (blue) to high (red) richness. LIG and current predictions of diversity reflect 25 summed SDM layers (one per taxon). LGM and future projections reflect 25 summed SDM layers, with each of the 25 taxa predicted by between 0 and 3 models (resulting in each pixel ranging from 0-75).

scheme to highlight variation in richness from low (cool blues) to high (warm reds) diversity.

National Monument (0.3M ha), Gates of the Arctic National Park and Preserve (3.4M ha), Kobuk Valley National Park (0.7M ha), and Noatak National Preserve (2.7M ha). The legislative purposes of each land-management unit varies, but in general include protecting and interpreting changes in habitats and populations of fish, wildlife and plants. The second federal management region is the Arctic National Wildlife Refuge (Arctic NWR), a 7.8M ha conservation area in northeastern Alaska that is administered by the U.S. Fish and Wildlife Service. Purposes of the Arctic NWR include conservation of fish and wildlife populations and habitats in their natural diversity, and providing continued opportunities for traditional harvest of those resources. Both the Arctic NWR and the ARCN span the boreal-tundra ecotone and include coastal plain lowlands as well as alpine and taiga environments. The third federal-management region is the National Petroleum Reserve–Alaska (NPR-A), administered by the U.S. Bureau of Land Management. The 9.2M ha NPR-A is distributed across the Arctic coastal plain and northern Brooks Range foothills and is legislatively managed as multi-use, including fish and wildlife conservation, and energy development. We evaluated changes in number of species through time from SDMs, and then predicted changes in the distribution of richness within each federal-management region based on climate predictions from the present and 2070s timeframes to produce a quantitative measure of relative diversity through time. Finally, we used Arc Calculator in ArcGIS 10.0 to generate total potential diversity maps by overlaying speciesspecific threshold-based predicted distributions for each period and for each ecological association. Resulting ensemble maps used a continuous color v www.esajournals.org

RESULTS Consistency among ecological and evolutionary analyses We performed a range of analyses to investigate species-level ecological and evolutionary responses to climate change through time. SDMs, genetic demographic summary statistics and ABC coalescent simulations represented both ecological and evolutionary methods for estimating population changes from archived specimens that were assessed both independently and also refined by reciprocally incorporating data into analyses. SDMs.—In general, SDM procedures performed well with AUC values consistently .0.70 with optimized b-regularizations applied (Appendix: Table A2). Spatial projections of the SDMs generally predicted population expansion for species associated with the boreal biome, but minimal change or predicted contraction of range size for tundra and widespread taxa that in most cases represent Beringian endemics (Fig. 2; Appendix: Figs. A2, A3). Increase in mean annual temperature was consistently the most important bioclimatic variable explaining range expansions in 15 of the 18 boreal and widespread taxa. However the variable of most importance was not consistent among tundra species (six important variables across seven species; Appendix: Table A1). Tundra species variably respond to increasing temperature (G. gulo, M. broweri ), diurnal temperature range (Lepus othus), highest temperature of the warmest month (Microtus miurus), isothermality (Sorex ugyunak), warmest quarter temperature (Dicrostonyx groenlandicus), or seasonality of precipitation (L. trimucronatus). 9


HOPE ET AL.

Fig. 2. Species distribution model predictions for boreal taxa (left column; blue) and tundra taxa (right column; red) based on (A) and (B) current (Now), (C) and (D) 2050s, and (E) and (F) 2070s climate projections. The color gradient reflects areas of low (light) to high (dark) diversity. Distribution of diversity reflects compilation of predictive maps for 12 boreal and 7 tundra/alpine associated species, each thresholded as presence/absence. Black lines provide a guideline to the extent of 50% of the diversity within boreal and tundra biomes respectively during each timeframe.

Current predictions for all species coincided closely with their known distributions based on specimens (Appendix: Fig. A2). Optimal Beta selection, AUC scores and threshold levels were independent of number of locality samples used v www.esajournals.org

for each taxon, suggesting no effect of sampling on robustness of SDM predictions (Appendix: Tables A1, A2). In addition, tests to assess the effect of varying sampling from 20 to 50 randomly sampled data points had little effect 10


HOPE ET AL. Table 2. Results of ABC scenario tests indicating the proportion of simulated datasets that supported each scenario based on the logistic regression method. Logistic regressions considered 30,000 closest datasets to empirical data reported as a single point-estimate where the highest values indicate the scenario that produced most similar simulations to the empirical data for each taxon. The first three competing ABC scenarios tested the relative strength of Bayesian skyline plot predictions of population size change through time (individual Bayesian skyline plots provided in Appendix: Fig. A4). BSP1 scenario reflects the empirical BSP; BSP2 scenario reflects an inverse of the empirical BSP; BSP3 scenario reflects no change from an average population size through time. The second three scenarios tested for the best-supported recent population size trends since the Last Glacial Maximum that may be indicative of near-future population trends. LGM1 scenario reflects growth (þ) since the LGM; LGM2 scenario reflects decline () since the LGM; LGM3 scenario reflects no change (¼) in population size since the LGM. ABC simulations Bayesian skyline scenarios of dNe

Post-LGM scenarios of dNe

Species

BSP 1

BSP 2

BSP 3

LGM 1 (þ)

LGM 2 ()

LGM 3 (¼)

Dicrostonyx groenlandicus Lemmus trimucronatus Microtus longicaudus Microtus miurus Microtus oeconomus Microtus pennsylvanicus Myodes rutilus Synaptomys borealis Zapus hudsonius Zapus princeps Marmota broweri Marmota caligata Urocitellus parryii Tamiasciurus hudsonicus (N) Tamiasciurus hudsonicus (S) Ochotona collaris Lepus americanus Lepus othus Gulo gulo Mustela erminea (H) Mustela erminea (C) Sorex cinereus Sorex hoyi Sorex minutissimus Sorex monticolus Sorex navigator Sorex tundrensis Sorex ugyunak

0.667 0.943 0.997 0.827 0.887 0.983 1.000 n/a 0.996 0.939 0.439 0.584 n/a 1.000 1.000 0.937 n/a n/a 0.599 0.905 0.707 1.000 0.987 0.999 0.999 1.000 0.998 0.999

0.081 0.002 0.000 0.001 0.045 0.000 0.000 n/a 0.000 0.000 0.047 0.133 n/a 0.000 0.000 0.000 n/a n/a 0.087 0.000 0.047 0.000 0.001 0.000 0.000 0.000 0.000 0.000

0.252 0.055 0.003 0.173 0.149 0.018 0.000 n/a 0.004 0.061 0.523 0.283 n/a 0.000 0.000 0.063 n/a n/a 0.315 0.095 0.247 0.000 0.012 0.001 0.001 0.000 0.002 0.001

0.443 0.406 0.383 0.347 0.226 0.534 0.381 n/a 0.639 0.440 0.265 0.400 n/a 0.530 0.578 0.408 n/a n/a 0.295 0.512 0.450 0.799 0.465 0.452 0.576 0.673 0.498 0.424

0.197 0.241 0.241 0.227 0.356 0.121 0.266 n/a 0.026 0.214 0.358 0.264 n/a 0.067 0.124 0.058 n/a n/a 0.354 0.146 0.220 0.001 0.205 0.246 0.118 0.027 0.172 0.239

0.360 0.353 0.375 0.427 0.419 0.345 0.353 n/a 0.335 0.346 0.377 0.336 n/a 0.403 0.299 0.535 n/a n/a 0.352 0.343 0.330 0.200 0.330 0.303 0.306 0.300 0.330 0.338

on model predictions projected across any timeframe, with the exception of the lowest sample size (n ¼ 20) for Lemmus (Appendix: Fig. A4). Genetic analyses.—Genetic signatures of most species exhibited significant values of Tajima’s D and R2 indicating recent demographic expansion—exceptions include some of the tundra and alpine species, notably M. miurus, Marmota broweri and Gulo gulo (Table 1). Timing of expansions inferred from BSPs that were calibrated to taxon-specific substitution rates generally indicated population growth coincident with climate change surrounding the LGM (Appenv www.esajournals.org

dix: Fig. A5). Population increases for the boreal species were consistent with post-LGM climate warming and associated northward shifts in boreal forests through North America. Population size changes of Beringian endemic species or those with a Palearctic evolutionary origin were more idiosyncratic, exhibiting gradual changes through time, earlier increases in population size, and/or accompanied by very recent population declines and some increases. Based on coalescent ABC analysis of demography according to scenarios set to test the reliability of BSP trends (Table 2, BSP scenarios), most empirical BSPs 11


HOPE ET AL.

were strongly supported by simulated data, indicating not only population growth in most species, but also recent decline in some species, particularly M. miurus (Table 2; Appendix: Fig. A5). Lowest support for BSPs was found for those species with conflict between results of expansion statistics, shape of BSP, and SDMs; namely M. broweri and G. gulo. Apparent recent population growth from BSPs for these two species is likely spurious, particularly given predicted severe range contractions from SDMs following the LGM (Appendix: Fig. A2). To further assess recent demographic changes, and to infer population trends into future time frames, we again used simulations of alternative scenarios (Table 2, LGM scenarios). While continued population increase received highest support for most species, population decline was best supported for M. broweri and G. gulo. Alternative scenarios to recent population growth were likewise supported for those species where BSPs indicated a recent decline (M. miurus, Microtus oeconomus, and Ochotona collaris), all representing endemic lineages within Beringia. Comparative ecological and evolutionary analyses.—When comparing the ecological and evolutionary evidence of population changes, species responses across all demographic tests were generally consistent, increasing support for observed trends (Table 1). Additionally, with only two exceptions (D. groenlandicus and S. ugyunak), observed genetic demographic results were consistent with increasing or decreasing predicted ranges from SDMs through time, considering a continued warming climate trend from the LGM into future timeframes (Table 1). Predictions for most boreal taxa reflected range expansions with concordant genetic signatures of significant recent demographic growth. Future SDMs of widespread and tundra species were predicted to either retain their present occupancy or contract in range, again reflected by genetic evidence. The most dramatic range contractions among highlatitude taxa were predicted for species associated with xeric (including alpine) areas and steppe tundra (e.g., D. groenlandicus, M. miurus, O. collaris, G. gulo, Urocitellus parryii; Table 1).

tire small and medium sized mammalian fauna under investigation reflected several salient findings (Fig. 1). During the previous interglacial warm period (i.e., LIG), distributional predictions were similar to those for the near future (2050s), although with less conservative prediction at the LIG, likely reflecting a more distant (broader confidence) climatic projection from the present (Fig. 1). Predictions for the LGM indicated lower but still moderate richness, suggesting potentially viable climatic conditions for most of these species in at least some regions of the study area, including the boreal taxa (Appendix: Fig. A6). The region of highest biodiversity projected to the LGM was coincident with central and south-central Beringia, including large areas currently submerged beneath the Bering Strait (Fig. 1; Appendix: Fig. A6). Western portions of the present-day study area, in particular the Seward Peninsula, were consistently predicted to accommodate high species richness across all historical, contemporary, and future timeframes. Future shifts in species richness through time (present through 2070s) indicate areas of high richness will shift from lower (e.g., valleys) to higher elevation upland areas, and from generally south and west to north and east within the study area. This shift in richness also reflected a shift from areas currently associated with boreal forest, taiga, and western mesic tundra ecoregions into areas currently dominated by both mesic and xeric tundra vegetation (Fig. 1). Future predictions of the distribution of richness between 2050s and 2070s exceed the northward extent of highest richness projected during the previous interglacial (LIG). Considering distinct communities, SDMs provided contrasting patterns of distributional change for tundra and boreal faunas (Fig. 2). Boreal species richness was predicted to shift northward and expand in total extent into the future, whereas tundra species richness was predicted to retract northward and decrease in total extent, although at a slower velocity of retreat than the simultaneous advance of the boreal biome, leading to an expanding zone of overlap between the two biomes through time. Width of overlap across the boreal-tundra ecotone was highest in western regions and decreased further east. Core refugial areas for highest predicted richness of tundra species in

Total biodiversity versus distinct community responses Temporally-predicted distributions of the env www.esajournals.org

12


HOPE ET AL.

Fig. 3. Plots of predicted small mammal community diversity in (A) the Arctic Network of National Parks (ARCN; squares), (B) the Arctic National Wildlife Refuge (Arctic NWR; circles), and (C) the National Petroleum Reserve–Alaska (NPR-A; triangles), for current (Now; left column) and 2070s (right column) timeframes, ranging from blue (low) to red (high) where colors reflect species distribution model predictions of total diversity (Fig. 1). Values reflect percentage area (y-axis) of each management region predicted to support a given proportion of richness (x-axis). For example, for current predictions within the NPR-A, almost half of total richness (12 taxa) are predicted to occur across roughly 45% of the management region, whereas peak richness (roughly two thirds of the taxa) is only predicted to occur through roughly 5% of this region. Into the 2070s timeframe, almost three quarters of the species are predicted to occur through over half of the total area of the NPR-A. Data points for the current timeframe indicate 0–25 taxa, whereas those for the 2070s timeframe indicate the 25 taxa predicted by between 0 and 3 model algorithms (resulting in each pixel ranging from 0-75).

future timeframes were coincident with areas of lowest predicted richness for boreal species during these same intervals (Fig. 2).

highest stability between these timeframes, reflecting persistence through time within this biome. Into future timeframes, total richness increased within all three land-use regions (Appendix: Fig. A6). However, by examining independent boreal and tundra biomes, major predicted increases were due to additional boreal species shifting northward into these areas. Tundra species richness was projected to remain

Variability among management regions Predicted number of species decreased dramatically between the LIG and LGM in all management regions (Fig. 1; Appendix: Figs. A6, A7). Tundra species richness exhibited v www.esajournals.org

13


HOPE ET AL.

constant or decrease through time. Quantification of the relative distribution of richness within major land-use regions (i.e., a measure of diversity; Fig. 3) indicated that even though the ARCN was predicted to accommodate the highest total richness at present, richness in this region was most unevenly distributed, consistent with vast areas of heterogeneous tundra and boreal habitats within the ARCN. The NPR-A, in contrast, is entirely within Arctic tundra habitats and this is reflected by the most evenly distributed richness. Richness in all regions was predicted to become more evenly distributed in the future (loss of b-diversity), although the ARCN and NPR-A were still predicted to have the most uneven and even distributions of richness, respectively (Fig. 3). Richness within the ARCN reflected a gradient from high in the west to low in the east, and within the Arctic NWR from high in the south (taiga) to low in the north (tundra; Fig. 1).

incorporating knowledge of the geographic limits of distinct evolutionary lineages, and demographic analyses by using change in distributional predictions to refine tests of competing scenarios. In particular, because understanding rates and empirical timing of changes is critical to managers, calibration of evolutionary rates (Hope et al. 2014) allowed for a comparative faunal evaluation across distinct communities. Further comparative analyses that integrate both ecotypic and genotypic characteristics will best increase our understanding of historical responses, and what we may expect in light of potential trajectories. SDM current predictions reflected known contemporary distributions based on museum specimens and are consistent with results using other predictive methods including machine learning (Baltensperger and Huettmann 2015). Importantly, predicted distributional changes through time were reflected in genetic demographic change, further supporting these environmental variables and methods as reasonable predictors for all distributions modeled for these species. These scenopoetic climate-based variables therefore provide accurate proxies for suitable conditions across an evolutionary timescale (Peterson et al. 2004). SDMs were also robust to random subsampling experiments, supporting the finding that thinning of records to reduce spatial bias can produce more robust predictions even at the cost of reduced data (Beck et al. 2014). Through careful optimization and refinement of variables on a per-species basis, and by utilizing genetically defined lineages for analysis, SDMs constitute conservative predictions that avoid over-fitting. However, this exploration of major community processes should be evaluated by investigating other taxonomic groups, and by refinement of decadal-scale changes using increasingly refined ecological (e.g., shrub cover, habitat associations, fire dynamics; Payette et al. 2001, Faleiro et al. 2013) and geophysical (e.g., hydrology, biogeochemistry, permafrost; Prowse et al. 2006, Raynolds et al. 2014) variables within predictive models. Other vertebrates, including resident and migratory birds and larger mammals, may exhibit responses to climate shifts disparate from these resident mammals. For example, birds may be limited by discrete resources for food and shelter, where a shift in species distribution may be more closely tied to community-level responses (Tulp and Schekkerman 2008, Virkkala et al. 2008,

DISCUSSION Comparative studies of resident high latitude species are providing key perspectives on community dynamics in the north (Eidesen et al. 2013, Reu et al. 2014). Our synthesis explored hypothesized responses to climate change among northern mammals and inferred the direction, tempo, and extent of responses based on statistical integration of contemporary distributions, fossil locations, abiotic climate variables, and genetic sequences. Assessment was possible for mammals due to their relatively widespread representation in museum archives and previous phylogeographic work that demonstrated sufficient levels of genetic variability for fine-scale comparison of spatiotemporal responses. Through multiple independent lines of evidence, we identified the potential for demographic responses to environmental change that have direct implications for conservation and management of high latitude biodiversity.

Consistency among ecological and evolutionary analyses Although both SDMs and genetic demographic results were largely congruent across all analyses, both inferences were improved by statistically combining their independent data; SDMs by v www.esajournals.org

14


HOPE ET AL.

Gilg et al. 2012). Often, carnivores and large herbivores exhibit low regional phylogeographic structure at Arctic and subarctic latitudes, as high mobility and migratory life-histories result in broad connectivity, buffering the localized effects of climate change (Teacher et al. 2011, Sommer et al. 2013). Groups such as invertebrates (Peterson et al. 2004, Hunter et al. 2014), plants, fungi, and lichens (e.g., Alsos et al. 2009, Frost and Epstein 2014) are more diverse, generally poor dispersers, and when coupled with phenological changes, these taxa may provide key insights into the complexity of local responses to climatic shifts. Fish and other aquatic components provide critical genetic information on connectivity across high-latitude landscapes and between freshwater and marine ecosystems (Reist et al. 2006, April et al. 2013). Collectively, these organisms are useful for assessing the impacts of environmental change and comparative approaches offer a rare opportunity to gain perspectives not offered by single-species studies. Although SDMs were generally robust with respect to predicting ranges of species and congruence with demographic trends, predictions occasionally deviated from expectations for the LGM timeframe where broad areas of potential niche were seemingly available for boreal species (Fig. 1; Appendix: Fig. A6). Few fossils support persistence of boreal taxa within Beringia during the last glacial period, nor is there evidence of persistence through the LGM from phylogeographic assessments, except possibly for T. hudsonicus (Chavez et al. 2014). Instead, molecular signatures indicate recent demographic expansions for boreal species now found in this region. The broad SDM predictions of boreal taxa in Beringia at the LGM, particularly in south-central areas, present a conundrum that may reflect our conservative methodological approach. While we were careful to avoid overfit distributions, the resulting models likely predict some non-realized niche space, indicating suitable conditions but unoccupied areas during this timeframe (Peterson et al. 2011). Our approach, while potentially overestimating actual paleodistributions, necessarily ensured that known fossil localities were predicted as suitable for species that persisted through the LGM. Similarly, our approach may be conservative with respect to future predictions. For instance, our results indicate large areas suitable for both v www.esajournals.org

tundra and boreal species into future timeframes, but competitive interactions across a broadening zone of contact between biomes may play a role in future community structure. Modeling that accounts for biotic interactions is complex (Araujo and Luoto 2007, Gutie´rrez et al. 2014). However if faster velocities of change among boreal taxa reflect competitive dominance then in the case of future timeframes, our predictions for tundra species distributions could reflect nonrealized niches, resulting in more pronounced range contractions than predicted among the tundra species. Although as yet unsubstantiated, these lines of investigation constitute important criteria for continued refinement of predictions into future decades by long-term field and analytical initiatives. An alternative explanation for the apparent absence of boreal species from the study area during the LGM is extreme climate events. Fossil evidence exists for some of these species within Beringia during early stages of the Wisconsinan (e.g., Lepus americanus, Microtus pennsylvanicus; Harington 2011), coupled with recent evidence indicating that forests persisted in isolated regions through the entire glacial phase (Brubaker et al. 2005, Zazula et al. 2006). The emerging view that central Beringia was characterized by more temperate, mesic conditions through much of the glacial phase (Elias and Crocker 2008, Lozhkin and Anderson 2011) supports Guthrie’s (2001) concept of a mesic buckle through this region. Our LGM predictions for both boreal and widespread mammals (Appendix: Fig. A5) lend further evidence to this assertion. Consequently, at least some boreal species likely persisted in Beringia through at least the early stages of the Wisconsinan glacial phase. These remnant populations would be near their ecological limits, and may subsequently have experienced stochastic extirpations by extreme climatic events of short duration throughout this period, of which many are recorded (Hofreiter and Stewart 2009, Brace et al. 2012). Considering still the relative paucity of sampling across this vast region, evidence of glacial phase persistence of other boreal mammals within Beringia may arise given more comprehensive field and genetic exploration. For instance, new evidence indicates persistence of both Microtus longicaudus and Peromyscus sp. 15


HOPE ET AL.

within eastern Beringia and the Alexander Archipelago during the LGM (Sawyer 2014). Concomitantly, because past interglacial stages also exhibited punctuated warm extremes, some tundra species may experience similar local or regional stochastic extirpations during warmer intervals in the near future. Our understanding of the intrinsic connection between ecological and historical biogeographic perspectives will undoubtedly improve through ongoing studies of genome- and range-wide genetic variation. For example, rapid changes in range may lead to adaptive genetic responses in small mammals (White et al. 2013). For high-latitude taxa experiencing rapidly changing environments, emerging evolutionary methods will allow us to assess how selective pressures are acting on genomes (Stapley et al. 2010).

tundra species (coincident with the Arctic-NWR) was also coincident with predicted low diversity areas for boreal species, reflecting regions where the tundra community may persist independent of boreal taxa. Importantly, the minimum predicted areas for tundra diversity (illustrated by 2070s predictions; Fig. 2), although fragmented across northern portions of the study area, were consistently recovered as suitable (often interconnected) habitat through all timeframes, indicating that these refugia are potentially critical to the future persistence of the tundra biome. In terms of individual taxa, species associated with mesic tundra (L. trimucronatus, S. ugyunak) and some widespread endemic species were not predicted to significantly decline in range size as reflected by demographic signatures of moderate growth or stability. The generally stable SDM predictions for S. ugyunak, contrast with significant genetic signatures of expansion, possibly reflecting in situ population growth without distributional change. However, species associated with xeric tundra or alpine habitats (e.g., M. miurus; O. collaris, U. parryii) primarily responded negatively to warming climate (reduction in range size and/or decline in population size), likely indicating adaptation to cold and dry polar environments that are in decline (Callaghan et al. 2004). Inconsistent results for D. groenlandicus are more difficult to explain. Genetic evidence suggests sustained but gradual population growth through the extent of the last glacial phase, reflecting high genetic diversity for this taxon (second highest for all study species after L. trimucronatus), buffering signals of predicted range contraction. This also highlights inherent idiosyncrasies among the coldadapted taxa and extended persistence in the study area (Cook et al. 2005).

Total biodiversity versus distinct community responses Our combined faunal data indicate an expected increase in total mammalian diversity for the northern high latitudes in the near-future. However, we demonstrated that distinct components of ecological communities exhibited independent, but generally predictable, responses to regional climate trends across time. Evolutionary signatures indicate a predicted northward advance across boreal species (Lessa et al. 2003), but greater heterogeneity among tundra species (earlier and/or more gradual population changes; Morris et al. 2011). The different velocities of change between boreal (fast) and tundra (slow) communities may constitute the most critical observation from our data, as this translates to increasing community overlap through time, although the consequences of community reorganization are poorly understood. Although most tundra-adapted mammal species are predicted to persist, their combined distributions will likely be reduced and fragmented, and tundra species may overlap more with expanding boreal species resulting in novel community associations. Substantial northward range expansions of boreal species are predicted to extend beyond the hypothesized maximum extent during the LIG, translating to rapid community turnover in these northern latitudes and a potentially tenuous future for cold-associated taxa. One of the two major refugia predicted for v www.esajournals.org

Variability among management regions This comparative assessment constitutes a basis for informed decisions about the future distribution and diversity of wildlife that may be applicable to management of Arctic resources (Hinzman et al. 2005), especially given predicted differences among major management regions. While predicted diversity is relatively evenly distributed through the NPR-A, reflecting a single tundra ecoregion and low topographic variability, diversity is unevenly distributed through the more ecologically variable ARCN and Arctic 16


HOPE ET AL.

NWR. Although all regions are predicted to increase in total richness through time, each region reflects unique attributes of this heterogeneous study region, with the most diverse mammal communities found in the western ARCN and southern Arctic NWR into 2070s. However, areas of highest conservation potential for the tundra biome may be centered where tundra diversity is highest relative to overall diversity. Although areas of refugial tundra diversity are predicted within all three of the management regions, northern portions of the Arctic NWR represent lowest overlap between tundra and boreal species predictions. Importantly, our results emphasize the importance of considering ecological, biogeographical, and political boundaries for future investigations through the Holarctic, considering that interactions among these layers of complexity are highly regionalized. Past and present SDM predictions for these regions were well supported by both fossil and genetic data, and from documented historic occurrences of these species. For example, the current east to west gradient of diversity predicted within the ARCN is consistent with recent survey efforts across this region (MacDonald and Cook 2009, Hope 2012). Recent sampling has also confirmed the occurrence of some boreal small mammal species north of subarctic treeline (e.g., Sorex; Hope 2012), as predicted by SDMs. SDMs projected through the next century also constitute distributional hypotheses of shifting ranges that may be tested into the future with continued field monitoring. Regular long-term specimenbased surveys are therefore vital for understanding species responses to changing and novel climate conditions (Hope et al. 2013b). There is growing concern over the persistence of Arctic biodiversity due to drastic declines in both distribution and diversity of some species (Hunter et al. 2010, Bellard et al. 2012). However other components of Arctic biodiversity, such as some migratory species, reflect recent population increases (e.g., Flint et al. 2008). Similarly, predictions for small and medium-sized mammals in Beringia include both declines and expansions that reflect broad ecological associations. Importantly, most taxa should persist into the next century although their distributions and abundances may decline. The few tundra species that reflect recent population declines are largely cov www.esajournals.org

distributed and their predicted refugial areas may be sufficient to maintain community associations. Evaluation of spatial changes in biodiversity through time also provides critical insight to comparative land-use between humans and wildlife (Faleiro et al. 2013, Kuemmerle et al. 2014). Focusing management on both total biodiversity and refugial areas associated with the tundra biome will therefore represent a balance of priorities. Designated land use for tundra refugial areas includes a mix of both protected wildlife refuges and parks, and regions of active resource exploration and development. Despite differences in land use however, all three of the federal management regions we studied have a conservation mandate. As the distribution and abundance of individual taxa shift and novel community associations inevitably form, land managers will be challenged to define and manage for ‘‘natural’’ diversity in the face of competing and potentially escalating demands for energy development in (NPR-A) or adjacent to (ARCN, Arctic NWR) these areas, as well as for other uses such as hunting and recreation (ARCN, Arctic NWR, NPR-A).

Local ecosystem to global biosphere A drawback of some meta-analyses is that independent studies often lack consistency across methods, particularly concerning geographic representation of samples, choice of genetic loci, and assumed evolutionary rates for accurate temporal inferences (Bellard et al. 2012). It is therefore unsurprising that few community-level assessments incorporate statistical phylogeographic approaches (e.g., as outlined by Richards et al. 2007, but see Dolman and Joseph 2012, Oaks et al. 2013). Beringia, however, offers an exceptional comparative system. The number of previous phylogeographic assessments available that span this region is comparable to the rich history of this discipline in Europe (Taberlet et al. 1998, Hewitt 2004) and temperate North America (Avise 2000, Soltis et al. 2006), but Beringia has the added complexity of combining evolutionary endemism with a dynamic ecological history of colonization and movement originating from multiple continents. Beringian phylogeography has also been complemented by a rich history of paleo-ecological (e.g., Hopkins 1982, Elias and Crocker 2008) and climatological investigations (e.g., Brigham17


HOPE ET AL.

Grette 2001, Miller et al. 2010). These multidisciplinary perspectives highlight the significant role of Beringia for intercontinental exchange and diversification of the northern faunas. Because this biogeographic crossroads between the northern continents is now expected to experience significant human modifications, comparative approaches may also provide a basis for management agencies with diverse mandates to assess how changing conditions and land uses will impact these highly productive environments. The implications of this assessment, however, transcend contrasting patterns of mammalian biodiversity and the observation of geographic overlap between localized human land-use and changing vertebrate communities. Significant future changes in community structure and increases in diversity, as indicated by our predictions, will likely be complicated by increased competition (Prowse et al. 2009), hybridization (Hewitt 2011, April et al. 2013), and introduction of parasites and disease agents. Pathogen lifecycles can abruptly shift across climate-driven thresholds, increasing morbidity and mortality of hosts, and leading to emerging zoonoses (Kutz et al. 2005, Hoberg et al. 2012). All of these dynamics across the Holarctic may strongly influence genetic diversity and viability of northern highlatitude species, many of which are integral to local human subsistence as well as the global economy. Our projected distributions for the previous interglacial warm period (LIG) closely resembled distributions predicted into the near future, but SDMs further into the coming century indicate a continued northward shift of biodiversity. These predictions reflect a warmer world than at any time within the evolutionary history of extant Arctic-adapted (cold-associated) species and emphasize the importance of locating and accommodating vital refugial areas for these taxa into future land-management scenarios (Cook et al. 2006, Woodruff 2010). Novel ecological conditions that drive community reorganization may lead to a state shift within the tundra biome (Barnosky et al. 2012) although such threshold changes are often difficult to predict (Hinzman et al. 2005). The SDMs and genetic signatures we have presented constitute ‘‘business-as-usual’’ trajectories following current environmental rates of change. Without considerable changes in strategies for reducing current environmental v www.esajournals.org

impacts, these forecasts may reflect conservative trends for the future of high latitude biodiversity (Ehrlich and Pringle 2008).

ACKNOWLEDGMENTS Support was provided by U.S. Geological Survey’s Science Support Program, Alaska Regional Executive DOI on the Landscape initiative and U.S. Geological Survey’s Changing Arctic Ecosystems Initiative, which is supported by funding from the Wildlife Program of the USGS Ecosystem Mission Area. Fieldwork producing the georeferenced specimens was supported by the National Park Service Arctic Network and National Science Foundation (NSF-DEB 0415668 and 1258010). Statistical analyses were facilitated by the University of Alaska, Fairbanks (UAF) Life Science Informatics Portal, http:// biotech.inbre.alaska.edu, supported by Grant Number RR016466 from the National Center for Research Resources (NCRR), a component of the National Institutes of Health (NIH). Thanks to J. R. Stewart, R. D. Guthrie, and J. Pearce for comments on an earlier draft. Any use of trade names is for descriptive purposes only and does not imply endorsement by the U.S. Government. The findings and conclusions in this manuscript are those of the authors and do not necessarily represent the views of the U. S. Fish and Wildlife Service.

LITERATURE CITED Abbott, R. J., and C. Brochmann. 2003. History and evolution of the Arctic flora: in the footsteps of Eric Hulteń. Molecular Ecology 12:299–313. Alsos, I. G., T. Alm, S. Normand, and C. Brochmann. 2009. Past and future range shifts and loss of diversity in dwarf willow (Salix herbacea L.) inferred from genetics, fossils and modeling. Global Ecology and Biogeography 18:223–239. Anderson, R. P., and I. Gonzalez, Jr. 2011. Speciesspecific tuning increases robustness to sampling bias in models of species distributions: an implementation with Maxent. Ecological Modelling 222:2796–2811. Anderson, R. P., and A. Raza. 2010. The effect of the extent of the study region on GIS models of species geographic distributions and estimates of niche evolution: preliminary tests with montane rodents (genus Nephelomys) in Venezuela. Journal of Biogeography 37:1378–1393. April, J., R. H. Hanner, A. Dion-C oˆ te´, and L. Bernatchez. 2013. Glacial cycles as an allopatric speciation pump in north-eastern American freshwater fishes. Molecular Ecology 22:409–422. Arau´jo, M. B., and M. Luoto. 2007. The importance of biotic interactions for modelling species distributions under climate change. Global Ecology and

18


HOPE ET AL. Biogeography 16:743–753. Avise, J. C. 2000. Phylogeography: the history and formation of species. Harvard University Press, Cambridge, Massachusetts, USA. Baltensperger, A. P., and F. Huettmann. 2015. Predictive spatial niche and biodiversity hotspot models for small mammal communities in Alaska: applying machine-learning to conservation planning. Landscape Ecology 30:681–697. Barnosky, A. D., et al. 2012. Approaching a state shift in Earth’s biosphere. Nature 486:52–58. Barve, N., V. Barve, A. Jimenez-Valverde, A. LiraNoriega, S. P. Maher, A. T. Peterson, J. Soberoń, and F. Villalobos. 2011. The crucial role of the accessible area in ecological niche modeling and species distribution modeling. Ecological Modelling 222:1810–1819. Beaumont, M. A., W. Zhang, and D. J. Balding. 2002. Approximate Bayesian computation in population genetics. Genetics 162:2025–2035. Beck, J., M. Bo¨ller, A. Erhardt, and W. Schwanghart. 2014. Spatial bias in the GBIF database and its effect on modeling species’ geographic distributions. Ecological Informatics 19:10–15. Bellard, C., C. Bertelsmeier, P. Leadley, W. Thuiller, and F. Courchamp. 2012. Impacts of climate change on the future of biodiversity. Ecology Letters 15:365–377. Brace, S., E. Palkopoulou, L. Daleń, A. M. Lister, R. Miller, M. Otte, M. Germonpre´, S. P. E. Blockley, J. R. Stewart, and I. Barnes. 2012. Serial population extinctions in a small mammal indicate Late Pleistocene ecosystem instability. Proceedings of the National Academy of Sciences USA 109:20532–20536. Brigham-Grette, J. 2001. New perspectives on Beringian Quaternary paleogeography, stratigraphy, and glacial history. Quaternary Science Reviews 20:15–24. Bromham, L. 2011. The genome as a life-history character: why rate of molecular evolution varies between mammal species. Philosophical Transactions of the Royal Society B 366:2503–2513. Brubaker, L. B., P. M. Anderson, M. E. Edwards, and A. V. Lozhkin. 2005. Beringia as a glacial refugium for boreal trees and shrubs: new perspectives from mapped pollen data. Journal of Biogeography 32:833–848. Brunhoff, C., K. E. Galbreath, V. B. Fedorov, J. A. Cook, and M. Jaarola. 2003. Holarctic phylogeography of the root vole (Microtus oeconomus): implications for late Quaternary biogeography of high latitudes. Molecular Ecology 12:957–968. Callaghan, T. V., et al. 2004. Biodiversity, distributions and adaptations of Arctic species in the context of environmental change. Ambio 33:404–417. Carnaval, A. C., M. J. Hickerson, C. F. B. Haddad, M. T. Rodrigues, and C. Moritz. 2009. Stability predicts genetic diversity in the Brazilian Atlantic Forest hotspot. Science 323:785–789.


Chan, L. M., J. L. Brown, and A. D. Yoder. 2011. Integrating statistical genetic and geospatial methods brings new power to phylogeography. Molecular Phylogenetics and Evolution 59:523–537. Chapin, F. S., A. L. Lovecraft, E. S. Zavaleta, J. Nelson, M. D. Robards, G. P. Kofinas, S. F. Trainor, G. D. Peterson, H. P. Huntington, and R. L. Naylor. 2006. Policy strategies to address sustainability of Alaskan boreal forests in response to a directionally changing climate. Proceedings of the National Academy of Sciences USA 103:16637–16643. Chavez, A. S., S. P. Maher, B. S. Arbogast, and G. J. Kenagy. 2014. Diversification and gene flow in nascent lineages of island and mainland North American tree squirrels (Tamiasciurus). Evolution 68:1094–1109. Cook, J. A., et al. 2005. Beringia: intercontinental exchange and diversification of high latitude mammals and their parasites during the Pliocene and Quaternary. Mammal Study 30:S33–S44. Cook, J. A., N. G. Dawson, and S. O. MacDonald. 2006. Conservation of highly fragmented systems: the north temperate Alexander Archipelago. Biological Conservation 133:1–15. Cornuet. J.-M. , P. Pudlo, J. Veyssier, A. Dehne-Garcia, M. Gautier, R. Leblois, J.-M. Marin, and A. Estoup. 2014. DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics 30:1187–1189. Davis, E. B., J. L. McGuire, and J. D. Orcutt. 2014. Ecological niche models of mammalian glacial refugia show consistent bias. Ecography 37:1133–1138. Diffenbaugh, N. S., and F. Giorgi. 2012. Climate change hotspots in the CMIP5 global climate model ensemble. Climate Change 114:813–822. Dolman, G., and L. Joseph. 2012. A species assemblage approach to comparative phylogeography of birds in southern Australia. Ecology and Evolution 2:354–369. Drummond, A. J., A. Rambaut, B. Shapiro, and O. G. Pybus. 2005. Bayesian coalescent inference of past population dynamics from molecular sequences. Molecular Biology and Evolution 22:1185–1192. Duerden, F. 2004. Translating climate change impacts at the community level. Arctic 57:204–212. Ehrlich, P. R., and R. M. Pringle. 2008. Where does biodiversity go from here? A grim business-asusual forecast and a hopeful portfolio of partial solutions. Proceedings of the National Academy of Sciences USA 105:11579–11586. Eidesen, P. B., D. Ehrich, V. Bakkestuen, I. G. Alsos, O. Gilg, P. Taberlet, and C. Brochmann. 2013. Genetic roadmap of the Arctic: plant dispersal highways, traffic barriers and capitals of diversity. New Phytologist 200:898–910.

19


HOPE ET AL. Elias, S. A., and B. Crocker. 2008. The Bering Land Bridge: a moisture barrier to the dispersal of steppe-tundra biota? Quaternary Science Reviews 27:2473–2483. Faleiro, F. V., R. B. Machado, and R. D. Loyola. 2013. Defining spatial conservation priorities in the face of land-use and climate change. Biological Conservation 158:248–257. Fedorov, V. B., and N. C. Stenseth. 2002. Multiple glacial refugia in the North American Arctic: inference from phylogeography of the collared lemming (Dicrostonyx groenlandicus). Proceedings of the Royal Society B 269:2071–2077. Flint, P. L., E. J. Mallek, R. J. King, J. A. Schmutz, K. S. Bollinger, and D. V. Derksen. 2008. Changes in abundance and spatial distribution of geese molting near Teshekpuk Lake, Alaska: interspecific competition or ecological change? Polar Biology 31:549–556. Ford, J. D., and B. Smit. 2004. A framework for assessing the vulnerability of communities in the Canadian Arctic to risks associated with climate change. Arctic 57:389–400. Frost, G. V., and H. E. Epstein. 2014. Tall shrub and tree expansion in Siberian tundra ecotones since the 1960s. Global Change Biology 20:1264–1277. Gilg, O., et al. 2012. Climate change and the ecology and evolution of Arctic vertebrates. Annals of the New York Academy of Sciences 1249:166–190. Guthrie, D. R. 2001. Origin and causes of the mammoth steppe: a story of cloud cover, woolly mammal tooth pits, buckles, and inside-out Beringia. Quaternary Science Reviews 20:549–574. Gutie´rrez, E. E., R. A. Boria, and R. P. Anderson. 2014. Can biotic interactions cause allopatry? Niche models, competition, and distributions of South American mouse opossums. Ecography 37:741–753. Harington, C. R. 2011. Pleistocene vertebrates of the Yukon Territory. Quaternary Science Reviews 30:2341–2354. Hewitt, G. M. 2004. Genetic consequences of climatic oscillations in the Quaternary. Philosophical Transactions of the Royal Society B 359:183–195. Hewitt, G. M. 2011. Quaternary phylogeography: the roots of hybrid zones. Genetica 139:617–638. Hickerson, M. J., B. C. Carstens, J. Cavender-Bares, K. A. Crandall, C. H. Graham, J. B. Johnson, L. Rissler, P. F. Victoriano, and A. D. Yoder. 2010. Phylogeography’s past, present, and future: 10 years after. Molecular Phylogenetics and Evolution 54:291–301. Hijmans, R. J., S. E. Cameron, J. L. Parra, P. G. Jones, and A. Jarvis. 2005. Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25:1965–1978. Hinzman, L. D., et al. 2005. Evidence and implications of recent climate change in northern Alaska and


other Arctic regions. Climatic Change 72:251–298. Ho, S. Y. W., U. Saarma, R. Barnett, J. Haile, and B. Shapiro. 2008. The effect of inappropriate calibration: three case studies in molecular ecology. PLoS ONE 3:e1615. Ho, S. Y. W., and B. Shapiro. 2011. Skyline-plot methods for estimating demographic history from nucleotide sequences. Molecular Ecology Resources 11:423–434. Hoberg, E. P., K. E. Galbreath, J. A. Cook, S. J. Kutz, and L. Polley. 2012. Northern host-parasite assemblages: history and biogeography on the borderlands of episodic climate and environmental transition. Advances in Parasitology 79:1–81. Hof, A. R., R. Jansson, and C. Nilsson. 2012. Future climate change will favour non-specialist mammals in the (sub)Arctics. PLoS ONE 7:e52574. Hofreiter, M., and I. Barnes. 2010. Diversity lost: are all Holarctic large mammal species just relict populations? BMC Biology 8:46. Hofreiter, M., and J. Stewart. 2009. Ecological change, range fluctuations and population dynamics during the Pleistocene. Current Biology 19:R584–R594. Hope, A. G. 2012. High shrew diversity on Alaska’s Seward Peninsula: community assembly and environmental change. Northwestern Naturalist 93:101–110. Hope, A. G., S. Y. W. Ho, J. L. Malaney, J. A. Cook, and S. L. Talbot. 2014. Accounting for rate variation among lineages in comparative demographic analyses. Evolution 68:2689–2700. Hope, A. G., N. Takebayashi, K. E. Galbreath, S. L. Talbot, and J. A. Cook. 2013a. Temporal, spatial and ecological dynamics of speciation among amphiBeringian small mammals. Journal of Biogeography 40:415–429. Hope, A. G., E. Waltari, D. C. Payer, J. A. Cook, and S. L. Talbot. 2013b. Future distribution of tundra refugia in northern Alaska. Nature Climate Change 3:931–938. Hopkins, D. M. 1982. Paleoecology of Beringia. Academic Press, New York, New York, USA. Hulteń, E. 1937. Outline of the history of Arctic and boreal biota during the Quaternary period. J. Cramer, New York, New York, USA. Hunter, C. M., H. Caswell, M. C. Runge, E. V. Regehr, S. C. Amstrup, and I. Stirling. 2010. Climate change threatens polar bear populations: a stochastic demographic analysis. Ecology 91:2883–2897. Hunter, M. D., M. V. Kozlov, J. Ita¨mies, E. Pulliainen, J. Ba¨ck, M. Kyro¨, and P. Niemela¨. 2014. Current temporal trends in moth abundance are counter to predicted effects of climate change in an assemblage of subarctic forest moths. Global Change Biology 20:1723–1737. Knowles, L. L., and D. F. Alvarado-Serrano. 2010. Exploring the population genetic consequences of

20


HOPE ET AL. the colonization process with spatio-temporally explicit models: insights from coupled ecological, demographic and genetic models in montane grasshoppers. Molecular Ecology 19:3727–3745. Kramer-Schadt, S., et al. 2013. The importance of correcting for sampling bias in MaxEnt species distribution models. Diversity and Distributions 19:1366–1379. Kuemmerle, T., L. Baskin, P. J. Leitão, V. Prishchepov, K. Thonicke, and V. C. Radeloff. 2014. Potential impacts of oil and gas development and climate change on migratory reindeer calving grounds across the Russian Arctic. Diversity and Distributions 20:416–429. Kutz, S. J., E. P. Hoberg, L. Polley, and E. J. Jenkins. 2005. Global warming is changing the dynamics of Arctic host-parasite systems. Proceedings of the Royal Society B 272:2571–2576. Lessa, E. P., J. A. Cook, and J. L. Patton. 2003. Genetic footprints of demographic expansion in North America, but not Amazonia, during the Late Quaternary. Proceedings of the National Academy of Sciences USA 100:10331–10334. Librado, P., and J. Rozas. 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25:1451–1452. Lorenzen, E. D., et al. 2011. Species-specific responses of Late Quaternary megafauna to climate and humans. Nature 479:359–364. Lozhkin, A. V., and P. M. Anderson. 2011. Forest or no forest: implications of the vegetation record for climatic stability in western Beringia during oxygen isotope stage 3. Quaternary Science Reviews 30:2160–2181. MacDonald, S. O., and J. A. Cook. 2009. Recent mammals of Alaska. University of Alaska Press, Fairbanks, Alaska, USA. Meltofte, H., editor. 2013. Arctic biodiversity assessment. Status and trends in Arctic biodiversity. Conservation of Arctic Flora and Fauna, Akureyri, Iceland. Merow, C., M. J. Smith, and J. A. Silander. 2013. A practical guide to MaxEnt for modeling species’ distributions: what it does, and why inputs and settings matter. Ecography 36:1058–1069. Miller, G. H., et al. 2010. Temperature and precipitation history of the Arctic. Quaternary Science Reviews 29:1679–1715. Morris, D. W., D. E. Moore, S. B. Ale, and A. Dupuch. 2011. Forecasting ecological and evolutionary strategies to global change: an example from habitat selection by lemmings. Global Change Biology 17:1266–1276. Oaks, J. R., J. Sukumaran, J. A. Esselstyn, C. W. Linkem, C. D. Siler, M. T. Holder, and R. M. Brown. 2013. Evidence for climate-driven diversification? A caution for interpreting ABC inferences of simultaneous historical events. Evolution 67:991–1010.


Otto-Bliesner, B. L., S. J. Marshall, J. T. Overpeck, G. H. Miller, and A. Hu. 2006. Simulating Arctic climate warmth and icefield retreat in the last interglaciation. Science 311:1751–1753. Overland, J. E., M. Wang, J. E. Walsh, and J. C. Stroeve. 2013. Future Arctic climate changes: adaptation and mitigation timescales. Earth’s Future 2:68–74. Palsbøll, P. J., M. Z. Peery, M. T. Olsen, S. R. Beissinger, and M. Be´rube´. 2013. Inferring recent historic abundance from current genetic diversity. Molecular Ecology 22:22–40. Pauls, S. U., C. Nowak, M. Ba´lint, and M. Pfenninger. 2013. The impact of global climate change on genetic diversity within populations and species. Molecular Ecology 22:925–946. Payette, S., M. J. Fortin, and I. Gamache. 2001. The subarctic forest-tundra: the structure of a biome in a changing climate. BioScience 51:709–718. Peterson, A. T., E. Martıńez-Meyer, C. Gonza´lez-Salazar, and P. W. Hall. 2004. Modeled climate change effects on distributions of Canadian butterfly species. Canadian Journal of Zoology 82:851–858. Peterson, A. T., J. Soberoń, R. G. Pearson, R. P. Anderson, E. Martıńez-Meyer, M. Nakamura, and M. B. Arau´jo. 2011. Ecological niches and geographic distributions. Princeton University Press, Princeton, New Jersey, USA. Post, E., et al. 2009. Ecological dynamics across the Arctic associated with recent climate change. Science 325:1355–1358. Prost, S., R. P. Guralnick, E. Waltari, V. B. Fedorov, E. Kuzmina, N. Smirnov, T. Kolfschoten, M. Hofreiter, and K. Vrieling. 2013. Losing ground: past history and future fate of Arctic small mammals in a changing climate. Global Change Biology 19:1854–1864. Prost, S., N. Smirnov, V. B. Fedorov, R. S. Sommer, M. Stiller, D. Nagel, M. Knapp, and M. Hofreiter. 2010. Influence of climate warming on arctic mammals? NewinsightsfromancientDNA studiesofthecollared lemming Dicrostonyx torquatus. PLoS One 5:e10447. Prowse, T. D., C. Furgal, F. J. Wrona, and J. D. Reist. 2009. Implications of climate change for northern Canada: freshwater, marine and terrestrial ecosystems. Ambio 38:282–289. Prowse, T. D., F. J. Wrona, J. D. Reist, J. J. Gibson, J. E. Hobbie, L. M. Le´vesque, and W. F. Vincent. 2006. Climate change effect on hydroecology of Arctic freshwater ecosystems. Ambio 35:347–358. Ramos-Onsins, S. E., and J. Rozas. 2002. Statistical properties of new neutrality tests against population growth. Molecular Biology and Evolution 19:2092–2100. Raynolds, M. K., D. A. Walker, K. J. Ambrosius, J. Brown, K. R. Everett, M. Kanevskiy, G. P. Kofinas, V. E. Romanovsky, Y. Shur, and P. J. Webber. 2014. Cumulative geoecological effects of 62 years of infrastructure and climate change in ice-rich

21


HOPE ET AL. permafrost landscapes, Prudhoe Bay Oilfield, Alaska. Global Change Biology 20:1211–1224. Reddy, S., and L. M. Da´valos. 2003. Geographical sampling bias and its implications for conservation priorities in Africa. Journal of Biogeography 30:1719–1727. Reist, J. D., F. J. Wrona, T. D. Prowse, M. Power, J. B. Dempson, R. J. Beamish, J. R. King, T. J. Carmichael, and C. D. Sawatzky. 2006. General effects of climate change on Arctic fishes and fish populations. Ambio 35:370–380. Reu, B., S. Zaehle, K. Bohn, R. Pavlick, S. Schmidtlein, J. W. Williams, and A. Kleidon. 2014. Future noanalogue vegetation produced by no-analogue combinations of temperature and insolation. Global Ecology and Biogeography 23:156–167. Richards, C. L., B. C. Carstens, and L. L. Knowles. 2007. Distribution modeling and statistical phylogeography: an integrative framework for generating and testing alternative biogeographical hypotheses. Journal of Biogeography 34:1833–1845. Sawyer, Y. E. 2014. Living on the edge: a comparative phylogeographic study of refugial and insular fragmentation. Dissertation. University of New Mexico., New Mexico, USA. Schloss, C. A., T. A. Nuñez, and J. J. Lawler. 2012. Dispersal will limit ability of mammals to track climate change in the Western Hemisphere. Proceedings of the National Academy of Sciences USA 109:8606–8611. Soltis, D. E., A. B. Morris, J. S. McLachlan, P. S. Manos, and P. S. Soltis. 2006. Comparative phylogeography of unglaciated eastern North America. Molecular Ecology 15:4261–4293. Sommer, R. S., J. Kalbe, J. Ekstro¨m, N. Benecke, and R. Liljegren. 2013. Range dynamics of the reindeer in Europe during the last 25,000 years. Journal of Biogeography 41:298–306. Sta¨dler, T., B. Haubold, C. Merino, W. Stephan, and P. Pfaffelhuber. 2009. The impact of sampling schemes on the site frequency spectrum in nonequilibrium subdivided populations. Genetics 182:205–216. Stapley, J., J. Reger, P. G. D. Feulner, C. Smadja, J. Galindo, R. Ekblom, C. Bennison, A. D. Ball, A. P. Beckerman, and J. Slate. 2010. Adaptation genomics: the next generation. Trends in Ecology and Evolution 25:705–712. Stewart, J. R., and A. M. Lister. 2001. Cryptic northern refugia and the origin of the modern biota. Trends in Ecology and Evolution 16:608–613. Taberlet, P., L. Fumagalli, A. G. Wust-Saucy, and J. F. Cosson. 1998. Comparative phylogeography and postglacial colonization routes in Europe. Molecular Ecology 7:453–464. Tajima, F. 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymor-


phism. Genetics 123:585–595. Teacher, A. G. F., J. A. Thomas, and I. Barnes. 2011. Modern and ancient red fox (Vulpes vulpes) in Europe show an unusual lack of geographical and temporal structuring, and differing responses within the carnivores to historical climatic change. BMC Evolutionary Biology 11:214. Tulp, I., and H. Schekkerman. 2008. Has prey availability for Arctic birds advanced with climate change? Hindcasting the abundance of tundra arthropods using weather and seasonal variation. Arctic 61:48–60. Usher, M. B., T. V. Callaghan, G. Gilchrist, B. Heal, G. P. Juday, H. Loeng, M. A. K. Muir, and P. Prestrud. 2005. Principles of conserving the Arctic’s biodiversity. Pages 539–596 in Arctic Climate Impact Assessment. ACIA Overview report, Cambridge University Press, Cambridge, UK. Veloz, S. D. 2009. Spatially autocorrelated sampling falsely inflates measures of accuracy for presenceonly niche models. Journal of Biogeography 36:2290–2299. Virkkala, R., R. K. Heikkinen, N. Leikola, and M. Luoto. 2008. Projected large-scale range reductions of northern-boreal land bird species due to climate change. Biological Conservation 141:1343–1353. Waltari, E., E. P. Hoberg, E. P. Lessa, and J. A. Cook. 2007. Eastward ho: phylogeographical perspectives on colonization of hosts and parasites across the Beringian nexus. Journal of Biogeography 34:561–574. Warren, D. L., R. E. Glor, and M. Turelli. 2010. ENMTools: a toolbox for comparative studies of environmental niche models. Ecography 33:607–611. Warren, D. L., and S. Seifert. 2011. Ecological niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria. Ecological Applications 21:335–342. White, T. A., S. E. Perkins, G. Heckel, and J. B. Searle. 2013. Adaptive evolution during an ongoing range expansion: the invasive bank vole (Myodes glareolus) in Ireland. Molecular Ecology 22:2971–2985. Wood, D. A., A. G. Vandergast, K. R. Barr, R. D. Inman, T. C. Esque, K. E. Nussear, and R. N. Fisher. 2013. Comparative phylogeography reveals deep lineages and regional evolutionary hotspots in the Mojave and Sonoran deserts. Diversity and Distributions 19:722–737. Woodruff, D. S. 2010. Biogeography and conservation in Southeast Asia: how 2.7 million years of repeated environmental fluctuations affect today’s patterns and the future of the remaining refugialphase biodiversity. Biodiversity and Conservation 19:919–941. Worth, J. R. P., G. J. Williamson, S. Sakuguchi, P. G. Nevill, and G. J. Jordan. 2014. Environmental niche modelling fails to predict Last Glacial Maximum refugia: niche shifts, microrefugia or incorrect

22


HOPE ET AL. Zazula, G. D., A. M. Telka, C. R. Harington, C. E. Schweger, and R. W. Matthews. 2006. New spruce (Picea spp.) macrofossils from Yukon Territory: implications for Late Pleistocene refugia in eastern Beringia. Arctic 59:391–400.

palaeoclimate estimates? Global Ecology and Biogeography 23:1186–1197. Yannic, G., et al. 2014. Genetic diversity in caribou linked to past and future climate change. Nature Climate Change 4:132–137.

SUPPLEMENTAL MATERIAL APPENDIX

Table A1. Parameters and output statistics for species distribution models. Species are arranged by major biome associations. Extent used in initial SDMs Species

Biome

n

W

E

S

N

Non-correlated rasters (,0.8)

Optimal betas

Dicrostonyx groenlandicus Lemmus trimucronatus Microtus miurus Marmota broweri Lepus othus Gulo gulo Sorex ugyunak Microtus longicaudus Microtus pennsylvanicus Synaptomys borealis Tamiasciurus hudsonicus Marmota caligata Zapus hudsonius Zapus princeps Lepus americanus Sorex cinereus Sorex hoyi Sorex monticolus Sorex navigator Microtus oeconomus Myodes rutilus Urocitellus parryii Ochotona collaris Sorex minutissimus Sorex tundrensis Mean no. occurrences Range

Tundra/alpine

64

179

20

50

85

1, 4, 10, 11, 12, 15, 16, 17

8 variables, beta9

Tundra/alpine

87

160

135

54

72

1, 2, 3, 4, 5, 12, 15

7 variables, beta7

Tundra/alpine Tundra/alpine Tundra/alpine Tundra/alpine Tundra/alpine Boreal

68 21 28 177 21 42

168 168 168 168 168 168

132 132 132 97 132 109.9

54 54 54 49 54 49.9

72 72 72 72 72 72

1, 1, 1, 1, 1, 1,

7 7 7 6 7 6

Boreal

139

168

109.9

49.9

72

1, 2, 4, 5, 12, 15

6 variables, beta5

Boreal

78

168

118

52

72

1, 2, 3, 4, 5, 12, 15

7 variables, beta7

Boreal

93

168

125

58

72

1, 2, 3, 4, 5, 12, 15

7 variables, beta3

Boreal Boreal Boreal Boreal Boreal Boreal Boreal Boreal Widespread

55 73 51 125 221 58 185 73 193

168 168 168 168 168 168 168 168 170

126 104 108 84 55 110 110 105 116

55 38 42 40 35 52 42 40 52

72 72 60 72 72 72 70 72 72

1, 1, 1, 1, 1, 1, 1, 1, 1,

2, 2, 2, 2, 2, 2, 2, 2, 2,

3, 4, 3, 4, 4, 4, 5, 4, 3,

4, 5, 12, 15 12, 15, 18 4, 5, 12, 15, 18 12, 15, 18 12, 15 5, 12, 15 7, 12, 15 12, 15, 18 4, 5, 12, 15

7 6 8 6 5 6 6 6 7

variables, variables, variables, variables, variables, variables, variables, variables, variables,

beta7 beta7 beta3 beta3 beta5 beta9 beta5 beta3 beta7

Widespread Widespread Widespread Widespread Widespread

157 135 45 31 77 102

179 140 168 130 168

80 96 126 140 132

55 50 55 54 54

83 72 72 72 72

1, 1, 1, 1, 1,

2, 2, 2, 2, 2,

4, 3, 3, 3, 3,

5, 4, 4, 5, 4,

6 7 7 6 7

variables, variables, variables, variables, variables,

beta11 beta1 beta5 beta1 beta5

2, 2, 2, 2, 2, 2,

3, 3, 3, 4, 3, 4,

4, 4, 4, 5, 4, 5,

5, 12, 15 5, 12, 15 5, 12, 15 12, 15 5, 12, 15 12, 15

12, 15 5, 12, 15 5, 12, 15 12, 15 5, 12, 15

variables, variables, variables, variables, variables, variables,

beta11 beta5 beta5 beta3 beta5 beta5

21–221


23


HOPE ET AL. Table A2. Additional parameters and output statistics for species distribution models. Most important variable (MIV; description)

Species Dicrostonyx groenlandicus Lemmus trimucronatus Microtus miurus Marmota broweri Lepus othus Gulo gulo Sorex ugyunak Microtus longicaudus Microtus pennsylvanicus Synaptomys borealis Tamiasciurus hudsonicus Marmota caligata Zapus hudsonius Zapus princeps Lepus americanus Sorex cinereus Sorex hoyi Sorex monticolus Sorex navigator Microtus oeconomus Myodes rutilus Urocitellus parryii Ochotona collaris Sorex minutissimus Sorex tundrensis Mean AUC Range

Mean temperature of warmest quarter Precipitation seasonality Maximum temperature of warmest month Mean annual temperature Mean diurnal range Mean annual temperature Isothermality Mean annual temperature Mean annual temperature Mean annual temperature Mean annual temperature Mean annual temperature Mean annual temperature Mean annual temperature Mean annual temperature Mean annual temperature Mean annual temperature Temperature annual range Temperature seasonality Mean annual temperature Mean annual temperature Maximum temperature of warmest month Mean annual temperature Mean annual temperature Mean annual temperature


MIV

AUC

LPT threshold

ES&S threshold

Fossil threshold

bio10

0.823

0.254

0.492

0.379

Beringia

bio15 bio5

0.712 0.680

0.294 0.360

0.459 0.493

0.510 0.283

Beringia Rangewide

bio1 bio2 bio1 bio3 bio1 bio1 bio1 bio1 bio1 bio1 bio1 bio1 bio1 bio1 bio7 bio4 bio1 bio1 bio5

0.849 0.967 0.739 0.849 0.621 0.688 0.679 0.782 0.759 0.786 0.838 0.795 0.721 0.695 0.709 0.812 0.727 0.775 0.839

0.314 0.255 0.136 0.349 0.349 0.190 0.204 0.056 0.262 0.082 0.175 0.094 0.129 0.380 0.252 0.129 0.235 0.184 0.085

0.502 0.525 0.448 0.472 0.504 0.494 0.519 0.468 0.482 0.491 0.395 0.398 0.487 0.509 0.482 0.420 0.443 0.465 0.394

... ... 0.521 ... ... 0.396 ... ... 0.610 ... ... ... 0.383 ... ... ... 0.426 0.461 0.404

Rangewide Rangewide Western Beringia Western Western Rangewide Alaska Rangewide Northern Boreal Western Western Western Western Western Beringia Beringia Rangewide

bio1 bio1 bio1

0.774 0.773 0.622 0.744 0.621–0.967

0.426 0.170 0.361

0.519 0.417 0.506

0.352 ... ...

Rangewide Beringia Nearctic

24

Lineage used


HOPE ET AL. Table A3. Parameters and taxon-specific values for approximate Bayesian computation scenario tests in DIYABC. Parameters include G ¼ generation time in years; l ¼ taxon-specific mutation rate (3 102 substitutions/site/ Myr); K ¼ transition/transversion ratio; PInv ¼ proportion of invariant sites; Ne1 ¼ effective population size (3 105 ) at time of maximum population decline; Ne2 ¼ maximum effective population size (3 105 ); Ne3 ¼ minimum effective population size (3 105 ); NeA ¼ average effective population size (3 105 ); NeI and NeD ¼ effective population sizes following the Last Glacial Maximum represented as an increase or decrease, respectively, from Ne2, where magnitude of population size change reflects SDM percentage change in predicted area from present to 2070s timeframes (values taken from Table 1); if no SDMs were available for a taxon, or percentage change was ,20%, then NeI and NeD reflected a standard 6 20% change in population size; t1 ¼ time (ka) of maximum population decline; t2 ¼ time (ka) of maximum population growth; t3 ¼ time (ka) of minimum population size; tR ¼ time (ka) of recent population change following the Last Glacial Maximum. All Ne and t values were calculated from data outputs of taxon-specific Bayesian skyline plots. Non-applicable spaces are denoted by n/a. Values correspond to two sets of three competing scenarios for population size change through time as illustrated in Appendix: Fig. A1. Species

G

là

K§

PInv§

Ne1}

Ne2}

Ne3

NeA

NeI

NeD

t1

t2

t3

tR#

Dicrostonyx groenlandicus Lemmus trimucronatus Microtus longicaudus Microtus miurus Microtus oeconomus Microtus pennsylvanicus Myodes rutilus Zapus hudsonius Zapus princeps Marmota broweri Marmota caligata Tamiasciurus hudsonicus (N) Tamiasciurus hudsonicus (S) Ochotona collaris Gulo gulo Mustela erminea (H) Mustela erminea (C) Sorex cinereus Sorex hoyi Sorex minutissimus Sorex monticolus Sorex navigator Sorex tundrensis Sorex ugyunak

1 1 1 1 1 1 1 1 1 2 2 1 1 1 2 2 2 1 1 1 1 1 1 1

5.60 6.51 8.11 6.73 6.48 7.72 4.98 9.70 5.90 9.15 5.05 5.15 4.77 6.17 5.77 5.17 5.18 5.54 5.11 7.40 6.57 6.86 6.53 6.29

n/a 15.6 5.0 6.0 22.3 25.5 7.7 4.1 n/a 11.1 n/a 21.1 5.4 9.1 n/a 10.6 n/a 14.1 5.7 9.0 4.0 10.1 18.5 n/a

0.81 0.81 0.78 0.76 0.79 0.88 0.80 0.88 0.81 n/a n/a n/a n/a 0.75 n/a n/a n/a 0.68 0.83 n/a 0.76 0.80 0.75 0.82

3.30 5.01 5.25 2.34 1.63 3.28 15.0 2.46 3.73 1.66 1.89 4.27 6.23 2.79 1.56 4.67 2.49 7.62 2.79 2.71 4.96 3.78 5.03 4.03

3.46 5.26 5.56 7.85 2.61 3.42 15.1 2.46 3.73 1.66 1.89 5.02 6.39 6.27 1.56 5.23 2.49 9.95 2.89 2.71 4.96 4.19 5.32 4.03

1.61 1.14 0.62 1.31 0.86 0.67 0.55 0.42 0.74 0.21 1.05 0.45 0.53 0.57 0.42 0.80 0.69 0.27 0.59 0.32 0.52 0.43 0.48 0.37

2.35 2.72 3.10 2.52 1.62 1.67 6.97 1.18 1.69 0.49 1.38 2.56 2.70 3.27 0.75 2.55 1.28 5.01 1.63 1.36 2.39 2.31 2.42 1.35

4.75 5.84 7.90 10.3 3.26 4.67 19.4 3.76 4.66 2.04 2.36 8.04 7.98 8.49 2.09 6.53 3.11 14.3 3.39 3.40 6.20 5.73 6.65 5.04

2.17 4.68 3.22 5.44 2.09 2.17 10.8 1.16 2.80 1.28 1.51 2.00 5.11 4.05 1.03 4.18 1.99 5.58 2.39 2.17 3.97 2.65 4.26 3.23

2.7 2.6 2.6 4.1 8.3 1.0 2.3 n/a n/a n/a n/a 2.9 1.2 6.7 n/a 1.2 n/a 2.3 1.8 n/a n/a 1.8 1.3 n/a

32.6 24.1 22.7 16.2 35.1 12.6 24.0 9.2 9.2 4.0 10.5 18.2 13.2 34.3 5.7 24.2 8.9 15.7 21.4 11.0 16.0 26.8 13.2 7.9

100.5 94.6 55.6 88.9 64.0 45.5 62.9 30.2 53.0 8.9 29.4 40.2 41.9 64.1 14.5 63.5 34.3 29.2 40.1 21.4 41.3 42.3 40.5 24.8

12.0 12.0 12.0 12.0 12.0 12.0 12.0 6.0 6.0 2.0 8.0 12.0 12.0 12.0 2.0 12.0 6.0 12.0 12.0 6.0 12.0 12.0 12.0 6.0

For taxa with generation times of two years, all values of Ne and t were halved for simulations. à Mutation rates were calculated previously (Hope et al. 2014) and calibrated to a reference rate of 5.6 3 102 substitutions/ site/Myr applied to Dicrostonyx (Prost et al. 2010). § Empirical values were recorded from MRMODELTEST results. For all models of evolution that did not incorporate rate heterogeneity or transition/transversion ratios in their outputs (n/a), K was set to 10 and Pinv was set to 0.75 (default mean values). } Inflection points (maximum population changes) from each Bayesian skyline plot were calculated from changes in median population size (Ne) based on outputs partitioned into 100 bins (where dNe ¼ Nebinx – Nebinxþ1), plotted against time, and retaining maximum positive and negative values for dNe. Where no recent population decline was detected from Bayesian skyline plots, Ne1 ¼ Ne2. # Recent population change was dated for most taxa to 12 ka, coincident with onset of major climate following the Last Glacial Maximum. More recent dates were necessary for taxa where Ne2 , 12 ka. All values for NeI and NeD were conceived to statistically compare different possible scenarios of recent population size changes as a proxy for possible future trajectories under a continued warming trend with the assumption that actual values encompass informative genetic change.


25


HOPE ET AL.

Fig. A1. All DIY-ABC scenarios, including parameters and scripts for defining models. (A) Bayesian skyline scenarios for change in population size tested for support of population size changes estimated previously from Bayesian skyline plots (BSPs) for each taxon (Appendix: Fig. A5). BSPs estimated in BEAST for all taxa conformed to one of two general shapes, either exhibiting growth followed by a recent decline (left curve) or growth with no subsequent decline up to the present (right curve). ABC tests under scenario BSP1 incorporated empirical parameter values calculated from BEAST Bayesian skyline plot outputs. Tests under scenario BSP2 reversed empirical population sizes for Ne2 and Ne3. Tests under scenario BSP3 assumed a constant average population size through time (substituted in NeA). (B) Post-Last Glacial Maximum (LGM) scenarios for population size change through time tested for support of either an increase (LGM1), decrease (LGM2), or no change (LGM3) in population size since the LGM as a potential indicator of contemporary trends that may continue into future timeframes. Results of ABC tests are reported in Table 2. Values for all parameters are reported in Appendix: Table A3.


26


HOPE ET AL.

Fig. A2. All species distribution models for individual taxa, based on (A) Last Interglacial (LIG), (B) Last Glacial Maximum (LGM), (C) current (Now), (D) 2050s, and (E) 2070s climate projections. Taxa shown in blue have a Nearctic evolutionary origin and are associated with the boreal biome. Taxa shown in red constitute Beringian lineages associated with tundra habitats. Taxa shown in green have Beringian or Palearctic evolutionary associations but are generally widespread within the study area. Future projections are shaded according to thresholds based on a predictive scale of 0–3 (light to dark) where darkest areas represent most conservative predictions.


27


HOPE ET AL.

Fig. A2. Continued.


28


HOPE ET AL.

Fig. A2. Continued.


29


HOPE ET AL.

Fig. A2. Continued.


30


HOPE ET AL.

Fig. A2. Continued.


31


HOPE ET AL.

Fig. A2. Continued.


32


HOPE ET AL.

Fig. A2. Continued.


33


HOPE ET AL.

Fig. A2. Continued.


34


HOPE ET AL.

Fig. A2. Continued.


35


HOPE ET AL.

Fig. A2. Continued.


36


HOPE ET AL.

Fig. A2. Continued.


37


HOPE ET AL.

Fig. A2. Continued.


38


HOPE ET AL.

Fig. A2. Continued.


39


HOPE ET AL.

Fig. A2. Continued.


40


HOPE ET AL.

Fig. A2. Continued.


41


HOPE ET AL.

Fig. A2. Continued.


42


HOPE ET AL.

Fig. A2. Continued.


43


HOPE ET AL.

Fig. A2. Continued.


44


HOPE ET AL.

Fig. A2. Continued.


45


HOPE ET AL.

Fig. A2. Continued.


46


HOPE ET AL.

Fig. A2. Continued.


47


HOPE ET AL.

Fig. A2. Continued.


48


HOPE ET AL.

Fig. A2. Continued.


49


HOPE ET AL.

Fig. A2. Continued.


50


HOPE ET AL.

Fig. A2. Continued.


51


HOPE ET AL.

Fig. A3. Species distribution model predictions for combined widespread taxa based on (A) current (Now), (B) 2050s, and (C) 2070s climate projections. The color gradient reflects areas of low (light) to high (dark) diversity. Distribution of diversity reflects compilation of predictive maps for 6 widespread species, each thresholded as presence/absence.


52


HOPE ET AL.

Fig. A4. Results from varying numbers of randomly sampled locality records on species distribution model predictions. (A–E) Predictions for Sorex cinereus (blue). (F–J) Predictions for Lemmus trimucronatus (red). For each species, four maps are shown for each of five projected timeframes including Last Interglacial (LIG), Last Glacial Maximum (LGM), current (Now), 2050s, and 2070s. The four SDM maps for each timeframe represent modeled distributions using 20, 30, 40, or 50 randomly sampled occurrence records, respectively, for each species. Future projections are shaded according to thresholds based on a predictive scale of 0-3 (light to dark) where darkest areas represent most conservative predictions. Summary data for all models are available from E. Waltari upon request.


53


HOPE ET AL.

Fig. A4. Continued.


54


HOPE ET AL.

Fig. A4. Continued.


55


HOPE ET AL.

Fig. A4. Continued.


56


HOPE ET AL.

Fig. A4. Continued.


57


HOPE ET AL.

Fig. A4. Continued.


58


HOPE ET AL.

Fig. A4. Continued.


59


HOPE ET AL.

Fig. A4. Continued.


60


HOPE ET AL.

Fig. A4. Continued.


61


HOPE ET AL.

Fig. A4. Continued.


62


HOPE ET AL.

Fig. A5. Bayesian skyline plots for all taxa incorporating taxon-specific substitution rates (adapted from Hope et al. 2014). Solid lines represent median effective population size change through time (past on the right to present on the left). Shaded interval represents the 95% highest probability distribution. The x-axis is time (thousands of years ago); the y-axis is log-transformed and represents s as a function of generation time and effective population size.


63


HOPE ET AL.

Fig. A5. Continued.


64


HOPE ET AL.

Fig. A5. Continued.


65


HOPE ET AL.

Fig. A6. Species distribution model predictions for (A) and (B) boreal taxa, (C) and (D) tundra taxa, and (E) and (F) widespread taxa for the Last Interglacial warm period (LIG; left column) and Last Glacial Maximum (LGM; right column). The color gradients reflect areas of low (light) to high (dark) diversity. Distribution of diversity reflects compilation of predictive maps for 12 boreal, 7 tundra/alpine, and 6 widespread species, each thresholded as presence/absence.


66


HOPE ET AL.

Fig. A7. Plot of species richness based on species distribution model predictions within the Arctic Network of National Parks (ARCN; squares), Arctic National Wildlife Refuge (Arctic NWR; circles), and National Petroleum Reserve–Alaska (NPR-A; triangles) for all timeframes including Last Interglacial (LIG), Last Glacial Maximum (LGM), present (Now), 2050s, and 2070s. Black lines consider all 25 taxa, blue lines consider only boreal taxa (Table 1, B), and red lines consider only tundra taxa (Table 1, T).


67
