rspb.royalsocietypublishing.org

Research Cite this article: Sultaire SM, Pauli JN, Martin KJ, Meyer MW, Notaro M, Zuckerberg B. 2016 Climate change surpasses land-use change in the contracting range boundary of a winteradapted mammal. Proc. R. Soc. B 283: 20153104. http://dx.doi.org/10.1098/rspb.2015.3104

Received: 28 December 2015 Accepted: 1 March 2016

Subject Areas: ecology, environmental science Keywords: forest fragmentation, global change, mammal distribution, range boundary, snow cover, trailing edge

Climate change surpasses land-use change in the contracting range boundary of a winter-adapted mammal Sean M. Sultaire1, Jonathan N. Pauli1, Karl J. Martin2, Michael W. Meyer3, Michael Notaro4 and Benjamin Zuckerberg1 1 Department of Forest and Wildlife Ecology, University of Wisconsin-Madison, 1630 Linden Dr. Madison, WI 53706-1598, USA 2 Community, Natural Resource, and Economic Development Program, University of Wisconsin-Extension, 432 N. Lake St Madison, WI 53706, USA 3 Wisconsin Department of Natural Resources, 107 Sutliff Ave, Rhinelander, WI 54501, USA 4 Center for Climactic Research, University of Wisconsin-Madison, 1225 West Dayton St Madison, WI 53706, USA

The effects of climate change on biodiversity have emerged as a dominant theme in conservation biology, possibly eclipsing concern over habitat loss in recent years. The extent to which this shifting focus has tracked the most eminent threats to biodiversity is not well documented. We investigated the mechanisms driving shifts in the southern range boundary of a forest and snow cover specialist, the snowshoe hare, to explore how its range boundary has responded to shifting rates of climate and land cover change over time. We found that although both forest and snow cover contributed to the historical range boundary, the current duration of snow cover best explains the most recent northward shift, while forest cover has declined in relative importance. In this respect, the southern range boundary of snowshoe hares has mirrored the focus of conservation research; first habitat loss and fragmentation was the stronger environmental constraint, but climate change has now become the main threat. Projections of future range shifts show that climate change, and associated snow cover loss, will continue to be the major driver of this species’ range loss into the future.

1. Introduction Author for correspondence: Benjamin Zuckerberg e-mail: [email protected]

Electronic supplementary material is available at http://dx.doi.org/10.1098/rspb.2015.3104 or via http://rspb.royalsocietypublishing.org.

Habitat loss and fragmentation have been a primary focus of conservation research for decades [1,2], and consequences such as reduced population size and connectivity leading to increased extinction probabilities have been studied extensively [3,4]. In recent decades, however, the emphasis on habitat loss as the primary threat to biodiversity has waned in the face of increasing concern over anthropogenic climate change [5]. In many regions characterized by high climate change velocity, the ecological impacts of habitat fragmentation appear to be exacerbated [6], and theory supports a potent synergy between these threats [7]. However, if habitat conversion is biased towards warmer areas on the landscape, then habitat loss may confound the effects of climate change on a species range [8]. Studies evaluating the relative effects of habitat loss and climate change so far indicate that habitat loss remains a more important population driver than climate change for many species [9,10], although this may be changing in some systems [11]. Both factors can lead to shifts in species ranges [12,13], but the response of species’ distributions to shifting contributions of these threats over broad timescales is less certain. Range boundaries offer unique opportunities to study the effects of climate change and habitat loss on species distributions as these populations are likely most sensitive to a changing climate [14], and habitat is often patchy [15]. Studies quantifying the distributional responses of species to climate change generally focus on two types of range boundaries: cold, leading edge boundaries, and warm, trailing edge boundaries [16]. Leading edge range boundaries are the result of population processes, such as colonization, that can be inhibited by

& 2016 The Author(s) Published by the Royal Society. All rights reserved.

later by a reduction in the snow cover season. Under this hypothesis, we predicted that snowshoe hare occupancy at the historical range boundary would coincide with the distribution of forest cover at that time period, and snow cover played a comparatively small role. Furthermore, we predicted that contemporary changes in occupancy (primarily extinction) to be driven by snow cover conditions along the boundary, while the contribution of forest cover has declined.

2. Material and methods The southern distribution of snowshoe hares in the US state of Wisconsin was first mapped using expert opinion in 1944 [31], and later more comprehensively via snow-tracking surveys in 1979–1980 [32]. For both these studies, the survey area was approximately 100 km north–south by approximately 400 km, and coincided with a transition zone between southern hardwoods and northern mixed forest [36]. Historically, the eastern portion of the snowshoe hare range boundary in Wisconsin was considerably north of this ecotone, attributed to high levels of agricultural conversion in the area. Topographic relief within the region is low, with all areas less than 600 m in elevation. We used these historical surveys to quantify the range shift of snowshoe hares in response to changing environmental conditions at the range boundary. For the more recent historical survey conducted in the winters of 1979–1980, the study area was stratified into 92.16 km2 blocks and at least one location was surveyed within each block [32]. Historic locations were provided at a 2.56 km2 resolution and we conducted surveys within or as close to these locations as possible (i.e. less than 1.60 km). As a primary goal of this study was to fully quantify the latitudinal shift of the range boundary, we surveyed additional sites north of historical sites. To select these sites, we used the same stratification scheme as the historical survey [32], and within each strata one section was randomly chosen from the five most forested 2.56 km2 sections.

