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.
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 . In many regions characterized by high climate change velocity, the ecological impacts of habitat fragmentation appear to be exacerbated , and theory supports a potent synergy between these threats . 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 . 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 . 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 , and habitat is often patchy . 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 . Leading edge range boundaries are the result of population processes, such as colonization, that can be inhibited by 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 . 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 , the phenological mismatch between seasonal coat colour change and background snow conditions is becoming more pronounced , and the mismatch has strong selective costs on survival . 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 . Comparisons of the range boundary between the mid- and late-twentieth century revealed a dynamic range boundary without a clear directional shift . 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 , and the most pronounced changes have occurred during winter months . During this time period, land cover in the region has remained fairly static and the conversion of forest has been minimal . 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 two-phase range shift—first driven by rapid land-use change and 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
(a) Study area and site selection
The southern distribution of snowshoe hares in the US state of Wisconsin was first mapped using expert opinion in 1944 , and later more comprehensively via snow-tracking surveys in 1979–1980 . 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 . 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 km² blocks and at least one location was surveyed within each block . Historic locations were provided at a 2.56 km² 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 , and within each strata one section was randomly chosen from the five most forested 2.56 km² 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 . 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 . 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) . 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 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  and approximates the maximum dispersal distance documented for hares in our study area . In addition, we calculated the number of houses within the 5 km buffer for 1980 (Houses1), 2010 (Houses2), and housing change (HouseΔ) using an existing housing density database  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 , 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 . 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.
(d) Climate data
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.1° (≈ 13 × 18 km) resolution. Snow cover for each 0.1° 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 (OnsetΔ,DurΔ, TermΔ). 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 . 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 .
(e) Occupancy modelling
Comparing current and historic survey data to infer range shifts is inherently problematic . One of the challenges in using 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 .
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  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 , 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 (ψ1) and extinction probability (ɛ). We did not seek to parametrize colonization dynamics (γ) 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 ψ1. 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 ψ 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 ψ1 and ɛ across Wisconsin at the resolution of our site selection (2.56 km2). Occupancy in the modern era (ψ2) was derived using values of ψ1, ɛ, and γ and modelled at the same resolution. We evaluated our distribution model for ψ2 by calculating the Area-Under-the-Curve (AUC)  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 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  in the R statistical platform .
(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.
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: −32.9, max: −6.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 time periods to be high (0.73, 95% CI: 0.63–0.81). The multi-season model with the strongest support (wi = 0.89), given the data, included the effect of Snowdur1, Forest1, and Houses1 on ψ1 (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 (β = 1.55, standard error (s.e.) = 0.44), longer duration of snow cover (β = 1.26, s.e. = 0.36), and lower housing density (β = −0.63, s.e. = 0.37) (electronic supplementary material, table S4). For every 19% (s.d.) increase in forest cover ψ1increased fourfold, while every 3.25 day (s.d.) increase in snow cover resulted in a threefold increase in ψ1. 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 (β = −0.78, s.e. = 0.31) and shorter snow cover duration (β = −1.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 (ΔAIC < 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 DurΔ as an extinction parameter, as opposed to Snowdur2 (electronic supplementary material, table S5). As predicted, the coefficient estimate for the effect of landscape-scale forest disturbance was negative (β = −0.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.
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 ψ2 included the variables Snowdur2 and Forest2 with a larger coefficient estimate for snow cover duration than forest cover (table 1), which is consistent with ɛ parameters in the top-ranked multi-season model. Comparing the ratio between ψ coefficients across time periods demonstrated a larger snow to forest effect on ψ2 (1.13:0.85) than ψ1 (1.26:1.55).
The AUC value for the spatial predictions of ψ2 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 for the current survey period (electronic supplementary material, figure S4) although the performance of the second-ranked 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).
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 ; 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 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 . 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 mid-twenty-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 . We do not know the distribution of young forest historically, but patterns of snowshoe hare extinction do not correlate with landscape-scale 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 . 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 . 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 . 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 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.
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.
We have no competing interests.
Wisconsin Department of Natural Resources and the University of Wisconsin-Madison.
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.
- Received December 28, 2015.
- Accepted March 1, 2016.
- © 2016 The Author(s)
Published by the Royal Society. All rights reserved.