Li. I «I. z. 5 kmwbmguh. hr 1 mummy; "‘34,. . any.» ..... J; » Pol .sz We. 5. :. {1... D t V1 134' I. t s A J. .413. . Whammy} 93W?“ ‘O a t .k “.129”?! ‘90 I.\ n55! 33“. .1. 3"..- \ 3..)V\\. .. . .11) a 35.3.. . .1: fifigfw tllnl FF.;"‘) :1 r3»; THESIS 1 EGI‘I.‘\|"H « ‘ ‘ ‘ ' {junta :4: lLblu’. Jive .II V‘LQLU I U. WrarqH DID . allIvVuALy This is to certify that the thesis entitled MODELING HABITAT SUITABILITY AND POPULATION DEMOGRAPHICS OF THE EASTERN MASSASAUGA RATTLESNAKE IN MANAGED LANDS IN SOUTHWESTERN MICHIGAN presented by Robyn Leah Bailey has been accepted towards fulfillment of the requirements for the MS. degree in Fisheries and Wildlife Mum“ km igajor Prhfessor’s Signature I 5 APT \ I LO! 0 Date MSU is an Affirmative Action/Equal Opportunity Employer PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE R0 "7; osoi‘ié 5f08 KzlProlechresICIRC/DaleDue.Indd MODELING HABITAT SUITABILITY AND POPULATION DEMOGRAPHICS OF THE EASTERN MASSASAUGA RATTLESNAKE IN MANAGED LANDS IN SOUTHWESTERN MICHIGAN By Robyn Leah Bailey A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Fisheries and Wildlife 2010 ABSTRACT MODELING HABITAT SUITABILITY AND POPULATION DEMOGRAPHICS OF THE EASTERN MASSASAUGA RATTLESNAKE IN MANAGED LANDS IN SOUTHWESTERN MICHIGAN By Robyn Leah Bailey The eastern massasauga rattlesnake (Sistrurus catenatus catenatus; hereafter, massasauga) is a candidate for federal listing as threatened or endangered, and a species of Special Concern in Michigan. Massasauga populations in Michigan are poorly understood, even though Michigan is considered the stronghold for the species. My objectives were to: 1) estimate adult survival for the activity season, 2) model population viability for the next 50 years, 3) quantify resource selection in a managed landscape, and 4) evaluate a habitat suitability index (HSI) model for massasaugas in southwestern Michigan. I captured (May—Aug), radiomarked, and monitored 27 adult massasaugas in 2008 and 2009 and pooled data for analyses; data on snakes from previous work in this area were included as necessary. I relocated snakes throughout the active season and estimated multiple demographic rates. I quantified the spatial and structural extent of habitat management practices at the study site by collecting stand- level vegetation data for used and random sites. My results indicate that massasaugas have high survival and will likely persist at the study site for 50 years. A habitat selection model ranked early and mid-successional wetlands and uplands as most- preferred, and a revised HSI model is presented. Biologists can use these results collectively to manage habitat for massasaugas in southwestern Michigan and perhaps elsewhere within the species’ range. For Mom. iii ACKNOWLEDGMENTS Funding and support for this project was provided by the Michigan Department of Natural Resources and Environment (MDNRE), Wildlife Division; Michigan State University (MSU); Potter Park Zoo (PPZ); and Pierce Cedar Creek Institute for Ecological Education (PCCI). I am very grateful to Dr. Tara Harrison and Deb Paperd (PPZ) for generously contributing their time, facilities, and expertise to the project. They served as invaluable reservoirs of knowledge for this project. I thank the PPZ docents for their generous gift. I am grateful to John Mitchell (Abaxis, Inc.) for donating supplies for diagnostic blood testing and to Chris Hencken for putting them to use. Michelle Skedgell and Matthew Dysktra (PCCI) deserve special thanks for their assistance and for their enthusiasm for massasauga conservation. I appreciate the cooperation of Mr. and Mrs. Walt Soya, Mr. and Mrs. Gene Willison, and Mr. and Mrs. Steve Wiersum for allowing access to their properties. Thank you Jon Pruitt for your work on the parasite investigation. I would further like to thank my committee members for their instrumental roles in this research as well as my professional development. First, I owe my major advisor, Dr. Rique Campa, much gratitude for his support and encouragement. Thank you for your guidance, humor, and patience. I owe Dr. Kelly Millenbah sincere thanks for introducing me to GIS, for always giving thoughtful insights, and for many helpfiIl remarks on study design. To Dr. Dennis Propst, thank you for lending your expertise in the translation and communication of research for the benefit of massasaugas. Your particular questions and approaches are greatly appreciated. iv I am also incredibly grateful to several individuals without whom this work would not have been possible. Thank you Kristin Bissell for your willingness to instruct and advise me in the ways of rattlesnakes, and for the use of your data for chapters 2 and 3. Your continued excitement for massasaugas is inspirational. My two field assistants, Wesley Anderson and Stephanie Zirnmer, were essential to the project. Thank you both for your hard work and willingness to roll up your sleeves. I appreciate the helpful comments of Dan Linden and Jason Karl in GIS support. Thanks to Christopher Hoving and Jen Howell for coming to bat when asked and for helping me sort things out. I am also very thankful for the handfiil of volunteers who took time out to help survey for massasaugas during crunch time. I am grateful for the statistical support from Wei Wang of the MSU Statistical Counseling Service and that of Dr. Scott Winterstein of this department. Bret Muter provided helpful reviews over the years. Finally, I must thank my family for their constant confidence, love, and faith in me. Mom and Dad, thank you for supporting me and this research in a myriad of ways. Thanks especially for helping me believe that I can do anything I set my mind on. Thanks Nanny for always saying you knew I could do it. Marie Clark, thanks for keeping me grounded and for always listening to me. And finally, a very special thank you to my best friend and partner, Paul Paradine, who endured the chaos and the many long drives just so he could stand by my side and cheer me on. I honestly could not have accomplished this without you. TABLE OF CONTENTS LIST OF TABLES ......................................................................................................... viii LIST OF FIGURES ........................................................................................................... x GENERAL INTRODUCTION .................................................................... 1 CHAPTER 1: Survival and cause-specific mortality of eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Michigan. ABSTRACT .......................................................................................................... 6 INTRODUCTION .......................................................................... 7 STUDY AREA ............................................................................. 9 METHODS ................................................................................ 10 Capture and Telemetry .......................................................... 10 Survival Estimates ............................................................... 11 RESULTS .................................................................................. 12 DISCUSSION ............................................................................. 16 ACKNOWLEDGMENTS ............................................................... l 8 CHAPTER 2: Population viability of the eastern massasauga rattlesnake (Sistrurus catenatus catenatus) in Michigan. ABSTRACT ................................................................. ' ............... 19 INTRODUCTION ........................................................................ 20 STUDY AREA ............................................................................ 21 METHODS ................................................................................ 22 Capture and Telemetry .......................................................... 22 Demographic Parameters ....................................................... 23 PVA and Sensitivity Analysis ................................................. 28 RESULTS .................................................................................. 33 Demographic Summary ......................................................... 33 Population Viability Model and Sensitivity Analysis ...................... 34 DISCUSSION ............................................................................. 38 ACKNOWLEDGMENTS ............................................................... 43 CHAPTER 3: Linking resource selection to habitat management for the eastern massasauga rattlesnake. ABSTRACT ................................................................................ 44 INTRODUCTION ........................................................................ 45 STUDY AREA ............................................................................ 47 METHODS ................................................................................ 49 Capture and Radiotelemetry ..................................................... 49 Resource Selection ............................................................... 51 Habitat Management ............................................................. 54 RESULTS .................................................................................. 57 Resource Selection .............................................................. 59 Habitat Management ............................................................ 59 DISCUSSION ............................................................................. 64 Management Implications ...................................................... 7O ACKNOWLEDGMENTS ............................................................... 70 CHAPTER 4: Testing a habitat suitability index model for the eastern massasauga rattlesnake. ABSTRACT ................................................................................ 72 INTRODUCTION ........................................................................ 73 STUDY AREA ............................................................................ 78 METHODS ................................................................................ 79 Capture and Radiotelemetry .................................................... 79 Vegetation Sampling ............................................................ 80 Model Evaluation ................................................................. 82 RESULTS .................................................................................. 89 DISCUSSION ............................................................................. 98 ManagementImplications.................... ................................ 105 ACKNOWLEDGMENTS ............................................................. 105 GENERAL CONCLUSIONS ................................................................... 106 APPENDICES APPENDD( A. Common plants associated with stands frequently used by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan ........................................ 109 APPENDIX B. Deterministic population growth rates for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan ........ 111 APPENDIX C. The snake community occurring in cover types used by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) ........................... 112 APPENDD( D. Metadata for radiolocations of 26 eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) tracked in 2008 and 2009 in Barry County, Michigan ........................................................................ 1 13 LITERATURE CITED ............................................................................................ 1 16 vii LIST OF TABLES Table 1.1. Periodic survival probabilities for radiotracked adult eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in southwestern lower Michigan estimated from the 168-day activity seasons (mid-May—Oct) of 2008 and 2009 .................... 14 Table 1.2. A compilation of range-wide survival estimates for adult eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) during the active season known from the literature and this study, listed in order of decreasing survival estimates. Data from this study were collected from May—October of 2008 and 2009 in southwestern lower Michigan ........................................................................................... 15 Table 2.1. Baseline parameters used for modeling eastern massasauga rattlesnake (Sistrurus catenatus catenatus) population viability in Barry County, Michigan, developed from field data collected 2004, 2005, 2008, and 2009 and data from the literature. Assumptions included no inbreeding depression, no concordance of reproduction and survival due to environmental variation, no change in carrying capacity, stable age distribution, and normal distribution of number of young/yr ...... 29 Table 2.2. Model parameter inputs were varied for sensitivity testing during a population viability analysis (PVA) for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan ............................................. 31 Table 2.3. Population viability model responses to variation in demographic parameters for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan, including probability of extinction in 50 years (PE50; n = 1,000), mean time to extinction (TE) rounded to nearest year (n = 1,000PE50), and stochastic growth rate (rs; n = 50,000) ..................................................................................... 35 Table 3.1. Herbicides applied from 2006—2009 for vegetation control at the study area in Barry County, Michigan, including species treated and year ............................ 48 Table 3.2. Integrated forest management analysis program/ gap analysis program (IFMAP; MDNR 2003) land cover categories found at the study area in Barry County, Michigan, reclassified according to floristic and structural characteristics. Cover types were considered wetland (W) or upland (U), coniferous (C) or deciduous (D), and qualified as either early to mid-successional (EM) or late-successional (L). All others were classified as DEVL ........................................................................ 53 Table 3.3. Stand-level vegetation variables collected at eastern massasauga rattlesnake (Sistrurus catenatus catenatus) locations in Barry County, Michigan from 2008 to 2009. .. ............................................................................................................................. 56 viii Table 3.4. Ranking matrix for selection of cover types used by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan from 2004— 2009 (n = 45). Larger ranks indicate higher preference .................................... 61 Table 4.1. Stand-level vegetation variables collected at eastern massasauga rattlesnake (Sistrurus catenatus catenatus) locations and random locations in Barry County, Michigan from 2008 to 2009. All variables were collected at snake presence locations, and LHC, DHC, SDTS, and ADT were also collected at random locations ............. 81 Table 4.2. Integrated forest management analysis prograrn/ gap analysis program (IF MAP; Michigan Department of Natural Resources 2003) vegetation type categories found at the study area in Barry County, Michigan, reclassified and quantified according to site-specific characteristics. Vegetation types were considered wetland (W) or upland (U), coniferous (C) or deciduous (D), and qualified as either early to mid-successional (E-M) or late-successional (L). All others were classified as DEVL ............................................................................................. 83 Table 4.3. Eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) that were radiotracked in Barry County, Michigan (2008 and 2009) were scored according to 3 fitness indices: breeding success, survival from release to hibernation, and home range size. Adding individual index scores yields the overall fitness score for any snake (possible scores range from 2—9). Breeding success and home range size point values were assigned according to whether a snake was within or beyond the mean i ISD for that value ......................................................................................... 86 Table A.1. Common plant species associated with 3 height strata in early to mid- successional deciduous uplands and wetlands inhabited by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) from 2008 to 2009 in Barry County, Michigan .......................................................................................... 109 Table B.l. Deterministic population growth rates calculated with Vortex 9.96 (Chicago Zoological Society, Brookfield, IL) for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan, derived fi'om a baseline scenario. Assumptions included no stochastic fluctuations, no inbreeding depression, no limitation of mates, no harvest, and no supplementation .................................... 111 Table CI. The snake community occurring in cover types used by the focal species, eastern massasauga rattlesnake (Sistrurus catenatus catenatus; n = 69), included 7 species of snake at the study area in Barry County, Michigan fi'om 2008 to 2009. Observations have been contributed to the statewide online public database, Michigan Herpetological Atlas (MDNRE 2010) ......................................................... 112 ix LIST OF FIGURES Figure 1.1 Shaded area represents the range of the eastern massasauga rattlesnake (Sistrurus catenatus catenatus) spanning the United States and southern Ontario, Canada (New York State Department of Environmental Conservation 2010).... . . . . . . . ...3 Figure 1.2. Eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) were studied on private property in Barry County, Michigan (represented by open star) in 2008 and 2009 ....................................................................................... 4 Figure 2.1. Over a simulated 50-year period for which population viability was modeled for eastern massasauga rattlesnakes (Sistrums catenatus catenatus) in Barry County, Michigan, probability of extinction, PE50, decreased as initial population size, No, was increased during a sensitivity analysis ................................................ 36 Figure 2.2. The sensitivity of a population viability analysis for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan, to increased mortality was assessed. Juvenile (<3 yr) and adult (23 yr) mortality were increased by 10 and 20% to measure effects on (A) stochastic growth rate, rs, (B) probability of extinction, PE50, and (C) final population size, N50 ........................................... 37 Figure 3.1. Home range sizes (:t 2 SE) for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan, by reproductive group: male (M), nongravid female (N GF), and gravid female (GF). Each reproductive group is summarized 3 ways: premanipulation (Pre, 2004—2005), postrnanipulation (Post, 2008— 2009), and pooled (Pooled, 2004—2009) ..................................................... 58 Figure 3.2. Percent of cover types available to eastern massasauga rattlesnakes (Sistmnrs catenatus catenatus) from 2004—2009 in Barry County, Michigan, at the study area level. Cover types were classified as wetland or upland (W or U), coniferous or deciduous (C or D), and designated as early-mid, or late successional (EM or L); all else was developed (DEVL) ..................................................................... 60 Figure 3.3. Mean percent of eastern massasauga rattlesnake (Sistrurus catenatus catenatus) home ranges (:L- 2 SE) occurring in each management type available at the study area in Barry County, Michigan, from 2004—2009: burning with shrub removal (B+SR), shrub removal only (SR), no shrub removal or burning (NM) .................. 63 Figure 4.1. Production functions for the 6 variables in a habitat suitability index (HSI) model (Bissell 2006, amended in 2010) for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in southwestern Michigan. Model variables include A) percent live herbaceous cover, B) percent dead herbaceous cover, C) stem density of trees and shrubs >3 m, D) absolute dominance of trees >3 m, E) percent area in early successional deciduous uplands, and F) percent area in early successional deciduous wetlands ...... 76 Figure 4.2. Percent live herbaceous cover in stands occupied by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) was different than in random stands in southwestern Michigan for one cover type and for all cover types pooled (a = 0.10). Frequently used cover types in 2008 and 2009 were early and mid-successional deciduous uplands (U DE and UDM, respectively) and early and mid-successional deciduous wetlands (WDE and WDM, respectively) ....................................... 90 Figure 4.3. Percent dead herbaceous cover in stands occupied by eastern massasauga rattlesnakes (Sistrums catenatus catenatus) was different than in random stands in southwestern Michigan for 2 cover types (a = 0.10). Frequently used cover types in 2008 and 2009 were early and mid-successional deciduous uplands (UDE and UDM, respectively) and early and mid-successional deciduous wetlands (WDE and WDM, respectively) ....................................................................................... 91 Figure 4.4. Stern density of trees and shrubs >3 m in stands occupied by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) was different than in random stands in southwestern Michigan for 2 cover types (a = 0.10). Frequently used cover types in 2008 and 2009 were early and mid-successional deciduous uplands (UDE and UDM, respectively) and early and mid-successional deciduous wetlands (WDE and WDM, respectively) .............................................................................. 92 Figure 4.5. Absolute dominance of trees >3 m in stands occupied by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) was different than in random stands in southwestern Michigan for one cover type and for all cover types pooled (a = 0.10). Frequently used cover types in 2008 and 2009 were early and mid-successional deciduous uplands (UDE and UDM, respectively) and early and mid-successional deciduous wetlands (WDE and WDM, respectively) ....................................... 93 Figure 4.6. I applied Bissell’s (2006) HSI model for eastern massasauga rattlesnake (Sistrurus catenatus catenatus) at 30 replication sites in the study area (Barry County, Michigan) in 2008 and 2009. Habitat suitability index (HSI) scores were assigned to each of 4 categories (poor: 0.00—0.25, marginal: 0.26—0.50, good: 0.51—0.75, and optimal: 0.76—1.00) to describe their distribution resulting from model application at 2 ha, 10 ha, and 20 ha ............................................................................... 96 Figure 4.7. Percent change in HSI scores resulting from a 10% change in each variable within Bissell’s (2006) habitat suitability index (HSI) model for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in southwestern Michigan. Sensitivity was tested for each of 6 variables individually across 30 replicates at 2 scales of model application (2 ha and 10 ha). Variables include percent live herbaceous cover (LHC), percent dead herbaceous cover (DHC), stem density of trees and shrubs >3 m (SDTS), absolute dominance of trees >3 m (ADT), percent area in early successional deciduous uplands (AEDU), and percent area in early successional deciduous wetlands (AEDW) ............................................................................................ 97 xi Figure 4.8. Proposed modifications to 4 production functions in Bissell’s (2006) habitat suitability index (HSI) model for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in southwestern Michigan. Dashed lines represent the modified maxima for A) percent live herbaceous cover, 62—100%, B) percent dead herbaceous cover, 51.5- 96%, C) percent area in early successional deciduous upland cover type, 0—57%, and D) percent area in early successional deciduous wetland cover type, 23—73% .................... 99 xii GENERAL INTRODUCTION Forest species composition is changing throughout the Great Lakes region, with early successional species giving way to late seral communities (Lorimer 2001). Charcoal and pollen records have shown that northern North America was historically subject to spatially heterogeneous and episodic disturbances such as fire and severe windthrow (Lorimer 2001). Many ecosystems in North America were shaped by these disturbance patterns. In Michigan, European settlement (18005) brought anthropogenic changes (e.g., agriculture, fire suppression), which dramatically altered Michigan’s landscape and caused native grassland species to decline markedly (W interstein et a1. 1995). At the same time, wetlands throughout the United States have been drained and converted to alternative land uses, practices often sanctioned by the government (Laubhan et a1. 2005). Although wetlands are now protected under several national regulations (e.g., Water Resources Development Act of 1990), impacts Of wetland loss on wildlife species have been severe (Laubhan et a1. 2005). Seldom found far from upland grasslands or wetlands, the eastern massasauga rattlesnake (Sistrurus catenatus catenatus; hereafter, massasauga) has also declined markedly throughout its range due to the loss of its early successional habitat, a problem exacerbated by human antagonism towards the species (Szyrnanski 1998, Parent and Weatherhead 2000). Although early observers noted population declines as early as the 1940s (Wright 1941), the massasauga was not widely recognized as imperiled until the 19708 (Szyrnanski 1998). The eastern massasauga is now classified as threatened or endangered in every state or province in which it occurs except Michigan, where it is a species of Special Concern. Massasaugas range from western New York and southern Ontario westward to Iowa and southward to Missouri (Figure 1.1). States and provinces comprising the range of the snake include: Illinois, Indiana, Iowa, Michigan, Minnesota, Missouri, New York, Ohio, Ontario, Pennsylvania, and Wisconsin (Johnson et al. 2000). Due to range contraction and population declines, by 1998 less than 45% of all states’ populations were considered secure (Szyrnanski 1998). Extirpation from many parts of its historical range (defined as no confirmed observation since 1967; Szyrnanski 1998) was a contributing factor when the massasauga was given candidate species status in 1999 by the US. Fish and Wildlife Service (U SFWS), granting it consideration for federal listing as threatened or endangered under the Endangered Species Act (U SFWS 1999). Because of their limited mobility, low levels of gene flow, and the high number of unique alleles in geographically distinct populations, each individual massasauga population is of high conservation value (Gibbs et al. 1997). Michigan is considered the stronghold of the species (Szymanski 1998), and yet long-terrn demographic data (i.e., >5 years) from massasaugas in Michigan is nonexistent in the primary literature. I address these data gaps by first quantifying adult survival using data collected in 2008 and 2009, so that survival results can be compared with studies of similar duration conducted elsewhere. I also combined my demographic data with that of a previous study (Bissell 2006) to model population viability for southwestern Michigan for the next 50 years. My study differs from previous research conducted in Barry County, Michigan (Figure 1.2), because the studied property has since undergone vegetation management to maintain open seral stages, including mowing, burning, and invasive species removal (mechanical and chemical). It is now Eastern Massasauga Figure 1.1. Shaded area represents the range of the eastern massasauga rattlesnake (Sistrurus catenatus catenatus) spanning the United States and southern Ontario, Canada (New York State Department of Environmental Conservation 2010). llL/ —L_L_L/ Figure 1.2. Eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) were studied on private property in Barry County, Michigan (represented by open star) in 2008 and 2009. possible to quantify snake habitat selection patterns in relation to active management and compare movement patterns to pre-manipulation conditions. Finally, I applied a proposed habitat model (Bissell 2006) to the study area to evaluate its performance at multiple spatial scales. A credible habitat model for massasaugas would give managers a reliable tool that can explicitly model the effects of proposed landscape changes on massasauga rattlesnakes. The specific objectives of my study were to: 1. Quantify survivorship and cause-specific mortality of massasaugas in southwestern Michigan, 2. Model population viability for the studied massasauga population in southwestern Michigan over the next 50 years, 3. Determine whether and to what extent resource selection by massasaugas is influenced by current habitat management practices applied in the study area, and 4. Evaluate whether a habitat suitability index (HSI) model for massasaugas accurately depicts habitat suitability in southwestern lower Michigan As a result of this research, state and federal wildlife biologists and researchers will have a more detailed understanding of the factors influencing massasauga demographic dynamics, population viability, habitat selection, and habitat quality. This information can be used to guide decision-making processes when managing land for eastern massasauga rattlesnakes. CHAPTER 1: Survival and cause-specific mortality of eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Michigan ABSTRACT The eastern massasauga rattlesnake, Sistrums catenatus catenatus, is a candidate for US. federal listing as threatened or endangered and is legally protected in every state or province in which it occurs. Habitat degradation and human persecution have contributed to range-wide population declines. Survival estimates are essential for a thorough understanding of population dynamics, yet are seldom reported for eastern massasaugas in the literature. There has been little research on massasauga survival in managed areas of Michigan, the state considered to be the species’ stronghold. My objectives were to: 1) quantify survival of adult eastern massasauga rattlesnakes during the active season (May—Oct), and 2) describe cause- specific mortality of massasaugas in southwestern lower Michigan. I captured, radiomarked, and monitored 27 adult massasaugas in 2008 and 2009 and pooled data for analyses. I observed snakes throughout the active season and estimated survival (Mayfield method) for that period (11 May—29 Oct; 168 days). Cause-specific mortality was investigated qualitatively. Estimated survival probability for the active season was 0.9472 (SD = 0.0509), higher than any estimate for similar studies. The ' single mortality event observed was due to predation. I discuss factors influencing this high survival estimate in a management context. These results suggest that adult massasauga survival is potentially geographically and temporally variable and can be high in areas not well investigated, such as southwestern lower Michigan. More research would yield long term survival trends for the studied population in Michigan. INTRODUCTION To effectively manage for a species of conservation concern, biologists must understand sources and rates Of mortality under field conditions. The eastern massasauga rattlesnake, Sistrurus catenatus catenatus, is a candidate for US. federal listing as threatened or endangered (U SFWS 1999) and a species of interest to managers trying to curb rapid population declines throughout the Midwestern United States and Canadian range. Primary reasons for population declines include habitat degradation and human persecution (Szyrnanski 1998, Johnson et al. 2000). Although population declines are well-documented, detailed survival data are lacking for massasaugas throughout their range (Szyrnanski 1998). Michigan has been described as the stronghold for massasaugas even though 33% of state historical populations have been extirpated (Szyrnanski 1998, Johnson et a1. 2000), making it critical that research be conducted to assist land and wildlife managers, biologists, and researchers in developing management plans for massasaugas. Understanding population demographics is essential for managing massasaugas, and yet few studies provide more than an estimate of apparent survival. Apparent survival is easy to calculate and requires few assumptions but is also an inefficient use of data and does not accommodate censoring, thus biasing the estimate if censoring occurs (W interstein et al. 2001). The peer-reviewed massasauga literature focuses primarily on spatial ecology, habitat modeling, or effects of disturbance, providing only cursory estimates of survival (King 1999, J ellen and Kowalski 2007). Published estimates of adult survival are typically specific to local phenomena, and I am aware of no peer-reviewed studies which have quantified survival in Michigan, where the massasauga is of Special Concern. Published survivorship data for free-ranging eastern massasauga populations are rare (Johnson et al. 2000), but some authors have provided informal estimates. In Wisconsin, King (1999) found apparent survival probability to be 0.33 for adults in a 3- year study (n = 12). J ellen and Kowalski (2007) reported apparent survival for 9 gravid females and 12 successfully tracked neonates in Pennsylvania as 0.56 and 0.92, respectively, based on one year of data. Harvey and Weatherhead (2006a) reported similar survival rates for winter and summer, 0.77 (SD = 0.06) and 0.79 respectively (n = 44), based on a 3-year study on Bruce Peninsula, Ontario. In southwestern lower Michigan, adult annual survival (n = 26) was estimated as 0.66 (SD = 0.22) from a 2- year study, the only available estimate for a Michigan population (Bissell 2006). With adult survival estimates ranging fiom 0.79 (Harvey and Weatherhead 2006a) to as low as 0.33 (King 1999), it is clear that more research on regional variability is needed to help model populations. Population viability analysis (PVA) for massasaugas has demonstrated population sensitivity to adult and neonate survival in Missouri and Ontario (Seigel and Sheil 1999, Miller 2006). I concur with Miller (2006) that effort should be directed toward more precise estimation of these parameters in the field. However, studies designed to document the effects of localized disturbances, both natural and experimental, could result in low survival estimates that apply to geographically or temporally restricted extents. For example, researchers have quantified the effects of mowing and summer burning (Durbian and Lenhoff 2004, Durbian 2006) and repatriation (King et al. 2004) on apparent survival, but these estimates may not be appropriate for populations elsewhere because of the disturbances studied. My Objectives were to: 1) quantify survival of adult eastern massasauga rattlesnakes during the active season (May—Oct), and 2) describe cause-specific mortality of massasaugas in southwestern lower Michigan. The results of my study can help inform other massasauga management and conservation approaches, such as population viability analysis or population modeling, repatriation efforts, and land management. STUDY AREA I conducted this study on a tract of contiguous non-profit private property totaling ~267 ha in Barry County, southwestern lower Michigan. Major plant community types included southern wet meadow, relict conifer swamp, mesic southern forest, dry southern forest, dry-mesic southern forest, mixed wetlands, and prairie fen (Slaughter and Skean 2003). Massasaugas at the site have been documented to primarily use early and mid-successional wetland and upland vegetation types (Bissell 2006, R. Bailey, Michigan State University, unpublished data), which are predominantly deciduous (Appendix A). Dominant soils included Penington and Marlette loam and Oshtemo sandy loam with some pockets of muck (Natural Resources Conservation Service 2000). Average summer and winter temperatures in nearby Hastings, Michigan (fi'om 1971—2000) were 20.4° C and -4.0° C, respectively (Michigan Department of Agriculture 2004). Average annual precipitation for the same time period was 90.9 cm (Michigan Department of Agriculture 2004). Southwestern lower Michigan is somewhat humid, with an average annual relative humidity of 62% (Michigan Department of Agriculture 2004). The terrain of central Barry County is rolling to hilly with many glacially-formed lakes (Michigan Department of Agriculture 2004). The property is bisected by an unpaved road that crosses Cedar Creek, which flows perpendicular to the road. Property managers used a combination of prescribed fire, herbicides, and mechanical control to promote early successional upland and wetland vegetation types and discourage invasive or encroaching woody species (e.g., multiflora rose, Rosa multrflora). METHODS Capture and Telemetry Radiotelemetry is a well-established and accepted method of estimating survival in free- ranging wildlife populations (W interstein et al. 2001). I used visual surveys to locate snakes in the spring and summer of 2008 and 2009, and all snakes were captured with snake tongs. Veterinary staff at Potter Park Zoo (Lansing, MI) surgically implanted, under sterile conditions, radiotransmitters in adult snakes weighing >100 g. Surgery generally followed the procedures described by Weatherhead and Anderka (1984) with modifications for this study. Modifications included anesthesia with Isoflurane, administration of pain relief and antibiotics post-operatively, introduction of a sterile 8 gauge polypropylene catheter (rather than a copper tube) to position antennas subcutaneously, and warming of snakes throughout surgery. Veterinary staff warmed snakes with a water-circulating heating pad during surgery and a 60-watt timed lamp during recovery. Beginning in June 2008, the Potter Park Zoo veterinarian (T. M. Harrison) began giving subcutaneous injections of an electrolyte solution post- operatively to enhance snake hydration. Depending on snake size, transmitters attached were either 0.8 g, 3.1 g, 8.0 g, 10.0 g, (Models A2410, R1680, R1520, R1520 respectively, Advanced Telemetry 10 Systems, Inc., Isanti, MN) or 5.0 g (Model SB-2, Holohil Systems Ltd., Ontario, Canada). I tagged all snakes with a Passive Integrated Transponder (PIT tag; Avid Identification Systems, Norco, CA) for permanent individual identification; transmitters and PIT tags weighed <5% of the individuals’ body weight to minimize marking effects on snakes (Parent and Weatherhead 2000, Withey et al. 2001, Fuller et al. 2005). Snakes were released at their capture site after being held individually for a 3—5 day recovery period at the Potter Park Zoo veterinary clinic. I relocated snakes in the field 6 days per week from May through August and at least once a week thereafter until overwintering began (Oct). I excluded 2 snakes fiom analysis that did not appear to be responding normally to surgery (i.e., lethargy preceded death at 5 and 13 days after surgery). Therefore, I only included data from snakes which responded normally in the results. Michigan State University’s Institutional Animal Care and Use Committee approved all capture, holding, and surgical procedures for this study (permit 05/08-067-00), and I conformed to regulations under the Michigan Department of Natural Resources scientific collector’s permits (issued 20 Nov 2007 and 12 Feb 2009). Survival Estimates I pooled data from 2008 and 2009 to increase sample size, assuming minimal yearly affects. I used Mayfield’s (1961) method adjusted for censoring (Bunck and Pollock 1993) to estimate daily and period survival probabilities. The Mayfield method is appropriate when n < 50 animals (W interstein et al. 2001). Assumptions of the Mayfield method include random sampling, independence of animals, independence Of observations, working radios are always located, transmitter attachment does not impact survival, random censoring, and constant survival for a period (W interstein et al. 2001). 11 Individuals tracked over both field seasons were treated as independent when they emerged successfully fi'om hibernation in spring 2009. For all calculations, entry time (to) was the day after a snake was returned to its capture location. Whenever the date of death or censoring was not known exactly, the midpoint between the last visual encounter and the last date tracked was substituted for the end time (ti) because I could not confirm that a snake was alive without a visual or aural observation. For snakes that survived to the end of the study, I defined ti as the midpoint between the last visual observation and the final tracking date, if they were not the same. Therefore, number of exposure days (E) was truncated slightly with negligible effects on results. The active season, defined as 11 May—29 October (168 days), corresponded with the earliest and latest days that researchers observed massasaugas active on the site during the study period. Daily, 30-day, and 90-day survival estimates are also reported to facilitate comparisons with existing or firture studies. I used Johnson’s (1979) procedure to calculate 95% confidence intervals for all estimates. Cause-specific mortality was investigated qualitatively. RESULTS I captured and radio-tagged 27 snakes from 11 May—l7 August in 2008 and 2009. Two snakes did not recover from surgery and were excluded fi'om subsequent analysis. Necropsy revealed these individuals were in poor health due to existing parasite burdens, causing adverse reactions to surgery. There were no signs of infection at the surgery sites. Of the 25 remaining snakes, 5 individuals were tracked both years. There were 10 males, 12 nongravid females, and 8 gravid females in the pooled sample. In 2008, no censoring occurred, and one snake died (n = 14). In 2009, I censored 7 snakes 12 but confirmed no mortalities (n = 16). I attributed censoring to transmitter expiration (n = 5) and snakes moving beyond the receiving range (n = 2) based on transmitter specifications. Total combined exposure days, E, was 3,104 during the 2 seasons of my study. The assumption of constant survival was appropriate because only one mortality was confirmed in 2 years. Low mortality precluded analysis of predictive variables (e.g., sex, reproductive condition). Daily survival probability, .5 d, was 0.9997 (95% CI = 0.9999—1.0003) and was similar to 30-day survival, 5* 30, Of 0.9904 (95% Cl: 0.9718—1.0093). Ninety-day survival, S 90, was 0.9714 (95% c1 = O.9177—1.0281) and 168-day survival, S 168, was 0.9472 (95% CI = 0.8518—1.0532) for the active season (Table 1.1). Circumstantial evidence suggests that the individual that died was most likely the result of mammalian predation. The death was discovered on 1 October 2008 (72 days after release); only the transmitter was found, which had been stripped, undamaged, of tissue. The depredated adult male weighed 298 g (snout-vent length = 61.4 cm) and was apparently healthy. I am also aware of at least 3 road-killed massasaugas during 2008 which were all untagged: l neonate observed, 1 adult reported by the facilities manager, and the desiccated remains of an adult known to property managers prior to the beginning of my study. I found no rattlesnake road mortalities in 2009, although other species were occasionally found dead on the road (e. g., Thamnophis sirtalis, Storeria dekayi, Coluber constrictor foxi, Heterodon platirhinos, and Ihamnophis sauritus septentrionalis). Road mortality is evidently a source of death as well as predation, even though it was not experienced by any tagged snakes. l3 Table 1.1. Periodic survival probabilities for radiotracked adult eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in southwestern lower Michigan estimated from the 168-day activity seasons (mid-May—Oct) of 2008 and 2009. Period 5‘ p SD 95% CI21 Daily 0.9997 0.0003 09999—10000 30-day 0.9904 0.0092 0.9718—1.0000 90-day 0.9714 0.0278 09177—10000 168-day 0.9472 0.0509 08518—10000 a Calculated using Johnson’s (1979) method and right-truncated. 14 Table 1.2. A compilation of range-wide survival estimates for adult eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) during the active season known fi'om the literature and this study, listed in order of decreasing survival estimates. Data from this study were collected from May—October of 2008 and 2009 in southwestern lower Michigan. Study 8’ p (SD) n Method Present studya 0.95 (0.05) 30 Mayfield Harvey and Weatherhead 2006ab 0.79 (0.06) 34 Apparent survival Bissell 2006al 0.66 (0.22) 26 Kaplan-Meier Jellen and Kowalski 2007c 0.56 (NA) 9 Apparent survival King 1999(1 0.33 (NA) 12 Apparent survival 3 Michigan. b Ontario, Canada. C . Pennsylvanla. d . . Wrsconsrn. 15 DISCUSSION My active season survival estimate, 0.9472, is the highest available estimate of similar studies of eastern massasauga (Table 1.2), but some caution is warranted in interpreting these results. Winterstein et al. (2001) noted that when censoring is >10—15% and n < 50, survival estimates may suffer fi'om decreased precision. In this study, 23.3% of the sample was censored, which could have decreased the precision of the estimate. Pooling data across 2 years may have resulted in lost information, but with low sample size and only one death, I would be unable to detect yearly differences. Small sample size and independence of animals are common tradeoffs among studies of eastern massasaugas. Because 5 snakes were treated independently across 2 years, the assumption of animal independence may not apply between years. Because the 2 snakes that died during the acclimation period were in ecdysis at the time of surgery, I recommend not surgically implanting transmitters in shedding snakes as the physical stress of shedding combined with an unknown parasite burden may additively affect snake recovery while obscuring the ability of observers to detect lethargy. Evidence suggests that direct sources of massasauga mortality in Barry County, Michigan, besides predation and road kills, include human persecution (R. Bailey, personal observation), factors which have long been recognized to contribute to massasauga population declines (Wright 1941, Johnson et al. 2000). Additionally, massasaugas have reportedly been killed accidentally by cyclists on public trails in Jackson County, Michigan (R. Grace, personal communication), whereas bicycles and recreational vehicles are prohibited at my study area. I attribute the high survival of this sample to the habitat quality within the study area. Snakes at the study area benefited 16 fiom a relatively undeveloped infi’astructure (i.e., few buildings and one unpaved road) under single ownership, bordered by some suitable habitat, as well as protection from persecution by humans. Vegetation management on the property is designed to increase open-canopy vegetation types. To compare, a study by Durbian et al. (2008) found as many road-killed massasaugas (at least 10) as live ones at one of their unmanaged study sites in Wisconsin (survival not estimated); the authors noted that snakes were limited by the availability and juxtaposition of open-canopy vegetation in the absence of active management. Attributes of habitat and cover must play an indirect role in massasauga survival, and further research needs to be done linking vegetation variables to survival. Lack of suitable habitat could concentrate snakes into unfavorable areas (i.e., railroad corridors, roads), leading to increased predation and road mortality (Durbian et al. 2008). Conversely, habitat fiagmentation and development could isolate individuals in patches of poor quality habitat, causing them to travel widely to fulfill life requisites. As this study was conducted on an actively managed property, I suggest habitat management activities such as prescribed fire and mechanical or chemical control of woody species can be compatible with high survival of eastern massasaugas if implemented in the appropriate season. Where undue mortality is imposed by roads, trails, or other infi’astructure on a property, habitat management for massasaugas could compensate for any incidental mortality caused by management activities. I encourage other authors to report survival results to facilitate regional comparisons associated with land ownership and habitat conditions. My survival estimates can help biologists develop biologically plausible minimum and maximum inputs for population viability l7 analysis and other population models. These results also have implications for guiding land management to conserve populations. Additionally, as King et al. (2004) discussed, baseline estimates of mortality are useful standards for comparison when evaluating the success of repatriation or reintroduction efforts. ACKNOWLEDGMENTS I thank W. Anderson, S. Zimmer, P. Paradine, T. Harrison, D. Paperd, and various volunteers for assistance with snake surveying, monitoring, and surgeries. Potter Park Zoo (Lansing, MI), Pierce Cedar Creek Institute (Hastings, MI), and Michigan State University Diagnostic Center for Population and Animal Health (East Lansing, MI) provided financial and logistical support. I also thank various landowners in the study area for allowing access to their property. I am gratefirl to S. Winterstein for analytical support. This study was funded by Michigan State University and the Michigan Department of Natural Resources and Environment (Wildlife Division). 18 CHAPTER 2: Population viability of the eastern massasauga rattlesnake (Sistrurus catenatus catenatus) in Michigan ABSTRACT The eastern massasauga rattlesnake (Sistrurus catenatus catenatus; hereafter, massasauga) is a candidate for federal listing as threatened or endangered, and is a species of Special Concern in Michigan. To help manage declining species, population viability analysis (PVA) has been used as a tool for predicting a species’ vulnerability to extinction. However, PVA has rarely been attempted for snakes. Even though Michigan is considered the stronghold for the eastern massasauga, little is known about their ecology in Michigan. My objectives were to: 1) quantify demographic parameters needed for PVA from field data, 2) predict whether massasaugas in southwestern lower Michigan will persist for the next 50 years, and 3) quantify how sensitive population responses are to uncertain components Of the developed model. I obtained data from 51 adult massasaugas captured and radiomarked (May—Aug) during 2004, 2005, 2008 and 2009 and pooled demographic data for analyses. I estimated 15 demographic parameters for PVA fiorn snakes Observed throughout the active season (May—Oct). I derived estimates from the literature when field data were insufficient. 'I used Vortex 9.96 to model population viability and to conduct sensitivity tests for some of my less certain estimates. Estimated adult mortality for the active season was low (20.8%), and mean litter size, 7.0, was near the lower end of range-wide estimates for massasaugas. Current demographic conditions resulted in a relatively low probability of extinction (0.002) for 1,000 simulated populations by 2060, if environmental and demographic conditions remain relatively l9 stable. Changes in the mortality of juvenile snakes influenced population persistence more than equivalent changes for adult snakes. Ecologists can use my results to help conserve massasaugas in Michigan and perhaps elsewhere. INTRODUCTION Understanding population demographics is essential for managing a species of conservation concern. The eastern massasauga rattlesnake (Sistrurus catenatus catenatus; hereafter, massasauga) is a candidate for federal listing as threatened or endangered, and a species of interest to biologists trying to curtail rapid population declines throughout its Midwestern and Canadian range. Biologists commonly use demographic models when managing threatened or endangered wildlife (Beissinger and Westphal 1998). One method, population viability analysis (PVA), is a stochastic modeling technique that estimates the probability of a population of a given size surviving for a specified time frame (White 2000), the results of which are used to estimate the vulnerability of a population to extinction (Miinzbergova and Ehrlén 2005). Yet, PVA for snakes has rarely been attempted (Seigel and Sheil 1999), likely because the extensive data required by the approach are generally lacking. When used as part of a suite of analytical techniques, PVA can serve as a useful tool to understand population characteristics and support conservation of the massasauga. Additionally, PVA has been used to determine the sensitivity of results to the variability associated with input parameters (Beissinger and Westphal 1998, Ellner et al. 2002, McLoughlin and Messier 2004). For example, factors such as error in initial population size (McLaughlin and Messier 2004), age-specific mortality rates (Seigel and Sheil 1999, Bissell 2006, Miller 2006), and reproductive traits (Miller 2006) have 20 been found to contribute heavily to projections made by models for snakes and other taxa. Sensitivity analysis is, therefore, an essential part of any PVA with high uncertainty in input parameters. Estimates of initial population size or density and juvenile mortality rates are especially lacking for massasaugas given their cryptic habits, low recapture rates, and conservation status (Johnson et a1. 2000). Even though Michigan is considered the eastern massasauga stronghold (Szyrnanski 1998), little is known about their ecology throughout Michigan, where they are listed as a species of Special Concern. Even with a range of plausible inputs for initial population size and subadult mortality rates, Bissell (2006) found that the likelihood of extinction was never <0.818 in 50 years for massasaugas in southwestern lower Michigan. My objectives were to: 1) quantify demographic parameters needed for PVA from field data, 2) predict whether eastern massasaugas in southwestem lower Michigan will persist for the next 50 years, and 3) quantify how sensitive population performance is to uncertain components of the developed model. My results should aid current conservation efforts for the massasauga, highlight which parameters contribute most to probability of population persistence, and identify which demographic parameters are most uncertain. STUDY AREA I conducted my study on a contiguous tract of non-profit private land in Barry County, southwestern lower Michigan in the Grand River watershed. Major plant community types on the area included southern wet meadow, relict conifer swamp, mesic southern forest, dry southern forest, dry-mesic southern forest, mixed wetlands, and prairie fen (Slaughter and Skean 2003). The study area is ~267 ha consisting of forest, wetland, 21 upland field, restored prairie, and Open water. Massasaugas in the area primarily use early and mid-successional wetland and upland vegetation types (Bissell 2006, R. Bailey, unpublished data). Average summer and winter temperatures in nearby Hastings, Michigan (from 1971—2000) were 20.4° C and -4.0° C, respectively (Michigan Department of Agriculture 2004). Average yearly precipitation for the same time period was 90.9 cm (Michigan Department of Agriculture 2004). Since 2003, the property has been managed using a combination of prescribed fire, herbicides, and mechanical control to promote early successional upland and wetland vegetation types and to discourage invasive or encroaching plant species (e.g., autumn olive, Elaeagnus umbellata; J. Howell, property land steward, personal communication). METHODS Capture and Telemetry I used visual surveys to locate rattlesnakes in the spring and summer of 2008 and 2009. I transported adult snakes weighing >100 g to Potter Park Zoo (Lansing, M1) for surgical implantation of radio transmitters by a veterinarian (T. M. Harrison). Surgery followed the procedures described by Weatherhead and Anderka (1984), with modifications for this study. Modifications included anesthesia with Isoflurane, administration of pain relief and antibiotics post-operatively, introduction of a sterile polypropylene catheter of 2.7 or 3.3 mm diameter (8 or 10 French, respectively) rather than a copper tube to position antennas subcutaneously, and warming of snakes throughout surgery and post-operatively. Radiotransrrritters weighed either 0.8 g, 3.1 g, 8 g, 10 g, (Models A2410, R1680, R1520, R1520 respectively, Advanced Telemetry Systems, Inc., Isanti, MN) or 5 g (Model SB-2, Holohil Systems Ltd., Ontario, Canada), 22 depending on snake size. I also tagged snakes with a Passive Integrated Transponder (PIT tag; Avid Identification Systems, Norco, CA) for permanent individual identification; transmitters and PIT tags weighed <5% of the individuals’ body weight to minimize marking effects on snakes (Parent and Weatherhead 2000, Withey et al. 2001, Fuller et al. 2005). After a 3—5 day post-surgery period, I released snakes at their capture site and relocated them in the field approximately 6 days/week fiom May—August and at least once a week thereafter until hibernation (mid to late Oct). The earliest and latest dates that I tracked snakes were 11 May and 29 October, respectively. I combined my data with that collected under similar protocol in 2004 and 2005 for a previous study (Bissell 2006). Michigan State University’s Institutional Animal Care and Use Committee approved all capture, holding, and surgical procedures for this study (permits 03/04- 040-00 and 05/08-067-00), which conformed to regulations under the Michigan Department of Natural Resources Scientific Collector’s Permit (dates issued: 3 Dec 2003, 22 Mar 2005, 24 Mar 2006, 20 Nov 2007, and 12 Feb 2009). Demographic Parameters Data necessary to develop my PVA included mortality rates, reproductive characteristics, catastrophe likelihood and its effects on reproduction and survival, mate monopolization, initial population size, carrying capacity, harvest, supplementation, and genetic management. I used the Kaplan-Meier method (Kaplan and Meier 1958) to estimate adult mortality. Start time, to, was considered the day after a snake was returned to its capture location. Whenever the date of death or censoring was not known exactly, then the midpoint between the last visual encounter and the last date 23 tracked was substituted for the end time, ti. For snakes that survived to the end of the study, ti was defined as the midpoint between the last visual observation and the final tracking date, if they were not the same. Assumptions Of the Kaplan-Meier method include random sampling, independence of animals, independence of observations, working radio transmitters are always located, random censoring is allowed, and radio transmitter attachment does not impact survival (Winterstein et a1. 2001). I estimated survival using Kaplan-Meier Survivorship Analysis Program 1.0 (Missouri Department of Conservation, Columbia, MO). Except where noted, I pooled data across years to increase sample size, assuming no yearly effects. For telemetered snakes, I estimated the average number of viable young by analyzing radiographs and sonograms Of gravid females (n = 17). Dr. Harrison imaged radiographs and sonograms at Potter Park Zoo at the time of transmitter implantation and again in mid-July using a portable ultrasound (MicroMaxx Ultrasound System, Sonosite, Inc., Bothell, WA) in the field for 2 reasons: 1) to minimize transporting gravid females, and 2) to confirm that embryos observed in early spring were developing normally. Dr. Harrison used the Doppler color feature to isolate embryos and detect neonate heartbeats. This procedure enabled me to exclude embryos lacking heartbeats fiom my estimate, as most gravid massasaugas were observed to retain 1 or 2 undeveloped young. In rare cases when parturition occurred before the second sonogram, I used the first estimate. When I encountered gravid eastern massasaugas incidentally, I used palpation to determine number Of young (n = 2). For an additional 4 females, litter size was only estimated from the number of neonates Observed at a parturition site. I calculated percent breeding females by averaging 4 yearly estimates 24 Of percent gravid females. This gives the same percent as simply pooling across years but yields a variability estimate. I excluded incidental captures from this variable because determining a snake was nongravid was more difficult than determining that it was gravid in the field. I used mortality data from the literature for snakes <3 years Old due to insufficient field data. Mortality estimates for neonates and juveniles are critically lacking in the literature and, when reported, are highly variable among populations (Johnson et al. 2000). For example, J ellen and Kowalski (2007) found that 6% of a total of 16 neonate eastern massasaugas radiotracked for 50 days were depredated in Pennsylvania, the lowest mortality estimate of neonates available. Miller (2006) estimated 65% mortality of massasaugas age 0—1 for snakes on Bruce Peninsula, Ontario. In Wisconsin, King (1999) found neonate massasaugas to have a mortality rate of a 78.1% (n = 32), but radio-transmitters were >5% of neonates’ weight, so I reduced this to 68.1% assuming that this could have increased mortality by 10% at most, following Bissell (2006). Keenlyne (1968) estimated 50% mortality for first year snakes in Wisconsin. Neonate mortality estimates reported for other late-maturing temperate viperids were 33% (Agkistrodon contortrllx; Fitch 1960), 60% (Crotalus viridis; Woodbury 1951), and 54% (C. viridl's; Macartrrey 1985), as presented by Parker and Plurnmer (1987). I averaged these 7 estimates to derive mortality for the 0—1 age class, 48.0% (SD = 21.8). The consequences of including 3 additional species are that mean mortality was <1% higher (massasaugas only: 47.3%), and variability associated with the estimate was decreased. 25 I report the same neonate and subadult mortality for each sex, presuming no mortality differences between sexes (Diller and Wallace 2002, Brown et al. 2007). I assumed subadult mortality was similar for age classes 1—2 and 2—3. King et al. (2004) reported a subadult mortality rate of repatriated massasaugas as 46.7%; however, since snakes were repatriated (King et al. 2004), I used data fi'om a study on another long- Iived temperate zone pit viper, the western rattlesnake (C. viridis oreganus) in Idaho, to supplement this estimate (Diller and Wallace 2002). These authors determined juvenile (1 - and 2-yr old) mortality to be 23.2%. The average mortality for snakes age 1—3 based on these 2 estimates was 34.9% (SD = 16.6). The consequences of using data from the more common western rattlesnake are that mortality was somewhat decreased and an estimate of variability could be generated. Variability is important for simulating environmental stochasticity. Age at first breeding, defined as the age at which first offspring are born, is 3 years for females and males (Keenlyne 1978, Seigel and Sheil 1999, King et al. 2004), and the mating system is polygamous. Miller (2006) reported that eastern massasaugas would likely have a maximum breeding age of 12 years in the wild based on field observations and captive breeding programs, but this is essentially a best guess given the lack of data. The eastern massasauga has a sex ratio at birth of 1:1 (Keenlyne and Beer 1973), consistent with many snake species (Parker and Plummer 1987, Diller and Wallace 2002, Brown et al. 2007). Mate monopolization, defined as the percent of males that are physically, physiologically, and socially capable of breeding (Miller and Lacy 2005), was set at 100% following Miller (2006). All males 22 years old are 26 assumed capable of breeding (King et a1. 2004); although some males probably fail to obtain mates, I do not consider mates to limit female breeding. The most likely catastrophe facing snakes in the study area is flooding. According to the Federal Emergency Management Agency (FEMA), the area immediately surrounding the major creek running through the property is subject to a 1% annual chance flood, or 100-year flood (FEMA 2009); many tracked snakes inhabited the emergent wetlands around this creek, but outlying high ground is available. Therefore, although a large number of snakes would be affected, a portion of the population not inhabiting the flood zone would remain secure. Seigel et al. (1998) documented the worst flood to affect a massasauga population in Missouri in 100 years, yielding insights into population responses. Because the Missouri population went fi'om 50% females breeding pre-flood to 22.7% post-flood (Seigel et al. 1998), the effect of flooding on reproduction factor is 0.454 (i.e., 50X0.454 = 22.7). Further, because the Missouri population exhibited changes in population structure (e. g., fewer immature snakes) and vigor (e.g., smaller clutch sizes) but abundance did not appear to change, the severity on survival factor is 1.0 (no effect). Only flood effects on percent breeding females and survival are considered in the model because there is no intrinsic way to model such changes in population structure or vigor due to catastrophe. I assumed the population had a stable age distribution at the beginning of each simulation because the precise age-class distribution was unknown (Miller and Lacy 2005). I also assumed no change in carrying capacity, K. Due to their long generation time and low reproductive rates, I did not expect snakes to be limited by carrying capacity within 50 years, and so used an arbitrary value of 30,000 following Bissell 27 (2006). I approximated initial population size (N 0) by using the average minimum and maximum density values, in snakes/ha, for 8 species of viperids over 12 studies (for details see Parker and Plummer 1987) and for 4 estimates for massasaugas (Maple 1968, Johnson 1995, Reinert 1978, Mauger in litt. 1996) as reported by Szymanski (1998). I rounded minimum density values of <1 snake/ha up because field observations suggest that at my study area, at least 1 snake/ha is reasonable (K. Bissell and R. Bailey, unpublished data). I used the resulting minimum of 3.3 snakes/ha and maximum of 5.5 snakes/ha as starting points for the calculation of initial population size based on the literature. Combining density estimates with data on area available to massasaugas at the study area, I created a lower and upper input for No. I interpreted total area available to massasaugas at my study area using radiotelemetry data and vegetation layers in ArcMap 9.2 and applied density estimates from the literature to the available area to generate 2 plausible initial population sizes. I considered approximately 217 ha available to snakes on the property based on habitat selection patterns (Bissell 2006; R. Bailey, unpublished data), yielding a lower bound of ~71 8 snakes and an upper bound of~l,185 snakes. Since harvest, supplementation, and genetic management do not affect massasaugas in the studied area, values for these parameters were not included in the model. PVA and Sensitivity Analysis I modeled population viability for massasaugas in southwestern lower Michigan using the individual-based simulation program Vortex 9.96 (Chicago Zoological Society, Brookfield, IL) customized with my baseline estimates (Table 2.1). Assumptions for 28 Table 2.1. Baseline parameters used for modeling eastern massasauga rattlesnake (Sistmrus catenatus catenatus) population viability in Barry County, Michigan, developed from field data collected 2004, 2005, 2008, and 2009 and data from the literature. Assumptions included no inbreeding depression, no concordance of reproduction and survival due to enviromnental variation, no change in carrying capacity, stable age distribution, and normal distribution of number of young/yr. Parameter Value (SD) Source Age of first Offspring, Q , (3‘ 3 , 3 Keenlyne 1978, Seigel and Sheil 1999, King et al. 2004 Max. age of reproduction 12 Miller 2006 Max. # broods/year 1 Field data Max. # progeny/year 12 Field data Sex ratio at birth 1:] Keenlyne and Beer 1973 % Adult females breeding 54 (14) Field data Mean litter size 7 (2.9) Field data Age 0—1 % mortality, 92 . 6” 48 (21.8) Keenlyne 1968, Parker and Plummer 1987, King 1999, Miller 2006, J ellen and Kowalski 2007 Age 1—3 % mortality, 92 , c? 34.9 (16.6) Diller and Wallace 2002, King et al. 2004 Age 3+ % mortality, SB , 6” 20.8 (13.6) Field data Annual % fi'equency of flood 1.0 FEMA 2009 29 Table 2.1 (cont’d). Severity on reproduction Severity on survival % Males in breeding pool Population size Carrying capacity 0.45 1.0 100 718 30,000 (0.001) Seigel et al. 1998 Seigel et al. 1998 Miller 2006 Field data, Parker and Plummer 1987, Szymanski 1998 Bissell 2006 30 Table 2.2. Model parameter inputs were varied for sensitivity testing during a population viability analysis (PVA) for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan. Level Model parameter Lower Baseline Upper Initial population size (snakes) 130 718 1,185 Carrying capacity (snakes) 1,436 and 3,000 30,000 Mortality rates +10% +20% Age 0—1 mortality (%) 48.0 52.8 57.6 Age 1—3 mortality (%) 34.9 38.4 41.9 + . Age 3 mortallty (%) 20.8 22.9 25.0 31 my model include no inbreeding depression, stable carrying capacity, stable age distribution, no concordance of reproduction and survival due to environmental variation, normal distribution of young per year, no immigration, and no emigration (Miller and Lacy 2005). Results are based on a stage-structured population and a 50- year simulation. Given the lack of field data for No, K, and neonate and juvenile mortality, I conducted sensitivity analyses using a range of inputs for these variables using Vortex 9.96 (Table 2.2). I compared 3 levels of initial population size, N0: 1) base model calculated from a minimum density estimate (718), 2) No calculated fi'om the maximum density (1,185), and 3) minimum number of snakes Observed (130). Vortex imposes additional mortality on simulations exceeding a designated K, thereby restricting population growth beyond environmental carrying capacity (Miller and Lacy 2005). To examine the effect of my choice of K in the absence of data, I simulated 3 levels Of carrying capacity: 1) 1,436 snakes (N0 is at 50% of K, following Miller 2006); 2) 3,000 snakes (No is at ~25% of K); and 3) 30,000 snakes (No is not approaching K; following Bissell 2006). I also compared 3 levels of neonate, juvenile, and adult mortality: base model, base + 10%, and base + 20%. Although I estimated adult mortality directly, I included it in sensitivity testing to assess impacts of different mortality levels on population vulnerability to extinction. For all variations, I ran models for 50 years because predictions become less reliable with increasing time frames (Mills et al. 2005). Response variables reported include stochastic growth rate prior to K truncation (rs), probability of extinction (PE50), time to extinction (TE), and final population size (N 50). 32 RESULTS Demographic Summary After I combined my field data with that of K. Bissell (unpublished data), a total of 51 adult massasaugas were radio-tagged and monitored in 2004, 2005, 2008, and 2009. For the survival analysis, I excluded all data for 3 snakes and one year of data for an additional snake because of insufficient observations (e. g., poor health, few relocations); sample size reflects pooling of 9 snakes across 2 years. Adult mortality was 20.8% (SD = 13.6; n = 57). Maximum number Of young/ year for snakes was 12. Average estimated viable young was 7.0 (SD = 2.9; n = 23). Average percent of breeding females/year was 54% (SD = 22), approximating a biennial breeding cycle. Biennial reproduction is clearly not strict, however, because for 7 females observed for 2 consecutive years, one was gravid both years, 4 were nongravid both years (although 2 mated), and 2 alternated between years. The net effect is biennial reproduction, given that the sample is not independent. Vortex typically uses SD to simulate environmental variability. However, the observed range of percent breeding females, 27—82%, represents more realistic boundaries on female reproductive success than would SD because my high year-to- year variation is almost certainly an artifact of sample size (4 yearly estimates). I substituted 14% for SD, thereby restricting the range of possible values for reproductive success to 26—81% in order to prevent Vortex from sampling from a range of unlikely extremes generated using i 28D (IO—98%). 33 Population Viability Model and Sensitivity Analysis For the baseline model, the probability of extinction, PE50, was 0.002 (SD = 1.41), with mean time to extinction (TE) of 27 years (SD = 9) (Table 2.3). Stochastic growth rate, rs, was positive at 0.072 (SD = 0.368). Viability of the population was relatively insensitive to No. As NO was varied, PE50 was highest for N0 = 130 (0.025; SD = 4.93) and zero for N0 = 1,185 (Figure 2.1), rS was unaffected, and TE did not vary consistently in one direction (Table 2.3). The effects of K on model estimates were negligible except for final population size, N50, which was minimally 1,017 snakes (SD = 424.7) at K = 1,436. I conclude that my choices of K did not greatly affect persistence because minimum K was still at least twice N0, leaving room for population growth. For all other sensitivity analyses, the most noticeable differences were between the baseline model and a 20% increase in the mortality rate, with baseline and 10% levels being similar (Table 2.3). Sensitivity analysis revealed that increased age 0—1 mortality or age 1—3 mortality had nearly identical affects on rs, PE50, and TE. When mortality for both groups increased by 20%, rS decreased by about 63% (Figure 2.2A). For both groups <3 yr, PE50 was maximally 0.013 (SD = 3.58) at the 20% level, 6.5 times the baseline estimate (Figure 2.28). The same was not true for the 20% increase in adult mortality, which caused only a 31% decrease in rS (0.049; SD = 0.372). The same 20% increase in adult mortality caused PE50 to increase only 4.5 times (Figure 2.23). For all sensitivity analyses, TE did not vary consistently in any direction. 34 Table 2.3. Population viability model responses to variation in demographic parameters for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan, including probability of extinction in 50 years (PE50; n = 1,000), mean time to extinction (TE) rounded to nearest year (17 = 1,000PE50), and stochastic growth rate (rs; n = 50,000). Model parameter PEso (SD) TE (SD) rS (SD) Baseline 0.002 (1.41) 27 (9) 0.072 (0.368) Initial population size 130 (snakes) 0.025 (4.93) 28 (12) 0.071 (0.375) 1,185 (snakes) 0.000 (0.00) NA 0.073 (0.366) Carrying capacity 1,436 (snakes) 0.001 (0.99) 31 (0) 0.072 (0.400) 3,000 (snakes) 0.001 (0.99) 23 (0) 0.070 (0.368) Age 0—1 mortality (%) + 10% 0.005 (2.23) 38 (15) 0.049 (0.378) + 20% 0.013 (3.58) 38 (8) 0.026 (0.393) Age 1—3 mortality (%) + 10% 0.005 (2.23) 38 (12) 0.048 (0.372) + 20% 0.013 (3.58) 35 (12) 0.027 (0.385) Age 3+ mortality (%) + 10% 0.000 (0.00) NA 0.061 (0.371) + 20% 0.009 (2.99) 34 (7) 0.049 (0.372) 35 0.030 — 0.025 — 0.020 r 0.015 n PE5° 0.010 ~ 130 718 1,185 No Figure 2.1. Over a simulated 50-year period for which population viability was modeled for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan, probability of extinction, PE50, decreased as initial population size, NO, was increased during a sensitivity analysis. 36 0.08 - A 0.014 , B 0.012 . 0.06 ‘ 8 0.010 . I_m III 0.008 - 0.04 ~ (I. 0.006 . 0.004 . 0'02 0.002 - o r . r 0 Base +10% +20% Base +10% +20% Mortality Level Mortality Level + < 3 yr Mortality -+- 3+ yr Mortality Base +10% +20% Mortality Level Figure 2.2. The sensitivity of a population viability analysis for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan, to increased mortality was assessed. Juvenile (<3 yr) and adult (_>_3 yr) mortality were increased by 10 and 20% to measure effects on (A) stochastic growth rate, rs, (B) probability Of extinction, PE50, and (C) final population size, N50. 37 For all scenarios, PE50 never exceeded 0.025, when N0 = 130 snakes. Among all scenarios examined, average rs was always small and positive (the same is true for deterministic growth rates). Differences between deterministic (Appendix B) and stochastic population growth rates (Table 2.3) for the baseline model can give an indication of the threat Of stochastic events on persistence (Miller and Lacy 2005). DISCUSSION To address my first objective, I calculated or assembled the minimum data necessary to conduct a PVA. Adult mortality (20.8%) was similar to the 23% annual post-juvenile mortality of viperids as reviewed by Parker and Plummer (1987). The median between my estimated neonate and adult mortality rates (34.4%) is similar to my literature-based estimate of subadult mortality (34.9%), meaning that my subadult mortality estimate from the literature is consistent with a linear reduction in mortality from neonates to adults. I believe the trend of adult snakes having lower mortality than juveniles or neonates to be accurate based on other studies (Parker and Plummer 1987, Diller and Wallace 2002, Brown et al. 2007). I nevertheless acknowledge that my inputs for neonate and subadult mortality have the most uncertainty of my base model parameters, having been derived from numerous methods not necessarily comparable. My PVA results should, therefore, be considered valuable based on the relative trends that are demonstrated and not absolute responses. The average litter size (7.0) was similar to an estimate of 6.55 for massasaugas in Pennsylvania, which Reinert (1981) obtained by combining his data with that of 2 other reports (Atkinson and Netting 1927, Swanson 1933), and data reported from Missouri by Seigel and Sheil (1999) of 6.1. This is markedly different from estimates 38 of 11.1 (Keenlyne 1978) in Wisconsin or 9.5 (Miller 2006) for Ontario. My estimate of percent females breeding (54%) was Similar to Reinert’s (1981) estimate of 52%. I emphasize the importance of using localized data for PVA when possible because massasauga productivity varies geographically. I posit that massasaugas will persist in my study area for the next 50 years and report several performance measures for characterizing persistence. Time to extinction responded unpredictably to sensitivity tests, making it a poor indicator of persistence. As White (2000) reports, probability of extinction for a specified time flame is a more realistic indicator of population viability than time to extinction. Therefore, time to extinction should not be considered especially realistic. Estimates of N50 should also be considered increasingly unrealistic through time, since uncertainty increases as the model is projected further into the firture (White 2000); its value lies in showing comparative trends and not in predicting absolute numbers of snakes (Figure 2C). Since certainty about PE50 actually increases with time, it is the best indicator of population performance (White 2000). To determine how sensitive these performance measures are to uncertainty in my inputs, I conducted several sensitivity analyses. Increasing juvenile mortality did not impact PE50 as severely as in Bissell’s model (2006); although PE50 was more than quadrupled after only a 20% increase in mortality, it still did not approach her estimates, although methodological differences exist. The trends of increased PE over time with increased neonate and adult mortality match those of Bissell (2006) and Miller (2006), although magnitude is dissimilar. I also found rS to decline disproportionately to 39 increased mortality of young snakes. Even though rS was always positive, it was still possible for populations to go extinct because of stochasticity (White 2000). My choice of a minimal K = 2N0 was arbitrary since data on carrying capacity for massasaugas are lacking. These results apply only if true N0 is no more than half of true K. However, Bissell (2006) used the same value of K (30,000) and her predictions were considerably more pessimistic than Miller’s (2006) simulations in which K was only 200. Population persistence did respond to decreased N0, and the population’s security may be questionable if N declines below 130. These sensitivity results indicate the quickest route to population recovery (i.e., maintaining or decreasing mortality of snakes <3 yr), but not why observed population declines initially occur (Beissinger and Westphal 1998) In addressing these objectives I have made many assumptions. This population is open, constituting a violation of one assumption. The degree to which this population is isolated (an unknown) will affect whether extinction risk can be mitigated by the occasional immigration Of snakes. I assumed that sampling variance is negligible due to my sample size and did not subtract it flom total variance; if this assumption is incorrect, the extinction probabilities presented would be biased high (White 2000, Miller and Lacy 2005, Mills et al. 2005). I ignored temporal variation in adult mortality, opting for improved sample size. These results should not be considered predictions for the future because many underlying processes (e. g., density dependence, inbreeding depression) simply cannot be modeled with existing data for massasaugas. The model further assumes that the system will not change in the next 50 years, which is unknown. This PVA did not duplicate the results of a previous PVA for the 40 massasauga in southwestern lower Michigan (Bissell 2006); even though I incorporated data used to build the prior model, the inclusion of additional data likely caused variations in results. A primary value of my model is the attention it draws to data gaps for the species and the trends in persistence it reveals. Many researchers have noted the lack of massasauga mortality estimates for immature age classes estimated flom a large sample size (Johnson et al. 2000, Bissell 2006, Miller 2006). Seigel and Sheil (1999) found that an eastern massasauga population in Missouri would need a neonate mortality rate of <80% to maintain a stable population. Even though I increased my baseline neonate mortality estimate (48.0%) by 20% (Table 2.2), it still met this criterion. Bonnet et al. (1999) showed that neonates of multiple genera (including viperids) experience high mortality, particularly immediately after birth when they disperse toward hibemacula. By providing suitable habitat for gestating females adjacent to hibemacula areas, managers can minimize dispersal distances neonates (and adult females) must move to find winter refugia and thereby facilitate population growth. In fact, most deaths in this study (56%) occurred in October, when snakes moved toward wintering areas. Bonnet et al. (1999) also Observed a peak in mortality of adult male snakes during the mating season, a consequence of dispersal to find receptive females. If it became necessary, imposing temporary road closures or reduced speed limits during this brief time period (peaking Jul—Sep) could be effective in reducing seasonal road mortality of adult males (Bonnet et al. 1999). Dispersing adult male snakes should be considered important for genetic diversity and colonization of underutilized habitat. I Observed female massasaugas 41 mating with >1 male per season, and. maintenance of a diverse breeding pool of males can increase population fitness. A large deficit of demographic data for massasaugas still exists. There is no reason to assume that unstudied populations are simpler than studied ones (White 2000), and massasauga population declines are well documented (Szymanski 1998). Given these well-documented declines, it is prudent to proceed as though the population were declining and to interpret my results as having ruled out the obvious causes (e.g., long generation time, high juvenile mortality). Other causes for decline could include disease, inbreeding depression, or changes in survival due to changing habitat or anthropogenic effects. My results indicate that massasaugas at my study area may be stable for the next 50 years if true demographic parameters do not vary more than were modeled (Table 2.3). The population would still be subject to potentially devastating stochastic events, and continued monitoring and management are needed. The landscape surrounding the study area is subject to increased inflastructural and economic development over the next 50 years, and it is imperative that the property remain a stable source population. This model can be used to assess relative risk and should stimulate discussion among biologists, land use planners, and private citizens about how to stabilize or enhance massasauga populations. Future work should focus on addressing data gaps described (e.g., neonate and subadult mortality, existence of inbreeding effects), evaluating effects of various habitat management scenarios on persistence, and documenting spatial or temporal variation in demographic parameters. 42 ACKNOWLEDGMENTS I thank M. Foster, K. Kucher, W. Anderson, S. Zimmer, and numerous volunteers for assistance with snake surveying and monitoring. I thank K. Bissell for the use of her data subset. Potter Park Zoo (Lansing, MI) and Pierce Cedar Creek Institute (Hastings, MI) provided financial and logistical support. I also thank various landowners in the study area for allowing access to their property and S. Winterstein for analytical support. I thank Michigan State University and the Michigan Department of Natural Resources and Environment (Wildlife Division) for providing firnding and permission to Obtain animals for this project. 43 CHAPTER 3: Linking resource selection to habitat management for the eastern massasauga rattlesnake ABSTRACT The eastern massasauga rattlesnake (Sistrurus catenatus catenatus; hereafter, massasauga) has experienced population declines throughout its range and is now a candidate for federal listing as threatened or endangered. However, little is known about massasauga habitat selection in Michigan, particularly in an actively managed landscape. The objectives of my study were to: 1) quantify whether massasaugas in southwestern lower Michigan are selecting certain vegetation types disproportionately to their availability, and 2) quantify the extent to which habitat manipulation influences snake resource selection. I obtained location data for 51 massasaugas telemetered in 2004, 2005, 2008, and 2009. I quantified second-order resource selection using compositional analysis and modeled the ability of 4 vegetation variables to classify habitat management practices using logistic regression. All snakes selected cover types disproportionately to their availability (P = 0.001), which did not differ by reproductive condition or year. A logistic regression model did not successfully classify manipulated and non-manipulated stands (P = 0.16), but snakes did use stands with habitat manipulation at least as much as those without. Habitat management for eastern massasaugas in southwestern Michigan should provide early to mid-successional conditions in deciduous upland and wetland vegetation types. My results can be used by biologists to guide the timing of prescribed fires in occupied massasauga habitat as well as assess the intensity of woody species removal needed to minimize habitat degradation for the massasauga. INTRODUCTION The eastern massasauga rattlesnake (Sistrurus catenatus catenatus; hereafter, massasauga) has experienced population declines throughout its range (Szyrnanski 1998) and is now a candidate for federal listing as threatened or endangered (U SFWS 1999). Ranging through the Great Lakes region of the Midwest and southern Ontario, the massasauga is now considered threatened or endangered in every state or province in which it occurs, except Michigan, where it is a species of Special Concern. There is a deficit of information on the massasauga’s habitat requirements in Michigan, the state considered to be the species’ stronghold (Szymanski 1998). Michigan has lost 33% of historical populations (Johnson et al. 2000), with the majority of those remaining confined to public land or nature preserves (Szymanski 1998). Reasons for the decline include habitat loss and degradation, urban encroachment, and human persecution (Wright 1941, Szymanski 1998, Johnson et al. 2000). Managers responsible for imperiled species face special challenges when dealing with small, disjunct populations because of limited time and vulnerability of small populations to extinction (Mills et al. 2005). To recover species such as the massasauga, an understanding of their habitat requirements and responses to management prescriptions is needed. Massasaugas generally use open wetland and upland vegetation types left by glacial retreat (Cook 1993, Johnson et al. 2000). Resource selection varies geographically, but several generalizations can be made. Massasaugas select stands that have a relatively open canopy compared to surrounding patches, a water table near the surface, and favorable juxtaposition of wetland and upland areas for seasonal use (Szymanski 1998, Johnson et al. 2000). Biologists managing land for massasaugas are, 45 therefore, flequently concerned with maintaining early seral stages (Johnson and Leopold 1998, Durbian 2006, Moore and Gillingham 2006) while avoiding harm to individual massasaugas incidentally (Seigel 1986, Durbian 2006). Vegetation management techniques proposed for the benefit of massasaugas include mowing, prescribed burning, herbicide application (Johnson and Leopold 1998, Durbian 2006), and browsing (Johnson and Leopold 1998). If early successional uplands and wetlands are limiting, then carefully-applied treatments (e. g., mowing, burning, woody species removal) should positively impact massasauga populations by creating or improving habitat (Johnson et al. 2000). Management activities can affect individual snakes (Seigel 1986, Durbian 2006, Moore and Gillingham 2006), with unknown impacts to populations. Managers will, therefore, need to use approaches that target site-specific threats to massasaugas, minimize accidental mortality, are implemented in locally-selected vegetation types, and are compatible with other management objectives for the property. My objectives were to: l) quantify whether massasaugas in southwestern lower Michigan are selecting certain vegetation types disproportionately to their availability, and 2) quantify the extent to which habitat manipulation influences snake resource selection. Linking population responses to habitat management over time has rarely been attempted for massasaugas (but see Johnson and Leopold 1998), and I provide a perspective flom the central portion of the species’ range. The results of my study can be used by biologists to guide habitat management activities and aid in efforts to minimize habitat degradation. 46 STUDY AREA I conducted this study on a contiguous tract of private property in Barry County, southwestern lower Michigan. Major plant community types in this 267 ha area include southern wet meadow, relict conifer swamp, mesic southern forest, dry southern forest, dry-mesic southern forest, mixed wetlands, and prairie fen (Slaughter and Skean 2003). Predominate soils included Penington and Marlette loam and Oshtemo sandy loam with some pockets of muck (Natural Resources Conservation Service 2000). Average summer and winter temperatures in nearby Hastings, Michigan (flom 1971—2000) were 20.4° C and —4.0° C, respectively (Michigan Department of Agriculture 2004). Average annual precipitation for the same time period was 90.9 cm (Michigan Department of Agriculture 2004). Southwestern lower Michigan has mild humidity, with an average annual relative humidity of 62% (Michigan Department of Agriculture 2004). The terrain of central Barry County is rolling to hilly with many glacially- forrned lakes (Michigan Department of Agriculture 2004). Property managers have used a combination of prescribed fire, prairie plantings, chemical control (Table 3.1), and mechanical control to promote early successional upland and wetland vegetation types and discourage invasive or encroaching species (e.g., autumn olive, Elaeagnus umbellata; multiflora rose, Rosa multiflora) concurrent with this study (J. Howell, property land steward, personal communication). Approximately 29.9 ha of fallow agricultural fields have been converted into tall grass prairie since 1999. Prairie and old field areas totaling approximately 30.6 ha were burned under a rotation flom 2003—2009 typically in March (excluding 2005), much of it more than once. Burned areas ranged flom about 0.27 ha to 11.57 ha in size. Tree 47 Table 3.1. Herbicides applied flom 2006—2009 for vegetation control at the study area in Barry County, Michigan, including species treated and year. Herbicide Active Species treated name ingredient Year Autumn olive (Elaeagnus umbellata) Crossbow triclopyr 2006 Rodeo glyphosate 2007, 2009 Pathfinder II triclopyr 2007—2009 Black locust (Robinia pseudoacacz'a) Pathfinder II triclopyr 2007—2009 Black walnut (Juglans nigra) Pathfinder II triclopyr 2007 Box elder (Acer negundo) Pathfinder II triclopyr 2007 Brome grass (Bromus sp.) Razor Pro glyphosate 2008 Common reed (Phragmites australis) Rodeo glyphosate 2006—2007 Dewberry (Rubus Sp.) Transline clopyralid 2009 Garlic mustard (A lliaria petiolata) Razor Pro glyphosate 2007 Misc. grasses Razor Pro glyphosate 2008 Misc. shrubs Pathfinder II triclopyr 2008 Mulberry (Morus rubra) Pathfinder II triclopyr 2007—2009 Multiflora rose (Rosa multiflora) Pathfinder II triclopyr 2007—2008 Poison sumac (Toxicodendron vernix) Rodeo glyphosate 2009 Purple loosestrife (Lythrum salicaria) Rodeo glyphosate 2006 Spotted knapweed (Centaurea stoebe) Transline clopyralid 2007—2009 48 and shrub removal were also conducted where early seral stages are desired. Invasive and undesirable native woody species were targeted with herbicides (Table 3.1) using a cut-stump herbicide application with an herbicide wand to limit runoff, and the invasive reed grass, Phragmites australis, is treated with a wipe-on application Of glyphosate to the stem or foliar spray (J. Howell, personal communication). Several herbaceous species (e.g., garlic mustard, Alliaria petiolata; spotted knapweed, Centaurea stoebe) are hand-pulled or sprayed with herbicide at select locations (J. Howell, personal communication). METHODS Capture and Radiotelemetry I used visual surveys to locate massasaugas in the spring and summer of 2008 and 2009. I transported adult snakes weighing >100 g to Potter Park Zoo (Lansing, M1) for surgical implantation of radio transmitters by the zoo veterinary staff. Before surgery I sexed snakes and used sonograms and radiographs to distinguish reproductive condition of females. Surgery generally followed the procedures described by Weatherhead and Anderka (1984), modified for this study. Modifications included anesthesia with Isoflurane, administration of pain relief and antibiotics post-operatively, introduction of a sterile catheter of 2.7 or 3.3 mm diameter (8 or 10 French, respectively) rather than a copper tube to position antennas subcutaneously, and warming of snakes throughout surgery. Transmitters attached were either 0.8 g, 3.1 g, 8 g, 10 g, (Models A2410, R1680, and R1520 respectively, Advanced Telemetry Systems, Inc., Isanti, MN) or 5 g (Model SB-2, Holohil Systems Ltd., ON, Canada). I also tagged snakes with a Passive Integrated Transponder (PIT tag; Avid Identification Systems, Norco, CA) for 49 permanent individual identification; transmitters and PIT tags weighed <5% of the individuals’ body weight to minimize marking effects on snakes (Parent and Weatherhead 2000, Withey et al. 2001 , Fuller et al. 2005). After a 3—5 day recovery period, I released snakes at their capture site and relocated them in the field approximately 6 days/week flom May—August and at least once a week thereafter until hibernation (mid to late Oct). I also included data collected under similar protocol in 2004 and 2005 for a previous study (Bissell 2006). Michigan State University’s Institutional Animal Care and Use Committee approved all capture, holding, and surgical procedures for this study (permits 03/04-040-00 and 05/08-067- 00), which conformed to regulations under the Michigan Department of Natural Resources Scientific Collector’s Permit (dates issued: 3 Dec 2003, 22 Mar 2005, 24 Mar 2006, 20 Nov 2007, and 12 Feb 2009). Seaman et al. (1999) generally recommend a minimum of 30 relocations to construct unbiased kernel home range estimates; in order to determine how many relocations are needed for unbiased estimation of massasauga home ranges, I constructed observation-area curves (Odum and Kuenzler 1955). After data collection ended, I selected a stratified random subset of snakes such that each year and each reproductive group (males, nongravid females, and gravid females) were represented because male massasauga home ranges are often larger than those of nongravid females, which in turn are larger than those of gravid females (W eatherhead and Prior 1992, Marshall et al. 2006), but see Reinert and Kodrich (1982) for an exception. I used ArcView 3.2 and Animal Movement Extension (P. N. Hooge, W. Eichenlaub, and E. Solomon, US. Geological Survey, Glacier Bay, AK) to construct 95% fixed kernel 50 home ranges by sequentially adding relocations in groups of 10 (or 5 if total relocations was <30) until all were included, using least squares cross-validation (LSCV) to select bandwidth (Seaman et al. 1999). I examined the points at which the curves became asymptotic for each reproductive group, flom which I developed a lower threshold of required numbers of relocations for each group. Massasaugas are viviparous with young being born flom August through September (Szymanski 1998), and I considered reproductive females gravid for an entire active season, including the postparturition and ingress period. Due to violations of parametric assumptions, I used the Kruskal-Wallis (1952) procedure to test whether home ranges were different by year and reproductive group. I used SAS 9.1 software (SAS Institute, Cary, NC) with a = 0.05 for home range comparisons between reproductive groups. Resource Selection I split relocation data into 2 groups for all further analyses, 2004—2005 and 2008-2009 for temporal replication. I tested whether snakes at the study area selected cover types in proportion to their availability. For snakes meeting relocation thresholds as determined by observation-area curves, I used individual home ranges (95% fixed kernels) to quantify use. I expected eastern massasaugas to select early and mid- successional wetlands and uplands disproportionately to their availability (Bissell 2006). To classify availability of resources, I used integrated forest management analysis program (IF MAP) coverages for my study area (MDNR 2003), reclassified based on structure, in ArcMap 9.2. I classified cover types as wetland (W) or upland (U), coniferous (C) or deciduous (D), and qualified them as either early and mid- 51 successional (EM) or late-successional (L); all else was bare ground (DEVL) (Table 3.2). The resulting 7 classes (UDEM, WDEM, UDL, UCL, WDL, WCL, and DEVL) were comparable to Bissell’s classification system (2006). The IF MAP data has a resolution of 30X30 m, and overall accuracy of the original IFMAP data generally exceeded 80% (Space Imaging 2004). Availability was defined at the population level, with resource use identified for each individual (i.e., Design 11 study; Erickson et al. 2001) I defined available resources according to 2 possible boundaries of Design 11 availability to assess any boundary effects caused by arbitrary definition of available resources (Aebischer et al. 1993, Erickson et al. 2001). The first. and smaller boundary of availability is a minimum convex polygon (MCP) of all known snake locations (excluding open water). The second, and larger, availability measure is an aggregation of stands known and suspected to be used by snakes. For this, I combined the MCP stands, adjacent stands with similar vegetation types, and stands where snakes had been reported to me (excluding open water, buildings, and parking areas). This second scale is akin to a study area definition of availability, though slightly smaller than the total property area. I hypothesized that selection order would not differ according to definition of availability. Because I could not spatially replicate the study, site effects will be confounded with selection results. Assumptions made when estimating resource selection included: marked animals are a random sample, telemetry locations are temporally independent, independence of animal resource use, and correct classification of used resources (Erickson et al. 2001). I used compositional analysis (Aebischer et a1. 1993) to quantify 52 Table 3.2. Integrated forest management analysis prograrn/ gap analysis program (IFMAP; MDNR 2003) land cover categories found at the study area in Barry County, Michigan, reclassified according to floristic and structural characteristics. Cover types were considered wetland (W) or upland (U), coniferous (C) or deciduous (D), and qualified as either early to mid-successional (EM) or late-successional (L). All others were classified as DEVL. Seral Code stage Description DEVL NA unpaved road; other bare, sparsely vegetated UD EM row crops; forage crops, non-tilled herbaceous; herbaceous openland; upland shrub, low-density trees UD L northern hardwood association; oak association; aspen association; other upland deciduous; mixed upland deciduous UC L other upland conifers; upland mixed forest WD L lowland deciduous forest WC L lowland coniferous forest; lowland mixed forest WD EM floating aquatic; lowland shrub; emergent wetland; mixed non-forest wetland 53 second-order selection (Johnson 1980) of cover types by massasaugas. I substituted null use values with 0.01% consistent with Aebischer et al.’s (1993) suggestion. I used PROC GLM to conduct a multivariate analysis of variance (MANOVA) in SAS 9.1, with the dependent variables equal to the differences in log-ratio transformed used and available components (Aebischer et al. 1993). Due to lack of multivariate normality, I used randomization to generate p-values for hypothesis testing; only randomized p-values are reported. I tested for effects of boundary (study area and MCP), reproductive group (male, nongravid female, and gravid female) and manipulation cohort (2004—2005 and 2008—2009). I used the BYCOMP.SAS program (P. Ott and F. Hovey, British Columbia Forest Service, BC, Canada) to test the omnibus hypothesis (H0: habitat use is proportional to availability), test main effects, and rank habitat components using randomization (or = 0.05). Habitat Management Habitat management interventions were limited in 2003 and 2004 (~10.7 ha burned, no shrub removal) so I only considered snakes telemetered in’2008 and 2009 to have responded to manipulations. I classified cover types at the study area as having one of 3 mutually exclusive manipulation regimes: 1) prescribed burning plus shrub removal, B+SR, 2) shrub removal only, SR, and 3) neither burning nor shrub removal-no management, NM. I excluded certain practices flom my analysis of management areas because scale of implementation was too small (hand pulling of weeds and mowing) or geographically restricted (native grass and forb plantings). For classification purposes, P. australis is included as a shrub removed due to its exceptional height and stem 54 density, forming thick colonies in Open wetlands. I quantified the percent of land in each management regime at the study area. Managers at the study area tended to select stands for burning or shrub and tree removal based on their structure and the objective to convert fallow agricultural land to tallgrass prairie, confounding a more straightforward test of selection of managed units with cover type because managed units were not randomly assigned. Instead, I examined home range data to determine the proportion of the 2008—2009 snakes’ home ranges that fell within each category of management intervention by year. If snakes were not using managed stands, I would expect the mean proportion Of home ranges within manipulated units to be less than that within non-manipulated units. I also analyzed stand-level vegetation data using logistic regression with SAS 9.1. Four vegetation variables thought to be important to massasaugas at the study area (Bissell 2006) were collected biweekly-weekly and averaged for each snake (Table 3.3). Each vegetation survey could have been flom either a managed (B+SR and SR) stand or not (NM), but I only used surveys flom similar cover types (early and mid-successional stages) to simplify the model. For each snake, I randomly assigned vegetation surveys flOm either managed or unmanaged stands to the study so that all samples were independent and approximately evenly distributed. Vegetation variables included percent live herbaceous cover (LHC), percent dead herbaceous cover (DHC), stem density of trees and shrubs >3 rn (SDTS), and absolute dominance of trees >3 m (ADT). I measured LHC and DHC using 2 20-m line intercepts (Canfield 1941) placed 5 m away flom snake locations. I quantified SDTS using 5X20 m plots adjacent to line intercepts; ADT was the average diameter at breast height (DBH) multiplied by the 55 Table 3.3. Stand-level vegetation variables collected at eastern massasauga rattlesnake (Sistrurus catenatus catenatus) locations in Barry County, Michigan flom 2008 to 2009. Variable Description Method LHC % live herbaceous cover 20-m line intercept DHC % dead herbaceous cover 20-m line intercept SDTS stem density of trees & shrubs >3 m tall 5X20 m plot, scaled by 100 ADT absolute dominance of trees >3 m tall mean DBH X number of trees/ha 56 number of trees/ha. I conducted a Hosmer and Lemeshow Goodness-of-Fit Test to evaluate whether the model fit the data well (a = 0.05). RESULTS From 2004—2009, 51 adult eastern massasaugas were radiotagged and monitored by K. Bissell and me. The earliest and latest snake release dates were 11 May and 9 August, and the latest tracking date was 29 October. I excluded all data for 3 snakes because of insufficient Observations (e.g., poor health, few observations). I had 2 years of data for 8 snakes, which I considered independent when they emerged flom hibernation (generally about 7 months) because they may have had a different reproductive status. Observation-area curves (n = 19) revealed that the minimum number of relocations needed to produce asymptotes for home ranges for males, nongravid females, and gravid females are 40, 50, and 40, respectively. On this basis, I excluded 11 home ranges flom analyses, resulting in 45 home ranges flom 39 snakes. Average home range areas were not different between 2004 and 2005 (x12 = 1.33, P = 0.248) or 2008 and 2009 ()02 = 0.136, P = 0.712) so home ranges were pooled into 2 groups, pre- and postrnanipulation, which were different (x12 = 17.74, P < 0.001). Premanipulation home ranges (3C = 1.04 ha, SD = 1.26, n = 22) were smaller than postrnanipulation ones (3': = 5.21 ha, SD = 4.28, n = 23), but more of the premanipulation sample was composed Of females. Pooled by year, male home ranges were larger than those of nongravid females ()(12 = 5.98, P = 0.015), which were larger than those of gravid females (ya2 = 8.63, P = 0.003). Average home range size for males was 7.42 ha (n = 10, SD = 5.05), for nongravid females was 3.15 ha (n = 18, SD = 2.58), and for gravid females was 0.71 ha (n = 17, SD = 0.58) (Figure 3.1). 57 14 . a Pre 12 1 I Post 10 . Pooled Home range size (ha) Reproductive group Figure 3.1. Home range sizes (:1: 2 SE) for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan, by reproductive group: male (M), nongravid female (N GF ), and gravid female (GF). Each reproductive group is summarized 3 ways: premanipulation (Pre, 2004—2005), postrnanipulation (Post, 2008- 2009), and pooled (Pooled, 2004—2009). 58 Resource Selection For resource selection analysis, I eliminated the UCL (Table 3.2) category because only one snake’s composition contained this category and to a very small extent (<0.5%). I combined WDL and WCL due to low overall usage of both cover types (<15% and <1 .4%, respectively). For the remaining 5 categories, I first tested for the existence of a boundary effect (other effects included), and finding none (Wilk’s A4, 80 = 0.938, P = 0.731), I conducted subsequent analyses using only reproductive group and manipulation cohort as factors at the study area level. I selected the study area level for subsequent analysis because this area corresponds well with the scale at which habitat management occurs. Reproductive group, manipulation cohort, and their interaction were all non-significant effects in the model (P = 0.154, P = 0.385, and P = 0.741, respectively). I rejected the null hypothesis of resource use proportional to availability at the study area level (Wilk’s A4, 41 = 0.486, P = 0.001, n = 45). A ranking matrix ordered the cover types, flom most used to least used, as WDEM > UDEM > DEVL >WDL-WCL > UDL (Table 3.4). The most abundant cover type was UDL (30.39%) and the least abundant was DEVL (2.77%) (Figure 3.2). The UDL cover type was used significantly less than all other categories (P _<_ 0.016). The WDEM category was used significantly more than the combined WDL-WCL category (P = 0.001). Although the UDEM class includes row crops and forage crops, no areas used by snakes were agricultural during this study. Habitat Management By 2008, 12.42% of available cover types were classified as B+SR, which increased to 15.66% by 2009 due to manipulation activities. Cover types classified as SR increased 59 Availability (%) N O ,' --'-". ' ".2 _1 DEVL UDEM UDL WDL-WCL WDEM Cover type Figure 3.2. Percent of cover types available to eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) flom 2004—2009 in Barry County, Michigan, at the study area level. Cover types were classified as wetland or upland (W or U), coniferous or deciduous (C or D), and designated as early-mid, or late successional (EM or L); all else was developed (DEVL). 60 Table 3.4. Ranking matrix for selection of cover types used by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan flom 2004— 2009 (n = 45). Larger ranks indicate higher preference. WDL- Cover typea WDEM UDEM DEVL WCL UDL Rankb WDEM 0.516 0.369 0.001 0.001 4 A UDEM 0.516 0.504 0.337 0.001 3 AB DEVL 0.369 0.504 0.698 0.016 2 ABC WDL-WCL 0.001 0.337 0 .698 0.001 1 BC UDL 0.001 0.001 0.016 0.001 0 D a Cover types classified as wetland or upland (W or U), coniferous or deciduous (C or D), and designated as early-mid, or late successional (EM or L); all else was developed (DEVL). b Ranks followed by the same letter are not significantly different (a = 0.05); smallest possible randomized p-value is 0.001. 61 flom 32.49% in 2008 to 35.18% in 2009. The NM portion of available cover types decreased flom 55.09% to 49.16% between years. On average, snakes (n = 23) had 61.21% (SD = 30.80) of their home ranges in SR areas and 32.24% (SD = 25.68) in NM areas (Figure 3.3). Although the average percent of home ranges in B+SR areas was 6.56% (SD = 17.85). the 95% confidence interval overlaps zero because Of the pattern of use (Figure 3.3). Just 7 snakes used B+SR areas, only 5 of which had >1% of their home range in this type, consisting of 2 gravid females, 2 nongravid females, and 1 male. For these 5 snakes, B+SR areas made up an average of 29.99% (SD = 29.60) of the home range. The same B+SR stands which were used in 2008 were used in 2009, all 1—3 years post burn; no stands >4 years post—bum or <1 year post-bum were used. Gravid female, nongravid female, and male home ranges were similarly composed of SR areas (.77? = 83.52%, SD = 40.62; J? = 66.83%, SD = 28.35; and 3c = 50.45%, SD = 28.34, respectively). The vegetation model included 13 vegetation samples flom managed stands and 10 vegetation samples flom unmanaged stands. A Hosmer and Lemeshow Goodness- of-Fit Test suggested that the model was adequate for my data ()(92 = 16.44, P > 0.05). Examination of Pearson correlation coefficients revealed that SDTS and ADT were positively correlated (P < 0.001), so I removed the variable ADT because it was the only category to have an outlier. My logistic regression model as a whole did not fit better than an empty model (likelihood ratio 39,2 = 5.17, P = 0.16, max-rescaled R2 = 0.27), with no individual variables having P < 0.13. None of the variables were significant predictors of management presence, with or without ADT included. 62 7o — 8 so - (=8 50 ~ 5 o 40 - E 30 .. I 20 . .\° 10 « I o W l r - B+SR SR Management type Figure 3.3. Mean percent of eastern massasauga rattlesnake (Sistrurus catenatus catenatus) home ranges (:t 2 SE) occurring in each management type available at the study area in Barry County, Michigan, flom 2004—2009: burning with shrub removal (B+SR), shrub removal only (SR), no shrub removal or burning (NM). 63 DISCUSSION The pattern of home range size and reproductive status evident in my study (males > nongravid females > gravid females) differs flom that reported by Marshall et a1. (2006) in Indiana (n = 26) and Johnson (1995) in New York (n = 14), who found males and nongravid females to differ flom gravid females but not flom each other. Weatherhead and Prior (1992) also found sex differences in home range size between males and females in Ontario (n = 12) but did not subdivide females. Reinert and Kodrich (1982) did not find differences in home range area in western Pennsylvania between sexes, but gravid females did have significantly shorter range lengths than nongravid females by 56 m (range length was defined as the distance between the 2 most distant location sites within an activity range). Sexual differences in home range size could be attributed to site effects, home range estimators, or the minimum number Of relocations used to estimate home range (Kemohan et al. 2001). Researchers should take into account how kernel home range estimators perform when fewer than 30 relocations per snake are available because it could affect other investigations in which the home range is interpreted as use, as in my study. Boundary effects on resource selection were not apparent in this study, likely because the 2 definitions chosen were not different enough to impact selection order. Neither was resource selection at the study area level shown to depend on reproductive group. For massasaugas in northeastern Indiana, Marshall et al. (2006) found that reproductive group had no significant effect on second-order resource selection when use was estimated with an MCP but became a significant effect when use was estimated with a 95% fixed kernel home range; Marshall et al. (2006) found gravid females (12 = 64 8) to differ flom nongravid females and males at this level, which my results do not corroborate even though my sample had over twice the number of gravid females (n = 17). I did not estimate use with an MCP because the fixed kernel home range is usually considered a more robust home range estimator than the MCP (Kemohan et al. 2001). I also found no effect of manipulation cohort, suggesting no temporal divergence in resource selection relative to habitat management. My results substantiate those of an earlier resource selection study at this area (Bissell 2006), despite somewhat differing definitions of availability and use due to different study designs. Bissell’s (2006) resource selection function showed WDEM and UDEM cover types were selected more than expected, while WCL and most types of UDL were selected less than expected. This is not surprising considering that the premanipulation cohort, which was half of my sample, included all but 4 of these same snakes. Several discrepancies did occur, however; my analysis documented that DEVL ranked third in overall selection and was not significantly different flom the top-ranking cover type (WDEM), whereas Bissell found it to be used less than expected (2006). I also found WDL-WCL to be used significantly less than the top ranking cover type ’ (ranked next to last), whereas Bissell found WDL to be used more than expected. Whether these discrepancies exist because of increased sample size, a shift flom a Design Ito a Design 11 study design, or changed resource selection, is unknown. My results show that early and mid-successional wetlands and uplands are equally important in southwestern Michigan, followed by roads and other bare ground features, all of which are preferred to heavily forested cover. 65 Massasaugas at my study area appear to select cover types similarly to snakes elsewhere, with a few differences. Marshall et al. (2006) found open wetland areas (e.g., shoreline, cattail marshes, sedge 'tussocks) to be preferentially selected over upland cover types (e.g., agriculture, old field) in northeastern Indiana; however, in my study, Open upland areas were also highly selected along with open wetlands. Harvey and Weatherhead (2006b) report preferential selection of open areas (e.g., abandoned field, crop, outcrops) over wetland areas (e. g., open bog, marsh, coniferous swamp) in Ontario, although the difference was not significant. Weatherhead and Prior (1992) reported a very minimal association of snake locations with open upland areas in Ontario; Moore and Gillingham (2006) also reported significantly higher disproportionate use of wetland cover types over upland areas (e.g., agriculture, bare ground, upland) in southeastern Michigan. Whether snakes prefer wetland or upland substrates varies by site and likely also by the relative availability and the quality of habitat (i.e., composition, structure) provided by the respective vegetation types in these areas. Range-wide discrepancies could also arise for reasons unrelated to habitat use, such as differing estimators, study designs, and sample attributes. I found that coniferous forests were seldom selected despite being readily available, contrary to their preferential selection in Ontario as reported by Weatherhead and Prior (1992) and Harvey and Weatherhead (2006b). Weatherhead and Prior (1992) stated that use of coniferous forest was often, though not always, associated with openings in the forest, which is consistent with my observations of the few snakes that used it. Some authors have found temporal shifts in habitat selection throughout an activity season (Reinert and Kodrich 1982, Seigel 1986, Bissell 2006, Harvey and 66 Weatherhead 2006b) while others have not (Weatherhead and Prior 1992). More research is needed on localized within-season variation in resource selection in order to coordinate the timing of management activities Of highly selected vegetation types to minimally disturb snakes. To link resource selection with habitat management, I investigated relationships between structural variables and management practices (e.g., shrub removal, no manipulation). None of the 4 vegetation variables I measured were sufficient to discriminate between managed and unmanaged stands. Structural variables similar to the ones I used have demonstrated importance for massasaugas (Johnson 1995) and effectively characterized differences in manipulated plots (Johnson and Leopold 1998), indicating that my variables were appropriate. Johnson (1995) found that massasaugas in New York, regardless of reproductive condition, were more likely to be found in areas with lower Shrub stem density, shrub height, shrub canopy coverage, and tree canopy closure. Johnson and Leopold (1998) were able to detect differences between management plots using similar structural variables (e. g., stem density, shrub height, basal area). One explanation of this result is that management activities at my study area have not been extensive enough to affect the variables that I measured throughout the manipulated stands. The number Of cut stems ranged flom 100—600 stems/ha in SR and B+SR areas (R. Bailey, Michigan State University, unpublished data). Alternatively, the variability associated with vegetation characteristics could have been too high to detect existing differences. Using averages for each of my vegetation variables (snake as sampling unit) rather than individual vegetation surveys (pooled presences as sampling unit) decreased the sample size of my regression analysis 67 (flom 304 to 23) but avoided using replicates flom a single sampling unit (Hurlbert 1984). I opted for this approach because I did not consider weekly-biweekly vegetation surveys for snakes spatially independent. Considering that snakes had on average 61 .21% of their home range in SR areas, I conclude that SR areas are at least as acceptable as their NM counterparts within the same cover types. Most snakes did not use B+SR areas, but the 5 that did had almost 30% of their home range in this treatment (at most 12—16% of what I considered available) suggesting disproportionate use by a subset of snakes. Johnson (1995) reported that 10.1% of above ground relocations for 9 massasaugas occurred inside a shrub removal treatment that constituted 2.5% of a wildfire burn area in which it was nested, flom which he concludes disproportionate use of burn plots with treated shrubs. Given that most home ranges in the 2008—2009 cohort (including those which were dropped flom the study due to insufficient relocations) used some combination of B+SR, SR, and NM treatments and few used any treatment exclusively, I suggest that habitat quality at the study area is generally sufficient in manipulated and unmanipulated areas. My inability to quantify the absolute number of snakes in stands treated complicates my ability to determine their responses to manipulated areas (Johnson and Leopold 1998). Because shrubs were removed flom the study area with a cut-stump herbicide application instead of broadcast application and all burns were conducted before April 15 (J. Howell, personal communication), I neither anticipated nor Observed negative impacts to snakes. Johnson and Leopold (1998) Observed significantly higher shrub mortality with broadcast-only herbicide application than with cut-wick application, and 68 these authors also did not explicitly observe any negative impacts to snakes. Cut-wick application does allow greater control than broadcast spraying, but it is more labor- intensive (Johnson and Leopold 1998). Implements such as brush hogs and wheeled diskers, although they can be effective at shrub removal, should only be used with care by persons familiar with snake activity patterns (Johnson et al. 2000). Burn regimes can affect herpetofaunal abundances through direct mortality flom fire, alteration of microhabitat structure, and alteration in abundance of prey items (W ilgers and Home 2006). Durbian (2006) found that pre-bum mowing (late June) and prescribed fire (early July), designed to increase massasauga habitat quality in Missouri, resulted in direct and indirect mortality (e.g., avian predation) of snakes. I concur with Seigel (1986) that negative impacts to massasauga populations flom habitat management activities may be mitigated by conducting them in an appropriate season, for example in the late fall after ingress occurs. Compelling evidence for population declines due primarily to habitat management activities is still lacking; in fact, a replicated experimental study by Jones et al. (2000) found that snakes in Kansas were most abundant in pastures treated with an herbicide and fire, compared to herbicide- only and control treatments. Management plans should, however, be periodically assessed for their effectiveness, and demographic data should be collected concurrently in an adaptive management flarnework. Where populations are particularly small, management and monitoring should begin as soon as possible because likelihood of extinction increases as management options become more and more constrained (Mills et al. 1995). 69 Management Implications Management for eastern massasaugas at my study area and perhaps elsewhere should provide early to mid-successional conditions in deciduous wetland and upland vegetation types. Stands with shrub removal were used extensively by massasaugas at this site, suggesting that SR stands are at least as suitable as their non-manipulated counterparts. Providing structure consistent with burning and shrub removal may attract a subset of snakes which prefer more open conditions. I propose that the maximum SDTS in an NM stand used by massasaugas included in my vegetation model (2,562.50 stems/ha) might be considered an upper threshold for the density of tall woody stems that are tolerated. The fact that burned stands were used 1—3 years post- bum by snakes with different reproductive conditions, including gravid females, has implications for planning the timing and rotation schedule of prescribed burns. Vegetation management implemented in occupied habitat should be monitored for long- terrn effectiveness (i.e., a measurable change in vegetation structure would indicate habitat management, as distinct from maintenance of existing conditions); demographic responses should be monitored concurrently (e.g., population survival, productivity). The considerable effort and resources involved in improving habitat for declining or threatened species justify an adaptive approach to habitat management, one in which responses are considered at the population level. ACKNOWLEDGMENTS I thank K. Kucher, M. Foster, W. Anderson, S. Zimmer, and numerous volunteers for assistance with snake surveying and monitoring. I thank K. Bissell for the use of her data subset. Potter Park ZOO and Pierce Cedar Creek Institute provided financial and 70 logistical support. I also thank various landowners in the study area for allowing access to their property. I am gratefill to D. Paperd and J. Howell for procedural and data assistance. I especially thank W. Wang of the Michigan State University Statistical Counseling Service and S. Winterstein for analytical support on this project. I thank Michigan State University and the Michigan Department of Natural Resources for funding this project. 71 CHAPTER 4: Testing a habitat suitability index model for the eastern massasauga rattlesnake ABSTRACT The eastern massasauga rattlesnake (Sistrurus catenatus catenatus; hereafter, massasauga) has experienced population declines throughout its range and is now a candidate for federal listing as threatened or endangered. A habitat suitability index (HSI) model was developed in 2006 to help managers evaluate habitat quality for massasaugas in southwestern Michigan, but this model has not been evaluated. My objectives were to: 1) quantify the appropriateness of the model variables, 2) quantify the relationship between habitat suitability and massasauga fitness, 3) determine an appropriate scale for model application, and 4) evaluate the sensitivity of the HSI model to variability in input values. I implanted radiotransmitters in 25 snakes flom May— August in 2008 and 2009 and collected demographic and habitat data. I used 4 performance measures to evaluate the HSI model: differences in vegetation of snake- selected and random sites, relationship of HSI scores and individual SI values to fitness surrogates, distribution of HSI scores, and sensitivity of model variables to changes in input data. The production functions of 5 variables did not sufficiently describe massasauga responses to habitat conditions, 4 of which I modified to better describe suitability. I suggest that the modified model be applied to 2-ha areas to assess relative, rather than absolute, habitat quality. Future applications Of the model would exploit the broader range Of snake responses suggested by the proposed adjustments. 72 INTRODUCTION The term habitat quality refers to the capacity of a species’ environment to supply conditions which sustain individual persistence and population growth (Hall et al. 1997, Garshellis 2000). Habitat suitability index (HSI) models are single-species or guild models originally developed by the US. Fish and Wildlife Service (U SFWS) to quantify habitat quality (USFWS 1981); the use of HSI models to assess habitat quality is widespread (Brooks 1997, Pauly and Christensen 2006). In general, HSI models and other rapid assessment indices are among the few methods that can be effective and inexpensive for managers to quantify habitat potential (Presser and Brooks 1998). Lower vertebrates (e. g., snakes, amphibians, fish) are relatively underrepresented in the literature pertaining to habitat relationships (Berry 1986), which could be because they are not as easily detected as birds and mammals. Nevertheless, managers responsible for such species need tools for quantifying habitat potential, especially when the species of interest is rare or threatened (Berry 1986). For example, the eastern massasauga rattlesnake (Sistrurus catenatus catenatus; hereafter, massasauga) has experienced population declines throughout its range at least since the 19708 (Szymanski 1998) and is now a candidate for federal listing as threatened or endangered (USFWS 1999). Ranging through the Great Lakes region of the Midwest and southern Ontario, the massasauga is now considered threatened or endangered in every state or province in which it occurs, except Michigan, where it is a species of Special Concern. Limited information is available in the peer-reviewed literature on the massasauga’s habitat requirements in Michigan, the state considered to be the species’ stronghold (Szymanski 1998). I concur with Johnson et al. (2000) that habitat 73 management for massasaugas should be based on specific Objectives, and they suggest that managers look to tested habitat models to monitor the efficacy of management strategies (Johnson et al. 2000). A tested HSI model would, therefore, be a valuable tool for massasauga conservation (Johnson et al. 2000, Bissell 2006). An HSI model for massasaugas was developed in southwestern lower Michigan by Bissell (2006) using 2 years of data flom telemetered snakes (n = 28), and one variable was amended in 2010 (K. Bissell, Michigan Department of Natural Resources and Environment, personal communication), hereafter called the reference model. The reference model (Figure 4.1) was designed for Michigan, where the massasauga primarily uses early successional upland and wetland vegetation types (Bissell 2006; R. Bailey, Michigan State University, unpublished data). Massasaugas throughout their range are thought to be dietary generalists, with factors related to other aspects of habitat quality limiting snakes more than any particular prey species (W eatherhead et al. 2009). Therefore, if therrnoregulation, gestation, and hibernation life requisites are satisfied, foraging opportunities should also be satisfied. Bissell’s (2006) model serves as a preliminary hypothesis of the massasauga-habitat relationship in southwestern Michigan, available for field-testing and potential modification. The adoption of HSI models that have not been tested with independent data has been justifiably criticized (Brooks 1997, Roloff and Kemohan 1999). Habitat models ideally should be evaluated before being implemented to describe model performance and model reliability, and to strengthen model performance where possible (U SF WS 1981, Scharnberger and O’Neil 1986). Model variables should be appropriate for the 74 Figure 4.1. Production functions for the 6 variables in a habitat suitability index (HSI) model (Bissell 2006, amended in 2010) for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in southwestern Michigan. Model variables include A) percent live herbaceous cover, B) percent dead herbaceous cover, C) stem density of trees and shrubs >3 m, D) absolute dominance of trees >3 m, E) percent area in early successional deciduous uplands, and F ) percent area in early successional deciduous wetlands. 75 14 A 0.8 X o ‘O .5 0.6 g 38 .g 0.4 m 0.2 0 7 r r r r 0% 20% 40% 60% 80% 100% F3 a: Suitability index 9 o A O) .0 N % Live herbaceous cover 0 1 0.8 ‘ Suitability index 0.2 0 0.6 ‘ 0.4 ‘ 1 1 0100200300400500600700800 Stem density of trees 8. shrubs >3 m (stems/ha) E 0% % Area in early successional 0.81 0.61 0.4‘ 0.2‘ 0 0% 1 . 0.8J 0.6’ 0.4’ 0.2i 20% 40% 60% 80% 100% % Dead herbaceous cover 0 O 0.8 ‘ 0.6 ‘ 0.4 ‘ 0.2 ‘ 0 20% 40% 60% 80% 100% 0% deciduous upland Figure 4.1 76 V 10 2o 30 40 Absolute dominance of trees >3m (mzlha) F 1 20% 40% 60% 80% 100% % Area in early successional deciduous wetland species, and the relationships between variables should be accurate (USFWS 1981). In general, HSI models assume that an animal’s presence in a particular habitat is positively correlated to fitness, an assumption which should be confirmed (Van Home 1983) Sensitivity analysis is also an important component of model evaluation, particularly if models will be used to evaluate changes in habitat quality after implementing land-management practices (Brooks 1997). Additionally, habitat models should be developed and applied at spatial scales which are biologically relevant to the focal species (Roloff and Kemohan 1999, Kroll and Haufler 2006). Brooks (1997) suggested that priority for model validation studies should be placed on species of special concern, like the massasauga. Because habitat relationships are modeled in Open natural systems, HSI models cannot be “validated” per se (Oreskes et al. 1994, Garshellis 2000); therefore, I use the terms “test” and “evaluate” to describe my review of model performance. Management decisions for snakes of conservation concern, like the massasauga, should be based upon data flom local populations (Moore and Gillingham 2006), and, as part Of a larger study of massasauga ecology, I tested Bissell’s (2006) model for southwestern lower Michigan. My objectives were to: 1) evaluate the appropriateness of the model variables, 2) quantify the relationship between habitat suitability and massasauga fitness, 3) determine an appropriate scale for model application, and 4) evaluate the sensitivity of the HSI model to variability in input values. Model diagnostics should identify strengths and weaknesses and aid managers in deciding 77 whether and how to apply this model when planning habitat management activities for I] the massasauga. STUDY AREA I conducted this study on a contiguous tract of non-profit private property in Barry County, southwestern lower Michigan. Major plant community types in this 267-ha area included southern wet meadow, relict conifer swamp, mesic southern forest, dry southern forest, dry-mesic southern forest, mixed wetlands, and prairie fen (Slaughter and Skean 2003). Predominate soils included Perrington and Marlette loam and Oshtemo sandy loam with some pockets of muck (Natural Resources Conservation Service 2000). Average summer and winter temperatures in nearby Hastings, Michigan (flom 1971—2000) were 20.4° C and —4.0° C, respectively (Michigan Department of Agriculture 2004). Average annual precipitation for the same time period was 90.9 cm (Michigan Department of Agriculture 2004). Southwestern lower Michigan is mildly humid, with an average annual relative humidity of 62% (Michigan Department of Agriculture 2004). The terrain of central Barry County is rolling to hilly with many glacially-forrned lakes (Michigan Department of Agriculture 2004). Property managers have used a combination of prescribed fire, prairie planting, herbicides (e. g., glyphosate), and mechanical control to promote early successional upland and wetland vegetation types and discourage invasive or encroaching species (e. g., autumn olive, Elaeagnus umbellata) concurrent with my study, and to a lesser extent, the model development period (2004 and 2005). 78 METHODS Capture and Radiotelemetry I used visual surveys to locate massasaugas in the spring and summer of 2008 and 2009. I transported adult snakes weighing >100 g to Potter Park Zoo (Lansing, M1) for surgical implantation of radio transmitters by the zoo veterinary staff. Before surgery, I sexed snakes and used sonograms and radiographs to distinguish reproductive condition of females. Surgery generally followed the procedures described by Weatherhead and Anderka (1984), modified for this study. Modifications included anesthesia with Isoflurane, administration of pain relief and antibiotics post-operatively, use of a sterile catheter of 2.7 or 3.3 mm diameter (8 or 10 French, respectively) rather than a copper tube to position antennas subcutaneously, and warrrring of snakes throughout surgery. We also shortened transmitter antennas so that they ended just caudal to the heart of the snake. Transmitters attached were either 0.8 g, 3.1 g, 8 g, 10 g, (Models A2410, R1680, and R1520 respectively, Advanced Telemetry Systems, Inc., Isanti, MN) or 5 g (Model SB-2, Holohil Systems Ltd., ON, Canada). I also tagged snakes with a passive integrated transponder (PIT tag; Avid Identification Systems, Norco, CA) for permanent individual identification; transmitters and PIT tags weighed <5% of the individuals’ body weight to minimize marking effects on snakes (Parent and Weatherhead 2000, Withey et al. 2001, Fuller et al. 2005). After a 3—5 day recovery period, I released snakes at their capture site and relocated them in the field approximately 6 days/week flom May—August and at least once a week thereafter until hibernation (mid to late Oct). Snakes were only relocated during daylight hours (typically between 8:00 and 17:00) because Bissell (2006) found 79 no evidence of temporal differences in massasauga resource selection patterns. Michigan State University’s Institutional Animal Care and Use Committee approved all capture, holding, and surgical procedures (05/08-067-00), which conformed to regulations under the Michigan Department of Natural Resources Scientific Collector’s Permit (dates issued: 20 Nov 2007 and 12 Feb 2009). Vegetation Sampling I measured the 6 variables in the HSI model (Bissell 2006; Figure 4.1) at the study area (Table 4.1): percent live herbaceous vertical cover (LHC), percent dead herbaceous vertical cover (DHC), stem density of trees and shrubs >3 m (SDTS), absolute dominance of trees >3 m (ADT), and area of early successional deciduous upland (AEDU) and early successional deciduous wetland (AEDW) within a spatially-explicit evaluation area. Once a week or whenever a snake moved to a new stand, the 4 vegetation variables (LHC, DHC, SDTS, and ADT) were measured at a snake location. After July 15, I sampled vegetation biweekly following peak plant productivity (Bissell 2006). To measure LHC and DHC, 2 20-m line intercept transects were established for each sampling event, one 5 m east and one 5 m west flom the snake location. I positioned intercepts perpendicular to a natural gradient (e.g., soil moisture, elevation), if one was apparent. I measured variables SDTS and ADT using 2 20X5 m plots adjacent to line intercepts. I measured the 4 vegetation variables the same way for 30 random points that were generated in ArcMap 9.2 in 2009; I chose 30 random locations to correspond with the number of snakes in my study. Variables AEDU and AEDW were remotely calculated flom land cover data in ArcMap 9.2. To quantify land cover area, I used integrated forest management analysis 80 Table 4.1. Stand-level vegetation variables collected at eastern massasauga rattlesnake (Sistrurus catenatus catenatus) locations and random locations in Barry County, Michigan flom 2008 to 2009. All variables were collected at snake presence locations, and LHC, DHC, SDTS, and ADT were also collected at random locations. Variable Description Method LHC % live herbaceous cover 20-m line intercept DHC % dead herbaceous cover same as LHC SDTS stem density of trees & shrubs >3 m 5X20 m plot ADT absolute dominance of trees >3 m same as SDTS AEDU area in early successional deciduous uplands 2, 10, and 20-ha . a crrcles AEDW area in early successional deciduous wetlands same as AEDU a Circles were centered over 95% fixed kernel home ranges of snakes. 81 program, IF MAP, data for my study area (Michigan Department of Natural Resources 2003), reclassified based on dominant over story species, hydrology, and successional attributes of vegetation types in the area (Table 4.2) following Bissell (2006). The IFMAP raster has a resolution of 30X30 m, and overall accuracy of the original IF MAP data generally exceeded 80% (Space Imaging 2004). I determined vegetation types (sensu Hall et al. 1997) for vegetation sampling in the field as in Table 4.2, with the exception that early and mid-successional stages were not collapsed. I compared the variable means against production functions in the reference model (Figure 4.1) to determine SI values to the nearest tenth. Habitat quality is depicted by HSI scores on a continuum flom 0.00 (poor) to 1.00 (optimal). The HSI score was calculated flom the equation (Bissell 2006): HSI = [(SILHC + SIDHC + SISDTS + SIADT)/4] x [(SIAEDU + SIAEDW)/2] (1) Model Evaluation I examined 4 measures of model performance. I first tested whether model vegetation variables were appropriate for this species. Following Bissell (2006), I excluded vegetation surveys flom cover types that were used least flequently (<10% of all samples) for both snake presence and random samples. Because an animal’s response to a habitat characteristic may depend on the vegetation type in which the variable is measured (Stauffer and Best 1986), I first tested whether LHC, DHC, SDTS, and ADT were different within each vegetation type for snake presence and random areas using Kruskal-Wallis tests (a = 0.10). Additionally, because Bissell (2006) pooled across cover types to construct production functions, I also similarly tested for differences with samples pooled across all flequently-used vegetation types. For all vegetation 82 Table 4.2. Integrated forest management analysis program/ gap analysis program (IFMAP; Michigan Department of Natural Resources 2003) vegetation type categories found at the study area in Barry County, Michigan, reclassified and quantified according to site-specific characteristics. Vegetation types were considered wetland (W) or upland (U), coniferous (C) or deciduous (D), and qualified as either early to mid-successional (E-M) or late-successional (L). All others were classified as DEVL. Vegetation Code structure Area (ha) Description DEVL NA 9.13 unpaved road; other bare, sparsely vegetated UD E-M 80.95 row crops; forage crops, non-tilled herbaceous; herbaceous openland; upland shrub, low-density trees UD L 77.90 northern hardwood association; oak association; aspen association; other upland deciduous; mixed upland deciduous UC L 6.21 other upland conifers; upland mixed forest WD L 24.11 lowland deciduous forest WC L 10.19 lowland coniferous forest; lowland mixed forest WD E-M 61.86 floating aquatic; lowland shrub; emergent wetland; mixed non-forest wetland 83 comparisons, the variables LHC, DHC, SDTS, and ADT were averaged individually by each snake to minimize any effects associated with individual snakes (e.g., sex, reproductive condition). Bissell (2006) recommended applying the model at a scale of 1—20 ha, preferably at >8.9 ha because that was the area of the 95% fixed kernel utilization distribution for her studied population; for my purposes, the scale at which the model is applied only affects SI values for AEDU and AEDW. I applied the model to 3 concentric circular areas centered over each snake’s home range. I estimated home ranges with a 95% fixed kernel method with bandwidth determined by least squares cross-validation (Seaman et al. 1999) using the Animal Movement extension (P. N. Hooge, W. Eichenlaub, and E. Solomon, US. Geological Survey, Glacier Bay, AK) in ArcMap 9.2. I converted each home range to a single centroid point in ArcMap 9.2, which was buffered to create 3 circles for model application, with total areas of 2 ha, 10 ha, and 20 ha. I chose these areas because 20 ha was approximately the largest area occupied by an individual massasauga during model development, and 10 ha exceeds the minimum recommended 8.9 ha spatial scale (Bissell 2006) and is half of the maximum area. Even though it was smaller than the preferred application scale, I also investigated 2 ha because this was the most flequent home range size (to the nearest integer) for massasaugas in my study (R. Bailey, Michigan State University, unpublished data). I then calculated HSI scores by combining the 4 vegetation variables with the GIS-derived estimates for AEDU and AEDW for each snake’s inhabited area (Equation 1). Density and relative abundance are not necessarily indicative of habitat quality 84 and should not be used to test an HSI model (Van Home 1983, Roloff and Kemohan 1999). To quantify the relationship between habitat productivity and habitat suitability, I derived several indices of fitness based on massasauga biology and available data (Table 4.3) which I presume to coarsely reflect individual fitness. Most female massasaugas in this area typically breed biennially (Bissell 2006), requiring that I separate snakes by reproductive categories (i.e., male, gravid female, and nongravid female) to assess fitness. I calculated a breeding success index for each snake according to its reproductive condition. For gravid females, 1 used the number Of viable offspring (estimated flom a radiograph and 1 or 2 sonograms, or in the case of one female, field observation) as an index of breeding success. I assigned points based on the number of young (Table 4.3), with an additional point given if females were Observed mating (annual breeders). For males and nongravid females, the breeding success index was the number of observed mates. Courting multiple females may benefit male snakes by increasing the chances of their Offspring surviving. Females of many species, including vipers, can also benefit flom having multiple mates (for review, see Shine 2003). Madsen et al. (1992) found that female adders (Vipera berus) that obtained multiple mates had enhanced offspring viability, presumably via sperm competition within the reproductive tract. Rattlesnakes in my study always mated with untagged snakes; therefore, mates were conservatively judged to be independent by considering the distances moved (in meters) and elapsed time (in days) between successive mating events, always with a period of solitude separating mating events (Table 4.3). This method makes no assumptions about how multiple copulations are achieved 85 Table 4.3. Eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) that were radiotracked in Barry County, Michigan (2008 and 2009) were scored according to 3 fitness indices: breeding success, survival flom release to hibernation, and home range size. Adding individual index scores yields the overall fitness score for any snake (possible scores range flom 2—9). Breeding success and home range size point values were assigned according to whether a snake was within or beyond the mean i ISD for that value. Score Index 1 2 3 Breeding success # matesa (males, nongravid females) 1 2 3 or 4 # young (gravid females) 1 or 2 3—8 29 Home range size (ha) gravid females 21.87 0.66—1.86 £0.65 nongravid females 26.55 1.98—6.54 £1.97 males 212.70 307—1269 $3.06 Survived season No Yes 3 Females who mated the same year they produced young (annual breeders) received one extra point. 86 (e. g., female selection, mate-searching success) but does assume that having multiple mates increases fitness in massasaugas and that independent mating events and number of young were correctly estimated. Massasaugas were never Observed mating during a 2-year study in southeastern Michigan (Moore and Gillingham 2006; n = 16), indicating that observable numbers of copulations cannot be expected everywhere. If movements are energetically costly and the extent of a snake’s movements depends on resource availability (Gregory et al. 1987), then smaller home ranges can be assumed to indicate good juxtaposition of necessary resources and lower energy expenditure by snakes. Garshellis (2000) states that home range size may be an appropriate surrogate for fitness and can contribute measurably to habitat quality assessments. I, therefore, used home range size as an additional index of fitness, with smaller home ranges receiving higher scores (Table 4.3). I used 3 broad categories to score home range size because I had 11-89 relocations per snake, which is likely insufficient to accurately describe space use (Seaman et al. 1999). I used reproductive condition to standardize scores because intersexual and reproductive effects influenced home range size (K. Bissell and R. Bailey, unpublished data). I also considered whether a snake survived to hibernation or not to be an important fitness metric (Johnson 1994), for which there is no evidence of a sex effect (K. Bissell and R. Bailey, unpublished data; Diller and Wallace 2002, Brown et al. 2007). For the breeding success and home range size indices, I developed 3 categories flom my data which were assigned point values as previously described: 1) mean :1: ISD, 2) values beyond ISD of the mean, and 3) values below ISD of the mean. I added individual index scores to calculate overall fitness score for each snake (Table 4.3). By 87 visually examining plotted data, I first confirmed that individual indices exhibited the :\ correct directional associations with HSI scores before combining indices into a fitness score. I tested for a linear relationship using linear regression in SAS 9.1 (SAS Institute, Inc., Cary, NC), with a = 0.10. I confirmed the assumption that the residuals were normally distributed by plotting residuals. 1 confirmed the assumption that variances of the errors were constant across the observations by plotting residuals and HSI scores. I repeated this process for a test of the relationship between fitness and individual SI values for all 3 scales. To determine whether any scale was more appropriate, I categorized replicates into 4 bins of HSI scores, poor (0.00—0.25), marginal (0.26—0.50), good (051-075), or optimal (0.76—1.00) following Negri (1995), to assess the distribution of scores. Scales producing wider distributions of scores and stronger relationships with the fitness index were considered more appropriate. I then conducted sensitivity analyses to analyze the effects of variable input data on model output, using only the 2 scales with the strongest relationships between HSI scores and fitness for sensitivity analyses. For each replicate (snake), I changed each input value by 10% in a random direction (holding the other 5 variables constant) and recorded the new HSI score resulting flom the change in the individual SI value. I then averaged the absolute value of the percent change in HSI score for the 30 replicates. I investigated sensitivity by inspecting whether percent change in HSI score was proportional to changes in input variables. 88 7-‘=—n - - ‘WMIMJLJ RESULTS .4 I captured 25 snakes: 16 females and 9 males. I tracked 5 of these snakes both years, and for all analyses, I treated snakes as independent when they emerged flom their hibemacula in 2009. The most flequently used cover types were early and mid- successional deciduous uplands (UDE and UDM) and wetlands (WDE and WDM). There were 374 vegetation transects used to quantify snake locations distributed among flequently used cover types. There were 46 (out of 60) random transects quantifying early and mid-successional vegetation types. Each Of the 4 vegetation variables at snake locations was significantly different . flom random locations in at least one of the commonly used vegetation types (Figures 4.2—4.5). I found differences most Often in early successional vegetation types, as could be expected flom the reference model; however, 2 of the model variables were also significant in mid-successional wetland cover. In WDE cover types (n = 24), LHC and DHC in snake-selected stands and random stands were significantly different ()(12 = 11.135, P < 0.001; x,2 = 10.240, P = 0.001, respectively; Figure 4.2, Figure 4.3); however, only the magnitude of the difference for DHC appears to be meaningful (presence stands had 36.3% less DHC than random stands). In WDM cover types (n = 32), presence sites had 19.3% more DHC than random sites ()(12 = 2.735, P = 0.098; Figure 4.3), and SDTS was lower in snake-selected stands by 589 stems/ha (x12 = 3.346, P = 0.067, Figure 4.4). In UDE cover types (11 = 32), SDT S was significantly lower for presence stands by only 9 stems/ha (x12 = 5.559, P = 0.018, Figure 4.4) and ADT in 89 Presence I Random 100 P< 0.001 P I 0.009 % Live herbaceous cover UDE UDM WDE WDM ALL Cover type Figure 4.2. Percent live herbaceous cover in stands occupied by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) was different than in random stands in southwestern Michigan for one cover type and for all cover types pooled (0. = 0.10). Frequently used cover types in 2008 and 2009 were early and mid-successional deciduous uplands (UDE and UDM, respectively) and early and mid-successional deciduous wetlands (WDE and WDM, respectively). 90 Presence 100* P=0.001 I Random % Dead herbaceous cover UDE UDM WDE WDM ALL Cover type Figure 4.3. Percent dead herbaceous cover in stands occupied by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) was different than in random stands in southwestern Michigan for 2 cover types (a. = 0.10). Frequently used cover types in 2008 and 2009 were early and mid-successional deciduous uplands (UDE and UDM, respectively) and early and mid-successional deciduous wetlands (WDE and WDM, respectively). 91 E Presence 1800 ‘ P ' 0067 I Random >3 m (stems/ha) Stem density of trees 8. shrubs UDE UDM WDE WDM ALL Cover type Figure 4.4. Stem density of trees and shrubs >3 m in stands occupied by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) was different than in random stands in southwestern Michigan for 2 cover types (a = 0.10). Frequently used cover types in 2008 and 2009 were early and mid-successional deciduous uplands (UDE and UDM, respectively) and early and mid-successional deciduous wetlands (WDE and WDM, respectively). 92 Presence I Random 4.3.3“) +0909 4.; meN Absolute dominance of trees >3 m (m2/ha) .b UDE UDM WDE WDM ALL Cover type Figure 4.5. Absolute dominance of trees >3 m in stands occupied by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) was different than in random stands in southwestern Michigan for one cover type and for all cover types pooled (a = 0.10). Frequently used cover types in 2008 and 2009 were early and mid-successional deciduous uplands (UDE and UDM, respectively) and early and mid-successional deciduous wetlands (WDE and WDM, respectively). 93 snake-selected stands was significantly higher than random, although only by 0.37 m2/ha()(12 = 4.441, P = 0.035, Figure 4.5). No variables were significant in the UDM cover type (n = 28). When all cover types were pooled (n = 76), only LHC and ADT were significantly different flom random ()(12 = 6.879, P = 0.009; X12 = 4.736, P = 0.030, respectively); however, the magnitude of differences was only 1.92% for LHC and 3.82 mZ/ha for ADT (Figure 4.2, Figure 4.5). Mean values for AEDU and AEDW (n = 30) were relatively stable (i.e., <5% different) with changes in scale. At the 2 ha scale, the mean for AEDU was 19.01% (SD = 31.15) and for AEDW was 50.78% (SD = 24.98). At the 10 ha scale, the mean for AEDU was 15.32% (SD = 23.02) and for AEDW was 51.08% (SD = 20.85). At the 20 ha scale, the mean for AEDU was 14.14% (SD = 17.59) and for AEDW was 48.39% (SD = 18.38). At the 2 ha scale, HSI scores showed the lowest average (7c = 0.26, SD = 0.19) but exhibited the widest distribution (0.00—0.77). When the model was applied at the 10 ha scale, the average HSI value was 0.28 (SD = 0.13; range = 009—060); the average HSI was highest at the 20 ha scale, but the range was narrowest (If = 0.31, SD = 0.08; range = 0.17-0.47). The fitness scores ranged from 3—7 (out of a possible range of 2—9), with the most common score being 4 (n = 30). Regression analysis revealed that fitness is not linearly related to HSI score at the 2 ha (F1 = 1.56, P = 0.222, R2 = 0.053), 10 ha (F1 = 0.61, P = 0.442, R2 = 0.021), or 20 ha scale (F1 = 0.02, P = 0.902, R2 = 0.001), and no other relationship was discemable flom the plotted data. The fitness-HSI association was weakly positive at all scales. Nor were SI values linearly 94 —'=, I i i 11. related to fitness scores at the 2 ha (F6 = 1.68, P = 0.170, R2 = 0.305), 10 ha (F6 = 0.85, ' P = 0.547, R2 = 0.181), or 20 ha scale (F6 = 0.76 P = 0.610, R2 = 0.165). None of the individual SI values were significantly related to fitness at 10 or 20 ha, but AEDU was a significant variable at the 2 ha scale (t1 = 2.43, P = 0.023). Intriguingly, only 1 or 2 of the SI variables were positively associated with fitness at these scales (AEDU and/or AEDW). Areas used by gravid females tended to have higher mean HSI scores, followed by those used by males and nongravid females, regardless of scale. The HSI scores for the 10 and 20 ha scales were narrowly distributed, with the 2 ha scale achieving the widest distribution (Figure 4.6). No scale achieved a range of HSI scores that correlated significantly with my fitness index, but the 2 ha and 10 ha scales exhibited stronger relationships. Therefore, I selected the 2 ha and 10 ha responses as the best candidates for sensitivity analysis. Model variables differed in their influence on the overall HSI equation output, but model sensitivity was relatively robust to scale of application. At the 2 ha and 10 ha scales, a 10% change in LHC caused roughly proportional changes in HSI scores (11.31% and 11.90%, respectively; Figure 4.7); this variable had the lowest coefficient of variation (CV) of all of the vegetation variables (11.9, pooled by cover type). However, the HSI formula was rather insensitive to 10% changes in DHC, SDTS, and ADT, each of which caused less than a 3% change in HSI scores at the 2 ha scale (AHSI = 2.28%, AHSI = 1.34%, and AHSI = 0.15%, respectively). Results at the 10 ha scale were similar for DHC and SDTS (AHSI = 2.38% and AHSI = 1.48%, respectively) and identical for ADT (Figure 4.7). The variables DHC, SDTS, and ADT had higher variability (CV = 52.6, 124.6, and 111.6, respectively) associated with their means. 95 I2ha I10ha D20ha No. of replicates —\ —i N N O 01 O 01 O 01 Poor Marginal Good Optimal HSI category Figure 4.6. I applied Bissell’s (2006) HSI model for eastern massasauga rattlesnake (Sistrurus catenatus catenatus) at 30 replication sites in the study area (Barry County, Michigan) in 2008 and 2009. Habitat suitability index (HSI) scores were assigned to each of 4 categories (poor: 0.00—0.25, marginal: 0.26—0.50, good: 0.51—0.75, and optimal: 0.76—1.00) to describe their distribution resulting flom model application at.2 ha, 10 ha, and 20 ha. 96 ......... " ’irtzv 11:1 :- M: ,au- ..... % Change in HSI LHC DHC sors ADT AEDU AEDW Model variables Figure 4.7. Percent change in HSI scores resulting flom a 10% change in each variable within Bissell’s (2006) habitat suitability index (HSI) model for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in southwestern Michigan. Sensitivity was tested for each of 6 variables individually across 30 replicates at 2 scales of model application (2 ha and 10 ha). Variables include percent live herbaceous cover (LHC), percent dead herbaceous cover (DHC), stem’density of trees and shrubs >3 m (SDTS), absolute dominance of trees >3 m (ADT), percent area in early successional deciduous uplands (AEDU), and percent area in early successional deciduous wetlands (AEDW). 97 Varying AEDU caused nearly proportional changes in HSI output at the 2 ha and 10 ha scales (AHSI = 7.05% and 7.44%, respectively) and AEDW caused the largest change in model output (AHSI = 14.02% and 11.14%; Figure 4.7). At 2 ha, AEDU and AEDW had relatively low variability (CV = 1.64 and 0.49, respectively), which was very similar at the 10 ha scale (CV = 1.50 and 0.41). Based on these results, I proposed modifications for 4 of the production functions (LHC, DHC, AEDU, and AEDW) using data flom snakes with fitness scores >5 points (11 = 12). I used fitness scores to adjust these maxima because, as Haufler et al. (2002) point out, critical habitat for a species is most likely to occur where the species is most productive rather than where it is most common. The proposed adjusted maxima are 62—100% for LHC, 51.5—96% for DHC, 0—57% for AEDU, and 23—73% for AEDW (Figure 4.8). The variable SDT S was problematic because snakes appeared to be productive even when stern densities exceeded 800 stems/ha; more research is needed to clarify the importance of this variable to massasaugas. I did not modify the production function for ADT, as it appeared to cover a wide range of potential snake responses. DISCUSSION Ideally one would test an HSI model for correlation with demographic data flom several populations in different locations (Van Home 1983, Lancia et al. 1986, Stauffer and Best 1986); however, such detailed demographic data are rarely available for a species of conservation concern. Instead, I attempted to test a proposed HSI model for eastern massasaugas in southwestern Michigan with data flom the same study area in which it was developed, several years later. Garshellis (2000) suggested than an effective design 98 /v rii 1.- 20M; .1? . .3.le Suitability Index .c o s: _o N. 0 0% Suitability index 5:: O\ .O N 0 °°. '0‘. an _o oo .o 05 .o Suitability index 4‘. .o N P °? P "t T 20% T 40% 60% 80% 100% % Live herbaceous cover '1 ----------- 0 C l _o co Suitability index O '05 53 n 0% 96 Area in early successional deciduous upland 20% 40% 60% r 1 0 80% 1009 .c 4:. 0% VT 40% 60% 80% 100% 96 Dead herbaceous cover 0% 20% 1 r-- L N ’ 20% 40% 60% 80% 100% ‘16 Area in early successional deciduous wetland Figure 4.8. Proposed modifications to 4 production functions in Bissell’s (2006) habitat suitability index (HSI) model for eastern massasauga rattlesnakes (Sistrurus catenatus. catenatus) in southwestern Michigan. Dashed lines represent the modified maxima for A) percent live herbaceous cover, 62—100%, B) percent dead herbaceous cover, 51 .5— 96%, C) percent area in early successional deciduous upland cover type, 0—57%, and D) percent area in early successional deciduous wetland cover type, 23—73%. 99 for studies conducted in a single site is to compare the habitat composition of an animal’s home range to the eventual reproduction and survival Of the occupant, similar to my study. One might expect that such intensive study of a species in a small geographic area would produce a model that performs well at that particular site (Stauffer and Best 1986); however, my performance measures only partially supported the reference model for the study area. Roloff and Kemohan (1999) assert that HSI evaluation studies should determine whether model variables have demonstrated importance for the biology Of the modeled species. A114 vegetation variables (LHC, DHC, SDTS, and ADT) were different in snake-selected stands than random stands for at least one of the flequently-used vegetation types. However, according to Oreskes et al. (1994), correspondence between predicted and observed responses merely supports, but does not validate, the appropriateness of these variables. Although LHC and ADT in snake-selected stands were different flom random stands when all vegetation types were pooled, the small magnitude of the differences (1.92% and 3.82 mZ/ha, respectively) suggests that factors other than the vegetation structure I measured could be influencing snake distribution (e. g., predation, foraging success; Morrison 2001). This means that random areas could readily produce high values for SILHC and SIADT. Average values for AEDU fell below the maxima, and values for AEDW were above, for both scales. I did not investigate microhabitat vegetation conditions which could function to regulate snake distribution within stands and perhaps better describe variability in habitat quality (Garshellis 2000). Perhaps microhabitat variables such as distance to nearest shrub, litter depth, or 100 maximum vegetation height (Harvey and Weatherhead 2006b, Moore and Gillingham 2006) could explain additional variation in fitness. According to mean HSI scores, overall habitat quality at the study area would appear to be biased low because no circle centered over a home range area had an HSI value over 0.77, regardless of the scale at which I applied the model. Areas used by gravid females tended to have HSI scores 012—0. 14 units higher than males and nongravid females at the 2 ha scale; this could be related to the over-representation of gravid females in the model development sample (50%). It is possible that Observed HSI ranges were narrow because I did not sample flom the entire range of habitat quality available on the property. My sample size (n = 30) was comparable to that of the model development phase (n = 28); this falls just short of a sample size criterion for modeling recommended by Roloff and Kemohan (1999), determined flom the formula given by Johnson (1981) which is 38 as applied here. Testing habitat-related fitness has been called “extraordinarily difficult” (Garshellis 2000: 144), but is recommended for HSI model evaluation by Van Home (1983) and Roloff and Kemohan (1999), among others. I attempted to test the reference model with surrogates for fitness, recognizing that even a demographic response design has its disadvantages (Garshellis 2000). I did not find a significant relationship between the HSI score of a circular area centered on snake home ranges and the fitness of the snake who occupied that area, although associations were always positive. Pauly and Christensen (2006) state that most habitat models can be expected to explain less than half of the variation in species density or abundance. The reference model explained only 5.3% of the variation in fitness, but even low correlations may provide insights 101 (Pauly and Christensen 2006). For example, hypothesized modifications to the reference model (or competing models) can be tested for a linear relationship with fitness surrogates to objectively compare their explanatory power. The fitness scoring process involved linking 3 fitness indicators with a simple additive relationship, rather than the more complex relationship between density, survival, and expected offspring described by Van Home (1983). I lacked density estimates needed to implement Van Home’s fitness per unit area measurement (1983) due to the elusive nature and cryptic habits of the snake. According to Haufler et al. (2002), animal productivity should be higher in good habitat than in poor habitat when habitat quality is variable. However, Crone (2001) suggests that for most long-lived species, annual adult survival is a better short-terrn fitness surrogate than productivity. High adult survival during the study period (5* = 95%; R. Bailey, unpublished data) necessitated the integration of additional surrogates that perhaps are less associated with habitat quality than survival. Furthermore, it may be infeasible for managers to consider separate models for individual demographic responses during the planning process (Rittenhouse et al. 2010), so a combined fitness index could be desirable flom a planning perspective. Because it is unlikely that a single demographic response would capture the full range of wildlife responses to habitat conditions (Rittenhouse et al. 2010), I assumed that multiple fitness surrogates (e. g., home range size, breeding success, survival) were adequate to link fitness to habitat quality. This assumption could be problematic given that at least 6 years of population data are recommended to demographically validate an HSI model for a biennially breeding species (Marcot et al. 1983, Roloff and Kemohan 1999). 102 Alternatively, the HSI model may not have been a significant predictor of snake \ fitness if habitat conditions of the inactive season (N ov—early May) affected fitness during the active season, which would have caused a decoupling of the HSI-fitness relationship (Rittenhouse et al. 2010). I collected other potentially useful fitness indicators of habitat quality, such as parasite burdens and birth weights of young, but these were unavailable for all snakes and, therefore, could not be used as fitness surrogates for this study. Future studies might assess the utility of these metrics as indicators of habitat quality, as they could elucidate limiting conditions experienced by massasaugas during the inactive season. The reference model showed poor responsiveness to changes in raw values of DHC, SDTS, and ADT, variables which are directly affected by habitat management for massasaugas; therefore, the model may not be well-suited to measuring the impacts of tree and shrub removal or prescribed fire on habitat quality. According to USFWS (1981), averaging functions like the one used here (Equation 1) may become insensitive to extreme values when more than 4 variables are used. Changes in HSI scores were minimally disproportionate when values for AEDW were varied. Variability in input data could potentially mask effects of habitat manipulations for AEDW if the model is applied at 2 ha, but this could probably be prevented with an adequate sampling design. Roloff and Kemohan (1999) suggest that minimum documented home range size may be an appropriate scale for model validation. I did not use the minimum documented home range size for my study (0.277 ha) for 2 reasons: 1) it is below the range recommended by Bissell (2006) for model application, and 2) it is too small for adequate spatial resolution of remotely—sensed variables. Such a small scale does 103 correspond to the minimum area burned on this property (0.27 ha), but the scale at which property managers have conducted vegetation manipulations is larger (maximum = 35.67 ha). Applying the model at 2 ha resulted in the widest range of HSI scores, with the majority of scores falling below 0.50. Based on my criteria, 2 ha would be the most appropriate scale at which to apply the model because it showed more variability in habitat quality, however, I acknowledge that snakes may react to their environment at other scales not studied or multiple scales (Garshellis 2000). General assumptions made during model development, all of which continue to appear valid, include: 1) the study area was within the range of the eastern massasauga, 2) snakes were unobstructed in their use of habitat, and 3) the population of snakes was unharvested. Only one road bisects the property, and though evidence for road-killed snakes exists (Bissell 2006, R. Bailey, personal Observation), there is no evidence that roads obstructed the movements of telemetered snakes to a meaningful extent; firrther, I observed no collecting or intentional killing. The study area was effectively positioned for modeling habitat suitability for massasaugas, given that the property owners’ objectives include conservation of rare species, and that massasaugas on the studied property are insulated flom neighboring land owners. Pauly and Christensen (2006) point out that HSI models are best viewed as testable hypotheses, the value of which lies in the documentation of a repeatable assessment procedure that allows for easier comparison among alternative management scenarios. While the reference model appears to have included appropriate variables and assumptions, I suggest that modified production functions (Figure 4.8) be used to better reflect the relationship between fitness and habitat quality as measured by HSI 104 I”: scores. This could improve the distribution of HSI scores for snakes and perhaps achieve a stronger relationship to fitness; however it may sacrifice some sensitivity. Management Implications Managers can use this model to identify appropriate vegetation types flom a landscape Of interest. At the macrohabitat scale, the production functions presented here can be used to identify ranges of vegetation conditions which appear to be suitable for snake productivity. Future work could identify microhabitat features (e. g., parturition sites) which might also influence fitness. I suggest that the modified model be applied to 2 ha areas to assess relative, rather than absolute, habitat quality. While the suitability relationship of SDTS remains problematic, I suggest that it might be resolved with more data flom different locations. I caution against using the model to predict absence Of massasaugas because their highly cryptic habits make it difficult to predict species absence with much confidence. Acknowledgments I thank W. Anderson, S. Zimmer, and numerous volunteers for assistance with snake surveying and monitoring. Potter Park Zoo staff and docents and Pierce Cedar Creek Institute provided financial and logistical support. I also thank various landowners in the study area for allowing access to their property. I am grateful to D. Paperd for assistance with surgeries. I thank Michigan State University and the Michigan _ Department of Natural Resources and Environment (Wildlife Division) for flmding this project. 105 GENERAL CONCLUSIONS Managers interested in conserving massasaugas can use my results in a variety of ways. For example, when comparing adult survival estimated flom a 2-year study (5’ = 94.72%) to that of a 4-year study (5’ = 79.28%), one might conclude that snake survival can be quite high over a short period, but that longer-term data may be needed to observe temporal variability. Furthermore, the mean number of viable young estimated florn sonograms and radiographs dropped from 9 (Bissell 2006) to 7 when the study duration increased flom 2 to 4 years. These small differences perhaps become meaningful when one considers the incongruity between 2 plausible PVA models, one Of which predicts low probability of population persistence (Bissell 2006) and another which predicts high probability of population persistence (present study). Having 2 additional years of data and more recent literature available to me made my results much more optimistic but still was insufficient to address several persistent data gaps. The importance of long-term population data for population viability analysis is clear in light of these examples; however, one consistent result is the sensitivity of the studied population to juvenile and neonate mortality, estimates of which are sparse and dissimilar for massasaugas. Low sample sizes might well be a continuing problem for massasauga studies given their cryptic and sedentary nature, and combining data flom multiple shorter-term studies within a region may be the best approach for developing reliable estimates. Population-level responses to habitat management were encouraging, yet inconclusive. Early and mid-successional wetlands and uplands are highly selected 106 cover types, whether manipulated or not. Providing these vegetation types in close proximity would minimize the amount of effort needed for snakes to utilize both cover types (i.e., smaller home ranges), which could result in reduced mortality risk. I found no difference in several structural characteristics between stands under different manipulation regimes, and stands with prescribed fire and shrub removal are used by productive snakes. It appears that habitat quality in unmanipulated stands is commensurate with that of manipulated ones. The fact that burned stands were used 1—3 years post-burn by snakes with different reproductive conditions has implications for planning the timing and return interval of prescribed burns. Impacts to massasaugas, and other species coinciding with them (Appendix C), can be minimized by conducting burns every 4 or 5 years in the early spring (late March or early April) before snakes firlly emerge. One unexpected yet alarming result was that areas classified as “developed” in my land cover database (i.e., mowed roadside right-Of-Ways; Table 3.2) were used often enough to rank third in an overall resource selection model. Even though this is likely influenced by the relative scarcity of roadsides at the study area, and no snakes were actually relocated in the road, it remains true that at least some rattlesnakes are killed trying to cross or bask in the road. Clark et al. (in press) document additional adverse road effects on rattlesnakes besides direct mortality (e.g., reduced dispersal and genetic connectivity). Managers can use this information to select appropriate cover types to adaptively manage as massasauga habitat, coordinate timing of management activities, and account for the detrimental effects of roads. This study might have been limited by the resolution of available land cover data; a vegetation map with a resolution finer than 107 N 30X30 m could elucidate snake selection of habitat features not apparent within 900-m parcels. Additionally, researchers can evaluate a modified HSI model’s performance at new locations in southern Michigan. I suggest that additional fitness surrogates (e.g., blood chemistry, parasite burdens, and neonatal birth weights) be investigated to corroborate the HSI model in the firture. Archived massasauga radiolocations (Appendix D) could prove useful in future analyses of geographic and temporal variation in habitat quality. If the modified model can demonstrate a relationship between HSI and fitness surrogates, it would be a valuable tool for managers whose Obj ectives include conserving the massasauga rattlesnake. 108 APPENDIX A. Common plants associated with stands flequently used by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan. Table A.1. Common plant species associated with 3 height strata in early to mid- successional deciduous uplands and wetlands inhabited by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) flom 2008 to 2009 in Barry County, Michigan. Symbol Latin name Common name Stratum 2FERN fern, unidentified fern herbaceous 2GP Poaceae perennial grasses herbaceous ACM12 Achillea millefolium yarrow herbaceous AN AT Angelica atropwpurea purplestem angelica herbaceous BEAL2 Betula alleghaniensis yellow birch tree or sapling CAREX Carex sp. sedge herbaceous CIMA2 Cicuta maculata water hemlock herbaceous CORA6 Camus racemosa gray dogwood shrub COST4 Camus stolonifera red-osier dogwood shrub DAUCU Daucus carota wild carrot herbaceous ELUM Elaeagnus umbellata autumn olive shrub EUMAB Eupatorium maculatum spotted joe-pye-weed herbaceous FRNI F raxinus nigra black ash tree or sapling flCA Impatiens capensis spotted touch-me-not herbaceous 109 Table A.l (cont’d). JUNI LALA LYSA2 MOFI ONSE PHAU7 PODE3 POTRl 0 PRSEZ ROMU RUID SAGIT SALIX SOLID SYFO THOC2 TOVE TYPHA SP. ULAM URDI Juglans nigra Larix’ larr'cina Lythrum salicaria Monardafistulosa Onoclea sensibilis Phragmites australis Populus deltoides Populus tremuloides Prunus serotina Rosa multiflora Rubus idaeus Sagittaria Sp. Salix sp. Solidago Sp. Symplocarpusfoetia'us Thuja occidentalis T oxicodendron vemz'x T ypha sp. Ulmus americana Urtica dior'ca 110 black walnut tamarack purple loosestrife wild bergamot sensitive fern phragrnites Eastern cottonwood quaking aspen black cherry multiflora rose raspberry common arrowhead willow Sp. goldenrod sp. skunk cabbage Northern white cedar poison sumac cattails American elm stinging nettle tree or sapling tree or sapling herbaceous herbaceous herbaceous herbaceous tree or sapling tree or sapling tree or sapling shrub shrub herbaceous shrub herbaceous herbaceous tree or sapling ' Shrub herbaceous tree or sapling herbaceous adu- APPENDIX B. Deterministic population growth rates for eastern massasauga § rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan. Table 3.1. Deterministic population growth rates calculated with Vortex 9.96 (Chicago Zoological Society, Brookfield, IL) for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan, derived flom a baseline scenario. Assumptions included no stochastic fluctuations, no inbreeding depression, no limitation of mates, no harvest, and no supplementation. Deterministic growth rate Projection Exponential rate of increase, r 0.108 Annual rate of change, A 1.115 Net replacement rate, R0 1.798 Mean generation time, (3‘ , S? 5.41 , 5.41 lll APPENDIX C. The snake community occurring in cover types used by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus). Table CI. The snake community occurring in cover types used by the focal species, eastern massasauga rattlesnake (Sistrurus catenatus catenatus; n = 69), included 7 species of snake at the study area in Barry County, Michigan flom 2008 to 2009. Observations have been contributed to the statewide online public database, Michigan Herpetological Atlas (MDNRE 2010). Common name Latin name No. processed Garter snake Thamnophis sirtalis 16 Northern ribbon snake Thamnophis sauritus septentrionalz‘s 6 Blue racer Coluber constrictor foxi 5 Brown snake Storeria dekayr' 3 Northern water snake Nerodr'a sipedon 2 Eastern hognose snake Heterodon platirhinos 1 Eastern milk snake Lampropeltis triangulum triangulum 1 Total 34 112 APPENDIX D. Metadata for radiolocations of 26 eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) tracked in 2008 and 2009 in Barry County, Michigan. Title: 08_09 pointsshp Format: ArcView 3.2 layer and all associated files (*shp, *.dbf, *.shx) Originator: Michigan State University Robyn L. Bailey, Graduate Research Assistant, Dept. of Fisheries and Wildlife, 13 Natural Resources, Michigan State University, East Lansing, MI 48824. Henry Campa, III, Professor, Dept. of Fisheries and Wildlife, 13 Natural Resources, Michigan State University, East Lansing, MI 48824. campa@msu.edu; 51 7/353-2042. Funding: Financial support for this researCh was provided by Michigan State University, Michigan Department of Natural Resources and Environment (Wildlife Division), Pierce Cedar Creek Institute (Hastings, MI), and Potter Park Zoo (Lansing, MI). 1 Date Created: 2009 Description: Digital point layer of all radiolocations recorded for eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) in Barry County, Michigan. Attributes include date, snake ID, waypoint name, coordinates, distance, bearing, estimate of precision error (EPE), dilution of precision (DOP), notes, new casting and northing (created flom distance and bearing whenever a snake was 3—30 m flom a previous waypoint), and source theme. Spatial Reference: Projected Coordinate System: NAD_1983__UTM_Zone_16N 113 Projection: Transverse_Mercator Geographic Coordinate System: GCS_North_American_l 983 Datum: D_North_American_l983 Scale Factor: 0.99960 Map Units: meters Contact: Technical questions regarding digital locations can be directed to: Robyn L. Bailey 13 Natural Resources East Lansing, MI 48824 or Methods: All snakes fitted with transmitters were relocated daily flom May through October of 2008 and 2009. Locations were recorded with a Garmin® Global Positioning System unit (Garrnin International Inc., Olathe, KS) in Universal Transverse Mercator (U TM) coordinates. Coordinates were recorded as northings and castings in a *.dbf file. In order to keep location error similar among nearby locations, subsequent locations were recorded as distances (m) and bearings (°) (measured using a meter tape and a compass) flom the former GPS-recorded location to the new locations. This was done only when the new location was 3—30 m flom the last GPS location. When snakes were 53 m flom a former location, they were recorded as not having moved. Locations were recorded using the GPS once the snake moved >30 m flom the previous GPS location, if there was too much vegetation between the previous GPS location and the new location for distance and bearing measurements, or if the snake was located in a 114 different vegetation type than the former location. The northing and casting coordinates recorded in the *.dbf files for each snake were added as themes and converted to shapefiles in ArcView 3.2. A script was written (K. Bissell, unpublished data) to generate coordinates flom the distance and bearing information collected in the field. The resulting coordinates were then converted to a point shapefile and transferred to ArcMap 9.2 for viewing. Shapefiles for GPS locations were then merged with shapefiles for distance-bearing locations and the resulting shapefile was projected in NAD_1983_UTM_Zone_l6N. Data Quality: Locations were compared visually to an IMAGINE image for the study area (which included building and road footprints) and found to be representative. The final layer file has 1,893 locations for 26 snakes. There is no data for snake 10R, and I recommend not using data flom snake 6R (poor health). Snakes 2R, 4R, 7R, 11R, and 16R were tracked in 2008 and 2009. Data Accuracy: All locations of rattlesnakes were recorded at the exact snake locations. Researchers generally had a visual of individuals at most locations and were able to pinpoint hidden individuals within 2 m using radiotelemetry. Estimated position error (EPE, a measurement of horizontal position error in meters based upon a variety of factors including dilution of precision and satellite signal quality) for GPS locations were never >8 m. Dilution of precision (DOP, a measure of the GPS receiver/satellite geometry, lower is better) for this data was never > 2.2. Therefore, human and equipment error was likely minimal. 115 Literature Cited Aebischer, N. J ., P. A. Robertson, and R. E. Kenward. 1993. Compositional analysis of habitat use flom animal radio-tracking data. Ecology 74:1313—1325. Atkinson, D. A., and M. G. Netting. 1927. The distribution and habits of the massasauga. Bulletin of Antivenin Institute of America 1:40—44. Beissinger, S. R., and M. I. Westphal. 1998. On the use of demographic models of population viability in endangered species management. Journal of Wildlife Management 622821—841 . Berry, K. H. 1986. Introduction: development, testing, and application of wildlife- habitat models. Pages 3—4 in J. Vemer, M. L. Monison, and C. J. Ralph, editors. Wildlife 2000: modeling habitat relationships of terrestrial vertebrates. University of Wisconsin Press, Madison, USA. Bissell, K. M. 2006. Modeling habitat ecology and population viability of the eastern massasauga rattlesnake in southwestern lower Michigan. Thesis, Michigan State University, East Lansing, USA. Bonnet, X., G. Naulleau, and R. Shine. 1999. The dangers of leaving home: dispersal and mortality in snakes. Biological Conservation 89:39—50. Brooks, R. P. 1997. Improving habitat suitability index models. Wildlife Society Bulletin 25:163—167. ' Brown, W. S., M. Kéry, and J. E. Hines. 2007. Survival of timber rattlesnakes (Crotalus horridus) estimated by capture-recapture models in relation to age, sex, color morph, time, and birthplace. Copeia 2007:656—67 1. Bunck, C. M., and K. H. Pollock. 1993. Estimating survival of radiotagged birds. Pages 51—63 in J. C. Lebreton and P. M. North, editors. Marked individuals in the study of bird populations. Birkhauser-Verlag, Basel, Switzerland. Canfield, R. H. 1941. Application of the line interception method in sampling range vegetation. Journal of Forestry 39:388-394. Clark, R. W., W. S. Brown, R. Stechert, and K. R. Zarnudio. In press. Roads, interrupted dispersal, and genetic diversity in timber rattlesnakes. Conservation Biology. Cook, F. R. 1993. After an ice age: zoogeography of the massasauga within a Canadian herpetofaunal perspective. Pages 19—25 in B. Johnson and V. Menzies, editors. Proceedings of the international symposium and workshop on the conservation of 116 the eastern massasauga rattlesnake Sistrurus catenatus catenatus. Toronto Zoo, 8—9 May 1992, Toronto, Ontario, Canada. Crone, E. E. 2001. Is survivorship a better fitness surrogate than fecundity? Evolution 55:2611—2614. Diller, L. V., and R. L. Wallace. 2002. Growth, reproduction, and survival in a population of Crotalus viridis oreganus in north central Idaho. Herpetological Monographs 2002:26—45. Durbian, F. E. 2006. Effects of mowing and summer burning on the massasauga (Sistrurus catenatus). American Midland Naturalist 155:329—334. Durbian, F. E., and L. Lenhoff. 2004. Potential effects of mowing prior to summer burning on the eastern massasauga (Sistrurus c. catenatus) at Squaw Creek National Wildlife Refuge, Holt County, Missouri, USA. Transactions of the Missouri Academy of Science 38:21—25. Durbian, F. E., R. S. King, T. Crabill, H. Larnbert-Doherty, and R. A. Seigel. 2008. Massasauga home range patterns in the Midwest. Journal of Wildlife Management 72:754—759. Ellner, S. P., J. Fieberg, D. Ludwig, and C. Wilcox. 2002. Precision of population viability analysis. Conservation Biology 16:258—261. Erickson, W. P., T. L. McDonald, K. G. Gerow, S. Howlin, and J. W. Kern. 2001. . Statistical issues in resource selection studies with radio-marked animals. Pages 209—242 in J. J. Millspaugh and J. H. Marzluff, editors. Radio tracking and animal populations. Academic Press, San Diego, California, USA. Federal Emergency Management Agency. 2009. Map Service Center. J essup, Maryland, USA. . Accessed 19 Nov 2009. Fitch, H. S. 1960. Autoecology of the copperhead. University of Kansas Publications of the Museum of Natural History 13:85—288. Fuller, M. R., J. J. Millspaugh, K. E. Church, and R. E. Kenward. 2005. Wildlife radiotelemetry. Pages 377—417 in CB. Braun, editor. Techniques for wildlife investigations and management. Sixth edition. The Wildlife Society, Bethesda, Maryland, USA. Garshellis, D. L. 2000. Delusions in habitat evaluation. Pages 111—164 in L. Boitani and T. K. Fuller, editors. Research techniques in animal ecology: controversies and consequences. Columbia University Press, New York, New York, USA. Gibbs, H. L., K. A. Prior, P. J. Weatherhead, and G. Johnson. 1997. Genetic structure 117 .4_. of populations of the threatened eastern massasauga rattlesnake, Sistrurus c. catenatus: evidence flom microsatellite DNA markers. Molecular Ecology 6:1123- 1 132. Gregory, P. T., J. M. Maccartney, and K. W. Larsen. 1987. Spatial patterns and movements. Pages 366—395 in R. A. Seigel, J. T. Collins, and S. S. Novak, editors. Snakes: ecology and evolutionary biology. Macmillan Publishing, New York, New York, USA. Hall, L. S., P. R. Krausman, and M. L. Morrison. 1997. The habitat concept and a plea for standard terminology. Wildlife Society Bulletin 25:173—182. Harvey, D. S., and P. J. Weatherhead. 2006a. Hibernation site selection by eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) near their northern range limit. Journal of Herpetology 40:66—73. Harvey, D. S., and P. J. Weatherhead. 2006b. A test Of the hierarchical model of habitat selection using eastern massasauga rattlesnakes (Sistrurus c. catenatus). Biological Conservation 130:206—216. Haufler, J. B., R. K. Baydack, H. Campa, III, B. J. Kemohan, C. Miller, L. J. O’Neil, and L. Waits. 2002. Performance measures for ecosystem management and ecological sustainability. Wildlife Society Technical Review 02-1. Hurlbert, S. H. 1984. Pseudoreplication and the design of ecological field experiments. Ecological Monographs 54:187—21 1. J ellen, B. C., and M. J. Kowalski. 2007. Movement and growth of neonate eastern massasaugas (Sistrurus catenatus). Copeia 2007:994-1000. Johnson, D. H. 1979. Estimating nest success: the Mayfield method and an alternative. Auk 96:651—661. Johnson, D. H. 1980. The comparison of usage and availability measurements for evaluating resource preference. Ecology 61 :65—71 . Johnson, D. H. 1994. Population analysis. Pages 419—444 in T. A. Bookhout, editor. Research and management techniques for wildlife and habitat. The Wildlife Society, Bethesda, Maryland, USA. Johnson, G. 1995. Spatial ecology, habitat preference, and habitat management of the eastern massasauga, Sistrurus c. catenatus in a New York weakly-minerotrophrc peatland. Dissertation, State University of New York, Syracuse, USA. Johnson, G. and D. J. Leopold. 1998. Habitat management for the eastern massasauga in a central New York peatland. Journal of Wildlife Management 62:84—97. 118 Johnson, G., B. Kingsbury, R. King, C. Parent, R. Seigel, and J. Szymanski. 2000. The eastern massasauga rattlesnake: a handbook for land managers. US. Fish and Wildlife Service, Fort Snelling, Minnesota, USA. Jones, 3., S. F. Fox, D. M. Leslie, Jr., D. M. Engle, and R. L. Lochmiller. 2000. Herpetofaunal responses to brush management with herbicide and fire. Journal of Range Management 53: 154—158. Kaplan, E. L., and P. Meier. 1958. Nonparametric estimation flom incomplete observations. Journal of the American Statistical Association 53:457—48 1. Keenlyne, K. D. 1968. The massasauga (Sistrurus catenatus catenatus) in Buffalo County, Wisconsin. Thesis, University of Minnesota, Minneapolis, USA. Keenlyne, K. D. 1978. Reproductive cycles in two species of rattlesnakes. American Midland Naturalist 100:368—3 75. Keenlyne, K. D. and J. R. Beer. 1973. Note on the size of Sistrurus catenatus catenatus at birth. Journal of Herpetology 34:381—382. Kemohan, B. J., R. A. Gitzen, and J. J. Millspaugh. 2001. Analysis of animal space use and movements. Pages 125—166 in J. J. Millspaugh and J. H. Marzluff, editors. Radio tracking and animal populations. Academic Press, San Diego, California, USA. King, R. S. 1999. Habitat use and movement patterns of the eastern massasauga rattlesnake in Wisconsin. Page 80 in B. Johnson and M. Wright, editors. Second international symposium and workshop on the conservation of the eastern massasauga rattlesnake, Sistrurus catenatus catenatus: population and habitat management issues in urban, bog, prairie, and forested ecosystems. Toronto Zoo, 2- 3 October 1998, Toronto, Ontario, Canada. King, R., C. Berg, and B. Hay. 2004. A repatriation study of the eastern massasauga (Sistrurus catenatus catenatus) in Wisconsin. Herpetologica 60:429—43 7. Kroll, A. J ., and J. B. Haufler. 2006. Development and evaluation of habitat models at multiple spatial scales: a case study with the dusky flycatcher. Forest Ecology and Management 229: 1 61—1 69. Kruskal, W. H., and W. A. Wallis. 1952. Use Of ranks in one-criterion variance analysis. Journal of the American Statistical Association 47 :583—621. Lander, R. A., D. A. Adams, and E. M. Lunk. 1986. Temporal and spatial aspects of Species-habitat models. Pages 65—69 in J. Vemer, M. L. Morrison, and C. J. Ralph, editors. Wildlife 2000: modeling habitat relationships of terrestrial vertebrates. University of Wisconsin Press, Madison, USA. 119 Laubhan, M. K., S. L. King, and L. H. Fredrickson. 2005. Managing inland wetlands for wildlife. Pages 797—838 in GE. Braun, editor. Techniques for wildlife investigations and management. Sixth edition. The Wildlife Society, Bethesda, Maryland, USA. Lorimer, C. G. 2001. Historical and ecological roles of disturbance in eastern North American forests: 9,000 years of change. Wildlife Society Bulletin 29:425—439. Macartney, J. M. 1985. The ecology of the northern pacific rattlesnake, Crotalus viridis oreganus, in British Columbia. Thesis, University of Victoria, British Columbia, Canada. Madsen, T., R. Shine, J. Loman, and T. Hz‘ikansson. 1992. Why do female adders copulate so flequently? Nature 335:440—441. Maple, W. T. 1968. The overwintering adaptations of Sistrurus c. catenatus in northeastern Ohio. Thesis, Kent State University, Kent, Ohio, USA. Marcot, B. G., M. G. Raphael, and K. H. Berry. 1983. Monitoring wildlife habitat and validation Of wildlife-habitat relationship models. Transactions of the North American Wildlife and Natural Resources Conference 48:315—329. Marshall, J. C., J. V. Manning, and B. A. Kingsbury. 2006. Movement and macrohabitat selection of the eastern massasauga in a fen habitat. Herpetologica 62: 141—150. Mayfield, H. 1961. Nesting success calculated flom exposure. Wilson Bulletin 73:255—261. McLoughlin, P. D. and F. Messier. 2004. Relative contributions of sampling error in initial population size and vital rates to outcomes of population viability analysis. Conservation Biology 18:1665—1669. Michigan Department of Agriculture. 2004. Climate of Hastings. Michigan Department of Agriculture, Climatology Program, East Lansing, USA. Michigan Department of Natural Resources [MDNR]. 2003. 2001 Integrated forest management analysis program/Gap analysis program land cover. Forest, Mineral, and Fire Management Department, Lansing, Michigan, USA. Michigan Department of Natural Resources and Environment [MDNRE]. 2010. Herp atlas. . Accessed 17 Mar 2010. Miller, P. 2006. Population viability assessment for the eastern massasauga rattlesnake (Sistrurus catenatus catenatus) on the Bruce Peninsula, Ontario, Canada. Third 120 international symposium and workshop on conservation of the eastern massasauga, Sistrurus catenatus: population viability and outreach. Massasauga Recovery Team, 12—14 October 2005, Toronto, Ontario, Canada. Miller, P. S., and R. C. Lacy. 2005. VORTEX: a stochastic simulation of the extinction process, version 9.50 user’s manual. Conservation Breeding Specialist Group, Apple Valley, Minnesota, USA. Mills, L. S., J. M. Scott, K. M. Strickler, and S. A. Temple. 2005. Ecology and management of small populations. Pages 691—713 in CB. Braun, editor. Techniques for wildlife investigations and management. Sixth edition. The Wildlife Society, Bethesda, Maryland. Moore, J. and J. C. Gillingham. 2006. Spatial ecology and multi—scale habitat selection by a threatened rattlesnake: the eastern massasauga (Sistrurus catenatus catenatus). Copeia 2006:742—75 1. Morrison, M. L. 2001. A proposed research emphasis to overcome the limits of wildlife-habitat relationship studies. Journal of Wildlife Management 65:613—623. Miinzbergova, Z. and J. Ehrlén. 2005. How best to collect demographic data for population viability analysis. Journal of Applied Ecology 42:1115—1120. Natural Resources Conservation Service. 2000. Soil survey geographic (SSURGO) database for Barry County, Michigan. US. Department of Agriculture, Fort Worth, Texas, USA. . Accessed 8 Mar 2010. Negri, S. J. 1995. Analysis of habitat suitability models for primary cavity-nesting birds in Michigan’s upper peninsula. Thesis, Michigan State University, East Lansing, USA. New York State Departurent of Enviromnental Conservation. 2010. Eastern massasauga fact sheet. Endangered Species Unit, Albany, New York, USA. . Accessed 17 Mar 2010. Odurn, E. P., and E. J. Kuenzler. 1955. Measurement of territory and home range size in birds. Auk 72:128—137. Oreskes, N., K. Shrader-Frechette, and K. Belitz. 1994. Verification, validation, and continuation of numerical models in the earth sciences. Science 263:641—646. Parent, C. and P. J. Weatherhead. 2000. Behavioral and life history responses of eastern massasauga rattlesnakes (Sistrurus catenatus catenatus) to human disturbance. Oecologica 125:170—178. 121 ./" Parker, W. S., and M. V. Plummer. 1987. Population ecology. Pages 253—301 in R. A. Seigel, J. T. Collins, and S. S. Novak, editors. Snakes: ecology and evolutionary biology. Macmillan Publishing, New York, New York, USA. Pauly, D., and V. Christensen. 2006. Modeling wildlife-habitat relationships. Pages 320—376 in M. L. Morrison, B. G. Marcot, and R. W. Manan, editors. Wildlife- habitat relationships: concepts and applications. Third edition. Island Press, Washington, DC, USA. Presser, D. J ., and R. P. Brooks. 1998. A verified habitat suitability index for the Louisiana watertlrrush. Journal of Field Ornithology 69:288—298. Reinert, H. K. 1978. The ecology and morphological variation of the massasauga rattlesnake (Sistrurus catenatus). Thesis, Clarion State College, Clarion, Pennsylvania, USA. Reinert, H. K. 1981. Reproduction by the massasauga (Sistrurus catenatus catenatus). American Midland Naturalist 105:393-395. Reinert, H. K. and W. R. Kodrich. 1982. Movements and habitat utilization by the massasauga, Sistrurus catenatus catenatus. Journal of Herpetology 16:162—1 71. Rittenhouse, C. D., F. R. Thompson, III, W. D. Dijak, J. J. Millspaugh, and R. L. Clawson. 2010. Evaluation of habitat suitability models for forest passerines using demographic data. Journal of Wildlife Management 74:411—422. Roloff, G. J., and B. J. Kemohan. 1999. Evaluating reliability of habitat suitability index models. Wildlife Society Bulletin 27:973—985. Schamberger, M. L. and L. J. O’Neil. 1986. Concepts and constraints of habitat-model testing. Pages 5—10 in J. Vemer, M. L. Monison, and C. J. Ralph, editors. Wildlife 2000: modeling habitat relationships of terrestrial vertebrates. University of Wisconsin Press, Madison, USA. Seaman, D. E., J. J. Millspaugh, B. J. Kemohan, G. C. Brundige, K. J. Raedeke, and R. A. Gitzen. 1999. Effects of sample size on kernel home range estimates. Journal of Wildlife Management 63:739—747. Seigel, R. A. 1986. Ecology and conservation of an endangered rattlesnake, Sistrurus catenatus, in Missouri, USA. Biological Conservation 35:333—346. Seigel, R. A., and C. A. Sheil. 1999. Population viability analysis: applications for the conservation of massasaugas. Pages 17—22 in B. Johnson and M. Wright, editors. Second international symposium and workshop on the conservation of the eastern massasauga rattlesnake, Sistrurus catenatus catenatus: population and habitat 122 management issues in urban, bog, prairie and forested ecosystems. Toronto ZOO, 2— 3 October 1998, Toronto, Ontario, Canada. Seigel, R. A., C. A. Sheil, and J. S. DOOdy. 1998. Changes in a population of an endangered rattlesnake Sistrurus catenatus following a severe flood. Biological Conservation 83:127—131. Shine, R. 2003. Reproductive strategies in snakes. Proceedings of the Royal Society of London, Biological Sciences 270:995—1004. Slaughter, B. S., and J. D. Skean, Jr. 2003. Annotated checklist of vascular plants in the vicinity of Cedar Creek and Brewster Lake, Pierce Cedar Creek Institute, Barry County, Michigan. Michigan Botanist 42:127—148. Space Imaging. 2004. Integrated forest monitoring assessment and prescription IFMAP: review of remote sensing technologies used in the IFMAP project. Report to Michigan Department of Natural Resources. Space Imaging, Ann Arbor, Michigan, USA. Stauffer, D. F ., and L. B. Best. 1986. Effects of habitat type and sample size on habitat suitability index models. Pages 71—77 in J. Vemer, M. L. Morrison, and C. J. Ralph, editors. Wildlife 2000: modeling habitat relationships of terrestrial vertebrates. University of Wisconsin Press, Madison, USA. Swanson, P. L. 1933. The size of Sistrurus c. catenatus at birth. Copeia 1933:3 7. Szymanski, J. 1998. Status assessment for the eastern massasauga (Sistrurus c. catenatus). US. Fish and Wildlife Service, Ft. Snelling, Minnesota, USA. US. Fish and Wildlife Service [USFWS]. 1981. Standards for the development of habitat suitability index models. US. Fish and Wildlife Service Release No. 1-81, 103 ESM. US. Fish and Wildlife Service [USFWS]. 1999. Endangered and threatened wildlife and plants; review of plant and animal taxa that are candidates or proposed for listing as endangered or threatened; annual notice of findings on recycled petitions; and annual description of progress on listing actions. Federal Register 64:57534— 57547. Van Horne, B. 1983. Density as a misleading indicator of habitat quality. J ournal of Wildlife Management 47:893—901. Weatherhead, P. J ., and F. W. Anderka. 1984. An improved radio transmitter and implantation technique for snakes. Journal of Herpetology 18:264—269. 123 Weatherhead, P. J ., and K. A. Prior. 1992. Preliminary Observations of habitat use and movements of the eastern massasauga rattlesnake (Sistrurus c. catenatus). Journal of Herpetology 26:447—452. Weatherhead, P. J ., J. M. Knox, D. S. Harvey, D. Wynn, J. Chiucchi, and H. L. Gibbs. 2009. Diet Of Sistrurus catenatus in Ontario and Ohio: effects of body size and habitat. Journal of Herpetology 43:693—697. White, G. C. 2000. Population viability analysis: data requirements and essential analyses. Pages 288—331 in L. Boitani and T.K. Fuller, editors. Research techniques in animal ecology: controversies and consequences. Columbia University Press, New York, New York, USA. Wilgers, D. J ., and E. A. Home. 2006. Effects of different burn regimes on tallgrass prairie herpetofaunal species diversity and community composition in the Flint Hills, Kansas. Journal of Herpetology 40:73—84. Winterstein, S. R., H. Campa, III, and K. F. Millenbah. 1995. Status and Potential of Michigan Natural Resources: Wildlife. Special report 75. Michigan Agricultural Experiment Station, Michigan State University, East Lansing, USA. Winterstein, S. R., Pollock, K. H., and C. M. Bunck. 2001. Analysis of survival data flom radiotelemetry studies. Pages 351—3 80 in J. J. Millspaugh and J. H. Marzluff, editors. Radio tracking and animal populations. Academic Press, San Diego, California, USA. Withey, J. C., T. D. Bloxton, and J. M. Marzluff, 2001. Effects of tagging and location error in wildlife radiotelemetry studies. Pages 43—75 in J. J. Millspaugh and J. H. Marzluff, editors. Radio tracking and animal populations. Academic Press, San Diego, California, USA. Woodbury, A. M. 1951. A snake den in Tooele County, Utah: introduction—a ten year study. Herpetologica 724—14. Wright, B. A. 1941. Habit and habitat studies of the massasauga rattlesnake (Sistrurus catenatus catenatus Raf.) in northeastern Illinois. American Midland Naturalist 25:659—672. 124