(b) Field methods We conducted snow track surveys between December and March during the winters of 2012–2013 and 2013 –2014. Historical surveys were conducted between January and March. Consistent with historical surveys, contemporary surveys consisted of a standardized 1 h search through forested habitat [32]. No information was available on the spatial configuration of historic surveys so we walked a series of 125 m transects spaced 75 m apart to maximize detectability [37]. During each survey period, we recorded the detection of each set of hare tracks encountered along transects, and pooled data across all transects for analysis. To account for imperfect detection, we used an occupancy modelling framework with repeat surveys at a subset of sites (see [e] occupancy modelling section) [38]. During each survey, we collected data on snow density and depth (cm) to quantify their effect on snowshoe hare detection probability. We measured snow surface density by dropping a 200 g cylindrical penetrometer from a height of 50 cm and measuring sinking depth (cm).

(c) Land cover analysis To characterize historical land cover conditions, we used LandSat Multi-spectral scanner imagery from 1977 to 1981 and created a 2-class map showing forested and non-forested areas during the historical survey period at a 60 m resolution. For contemporary land cover, we used National Landcover Database 2011 at a 30 m resolution and combined all forest classes into a single forest class to represent total forest and all other classes as non-forest. To capture landscape-scale forest cover, we calculated the total area of

Proc. R. Soc. B 283: 20153104

(a) Study area and site selection

2

rspb.royalsocietypublishing.org

habitat fragmentation, preventing a species from tracking shifting and potentially more favourable climate conditions [13,17]. Alternatively, the dynamics of trailing edge range boundaries are the culmination of multiple localized extirpations, as areas at the edge of a species’ climate niche become unsuitable [18,19]. If a species’ trailing edge range boundary has shifted north as a result of land-use change, then these populations may be less sensitive to more subtle changes in climate [8,20]. With increased warming, however, the mechanisms driving these trailing edge range boundaries may shift away from habitat availability as the direct impacts of climate change become increasingly pronounced. Snowshoe hares (Lepus americanus) are a dominant species of North America’s northern forests, and their population dynamics and habitat associations have been intensively studied for nearly a century [21]. This wealth of research has shown that snowshoe hares are more likely to occupy, and reach greater abundances, in forest stands characterized by high levels of vegetative cover, embedded within landscapes that are primarily forested [22–24]. Historically, the availability of forested habitats was considered an important component of this species’ southerly distribution and central in mediating predator-induced mortality [25,26]. In addition to their association with complex forests, snowshoe hares are adapted to cold, snow-dominated environments. Notably, snowshoe hares molt to white in the winter to match background snow conditions, and this molt is a fairly rigid phenotypic trait initiated by photoperiod rather than the presence of snow cover [27,28]. As the Northern Hemispheric snow season attenuates [29], the phenological mismatch between seasonal coat colour change and background snow conditions is becoming more pronounced [27], and the mismatch has strong selective costs on survival [30]. The extent to which long-term changes in snow cover has influenced the distribution of snowshoe hares has not been quantified. The central portion of snowshoe hare range reaches its southern extent in the Great Lakes region of the USA. Historically associated with a forest ecotone between southern deciduous and northern mixed forest, snowshoe hare distribution receded north by the early twentieth century, mainly attributed to the rapid conversion of forest habitat to agriculture during this time period [31]. Comparisons of the range boundary between the mid- and late-twentieth century revealed a dynamic range boundary without a clear directional shift [32]. Just north of this boundary, snowshoe hare populations appear to be in a long-term decline since the late 1970s, with population cycles dampening or collapsing entirely (Wisconsin Department of Natural Resources, electronic supplementary material, figure S1). Over the last 50 years, the Great Lakes region has experienced a relatively high magnitude of climate change [33], and the most pronounced changes have occurred during winter months [34]. During this time period, land cover in the region has remained fairly static and the conversion of forest has been minimal [35]. The availability of detailed historical datasets on the southern distribution of this forest and winter climate specialist offers a unique opportunity to track this species’ response over decades of environmental change. We revisited historical locations (ca 1980) along the southern extent of snowshoe hare distribution to quantify the response of their trailing edge to changing climate and land cover conditions. We tested the hypothesis that the southern range boundary of snowshoe hares has experienced a twophase range shift—first driven by rapid land-use change and

To characterize winter snow cover for each survey period, we used National Climatic Data Center data from 692 weather stations in the region to create a daily snow depth map (mm) from 1950 to 2013 at 0.18 ( 13  18 km) resolution. Snow cover for each 0.18 grid cell was computed from the 10 closest weather stations using inverse distance weighting. We investigated the duration of persistent snow cover as the major climatic driver on the range boundary due to the potential effects of coat colour mismatch [27,28,30]. For each site and year surveyed, we calculated the date of persistent snow cover onset (Snowonset), date of persistent snow cover termination (Snowterm), and duration of persistent snow cover (Snowdur). We defined onset of persistent snow cover as the first day of the first 14 day sequence with snow cover more than 1 cm, and termination as the final day of the last such 14 day period. Persistent snow cover duration was calculated as the time interval between these two dates. We averaged these variables over the five years prior to the survey of a site for each time period. Although we had no strong a priori reason to conclude that this temporal scale is most ecologically appropriate, average snow cover duration was highly correlated at a range of temporal scales (electronic supplementary material, figure S2), and a five-year period best reflected the scale this short-lived species likely responds to its environment. We further calculated the relative change in each variable between time periods (OnsetD,DurD, TermD). We calculated future snow cover (2040–2059) based on forcing an operational snow model with statistically downscaled projections of air temperature and precipitation for eastern North America [29]. For future range shift projections, we used the ensemble median snow cover values from nine different general circulation models for both the A1B and A2 emissions scenarios at mid-twenty-first century. We chose these emissions scenarios as the current global economic path makes reduced emission scenarios increasingly unlikely [41].

