Northeastern Naturalist Vol. 23, No. 2 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 229 2016 NORTHEASTERN NATURALIST 23(2):229–248 Evidence for Range Contraction of Snowshoe Hare in Pennsylvania Duane R. Diefenbach1,*, Stephen L. Rathbun2,3, Justin K. Vreeland4,5, Deborah Grove6, and William J. Kanapaux4 Abstract - In Pennsylvania, Lepus americanus (Snowshoe Hare) is near the southern limits of its range and at risk of range contraction because of loss of early-successional forest and impacts of climate change. We used hunter-harvest data to investigate changes in the distribution of Snowshoe Hare in Pennsylvania (1983–2011), forest inventory and land-use data to assess changes in amount and distribution of early-successional forest (1988–2011), and occupancy modeling (2004) to identify habitat and climate variables that explain the current distribution of Snowshoe Hare. We determined presence of Snowshoe Hare based on visual sightings, observations of tracks, and DNA analysis of fecal pellets, and used repeated visits to sampling sites and occupancy models to estimate occupancy rates (Ψ). Hunter-harvest data indicated the range of Snowshoe Hare in Pennsylvania contracted towards northwestern and northeastern portions of the state. Based on occupancy modeling, Snowshoe Hare were most likely to occupy early-successional and mixed deciduous–coniferous forest types and areas with colder winter temperatures, which coincided with the distribution of hunter harvests. Among the 4 forest types, we estimated Ψ = 0.52–0.79 and Ψ = 0.10–0.32 where winter temperatures were coldest and warmest, respectively. Total forest loss was less than 1% during 1988–2011, and the loss of early-successional forest in the current and former range of Snowshoe Hares was similar as were mean patch size and a fragmentation metric of early-successional habitat. Thus, changes in forest characteristics did not explain the range contraction we observed. We used climate-model predictions and our occupancy model to predict that average occupancy probability across northern Pennsylvania may decline from 0.27 in 2004 to 0.10–0.18 by 2050–2059, depending on the climate model. The range of Snowshoe Hare in Pennsylvania has contracted to regions of Pennsylvania with the coldest winter temperatures and most persistent snowpack, and based on projected climate change, our results suggest further range contraction of Snowshoe Hare in Pennsylvania. Introduction Lepus americanus Erxleben (Snowshoe Hare), is considered secure and abundant in North America (NatureServe 2012) as a whole, but ranges from secure to extirpated in the mid-Atlantic region of the eastern US. The Snowshoe Hare is ranked as 1US Geological Survey, Pennsylvania Cooperative Fish and Wildlife Research Unit, Pennsylvania State University, University Park, PA 16802. 2Pennsylvania State University, Department of Statistics, University Park, PA 16802. 3Current address - University of Georgia, Health Sciences Campus, B.S. Miller Hall, Athens, GA 30602. 4Pennsylvania Cooperative Fish and Wildlife Research Unit, Pennsylvania State University, University Park, PA 16802. 5Current address - Pennsylvania Game Commission, Huntingdon, PA 16652. 6Pennsylvania State University, Huck Institutes of the Life Sciences, Genomics Core Facility, University Park, PA 16802. *Corresponding author - ddiefenbach@psu.edu. Manuscript Editor: Fred Servello Northeastern Naturalist 230 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 Vol. 23, No. 2 secure in New York, vulnerable and apparently secure in Pennsylvania, vulnerable in West Virginia, critically imperiled in Virginia, possibly extirpated in Maryland, and presumed extirpated in Ohio (NatureServe 2012). The southern-range boundary of the species along the Appalachian Mountains is situated within those states reporting vulnerable or extirpated populations, which raised concern that climate change may be affecting these populations (Mills et al. 2013; Zimova et al. 2014, 2016); however, loss of habitat also could adversely affect their population status. Much of the land in the eastern US is privately owned, land parcels are becoming smaller (Brooks 2003, Litvaitis 2001), and early-successional habitats are likely to become increasingly fragmented (Brooks 2003), which has been identified as a problem for Snowshoe Hare in Pennsylvania (Brown 1984, Scott and Yahner 1989). Climate change is thought to be a concern for Snowshoe Hare because it could lead to less snowfall and fewer days with snow on the ground. Consequently, Snowshoe Hares may experience greater predation risk because their white pelage in winter would not match environmental conditions (Mills et al. 2013; Zimova et al. 2014, 2016). Annual temperature in Pennsylvania is predicted to increase 2–7 °C during the 21st century depending on human choices made regarding development, technological progress in reducing carbon emissions, and the type of climate model used (IPCC 2007a). Predictions regarding changes in precipitation are less certain, but overall precipitation may increase 7%, and whether increased precipitation will occur as snow or rain is unclear (Burakowski et al. 2008, IPCC 2007a, Kunkel et al. 2002). We expect that without compensatory changes in survival in other seasons or increased reproduction, Showshoe Hare populations will decline under current climate change scenarios, and the species’ range will likely contract towards regions with colder winters with more snow. Along with the potential for adverse effects of climate change (e.g., Zimova et al. 2016), recent trends in the composition and structure of eastern US forests may be detrimental to the Snowshoe Hare. Most research on habitat use of Snowshoe Hare, however, has occurred where populations are considered secure (Natureserve 2012). In Maine, Litvaitis et al. (1985), Monthey (1986), and Fuller and Harrison (2013) reported that Snowshoe Hare are associated with high understory densities of hardwood and softwood tree species, although softwood understories provide more visual cover from predators and greater thermal protection. Conroy et al. (1979) reported that Snowshoe Hare in Michigan were unlikely to use habitats >200–400 m from lowlands dominated by Thuja occidentalis L. (Northern White Cedar) and Abies balsamea L. (Balsam Fir). In contrast to more northern populations, Snowshoe Hare in Pennsylvania are generally limited to forested habitats at higher elevations (>450 m) and scrub–shrub-type wetlands (Brown 1984, Doutt et al. 1996, Glazer 1959). Coniferous tree species are considered an important preferred habitat component of northern Snowshoe Hare populations (e.g., Litvaitis et al. 1985); however, Pennsylvania has little of this forest type (McWilliams et al. 2007). Balsam Fir, Picea rubens L. (Red Spruce), and Northern White Cedar are tree species commonly associated with preferred Snowshoe Hare habitat (Conroy et al. 1979, Fuller and Harrison 2013, Litvaitis et al. 1985); however, these species Northeastern Naturalist Vol. 23, No. 2 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 231 are rare and only occur at high elevations or near bogs in Pennsylvania (Rhoads and Block 2007). Brown (1984) and Scott and Yahner (1989) reported that stem densities in regenerating hardwood stands in Pennsylvania (>10,000 stems/ha) provided sufficient browse and protective cover, and Brown (1984) noted that areas with Kalmia latifolia L. (Mountain Laurel) or Tsuga canadensis (L.) Carr. (Eastern Hemlock) were used less than 5–15-year-old clearcuts when all 3 habitat types were available in close proximity. Scott and Yahner (1989) concluded that proximity (less than 0.5 km) of other clearcut stands was important to Snowshoe Hares because these nearby habitats provided alternate sources of food and cover, and that habitat availability may be critical to the viability of Snowshoe Hare populations in Pennsylvania. Consequently, the vulnerability of Snowshoe Hare populations could be related to changes in the amount and dispersion of early-successional habitat. Our objectives were to determine if (1) there were changes in the distribution of Snowshoe Hare in Pennsylvania during 1983–2011 as reflected by the distribution of hunter harvests, (2) occupancy models would indicate that certain habitat and environmental variables plausibly explain the current distribution of Snowshoe Hare in Pennsylvania, and (3) climate variables could explain current distribution, thus enabling us to use models of climate change to predict the future distribution of Snowshoe Hare in Pennsylvania. Field-Site Description The study area included Pennsylvania north of Interstate 80 from Warren, Forest, and Clarion counties in the west to the eastern border of Pennsylvania with New York and New Jersey (39,516 km2). We used a geographic information system (GIS) with vegetative data from the Pennsylvania GAP Analysis Project (http:// www.orser.psu.edu/pagap) to classify forest vegetation as 1 of 4 types used in this study—early-successional forest (3076 km2; 5–40% woody plant foliage, shrubland or forest regeneration), deciduous forest (24,896 km2; ≤30% of tree canopy evergreen), mixed deciduous–coniferous forest (4112 km2; with deciduous trees present and evergreen trees comprising >30% of canopy cover), and coniferous forest (776 km2; ≤30% of tree canopy cover deciduous). We did not expect other vegetation types, including annual or perennial herbaceous vegetation, water, wetlands, or anthropogenic land-use types (6656 km2) to be occupied by Snowshoe Hares and excluded these areas from our study. Pennsylvania is currently in the transition zone between oak–hickory forests to the south and northern hardwood forests to the north (Cuff et al. 1989). In southern Pennsylvania, common tree species include Quercus alba L. (White Oak), Quercus prinus L. (Chestnut Oak), Quercus rubra L. (Red Oak), Liriodendron tulipifera L. (Yellow Poplar), and Acer rubrum L. (Red Maple). Species composition of understories was diverse, except on drier sites with shallow soils where the understory is commonly dominated by Mountain Laurel, Gaylusacia spp. (huckleberries spp.) and Vaccinium spp. (blueberries). The transition between the oak–hickory and northern hardwood forests occurrs at mid-latitudes in Pennsylvania, although the Northeastern Naturalist 232 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 Vol. 23, No. 2 oak–hickory forests extend farther north than northern hardwoods extend south. Common tree species in northern hardwood stands were Prunus serotina Ehrh. (Black Cherry), Red Maple, Betula alleghaniensis Britton (Yellow Birch), Red Oak, and Acer saccharum Marsh. (Sugar Maple). Temperatures in Pennsylvania are warmest in southeastern Pennsylvania (average of less than 20 days with daytime high less than 0 °C) and coolest in northcentral Pennsylvania (average of ≥40 days with daytime high temperatures less than 0 °C) (Cuff et al. 1989). Winter precipitation is greatest in northwestern Pennsylvania, in the south through the Laurel Highlands, and in eastern Pennsylvania along the Delaware River. Average annual snowfall is greatest in northwestern Pennsylvania (152–229 cm) because of lake-effect snowfall from Lake Erie, but is also relatively high in northeastern Pennsylvania (102–127 cm) because of cooler temperatures and greater amounts of winter precipitation. Recent trends (1980–2010) indicated increasing temperatures across Pennsylvania (≤0.75 °C/decade; http://climate.psu.edu/features/changing_climate/). Methods Snowshoe Hare distribution based on harvest We used data from the Pennsylvania Game Commission’s (PGC) annual gametake survey for 1983–2011 to determine the county or wildlife management unit (WMU) where hunters reported harvesting Snowshoe Hare. The PGC survey is a statistically based program for monitoring hunter effort and harvest. The Snowshoe Hare hunting seasons in Pennsylvania were short (less than 3 weeks), and the annual statewide harvest was 510–14,749 hares during 1983–2011. Each year, the PGC randomly surveyed 2% of hunting-license buyers and received ~10,000 responses (response rate of >60%). During the 1983–2002 period, harvest was reported by county, and during 2003–2011 by WMU. We summarized harvest data by 7-year intervals for analysis: 1983–1989, 1990–1996, 1997–2002, and 2005–2011. We excluded the 2003 survey data because it was the first year that harvest data were collected by WMU and the PGC did not conduct a survey in 2004. We assumed that hunter harvest was related to the relative abundance and range distribution of Snowshoe Hare in Pennsylvania, although harvest counts reflect the combination of hunter effort and hare abundance. Therefore, we used the estimated number of Snowshoe Hares harvested by county or WMU per 5000 hunters to adjust for hunter effort and map the distribution and relative size of harvest. The number of Snowshoe Hare hunters changed over time and declined by 1795 hunters/year between 1985–1989 (28,960 in 1983 to 17,568 in 1989) and declined by 170 hunters/year during 1990–2011 (7831 in 1990 to 4039 in 2011). We excluded data when harvests were reported in counties where it was unlikely that Snowshoe Hare existed according to results of mammal surveys conducted by the PGC in the 1940s or were unlikely to exist based on land use. Only 1–2 hunters per 7-y period per county or WMU reported any harvest in counties we classified as unlikely to have Snowshoe Hare. Snowshoe Hare distribution based on field surveys We selected 240 sampling sites and allocated sites equally (n = 60) among the 4 forest types under the constraint that no pair of sampling sites was located ≤10 Northeastern Naturalist Vol. 23, No. 2 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 233 km of one another. We included 6 extra sites per stratum as alternates intended to replace sampling sites when access was not possible. Snowshoe Hares are thought to favor early-successional forest, such that the distribution of Snowshoe Hares in this forest type may yield better information regarding the distribution of Snowshoe Hares in northern Pennsylvania. Therefore, when the primary forest type sampled was coniferous, mixed, or deciduous forest, we attempted to locate and sample an additional site with early-successional forest habitat within 5 km of the original non-early-successional sample site. We sampled sites during January–April 2004 and visited them in clusters in such a way that the date of visitation was not spatially correlated across the study area. This approach minimized travel distance between sites and also minimized any non-random ordering of date of visit to sample sites. At each sampling site, we surveyed a 1000-m transect; we selected the shape and orientation of transects in the field to ensure the entire transect fit within a single patch of the desired forest type. If the forest type at the sampling site differed from that determined by GIS, we identified a new sampling site of the prescribed forest type within 5 km of the randomly selected site. On each transect, we searched for lagomorph sign, which included fecal pellets, tracks, and visual observations of Snowshoe Hare. Investigation of probable tracks detected on transects to confirm species presence was permitted, but wandering off transects was not. To maximize the probability of encountering tracks or fecal pellets, we tried to visit sampling sites 12–48 hours after a snow event; otherwise, we conducted surveys any time regardless of snow cover. To be able to estimate the probability of detecting Snowshoe Hare, we revisited 24 sites where we previously detected hares (visually observing Snowshoe Hare or their tracks in snow). Also, we recorded if no or partial snow cover was present (= 0), or complete snow cover was present (= 1) to be used as a covariate to model detection-probability as a function of snow cover. We collected fecal pellets so that we could employ DNA techniques (Kovach et al. 2003) to distinguish fecal pellets of Snowshoe Hare, Silvilagus floridanus J.A. Allen (Eastern Cottontail), and Silvilagus obscurus Chapman (Appalachian Cottontail). When collecting pellets, technicians wore surgical gloves and stored pellets in Whirl-pak (Nasco, Fort Atkinson, WI) bags. We sometimes placed multiple pellets in a single Whirl-pak, but only pellets collected less than 3 m from each other. We kept pellet samples cool or frozen until they were brought to the PSU nucleic acid facility where they were stored at -20 °C until DNA extraction. Fecal DNA analysis We used a modified procedure to extract DNA with the QIAmp® DNA stool Mini kit (Qiagen, Valencia CA). We used wood applicators to break up the pellets into small pieces in 15-ml centrifuge tubes. We modified the stool-kit procedures by (1) adding 1.8 ml of buffer ASL supplied in the kit in 3 aliquots to the macerated pellet and then spinning down the particulate matter with each aliquot as in the kit protocol, and (2) centrifuging the lysate a second time after the 70 oC incubation Northeastern Naturalist 234 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 Vol. 23, No. 2 period and before the addition of ethanol. Amounts of DNA extracted ranged from 200 ng to 2000 ng per pellet. Primer/probe sets for the 3 lagomorphs were from consensus sequences for the mitochondrial (mt) DNA gene coding the proline tRNA gene and the conserved Dloop. We used Primer Express version 1.0 (Applied Biosystems, Foster City, CA) to select the primers/probe sets. We used sequences from Genbank (http://www. ncbi.nlm.nih.gov/genbank/) as identified by the accession numbers listed with each species below. The forward and reverse primers for Eastern Cottontail were 5'–TTC CCC ATG CAT ATA AGC TAG T–3' and 5'–AAA AGT ATA TGT GGA GTT AGG GTT AAG–3', respectively (Accession No. AF497543). For Appalachian Cottontail, the forward and reverse primers were 5'–TTA ACA AAT TTT TCC ACA ACC CTA TG 3' and 5' CCC ATG TTG GTT ATG GAA TTA TTG TAC–3' (Accession No. AF002244). For Snowshoe Hare, the forward and reverse primers were 5'– CGA AAA CCC TCT TCG TGC TAT G–3' and 5'–TAT GCA TGG GGC AGA ACT TTA–3' (Accession No. AF497544). The probe sequence for Eastern Cottontail was 5'–FAM-CAT TCC TGC TTT ATC GGA CAT AGA CCA–3'-BHQ1, and the shared Appalachian Cottontail and Snowshoe Hare probe was 5'–VIC-AAT TCG GGC ATT ACT GCT TTT CCC CA–3'-TAMRA. We synthesized primers using phosphoramidite chemistry in MerMade 12 (Bioautomation, Plano, TX) at the Penn State University Genomics Core Facility. The Eastern Cottontail probe was synthesized by Biosearch Technologies (Novato, CA), and the Snowshoe Hare and Appalachian Cottontail probe was synthesized by Applied Biosystems. We added DNA (1–50 ng) to a master mix containing 10X Taqman® Buffer (Applied Biosystems), 4-mM MgCl2, 10 nmol each of deoxynucleotidetriphosphates (dNTPs), 10 pmol of primer, 500 pmol of probe, and 2.5 units of Taqman® Gold polymerase. (Applied Biosystems). We assayed samples in duplicate in a 96-well thin-walled PCR plate in an Applied Biosystems 7300 Sequence Detector. The cycling protocol was 2 min at 50 oC, 10 min at 95 oC, and 45 cycles of 15 sec at 95 oC and 75 sec at 60 oC. We multiplexed Eastern Cottontail and Appalachian Cottontail primers/probes simultaneously and assayed Snowshoe Hare samples separately. Climate data We obtained daily minimum temperature, daily snowfall, and daily snow depth from 232 weather stations in Pennsylvania and surrounding states using daily Global Historical Climatology Network data from the National Oceanic and Atmospheric Administration’s (NOAA) National Climatic Data Center for December–March from 1995 to 2005. We selected weather stations that had data for at least 10 of the 11 winter seasons and at least 75% of days in each December–March period. We used these data to create 3 statistics to characterize winter weather conditions: sum of minimum temperature (degree days), number of days with snow on the ground (snow days), and total snowfall (snowfall) for December–March. To account for days with missing data at individual weather stations, we calculated the average daily value for degree days, snowfall, and snow days (and multiplied the average by the number of days in Northeastern Naturalist Vol. 23, No. 2 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 235 the December–March period. We were able to use data from 122 stations for snow days, 144 stations for snowfall, and 98 stations for degree days. We projected all spatial data used in this study in NAD 1983 State Plane Pennsylvania North FIPS 3701. We created shapefiles in ArcGIS 10.1 for degree days, snow days, and snowfall at each weather station’s geographic location. We used the data in these shapefiles to interpolate weather characteristics across the state at a resolution of 30 m × 30 m by kriging with a spherical semivariogram model. Using data from weather stations beyond the Pennsylvania border ensured that the interpolations created climate data for the entire spatial extent of the study area. We used the spatial coordinates for each of our sampling sites to extract interpolated data for the 3 weather statistics at each sampling location. To create spatial data to investigate effects of future climate change on occupancy probability, we used dynamically downscaled simulations of present and future climate over eastern North America developed by the US Geological Survey and researchers at the College of Oceanic and Atmospheric Sciences, Oregon State University, Corvallis, OR. These downscaled simulations used the regional climate model RegCM3 to simulate regional climate conditions from the outputs of 3 general circulation models: GFDL CM 2.0, from NOAA’s Geophysical Fluid Dynamics Laboratory, Princeton, NJ; MPI ECHAM5, from the Max Planck Institute for Meteorology, Hamburg, Germany; and GENMOM, a coupled atmosphere–ocean model (Hostetler et al. 2011). Dynamical downscaling accounts for the effects of terrain on coarser general-circulation models, and we used downscaled simulations of the study area’s climate for 2050–2059 described by Hostetler et al. (2011) based on the IPCC A2 scenario (2007b), which provides a best estimate of a 3.4-ºC increase in global temperature in the period 2090–2099 relative to 1980–1999, and a likely increase ranging from 2.0–5.4 ºC. We created maps of projected degree days using the 3 downscaled simulations, which are calculated for 15-km-wide grids, and resampled them to match the 30 m × 30 m grid used for the forest-type data. Occupancy modeling We used occupancy models (MacKenzie et al. 2006) to estimate the probabilities that Snowshoe Hares were detected (p) and the proportion of sampling sites occupied by Snowshoe Hares (Ψ). We used program MARK (White and Burnham 1999) and Akaike’s information criterion adjusted for sample size (AICc) to select the most parsimonious model. In all models, we used a logit link to model occupancy and detection probability as a linear function of predictive variables. We modeled detection probability as a function of snow cover when the sampling site was visited and assumed that detection probability did not differ among forest types. We developed a model that estimated occupancy by forest type as well as models that included climate variables. We created both additive models and crossclassified models. Additive models indicated that the climate variable had the same relationship (same slope but different intercept terms) across forest types, whereas cross-classified models allowed the climate effect (intercept and slope) to differ for each forest type. The climate variables we investigated were correlated with each Northeastern Naturalist 236 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 Vol. 23, No. 2 other; thus, we did not include more than 1 variable in a given model. (Pearson’s |r| ranging from 0.67 to 0.83, n = 750–1265, P < 0.001). Based on the selected model, we used the coefficients in a GIS to create maps of the probability of occupancy across the range of Snowshoe Hares in Pennsylvania. We used the forest type and climate variables associated with each 30 m × 30 m cell in the GIS to calculate predicted occupancy based on our best model. We created maps of future occupancy by Snowshoe Hares to reflect potential changes as a result of climate change. Because the scale at which we estimated occupancy (1000-m transects) differed from the scale at which we mapped occupancy (30 m × 30 m cells), the maps overestimate occupancy probabilities for each cell, but provide relative measures of occupancy among cells. Changes in amount and distribution of early successional habitat We obtained forest statistics for Pennsylvania for 1987–1988 (Alerich 1989) and 2007–2011 (US Department of Agriculture Forest Service, Forest Inventory and Analysis, Newtown Square, PA, unpubl. data) by county. We used the estimated area of small-diameter forest—defined as areas where >50% of live trees were seedlings (<2.5 cm diameter but >30.5-cm tall) and saplings (2.54–12.5 cm dbh)—as a measure of the area of available habitat for Snowshoe Hares. We placed counties into 2 groups: (1) core counties that represented the 15 counties where Snowshoe Hares were harvested during 1997–2002 (Bradford, Carbon, Clearfield, Elk, Forest, Jefferson, Lackawanna, Luzerne, McKean, Monroe, Pike, Sullivan, Susquehanna, Warren, and Wayne counties), and (2) 19 peripheral counties where Snowshoe Hares were harvested during 1983–1989 but not during 1997–2002. For each group of counties, we estimated the change in area and percent change in area between 1988 and 2011. The percent sampling error for these groups of counties was less than 5%. In addition to changes in the amount of early-successional forest, we investigated if the spatial configuration of this forest-vegetation type differed among 3 regions: the 15 core counties where Snowshoe Hares were most recently harvested, 10 counties in the Laurel Highland and Ridge and Valley region where Snowshoe Hares were no longer harvested (Bedford, Cambria, Columbia, Huntingdon, Indiana, Northampton, Northumberland, Schuylkill, Somerset, and Westmoreland counties), and 4 counties in north-central Pennsylvania where Snowshoe Hares were no longer harvested (Centre, Lycoming, Potter, and Tioga counties). We used a GIS with vegetative data from 2000 (Pennsylvania State University, Pennsylvania Spatial Data Access, http://www.pasda.psu.edu) and the SDMTools package in the statistical programming language R (R Core Team 2012) to calculate mean patch size and fragmentation (proportion of grid cells of early-successional forest that were adjacent to another grid cell of the same forest type) for each county. We hypothesized that if loss of habitat was the reason for range contraction, loss of early-successional forest would be greater and patches of early-successional forest would be smaller and more fragmented in counties where Snowshoe Hares were no longer harvested than in the counties where Snowshow Hares were still observed. Northeastern Naturalist Vol. 23, No. 2 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 237 Results Harvest records indicated a range contraction of Snowshoe Hares towards northeastern and northwestern Pennsylvania during 1983–2011 (Figs. 1, 2). During 1983–1989, Snowshoe Hare harvests occurred across the northern counties of Figure 1. Average harvest of Snowshoe Hare per 5000 hunters in Pennsylvania for (A) 1983–1989 and (B) 1990–1996. Based on data from PGC annual surveys. Northeastern Naturalist 238 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 Vol. 23, No. 2 Pennsylvania and extended from the northeastern part of the state south to Schuylkill County. From the northwest, harvests were reported south to the Maryland border in counties that encompassed the Laurel Highlands. Counties with the largest harvests were in the northeast and in the northwest, south through the Laurel Highlands (Fig. 1A). During 1990–1996, there was evidence of range contraction and limited Figure 2. Average harvest of Snowshoe Hare per 5000 hunters in Pennsylvania for (A) 1997– 2002 by county and (B) 2005–2011 by WMU. Based on data from PGC annual surveys. Northeastern Naturalist Vol. 23, No. 2 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 239 evidence of Snowshoe Hares being harvested throughout the Laurel Highlands (Fig. 1B). During 1997–2002, harvests were restricted to northeastern and some northwestern counties (Fig. 2A). The distribution of harvest did not seem to change for 2005–2011. Harvest records for this period were calculated by WMUs; thus, it was difficult to assess whether coarser-scale data collection was masking results and further range contraction had actually occurred or if there was some range expansion into the north-central region of PA. (Fig. 2B). We sampled 238 sites (43 conifer, 57 deciduous, 53 mixed, and 85 early successional) during January–April 2004 and detected Snowshoe Hares at 69 sites. Not accounting for detection probability, 8 of 57 deciduous sites were occupied (14%), 7 of 43 conifer sites were occupied (16%), 20 of 53 mixed deciduous–conifer sites were occupied (38%), and 34 of 85 early-successional sites were occupied (40%). Presence of Snowshoe Hares was detected via DNA analysis of pellets on 50 of 87 visits, via tracks on 14 visits, and the remainder via some combination of 2 or more types of sign. We visually detected Snowshoe Hares at 2 sites but also detected them at these same sites via tracks or DNA. We found that detection probability was greater when snow cover was present at the time we visited the sampling site, so we modeled detection probability as a function of snow cover in all models. The model with the lowest AICc value included different intercepts and slopes for the occupancy relationship of snowfall with each forest type, but not all parameters were estimable (Table 1). This model indicated that the deciduous forest type, had the lowest occupancy rates and earlysuccessional and mixed deciduous–conifer forest types had greater occupancy rates; no relationship was estimable for the conifer forest type. We selected the model that was ranked second by AICc because all parameters were estimable, the occupancy rates among forest types were similar to the best model, and the covariate that explained occupancy (degree days) was correlated Table 1. Models of presence of Snowshoe Hare in northern Pennsylvania, 2004, based on detection probability (p) as a function of snow cover (0 = no snow, 1 = full or partial snow cover) where occupancy (Ψ) was a function of 4 habitat types, as well as the number of days with snow on the ground (snow days), sum of minimum temperatures (degree days in C°), and total snowfall (snowfall) for the period December–March (1995–2005). K = number of parameters in the model. -2×log(L) = loglikelihood multiplied by -2. Models with a + indicate additive models where the climate variable had different intercepts but the same slope by habitat type. Models wi th a × indicate cross-classified models in which intercept and slope differed for each habitat type. Model Δ AICc AICc weight Model likelihood K -2×log(L) Ψ(habitat × snowfall)A 0.00 0.91 1.00 10 292.7 Ψ(habitat + degree days) 5.68 0.05 0.06 7 304.8 Ψ(habitat) 8.89 0.01 0.01 6 310.2 Ψ(habitat + snow days) 8.91 0.01 0.01 7 308.1 Ψ(habitat + snowfall) 9.22 0.01 0.01 7 308.4 Ψ(habitat × degree days) 10.25 0.01 0.01 10 302.9 Ψ(habitat × snow days) 14.78 less than 0.01 less than 0.01 10 307.4 AThis model was excluded from consideration because not all parameters were estimable and snowfall was correlated with degree days (second-ranked model). Northeastern Naturalist 240 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 Vol. 23, No. 2 with the snowfall covariate in the best model (Pearson’s r = -0.67, n = 917, P < 0.001). The selected model was additive; occupancy declined as degree days increased for all forest types, but the intercept differed among forest types (Table 2, Fig. 3). Estimated probability of detection was 0.71 (95% CI = 0.42–0.89) when snow covered the ground and 0.49 (95% CI = 0.32–0.67) when snow cover was partial or lacking. The estimated probability of occupancy ranged 0.52–0.79 among forest types when the value for degree days was -1085 Cº and 0.10–0.32 when the degree days value was -635 Cº (Fig. 3). Figure 3. Lines depicting predicted probability of occupancy for Snowshoe Hare in Pennsylvania for earlysucces s ional ( s q u a r e s ) , mixed deciduous- conif er ( t r i a n g l e s ) , conifer (diamonds), and d e c i d u o u s (x’s) forest as a function of the sum of December– March minimum temperatures (degree days in Cº), 2004. Table 2. Coefficient estimates to calculate the probability of detection (p) as a linear function of snow cover (1 = full; 0 = partial or no snow) and probability of occupancy (Ψ) as a linear function of habitat type (conifer, early successional, mixed conifer–deciduous; deciduous is the reference-habitat type) and degree days (sum of minimum temperatures December–March, Cº ) using a logit link function. Parameter Coefficient Estimate SE 85% CI Odds ratio (85% CI) p Intercept -0.022 0.371 -0.557–0.513 p Snow cover 0.916 0.603 0.048–1.785 2.50 (1.05–5.96) Ψ Intercept -5.127 2.014 -8.026– -2.227 Ψ Conifer 0.275 0.646 -0.655–1.205 1.32 (0.52–3.34) Ψ Early successional 1.481 0.602 0.615–2.348 4.40 (1.85–10.47) Ψ Mixed 1.319 0.600 0.455–2.184 3.74 (1.58–8.88) Ψ Degree days -0.005 0.002 -0.008– -0.001 1.58 (1.12–2.23)A AOdds of a habitat more likely to be occupied for every 100 degree-day (C°) decrease in December– March daily minimum temperature. Northeastern Naturalist Vol. 23, No. 2 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 241 Mapping the selected model of occupancy across the study area indicated that forests in the Poconos in northeastern Pennsylvania and Elk, Forest, McKean, Potter, and Warren counties in the northwest, which were associated with the coolest temperatures and greatest snowfall and snow days (Fig. 4B, C, D), had the greatest predicted occupancy (Fig. 5A). Using the downscaled climate simulations, mean occupancy probabilities within the study area indicated a decline from 0.27 for 2004 to a range of 0.10–0.18 in 2050–2059 (Table 3; Fig. 5B, C, D). The amount of forested area in Pennsylvania declined by less than 1% between 1988 and 2011. However, statewide, the area of small-diameter stands decreased by 33%. The area of small-diameter stands declined by 28.5% within the core counties (counties where Snowshoe Hares were harvested during 1997–2002; Fig. 2A), and declined 31.6% in peripheral counties (counties where Snowshoe Hares were harvested during 1983–1989 but not during 1997–2002; Fig. 1A). Mean patch size and fragmentation did not differ among core counties, areas in the Laurel Highland and Ridge and Valley Region where Snowshoe Hares used to be harvested, or areas of north-central Pennsylvania where they used to be harvested (Fig. 6). Discussion Our evaluation of hunter-harvest data indicated that range contraction for Snowshoe Hare has occurred in Pennsylvania. Furthermore, the most recent harvests of Snowshoe Hare (Fig. 2) occurred in 2 regions of Pennsylvania with the coldest temperatures and most days with snow on the ground (Fig. 4D). These findings were supported by our field survey of the occurrence of Snowshoe Hare, which indicated that probability of occupancy was highest in these same 2 regions of Pennsylvania (Fig. 5A). This range contraction was similar to the one in Wisonsin recently reported by Sultair et al. (2016) . If climate conditions are an important driver of the distribution of Snowshoe Hare in Pennsylvania, then predictions from climate change models indicate that further range contraction is likely (Fig. 5B, C, D). Based on work by the Intergovernmental Panel on Climate Change (IPCC 2007b), the mean projected temperature increase in Pennsylvania is 3.5 °C during this century (>400 degree days increase for December–March). Moreover, Table 3. A comparison of minimum, maximum, and mean occupancy probability for Snowshoe Hare in all habitat types in the study area in Pennsylvania based on the average sum of minimum temperatures for December–March, from 1995 to 2005 (degree days) and projected degree days based on 3 dynamically downscaled climate-change projections for the study area in 2050–2059. MPI ECHAM5 = general circulation model from the Max Planck Institute for Meteorology, GENMOM = a coupled atmosphere– ocean model (Hostetler et al. 2011), and GFDL CM 2.0 = a general circulation model from the Geophysical Fluid Dynamics Laboratory at the US National Oceanic and Atmospheric Administration. Probability of occupancy Scenario Degree days (C°) Maximum Minimum Mean SD Baseline (1995–2005) -854.6 0.82 0.08 0.27 0.13 MPI ECHAM5 -736.5 0.63 0.05 0.18 0.09 GENMOM -697.2 0.55 0.04 0.15 0.08 GFDL CM 2.0 -570.9 0.38 0.03 0.10 0.05 Northeastern Naturalist 242 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 Vol. 23, No. 2 Figure 4. Maps of variables used to estimate occ u r r e n c e o f Snowshoe Hare. (A) Distribution of 4 forest types. (B) Interpolated values for degree days. (C) Interpolated values for annual snowfall. (D) Interpolated values for number of days with snow on the ground. Degree days (minimum temperature, C°), snowfall, and snow days (days with snow on the ground) calculated as mean of daily values for December– March multiplied by number of days December– March and averaged for 1995– 2005. Northeastern Naturalist Vol. 23, No. 2 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 243 Figure 5. Probability of occupancy of Snowshoe Hare across northern Pennsylvania: (A) Estimated probability based on sum of December–March minimum temperatures (degree days in Cº), 2004. Predicted probability if sum of minimum temperatures (degree days) for December–March, 2050–2059 increased as projected by dynamically downscaled climate models from (B) GENMOM, (C) MPI ECHAM5, and(D) GFDL CM 2.0. Areas in white were non-forested habitat. Northeastern Naturalist 244 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 Vol. 23, No. 2 the IPCC (2007a) predicts that minimum winter temperatures may exhibit greater increases or greater variability than the average temperature. Consequently, predictive climate models may be conservative, and the potential decline in occupancy of Snowshoe Hares could be greater than we have modeled, with a mean probability of occupancy of 0.10–0.18 in northern Pennsylvania by 2050–2059 (Table 3; Fig. 5B, C, D). Similarly, Sultaire et al. (2016) suggested that climate change will continue to be the major driver of range contraction of Snowshoe Hares in Wisconsin. Climate change has already reduced winter duration in North America (Magnuson et al. 2000). Although models of climate change that include precipitation Figure 6. Mean patch-size (circles) and fragmentation (triangles) of early-successional forest (with 95% confidence intervals) in Pennsylvania, 2000. Patch size measured in ha and fragmentation represented as the proportion of 30 m × 30 m grid cells of early-successional forest adjacent to the same forest type. Statistics were calculated for 15 counties where Snowshoe Hares were harvested during 1997–2002 (Core), 10 counties in the Ridge and Valley and Laurel Highland regions where Snowshoe Hares were not harvested during 1997–2002 (Ridge and Valley), and 4 north-central counties (North-Central) where Snowshoe Hares were not harvested during 1997–2002 in Pennsylvania. Northeastern Naturalist Vol. 23, No. 2 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 245 and snowfall are not as well developed as temperature models, weather data for our study area in Pennsylvania indicated that the number of days with snow on the ground (snow days) was negatively correlated with degree days (r = -0.83, n = 750, P < 0.001). Therefore, predicted milder and shorter winters likely will provide fewer days of snow on the ground as winter temperatures warm. Only 2 harvests of Snowshoe Hare in Pennsylvania have been reported in WMUs south of Interstate 80 since 2006. Although the presence of Snowshoe Hare has been detected in Huntingdon and Westmoreland counties (E. Boyd, Pennsylvania Game Commission, Harrisburg, PA, pers. observ.), our results suggest that connectivity of Snowshoe Hare populations in the higher elevations of Maryland and West Virginia with northern Pennsylvania may be limited. The mechanism by which climate change is thought to adversely affect Snowshoe Hare is through mismatched pelage coloration and environmental conditions (Mills et al. 2013; Zimova et al. 2014, 2016). Shorter winters could increase the length of time in which the white camouflage coloration of Snowshoe Hare is mismatched with lack of snow, and this mismatch could put Snowshoe Hares at greater risk of predation. However, some subspecies of Snowshoe Hare in the southern range have evolved to forego pelage coloration change (Dalquest 1942, Nagorsen 1983), and recent research in Pennsylvania has found that some Snowshoe Hare in northeastern Pennsylvania do not develop a white pelage in winter or become only partially white (Gigliotti 2016). Several potential sources of error could have affected our results, including errors in the forest-type map and the kriging of weather data. Forest-type error did affect the map of occupancy based on forest type but did not affect occupancy modeling (Tables 1, 2) because forest type was determined in the field and not from the map. The kriging we conducted to interpolate weather data could have large and spatially autocorrelated interpolation errors among local grid-cells but less bias at broader scales (Holmes et al. 2000). Interpolation errors with a strong spatial pattern could result in local prediction errors, but if the error structure were consistent between the predicted values and the modeling data, the model should reflect the errors, and the predictions should be consistent with the data (Barry and Elith 2006). Measurement error in a covariate when using logistic regression usually attenuates the estimates of occupancy (overestimates low probabilities and underestimates high probabilities; Stefanski and Carroll 1985). Our sample locations may not have been completely representative due to constraints imposed by logistical challenges and random sampling. We had limited a priori information available on the expected occupancy of habitats except that we expected the early-successional forest type to have the greatest occupancy rates. Also, we were limited by the accuracy of the forest-type map and the relative rarity of early-successional habitat, which led us to direct field technicians to sample additional early-successional sites near the conifer, mixed deciduous–conifer, and deciduous forest-type sampling sites. However, this approach would not have affected our estimates of occupancy by forest type unless these additional sampling sites were not representative. Northeastern Naturalist 246 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 Vol. 23, No. 2 Despite the relationship we discovered between climate and distribution of Snowshoe Hare in Pennsylvania, the amount and distribution of suitable habitat is essential to the viability of the species. Snowshoe Hare are associated with early-successional forest types, and from 1978 to 2002 the total acreage of Pennsylvania forestland remained stable but the proportion in early-successional stages declined from 21% to 12% because succession of forest vegetation outpaced timber-harvest rates (Mc- Williams et al. 2007). Wildlife managers may be able to mitigate effects of climate change by employing strategies to increase resilience of ecosystems and species (Mawdsley et al. 2009), which for the Snowshoe Hare would include increasing the amount of early-successional forest and improving habitat connectivity. We found the greatest occupancy rates in early-successional forest, and other researchers have noted the association of Snowshoe Hares in Pennsylvania with young forest age-classes (Brown 1984, Scott and Yahner 1989), similar to populations studied in Michigan and Maine (Conroy et al. 1979, Fuller and Harrison 2013, Litvaitis et al. 1985). Consequently, the loss of early-successional habitat in Pennsylvania likely impacts the populations of the species. Despite the loss of habitat, however, we found no evidence that loss of early-successional habitat explained the pattern of range contraction we observed. Loss of habitat in Pennsylvania was similar in areas where Snowshoe Hare were no longer harvested compared to where harvests still occurred. Moreover, patch size and fragmentation of early-successional habitat was similar, with no evidence for differences among counties in the Ridge and Valley and Laurel Highland regions (southern Pennsylvania) or north-central Pennsylvania, where Snowshoe Hares have not been recently harvested compared to those counties where harvests still occurred (Fig. 7). We believe climate change will be a challenge for maintaining the viability of Snowshoe Hare in Pennsylvania. Although the species’ conservation status is currently defined as vulnerable (S3) and apparently secure (S4) in Pennsylvania (Natureserve 2012), our data suggest the status of Snowshoe Hare is not stable and becoming less secure. Better understanding of the effect of climate change on this species, and how that factor may interact with loss of early-successional forest, will be important to developing effective conservation strategies for the species. Acknowledgments This project could not have been completed without the outstanding field efforts of B. Carlson, E. Rogan, and M. Schrecengost. We thank K. Sitch and R. Myers for their assistance with the genetics analyses. We thank W.H. McWilliams and the USDA Forest Service’s Forest Inventory and Analysis for access to forest data for Pennsylvania. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US government. Literature Cited Alerich, C.L. 1993. Forest Statistics for Pennsylvania: 1978 and 1989. Resource Bulletin NE-126. US Department of Agriculture Forest Service, Northeastern Forest Experiment Station, Radnor, PA. 244 pp. Barry, S., and J. Elith. 2006. Error and uncertainty in habitat models. Journal of Applied Ecology 43:413–423. Northeastern Naturalist Vol. 23, No. 2 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 247 Brooks, R.T. 2003. Abundance, distribution, trends, and ownership patterns of earlysuccessional forests in the northeastern United States. Forest Ecology and Management 185:65–74. Brown, D.F. 1984. Snowshoe Hare populations, habitat, and management in northern hardwood forest regeneration areas. M.Sc. Thesis. The Pennsylvania State University, University Park, PA. 100 pp. Burakowski, E.A., C.P. Wake, B. Braswell, and D.P. Brown. 2008. Trends in wintertime climate in the northeastern United States: 1965–2005. Journal of Geophysical Research 113:D20114. Available online at doi:10.1029/2008JD009870. Last accessed 18 May 2016. Conroy, M.J., L.W. Gysel, and G.R. Dudderar. 1979. Habitat components of clear-cut areas for Snowshoe Hares in Michigan. Journal of Wildlife Management 43:680-690. Cuff, D.J., W.J. Young, E.K. Muller, W. Zelinsky, and R.F. Abler (Eds.). 1989. The Atlas of Pennsylvania. Temple University Press, Philadelphia, PA. 288 pp. Dalquest, W.W. 1942. Geographic variation in northwestern Snowshoe Hares. Journal of Mammalogy 23:166–183. Doutt, J.K., C.A. Heppenstall, and J.E. Guilday. 1966. Mammals of Pennsylvania. Pennsylvania Game Commission, Harrisburg, PA. 283 pp. Fuller, A.K., and D.J. Harrison. 2013. Modeling the influence of forest structure on microsite- habitat use by Snowshoe Hares. International Journal of Forestry Research 2013. Available online at http://dx.doi.org/10.1155/2013/892327. Last accessed 18 May 2016. Gigliotti, L.C. 2016. Ecology, habitat use, and winter thermal dynamics of Snowshoe Hares in Pennsylvania. M.Sc. Thesis. The Pennsylvania State University, University Park, PA. Glazer, R.B. 1959. An evaluation of a Snowshoe Hare restocking program in Centre County, Pennsylvania. M.Sc. Thesis. The Pennsylvania State University, University Park, PA. 137 pp. Holmes, K.W., O.A. Chadwick, and P.C. Kyriakidis. 2000. Error in a USGS 30-m digital elevation model and its impact on terrain modeling. Journal of Hydrology 233:154–173. Hostetler, S.W., J.R. Alder, and A.M. Allan. 2011. Dynamically downscaled climate simulations over North America: Methods, evaluation, and supporting documentation for users. US Geological Survey Open-file Report 2011-1238. 64 pp. Available online at http:// pubs.usgs.gov/of/2011/1238/. Last accessed 18 May 2016. Intergovernmental Panel on Climate Change (IPCC). 2007a. Climate Change 2007–Working Group II: Impacts, Adaptation, and Vulnerability. Cambridge University Press, Cambridge, UK. 987 pp. IPCC. 2007b. Climate Change 2007: The Physical Science Basis: Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK. 996 pp. Kovach, A.I., M.K. Litvaitis, and J.A. Livaitis. 2003. Evaluation of fecal mtDNA analysis as a method to determine the geographic distribution of a rare lagomorph. Wildlife Society Bulletin 31:1061–1065. Kunkel, K.E., N.E. Westcott, and D.A.R. Krtistovich. 2002. Assessment of potential effects of climate change on heavy lake-effect snowstorms near Lake Erie. Journal of Great Lakes Research 28:521–536. Litvaitis, J.A. 2001. Importance of early-successional habitats to mammals in eastern forests. Wildlife Society Bulletin 29:466–473. Litvaitis, J.A., J.A. Sherburne, and J.A. Bissonette. 1985. Influence of understory characteristics on Snowshoe Hare habitat use and density. Journal of Wildlife Management 49:866–873. Northeastern Naturalist 248 D.R. Diefenbach, S.L. Rathbun, J.K. Vreeland, D. Grove, and W.J. Kanapaux 2016 Vol. 23, No. 2 MacKenzie, D.I., J.D. Nichols, J.A. Royle, K.H. Pollock, L.L. Bailey, and J.E. Hines. 2006. Occupancy Estimation and Modeling: Inferring Patterns and Dynamics of Species Occurrence. Academic Press, San Diego, CA. 324 pp. Magnuson, J.J., D.M. Robertson, B.J. Benson, R.H. Wynne, D.M. Livingstone, T. Arai, R.A. Assel, R.G. Barry, V. Card, E. Kuusisto, N.G. Granin, T.D. Prowse, K.M. Stewart, and V.S. Vuglinski. 2000. Historical trends in lake and river ice-cover in the northern hemisphere. Science 289:1743–1746. Mawdsley, J.R., R. O’Malley, and D.S. Ojima. 2009. A review of climate-change adaptation strategies for wildlife management and biodiversity conservation. Conservation Biology 23:1080–1089. McWilliams, W.H., S.P. Cassell, C.L. Alerich, B.J. Butler, M.L. Hoppus, S.B. Horsley, T.W. Lister, R.S. Morin, C.H. Perry, J.A. Westfall, E.H. Wharton, and C.W. Woodall. 2007. Pennsylvania’s Forest 2004. Resource Bulletin NRS-20. US Department of Agriculture Forest Service, Northern Research Station, Newtown Square, PA. 86 pp. Mills, L.S., M. Zimova, J. Oyler, S. Running, J.T. Abatzoglou, and P.M. Lukacs. 2013. Camouflage mismatch in seasonal coat color due to decreased snow duration. Proceedings of the National Academy of Sciences 110:7360–7365. Monthey, R.W. 1986. Responses of Snowshoe Hares, Lepus americanus, to timber harvesting in northern Maine. Canadian Field-Naturalist 100: 568–570. Nagorsen, D.W. 1983. Winter pelage color in Snowshoe Hares (Lepus americanus) from the Pacific Northwest. Canadian Journal of Zoology 61:2313–2318. NatureServe. 2012. NatureServe Explorer: An online encyclopedia of life [web application]. Version 7.1. Available online at http://www.natureserve.org/explorer. Accessed 27 February 2016. R Core Team. 2012. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Available online at http://www.R-project. org. Last accessed 18 May 2016. Rhoads, A., and T. Block. 2007. The Plants of Pennsylvania. 2nd Edition. University of Pennsylvania Press, Philadelphia, PA. 1056 pp. Scott, D.P., and R.H. Yahner. 1989. Winter habitat and browse use by Snowshoe Hares, Lepus americanus, in a marginal habitat in Pennsylvania. Canadian Field-Naturalist 103:560–563. Stefanski, L.A., and R.J. Carroll. 1985. Covariate measurement error in logistic regression. Annals of Statistics 13:1335–1351. Sultaire, S.M., J.N. Pauli, K.J. Martin, M.W. Meyer, M. Notaro, and B. Zuckerberg. 2016. Climate change surpasses land-use change in the contracting range boundary of a winter- adapted mammal. Proceedings of the Royal Society B 283:20153104. White, G.C., and K.P. Burnham. 1999. Program MARK: Survival estimation from populations of marked animals. Bird Study 46 (Supplement):120–138. Zimova, M., L.S. Mills, P.M. Lukacs, and M.S. Mitchell. 2014. Snowshoe Hares display limited phenotypic plasticity to mismatch in seasonal camouflage. Proceedings of the Royal Society B. Available online at doi:10.1098/rspb.2014.0029. Last accessed 18 May 2016. Zimova, M., L.S. Mills, and J.J. Novak. 2016. High fitness-costs of climate change-induced camouflage mismatch. Ecology Letters 19(3):299–307.