(e) Occupancy modelling Comparing current and historic survey data to infer range shifts is inherently problematic [42]. One of the challenges in using

3

Proc. R. Soc. B 283: 20153104

(d) Climate data

historical data is lack of information on species’ detection probability for historical surveys, and failing to account for this heterogeneity can lead to biases in range boundary shift estimates. To account for this potential source of error, we used a multi-season occupancy model to compare survey data between time periods, account for imperfect detection, and estimate probabilities of colonization and extinction [43]. Although no repeat surveys were conducted for the historical time period, we used the current detection probability– covariate relationships to account for imperfect detection during historical surveys. Our contemporary survey dataset included 199 sites, of which 48 sites were surveyed three times, 21 sites twice, and the remainder surveyed once, and 53 sites detected snowshoe hares. In total, 148 of these sites were used in our multi-season model, and the remaining sites were used in model validation (see below). Analysing the contemporary dataset in a single-season occupancy modelling framework [38] indicated that surface snow density was the most important variable driving the detection of snowshoe hares (electronic supplementary material, table S1). Although snow cover characteristics were not collected during historic surveys, we found that snow depth (ObsDepth), Julian date (ObsDate), and maximum temperature since the last snowfall of more than 1 cm (ObsTmax) explained 60% of the variation in snow density measurements collected during current surveys. We used the Julian date that historical surveys were conducted and weather station data to calculate these values for historical surveys and fit them as model covariates for historical surveys. Because other sources of variation between time periods can influence detection probability [42], we also evaluated the robustness of our inference to changes in detection probability between time periods. If detection probability was higher during the historical survey era, applying current detection relationships to past surveys could result in an overestimation of the range shift. To account for this potential Type 1 error, we ran occupancy models with historical detection probability fixed at 1 (perfect detection), and a model that assumed the historical detection probability to be 0.5 (lower than contemporary detection), as well as applying contemporary detection–covariate relationships to historical surveys. We limited our multi-season comparison to quantifying states of historical occupancy probability (c1) and extinction probability (1). We did not seek to parametrize colonization dynamics ( g) due to a lack of colonization events (see Results). We fitted models testing the effects of historical snow cover, forest cover, and housing variables on c1. Extinction was parametrized with the snow cover change covariates, contemporary snow cover conditions, contemporary land cover, forest disturbance, and housing change variables. We also fitted a single-season occupancy model to the contemporary survey data with the same covariates. This enabled us to make a better comparison of the relative importance of each parameter on c for each time period, and test our hypothesis of shifting mechanisms operating on the range boundary. Prior to fitting any models, all covariates were standardized with a mean of zero and standard deviation of one. We used Akaike’s Information Criterion (AIC) to perform model selection and rank models based on AIC weights. We used the top-ranked multi-season model and covariate information to make spatial predictions of c1 and 1 across Wisconsin at the resolution of our site selection (2.56 km2). Occupancy in the modern era (c2) was derived using values of c1, 1, and g and modelled at the same resolution. We evaluated our distribution model for c2 by calculating the Area-Under-theCurve (AUC) [44] using the 65 supplementary survey sites north of historical sites as a validation dataset. We made future projections under a snow cover change model ( parametrized with current forest cover and A2 snow cover projections) and a forest cover change model ( parametrized with current snow cover and baseline forest cover projections) to examine which variable is most likely to drive the future range of snowshoe

rspb.royalsocietypublishing.org

forest within a 5 km buffer around survey sites for each time period (Forest1, Forest2); this landscape scale is important to snowshoe hare occupancy dynamics elsewhere in their southern range [22] and approximates the maximum dispersal distance documented for hares in our study area [26]. In addition, we calculated the number of houses within the 5 km buffer for 1980 (Houses1), 2010 (Houses2), and housing change (HouseD) using an existing housing density database [39] to model the effect of exurban development on occupancy and range boundary dynamics. For future occupancy projections, future land cover (2050) was parametrized with economic-based land-use change projections [40], which map the probability that a 1 ha cell will transition to one of five land cover types in 2050: forest, cropland, urban, pasture, or range. To summarize these data into future forest cover, we summed the probability of each 1 ha cell being forest within a 5 km buffer round each occupancy grid cell (2.56 km2), and divided by the buffer area, yielding the per cent forest cover for each cell in 2050. As snowshoe hares most often occur in forests characterized with understory cover, changes not only in the amount of forest over time, but also in forest age and structure could affect hare occupancy. To account for the influence of forest disturbance and succession on range dynamics, we used a forest disturbance dataset that maps all canopy disturbances that regenerated to forest in the region since 1990 [39]. Using this dataset, we summed the area of all disturbed forest in the 5 km buffer around survey sites, allowing us to quantify the effect of forest succession dynamics on snowshoe hare extinction. These values were then log-transformed as the distribution was highly right skewed.

(a)

(b)

4

rspb.royalsocietypublishing.org

0

(c)

45

90 km

(d)

Figure 1. The southern range boundary of the snowshoe hare terminates in Wisconsin, USA (a). Leopold [31] conducted anecdotal surveys of snowshoe hare presences (b). These surveys were followed by more rigorous snow-tracking surveys by Buehler & Keith [32] who documented detections (black) and non-detections (white) in the same region (c). Contemporary survey locations documenting sites where hares were detected (grey), detected in historical but not current surveys (black), and sites where hares were not detected in either current or historical surveys (white). The background represents forest cover (credit J. M. Rhemtulla and D. J. Mladenoff for 1940s era land cover).

hare in Wisconsin. We compared the magnitude of these future projections by summing occupancy probability across Wisconsin for each future projection and subtracting the sum of the snow cover change projection from that of the forest change model. Occupancy modelling was performed using the package unmarked [45] in the R statistical platform [46].

(f ) Range shift We used pairings of historical and contemporary snowshoe hare detections within the same longitudinal band of survey strata to quantify the magnitude of the range shift between time periods. For each of these longitudinal bands, we averaged the latitude (metres) of the most southerly detection for each survey period, and calculated the difference between these two values as the magnitude of the range shift for both the entire inter-survey period and decadal scales. We excluded sites in longitudinal bands that did not have detections in both time periods for this calculation.

3. Results We revisited 148 of 249 historic survey sites: 126 of the 151 historical detections, and 22 of the 98 historical non-detections,

with repeat surveys at 38 sites total. At the 126 historic survey sites that detected snowshoe hares, we detected hares at 28 sites (22%) and failed to detect hares at 98 sites (78%) (figure 1). We only detected hares at one site where hares were not detected historically. Between the two survey periods, snow cover duration decreased across the study area (min: 232.9, max: 26.4 days) (electronic supplementary material, table S2), while between-site variation in snow cover doubled (electronic supplementary material, table S2). The cross tabulation used to validate the historical land cover map showed an 85% concordance between the stratified random sample and our ground-truth points (95% CI: 79–90%). Using this map, we found that forest cover was static over the inter-survey time period with no observable change in the overall mean forest cover while housing density increased (electronic supplementary material, table S2), supporting findings from previous studies [35,39]. Including the detection covariates ObsTmax, ObsDepth, and ObsDate increased the model fit over the base detection model and was used in all subsequent models (electronic supplementary material, table S3). The base (unparametrized) multi-season model estimated extinction probability between

Proc. R. Soc. B 283: 20153104

N

(a)

(b)

5

rspb.royalsocietypublishing.org

historical contemporary

1.0

Proc. R. Soc. B 283: 20153104

probability

0.8 0.6 0.4 0.2 0 90

100 110 120 130 snow cover duration (days)

0.2

0.4 0.6 0.8 proportion of forest

Figure 2. Relationships between occupancy dynamics and snow and forest cover. We present the predicted relationships between historical occupancy (solid line) and extinction probability between the time periods (dashed line) as a function of snow cover duration (a) and forest cover (b). Histograms display the distribution of snow and forest cover at study sites for each survey period. time periods to be high (0.73, 95% CI: 0.63–0.81). The multiseason model with the strongest support (wi ¼ 0.89), given the data, included the effect of Snowdur1, Forest1, and Houses1 on c1 (electronic supplementary material, table S4). We found that, during the historical time period, hares were more likely to occupy sites in landscapes that had high forest cover (b ¼ 1.55, standard error (s.e.) ¼ 0.44), longer duration of snow cover (b ¼ 1.26, s.e. ¼ 0.36), and lower housing density (b ¼ 20.63, s.e. ¼ 0.37) (electronic supplementary material, table S4). For every 19% (s.d.) increase in forest cover c1increased fourfold, while every 3.25 day (s.d.) increase in snow cover resulted in a threefold increase in c1. Models with the individual effects of forest and snow cover on initial occupancy had no support (electronic supplementary material, table S4). The top-ranked model contained the effects of Snowdur2 and Forest2 on extinction probability. Between survey periods, we found that snowshoe hares were more likely to become extirpated in landscapes with low forest cover (b ¼ 20.78, s.e. ¼ 0.31) and shorter snow cover duration (b ¼ 21.34, s.e. ¼ 0.47; electronic supplementary material, table S4; figure 2). Snowshoe hares were four times more likely to persist for every 7.41 day (s.d.) increase in snow cover duration and twice as likely to persist for every 23% (s.d.) increase in forest cover. For both initial occupancy and extinction, we found stronger support (DAIC , 2) for models containing the effect of Snowdur than either models containing Snowterm or Snowonset (electronic supplementary material, table S4). The second-ranked model was identical to the top-ranked model but contained DurD as an extinction parameter, as opposed to Snowdur2 (electronic supplementary material, table S5). As predicted, the coefficient estimate for the effect of landscapescale forest disturbance was negative (b ¼ 20.67, s.e. ¼ 0.30), but it did not outcompete models that contained only forest cover as an extinction covariate, and so was not retained in further models.

Table 1. Coefficient estimates (and s.e.) for forest cover and snow cover duration in the top-ranked multi-season occupancy model (c1, 1), and the top-ranked single-season model (c2). parameter

forest cover

snow cover duration

c1 1

1.55 (0.44) 20.78 (0.31)

1.26 (0.36) 21.34 (0.47)

c2

0.85 (0.28)

1.13 (0.39)

We found that modelling perfect detectability for historical surveys had no influence on model selection and did not change the relative rank between model coefficients (electronic supplementary material, table S6), suggesting our results were robust to the possibility of higher historical detection probability. Conditional predictions of historical occupancy for the higher ( p ¼ 1) and lower ( p ¼ 0.5) historical detection models differed as expected, but predictions of extinction were not sensitive to varying historical detection probability (electronic supplementary material, figure S3). We parametrized our single-season occupancy model with the detection/non-detection data from the 148 resurveyed sites. The single-season model with strongest support for c2 included the variables Snowdur2 and Forest2 with a larger coefficient estimate for snow cover duration than forest cover (table 1), which is consistent with 1 parameters in the topranked multi-season model. Comparing the ratio between c coefficients across time periods demonstrated a larger snow to forest effect on c2 (1.13:0.85) than c1 (1.26:1.55). The AUC value for the spatial predictions of c2 from the top-ranked model and the independent validation data was 0.79 (s.d. ¼ 0.06) compared to an AUC of 0.77 (s.d. ¼ 0.06) for the second-ranked model. This indicates our top model had good explanatory power for determining hare presence

(a)

(b)

6

rspb.royalsocietypublishing.org

1.0

Proc. R. Soc. B 283: 20153104

N 0.5 0

50 100 km

0

(c)

(d)

Figure 3. Spatial predictions of the probability of snowshoe hare occupancy in historical (a) and contemporary (b) time periods. Using the baseline land-use change scenario, we projected future snowshoe hare distribution to the mid-century assuming no change in snow cover dynamics (c). We also projected future distributions assuming a static land cover (no land-use change), under the A2 emissions scenario (d ). We found a more pronounced northward shift under dynamic snow cover projections compared to the dynamic forest cover projections, indicating that snow cover will remain the primary driver of snowshoe hare distribution shifts into the future. (Online version in colour.)

for the current survey period (electronic supplementary material, figure S4) although the performance of the secondranked model was comparable. We observed an average northward shift of 29.6 km (s.d. ¼ 14.3 km) between 1980 and 2014, representing an 8.7 km per decade northward shift. Comparisons between the snow cover only change projections and the forest cover only change projections indicate an 11% higher reduction in mid-century-occupancy in response to future changes in snow cover as opposed to forest cover (figure 3).

4. Discussion The trailing range boundary of snowshoe hares in Wisconsin has continued to retract northward during the past three decades, and our results indicate that the most recent shift is primarily driven by a reduction in snow cover. Our findings lend support to the survival consequences of coat colour mismatch documented in southerly populations of snowshoe hares elsewhere [30]; and suggest that phenotypic mismatch associated with an altered winter climate is an

important population driver across the southern distribution of a winter-adapted mammal. As the current level of snow cover duration predicts snowshoe hare extinction better than change in snow cover duration between time periods, our results also support the impact of an increasingly novel and reduced snow cover regime along the range boundary. The current season of persistent snow cover in Wisconsin is much shorter than that described for Montana snowshoe hares [27,28], as might be expected at a southern range boundary. To the extent that our observed range contraction is driven by mortality due to coat colour mismatch, it appears that local adaptation of colour molt phenology is not keeping pace with recent and prominent declines in snow cover. Contrary to our predictions of a limited role of snow cover historically, our models had approximately equal support for snow cover and forest cover in driving historical occupancy, and a strong positive relationship between snow cover and historical occupancy (figure 2). However, only a few of the historical sites had snow cover duration less than 110 days—an apparent threshold above which extinction probability decreased markedly (figure 2a). This suggests that while low snow cover may have constrained the

(b)

7

0.1

rspb.royalsocietypublishing.org

snow cover duration anomaly

(a)

0 –0.1 –0.2

1980 year

2000

1960

1980 year

2000

Figure 4. Five-year (a) and 10-year (b) average proportional change in days with snow cover relative to a 1950– 1980 baseline, prior to historical surveys (1980) and the contemporary survey period, for our study area. Dashed lines represent least-squares regression for each time period. We found a dramatic reduction in the duration of snow cover prior to the contemporary survey period. (Online version in colour.) historical occupancy of hares in a few localities, most sites had snow cover duration at levels where predicted snowshoe hare occupancy was high (figure 2a), and hence, not limiting their occurrence. This conclusion is fairly robust to the possibility of lower or higher historical detection probability, particularly with respect to extinction (electronic supplementary material, figure S4). The predicted relationship between forest cover and historical occupancy show a more consistent pattern of increasing occupancy along a gradient of forest cover, and a shallower extinction curve between time periods (figure 2), supporting our hypothesis of a declining importance of forest cover to the range boundary. During the decades leading up to the historical survey, snow cover generally increased across the study area (figure 4), and although snowshoe hares did expand southward in some areas, they also lost range during this time period [32]. Conversely, average snow cover has decreased considerably since the more recent historical survey and is at its lowest levels since 1950. Consistent with this region-wide reduction in snow cover, we identified no regions of extensive range expansion since 1980. Future projections of snowshoe hare distributions under different land cover and climate change scenarios suggests this pattern will persist, and midtwenty-first century projections for snowshoe hare distribution in Wisconsin predict a continued northward shift in the range boundary, which was qualitatively similar between emissions scenarios (electronic supplementary material, figure S5). In the future, it is likely that decreasing snow cover will continue to drive a northward shift in this range boundary, while the role of changing land use will diminish (figure 3). Snow cover dynamics are rarely integrated into species distribution models and vulnerability assessments, but are likely a critical determinant of future distributions for winter-adapted species. Although we found a reduced role of forest cover in shaping current occupancy, this conclusion assumes that landscape-scale forest estimates accurately captured habitat availability. Snowshoe hares require forested or shrubland habitat, but they are most often associated with young forests that provide dense cover from predators [24]. We do not know the distribution of young forest historically, but patterns of snowshoe hare extinction do not correlate with landscapescale patterns of forest disturbance and regeneration in the 20 years preceding our survey. One potentially important factor in our study area are changes in the distribution and abundance of predator species. Fisher (Pekania penannti), which were rare during historical surveys, are now quite common

in the region. Predator communities can drive shifts in prey species’ ranges, particularly when in conjunction with environmental perturbations like climate change [47]. Although predation is likely central to this range shift, the relationship between extinction and snow cover duration would suggest that the conditions that mediate predation are more important than predator community composition. Another potentially confounding process is the highly cyclic nature of snowshoe hare populations that could influence the magnitude of range contraction [21]. There is evidence that snowshoe hares were particularly abundant in Wisconsin during historical surveys, but statewide monitoring surveys indicate that dampened cycling has occurred and population counts have steadily declined since the historical survey (electronic supplementary material, figure S1). As such, it is unlikely that range contraction is part of a broader pattern of population cycling (i.e. contemporary sampling concurrent with a low point in the cycle), but rather that dampened cycling is a prelude to chronic population declines and subsequent range retraction [48]. If cyclic decline is responsible for some degree of the observed range contraction, the strong relationship between extinction dynamics and declining snow cover conditions would suggest a continued erosion of cycling in the future. The duration and variability of snow cover along the southern range limit of snowshoe hares have changed dramatically, and although populations of snowshoe hares are not tracking these changes lockstep, this species’ distribution has tracked the emerging paradigm of conservation biology. For decades, conservation biologists focused on habitat loss, which was indeed the primary factor driving snowshoe hare range in the Great Lakes region during this time period. In recent decades, however, the regional distribution of snowshoe hares has become increasingly limited by winter climate, and at the same time, climate change has moved into the spotlight of conservation research. Although snowshoe hares are likely a highly sensitive species to a changing winter climate, our results support a shifting conservation focus. However, forces of climate and land-use change do not exist in separate vacuums, and our findings highlight this synergistic relationship; extinction was not only the result of decreased snow cover, but was also mediated by forest cover. Three decades ago, when the focus of conservation was squarely trained on the impacts of habitat loss, few likely anticipated the threat of climate change looming on the horizon, and it is unlikely that the conservation issues we are currently facing will be

Proc. R. Soc. B 283: 20153104

1960

Authors’ contributions. S.M.S. constructed the study design, conducted the fieldwork, carried out the statistical analysis, and drafted the manuscript; B.Z. and J.N.P. secured funding, conceived the study, coordinated the study, participated in study design, assisted with analysis, and drafted the manuscript; K.J.M. and M.W.M. provided logistical support and input on the manuscript. M.N. conducted the snow cover analysis and provided input on the manuscript. All authors gave final approval for publication.

Competing interests. We have no competing interests. Funding. Wisconsin Department of Natural Resources and the

8

University of Wisconsin-Madison.

rspb.royalsocietypublishing.org

the same in the next 30 years. Our results highlight the dynamic nature of the field of conservation biology and the need to constantly adjust our approach in this time of rapid global change.

Acknowledgements. We would like to thank D.A. Buehler and L.B. Keith, for providing the historical data and making this research possible. We thank the Wisconsin Department of Natural Resources for funding and providing logistical support with data collection. We thank M. Ozdogan for lending his expertise in developing the historical land cover map, the SILVIS Lab for providing housing and future land cover datasets, and D.J. Mladenoff and J.M. Rhemtulla for historical land cover data. The hard work of field technicians D. Laufenberg and C. Peterson is greatly appreciated. The manuscript benefited from the input of two anonymous reviewers.

1.

Diamond JM. 1975 The island dilemma: lessons of modern biogeographic studies for the design of natural reserves. Biol. Conserv. 7, 129– 146. (doi:10. 1016/0006-3207(75)90052-X) 2. Wiens JA. 1976 Population responses to patchy environments. Annu. Rev. Ecol. Evol. Syst. 7, 81 –120. (doi:10.1146/annurev.es.07.110176. 000501) 3. Fahrig L. 2003 Effects of habitat fragmentation on biodiversity. Annu. Rev. Ecol. Evol. Syst. 34, 487–515. (doi:10.1146/annurev.ecolsys.34.011802. 132419) 4. Bender DJ, Contreras TA, Fahrig L. 1998 Habitat loss and population decline: a meta-analysis of the patch size effect. Ecology 79, 517–533. (doi:10. 1890/0012-9658(1998)079[0517:HLAPDA]2.0.CO;2) 5. Parmesan C. 2006 Ecological and evolutionary responses to recent climate change. Annu. Rev. Ecol. Evol. Syst. 37, 637–669. (doi:10.1146/annurev. ecolsys.37.091305.110100) 6. Mantyka-pringle CS, Martin TG, Rhodes JR. 2012 Interactions between climate and habitat loss effects on biodiversity: a systematic review and meta-analysis. Glob. Chang. Biol. 18, 1239 –1252. (doi:10.1111/j.1365-2486.2011.02593.x) 7. Travis JMJ. 2003 Climate change and habitat destruction: a deadly anthropogenic cocktail. Proc. R. Soc. Lond. B 270, 467–473. (doi:10.1098/ rspb.2002.2246) 8. Pyke CR. 2004 Habitat loss confounds climate change impacts. Front. Ecol. Environ. 2, 178–182. (doi:10. 1890/1540-9295(2004)002[0178:HLCCCI]2.0.CO;2) 9. Warren MS, Hill JK, Thomas JA, Asher J, Fox R, Huntley B. 2001 Rapid responses of British butterflies to opposing forces of climate and habitat change. Nature 414, 65 –69. (doi:10.1038/ 35102054) 10. Reid BN, Peery MZ. 2014 Land use patterns skew sex ratios, decrease genetic diversity and trump the effects of recent climate change in an endangered turtle. Divers. Distrib. 20, 1425 –1437. (doi:10.1111/ ddi.12243) 11. Lemoine N, Bauer H-G, Peintinger M, Bo¨hningGaese K. 2007 Effects of climate and land-use change on species abundance in a Central European bird community. Conserv. Biol. 21, 495–503. (doi:10.1111/j.1523-1739.2006.00633.x)

12. Chen I-C, Hill JK, Ohlemu¨ller R, Roy DB, Thomas CD. 2011 Rapid range shifts of species associated with high levels of climate warming. Science 333, 1024 –1026. (doi:10.1126/science.1206432) 13. Opdam P, Wascher D. 2004 Climate change meets habitat fragmentation: linking landscape and biogeographical scale levels in research and conservation. Biol. Conserv. 117, 285–297. (doi:10. 1016/j.biocon.2003.12.008) 14. Sexton JP, McIntyre PJ, Angert AL, Rice KJ. 2009 Evolution and ecology of species range limits. Annu. Rev. Ecol. Evol. Syst. 40, 415–436. (doi:10.1146/ annurev.ecolsys.110308.120317) 15. Holt RD, Barfield M. 2011 Theoretical perspectives on the statics and dynamics of species’ borders in patchy environments. Am. Nat. 178, S6–25. (doi:10.1086/661784) 16. Thomas CD. 2010 Climate, climate change and range boundaries. Divers. Distrib. 16, 488– 495. (doi:10.1111/j.1472-4642.2010.00642.x) 17. Mair L, Hill JK, Fox R, Botham M, Brereton T, Thomas CD. 2014 Abundance changes and habitat availability drive species’ responses to climate change. Nat. Clim. Chang. 4, 1– 5. (doi:10.1038/ nclimate2086) 18. Anderson BJ, Akc¸akaya HR, Arau´jo MB, Fordham DA, Martinez-Meyer E, Thuiller W, Brook BW. 2009 Dynamics of range margins for metapopulations under climate change. Proc. R. Soc. B 276, 1415 –1420. (doi:10.1098/rspb.2008.1681) 19. Thomas CD, Franco AMA, Hill JK. 2006 Range retractions and extinction in the face of climate warming. Trends Ecol. Evol. 21, 415 –416. (doi:10. 1016/j.tree.2006.05.012) 20. Fox R, Oliver TH, Harrower C, Parsons MS, Thomas CD, Roy DB. 2014 Long-term changes to the frequency of occurrence of British moths are consistent with opposing and synergistic effects of climate and land-use changes. J. Appl. Ecol. 51, 949–957. (doi:10.1111/13652664.12256) 21. Krebs CJ. 2011 Of lemmings and snowshoe hares: the ecology of northern Canada. Proc. R. Soc. B 278, 481 –489. (doi:10.1098/rspb.2010.1992) 22. Thornton DH, Wirsing AJ, Roth JD, Murray DL. 2012 Habitat quality and population density drive occupancy dynamics of snowshoe hare in variegated

23.

24.

25.

26.

27.

28.

29.

30.

31. 32.

33.

landscapes. Ecography 36, 610 –621. (doi:10.1111/j. 1600-0587.2012.07737.x) Lewis CW, Hodges KE, Koehler GM, Mills LS. 2011 Influence of stand and landscape features on snowshoe hare abundance in fragmented forests. J. Mammal. 92, 561 –567. (doi:10.1644/10-MAMMA-095.1) Litvaitis JA, Sherburne JA, Bissonette JA. 1985 Influence of understory characteristics on snowshoe hare habitat use and density. J. Wildl. Manage. 49, 866–873. (doi:10.2307/3801359) Sievert PR, Keith LB. 1985 Survival of hares at a geographic range boundary. J. Wildl. Manage. 49, 854–866. (doi:10.2307/3801358) Keith LB, Bloomer SEM. 1993 Dynamics of a snowshoe hare population in fragmented habitat. Can. J. Zool. 71, 1385– 1392. (doi:10.1139/ z93-191) Mills LS, Zimova M, Oyler J, Running S, Abatzoglou JTM, Lukacs P. 2013 Camouflage mismatch in seasonal coat color due to decreased snow duration. Proc. Natl Acad. Sci. USA 110, 7360–7365. (doi:10. 1073/pnas.1222724110) Zimova M, Mills LS, Lukacs PM, Mitchell MS. 2014 Snowshoe hares display limited phenotypic plasticity to mismatch in seasonal camouflage. Proc. R. Soc. B 281, 20140029. (doi:10.1098/rspb. 2014.0029) Notaro M, Lorenz D, Hoving C, Schummer M. 2014 Twenty-first-century projections of snowfall and winter severity across Central-Eastern North America. J. Clim. 27, 6526–6550. (doi:10.1175/ JCLI-D-13-00520.1) Zimova M, Mills LS, Nowak JJ. 2014 High fitness costs of climate change induced camouflage mismatch in a seasonally colour moulting mammal. Ecol. Lett. 19, 299 –307. (doi:10.1111/ele.12568) Leopold A. 1945 The distribution of Wisconsin hares. Wisconsin Acad. Sci. Arts Lett. 37, 1–14. Buehler DA, Keith LB. 1982 Snowshoe hare distribution and habitat use in Wisconsin. Can. Field Nat. 96, 19 –29. Dobrowski SZ, Abatzoglou J, Swanson AK, Greenberg JA, Mynsberge AR, Holden ZA, Schwartz MK. 2013 The climate velocity of the contiguous United States during the 20th century. Glob. Chang. Biol. 19, 241–251. (doi:10.1111/gcb.12026)

Proc. R. Soc. B 283: 20153104

References

44.

45.

46.

47.

48.

detected imperfectly. Ecology 84, 2200 –2207. (doi:10.1890/02-3090) Fielding AH, Bell JF. 1997 A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 24, 38– 49. (doi:10.1017/S0376892997000088) Fiske IJ, Chandler RB. 2011 Unmarked?: an R package for fitting hierarchical models of wildlife occurrence and abundance. J. Stat. Softw. 43, 1–23. (doi:10.18637/jss.v043.i10) Core Team R. 2015 R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Pokallus JW, Pauli JN. 2015 Population dynamics of a northern-adapted mammal: disentangling the influence of predation and climate change. Ecol. Appl. 25, 1546– 1556. (doi:10.1890/ 14-2214.1) Ims R, Henden J, Killengreen S. 2008 Collapsing population cycles. Trends Ecol. Evol. 23, 79 –86. (doi:10.1016/j.tree.2007.10.010)

9

Proc. R. Soc. B 283: 20153104

39. Hammer RB, Stewart SI, Winkler RL, Radeloff VC, Voss PR. 2004 Characterizing dynamic spatial and temporal residential density patterns from 1940– 1990 across the North Central United States. Landsc. Urban Plan. 69, 183–199. (doi:10.1016/j.land urbplan.2003.08.011) 40. Radeloff VC et al. 2012 Economic-based projections of future land use in the conterminous United States under alternative policy scenarios. Ecol. Appl. 22, 1036–1049. (doi:10.1890/11-0306.1) 41. Peters GP, Andrew RM, Boden T, Canadell JG, Ciais P, Le Que´re´ C, Marland G, Raupach MR, Wilson C. 2012 The challenge to keep global warming below 28C. Nat. Clim. Chang. 3, 4–6. (doi:10.1038/nclimate1783) 42. Tingley MW, Beissinger SR. 2009 Detecting range shifts from historical species occurrences: new perspectives on old data. Trends Ecol. Evol. 24, 625 –633. (doi:10.1016/j.tree.2009.05.009) 43. Mackenzie DI, Nichols JD, Hines JE, Knutson MG, Franklin AB. 2003 Estimating site occupancy, colonization, and local extinction when a species is

rspb.royalsocietypublishing.org

34. Kucharik CJ, Serbin SP, Vavrus S, Hopkins EJ, Motew MM. 2010 Patterns of climate change across Wisconsin from 1950 to 2006. Phys. Geogr. 31, 1– 28. (doi:10.2747/0272-3646.31.1.1) 35. Rhemtulla JM, Mladenoff DJ, Clayton MK. 2012 Legacies of historical land use on regional forest composition and structure in Wisconsin, USA (mid1800s-1930s-2000s). Ecol. Appl 19, 1061 –1078. (doi:10.1890/08-1453.1) 36. Curtis JT, McIntosh RP. 1951 An upland forest continuum in the prairie-forest border region of Wisconsin. Ecology 32, 476 –496. (doi:10.2307/ 1931725) 37. Stevens Jr DL, Olsen AR. 2004 Spatially balanced sampling of natural resources. J. Am. Stat. Assoc. 99, 262–278. (doi:10.1198/016214504000000250) 38. Mackenzie DI, Nichols JD, Lachman GB, Droege S, Andrew J, Langtimm CA. 2002 Estimating site occupancy rates when detection probabilities are less than one. Ecology 83, 2248 –2255. (doi:10. 1890/0012-9658(2002)083[2248:ESORWD]2.0.CO;2)