LANDSCAPE GENETICS OF BLACK BEARS (URSUS AMERICANUS) IN THE NORTHERN LOWER PENINSULA (NLP) OF MICHIGAN, USA By Hope M. Draheim A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for a degree of Zoology—Doctor of Philosophy Ecology, Evolutionary Biology, and Behavior—Dual Major 2015 ABSTRACT LANDSCAPE GENETICS OF BLACK BEARS (URSUS AMERICANUS) IN THE NORTHERN LOWER PENINSULA (NLP) OF MICHIGAN, USA By Hope M. Draheim My dissertation combines population genetics, landscape ecology and spatial statistics to examine connectivity and degree of gene flow in black bears (Ursus americanus). More specifically, I used landscape genetic approaches to address an array of questions, including identifying specific landscape features that impede or promote dispersal, identify source-sink dynamics, and to infer the effect of landscape change on connectivity of a large, highly harvested, isolated black bear population in the Northern Lower Peninsula (NLP) of Michigan. My dissertation consists of four chapters following the introductory chapter (Chapter 1). Chapter 2 evaluates the effects of timing and design in collection methodology on our ability to draw inferences about the extent of gene flow in NLP black bears. This chapter shows black bears exhibit significant positive genetic spatial autocorrelation. Results were concordant across seasons and years indicating the overall pattern of spatial genetic structure would be detected regardless of season, year, or collection methodology. Chapter 3 results indicate that two genetically distinct groups exist in the NLP, defined as a western and an eastern genetic cluster. In addition, I found land cover was the landscape feature most strongly correlated with genetic distance. Chapter 4 investigates the presence of source-sink dynamics in NLP black bear population. I found asymmetric emigration and immigration, revealing habitat for black bears in the NLP is composed of source and sink areas. Further, I showed source strength was associated with black bear local harvest density (a proxy for bear density). Chapter 5 examines the influence of landscape change on spatial genetic structure over time. My times series analyses shows land cover was significantly correlated with gene flow; however, comparisons among sampling years revealed temporal variability in the predictive power and performance of landscape resistance models. I found the best model to explain temporal inconsistencies was land cover change over time. These data represents one of most comprehensive landscape genetic studies performed on a single black bear population. Our data enables managers to target regions or habitat types that are important for maintaining connectivity across anthropogenically-altered habitats, such as the NLP. In addition, here I present two novel extensions of current landscape genetics analytical approaches that provide flexible frameworks for understanding connectivity and assessing impacts of future landscape alternation which can be widely integrated into landscape genetics research and conservation planning at multiple spatial scales. Copyright by HOPE M. DRAHEIM 2015 ACKNOWLEDGEMENTS Science is a collaborative endeavor and completion of this dissertation would not have been possible without the assistance and support of many people and organizations. First and foremost, I thank Dr. Kim Scribner, my graduate advisor, who has been an immense and constant source of inspiration, encouragement, and guidance over the years. I am very fortunate to have been his student and have benefitted from his expertise. The success of this project would not have been possible without his foresight and ability to recognize the scientific merit of such a long term data set. In addition, my advisory committee: Scott Winterstein, Joseph Messina, Dwayne Etter, and Kay Holekamp, have all been extremely supportive and have made me a stronger scientist. Scott has gone above and beyond to ensure my best interests and provided me with a strong example of collegiality. I am especially thankful for Scott filling in as an “on campus” advisor during Kim’s sabbatical. I could always count on Joe to really push me to consider difficult questions related to my science, that ultimately strengthen my understanding of my field and confidence in my work. Dwayne really pushed me to consider the management applications of my work and helped me navigate between academia and management. Kay is an ideal role model for any women in science. Her passion and knowledge is inspirational and I consider it a privilege to have learned from her. I would like to thank the funding sources that have enabled this work. Funding for my research and data was provided by the Michigan Department of Natural Resources (MDNR) through the Wildlife and Sportfish Restoration Program F11AF00640 and the Department of Fisheries at Michigan State University. I want to thank the within Michigan State University: the College of Natural Sciences, the Department of Zoology, the Department Fisheries and Wildlife, v the Ecology, Evolutionary Biology and Behavior program, and the Graduate School for funding and providing logistical support to me as a graduate student. Also, I want to thank the Safari Club and NASA that have supported me as a student and a researcher. In addition, the large scale nature of sample collection for this project would not have been possible if it were not for the samples provided by Michigan bear hunters and collected by the many MDNR staff at registration stations. I want to extend special recognition to Jen Moore whom is responsible for generating a great deal of the genetic data and is an endless source of support both personally and professionally. Jen your scientific knowledge coupled with your infinite patience turned abstract ideas into realized scientific approaches. Your council is a gift I am constantly grateful for. I am particularly thankful to Kari Dammerman and Jeannette Kanefsky, not only for their guidance in the lab and with analyses, but for their friendship and laughter that has made my graduate experience so memorable. Additionally, I would like to thank the all members of the“Scribner” lab, especially Jeanette McGuire, Shairah Razak, Katy Jay, John Bauman, Lisette Delgado, Yen Duong, an Jan Hassenauer who have been a scientific family for me. Lastly, I want to thank my friends and family for their unwavering support for my aspirations and career goals. vi TABLE OF CONTENTS LIST OF TABLES ........................................................................................................................ ix LIST OF FIGURES ...................................................................................................................... xi CHAPTER 1 - GENERAL INTRODUCTION ..............................................................................1 RFERENCES ..................................................................................................................................6 CHAPTER 2 - EFFECTS OF SPATIAL AND TEMPORAL SAMPLING SCALE ON SAMPLE DISPERSION AND SPATIAL GENETIC STRUCTURE OF A BLACK BEAR POPULATION IN MICHIGAN ...............................................................................................................................9 ABSTRACT .................................................................................................................................9 INTRODUCTION .....................................................................................................................10 METHODS ................................................................................................................................11 Study Area ..............................................................................................................................11 Sampling .................................................................................................................................13 Individual Identification and Quality Control .........................................................................14 Point Pattern Analysis of Samples ..........................................................................................15 Genetic Diversity, Hardy-Weinberg, Linkage Disequilibrium ...............................................16 Spatial Genetic Analysis .........................................................................................................16 RESULTS ..................................................................................................................................18 Number of bears and genetic diversity ...................................................................................18 Point Pattern Analysis .............................................................................................................18 Spatial Genetic Structure ........................................................................................................22 DISCUSSION ............................................................................................................................25 MANAGEMENT IMPLICATIONS .........................................................................................28 ACKNOWLEDGEMENTS .......................................................................................................28 APPENDIX ..................................................................................................................................30 REFERENCES .............................................................................................................................34 CHAPTER 3 - SPATIAL GENETIC STRUCTURE AND LANDSCAPE CONNECTIVITY IN BLACK BEARS: INVESTIGATING THE SIGNIFICANCE OF DIFFERENT LAND COVER DATA SETS IN LANDSCAPE GENETICS ANALYSES .........................................................40 ABSTRACT ...............................................................................................................................40 INTRODUCTION .....................................................................................................................41 METHODS ................................................................................................................................43 Sampling, DNA extraction, PCR amplification ......................................................................43 Population genetic analysis .....................................................................................................45 Landscape Genetic Analysis ...................................................................................................46 RESULTS ..................................................................................................................................51 Genetic Diversity and Spatial Genetic Structure ....................................................................51 Landscape Genetic Analysis ...................................................................................................53 DISCUSSION ............................................................................................................................55 Black bear spatial genetic structure ........................................................................................55 vii Effect of landscape features on black bear connectivity .........................................................56 Land cover dataset comparison ...............................................................................................57 CONCLUSIONS ........................................................................................................................60 ACKNOWLEDGMENTS .........................................................................................................61 APPENDIX ..................................................................................................................................62 REFERENCES .............................................................................................................................65 CHAPTER 4 – DETECTING BLACK BEAR SOURCE-SINK DYNAMICS USING INDIVIDUAL-BASED GENETIC GRAPHS .............................................................................71 ABSTRACT ...............................................................................................................................71 INTRODUCTION .....................................................................................................................72 METHODS ................................................................................................................................75 Study Area ..............................................................................................................................75 Laboratory Analysis ................................................................................................................75 Defining the Nodes .................................................................................................................77 Graph construction ..................................................................................................................78 Node importance and source-sink analysis .............................................................................82 RESULTS ..................................................................................................................................84 Node removal ..........................................................................................................................86 DISCUSSION ............................................................................................................................86 CONCLUSIONS ........................................................................................................................90 ACKNOWLEDGMENTS .........................................................................................................91 APPENDIX ..................................................................................................................................92 REFERENCES .............................................................................................................................99 CHAPTER 5 – BEYOND THE SNAPSHOT: LANDSCAPE GENETIC ANALYSIS OF TIME SERIES DATA REVEALS RESPONSES OF AMERICAN BLACK BEARS TO LANDSCAPE CHANGE .....................................................................................................................................105 ABSTRACT .............................................................................................................................105 INTRODUCTION ...................................................................................................................106 METHODS ..............................................................................................................................109 Study Area and Sample Collection .......................................................................................109 Microsatellite Genotyping ....................................................................................................111 Defining the Temporal Sampling Locations .........................................................................111 Landscape Genetic Analysis .................................................................................................112 Genetic and Landscape Change Analysis .............................................................................114 RESULTS ................................................................................................................................117 Spatial Genetic Structure ......................................................................................................116 Landscape Genetic Analyses ................................................................................................121 Genetic Change and Landscape Change ...............................................................................121 DISCUSSION ..........................................................................................................................125 ACKNOWLEDGMENTS .......................................................................................................130 APPENDIX ................................................................................................................................131 REFERENCES ...........................................................................................................................136 viii LIST OF TABLES Table 2.1. Sample size (N), number of resampled individuals (NRES) removed, number of unique individuals (NUNIQ), probability of identity (PID), probability of identity for siblings (PSIB), mean number of alleles (Na), mean observed (HO) and expected heterozygosity (HE) over all loci, and mean relatedness (rxy) for each sample period for the Northern Lower Peninsula, Michigan black bear population. ....................................................19 Table 3.1. Mantel (r, IBD) and partial Mantel correlations (r, all resistance models) between geographic and genetic pairwise distances among 365 individual black bear in the NLP. Bold indicates significant correlations. .............................................................50 Table B1. Side by side comparisons of land cover types according to IFMAP and CCAP classification schemes IFMAP and CCAP data the class number and descriptions pulled from original data sets. Class numbers and descriptions for the black bear ecological model were originated from Carter (2007). ................................................63 Table 4.2. Spearman's rank correlations between genetic and ecological measures of node importance (source strength). The three graph metrics compared are net flux (Fout-Fin), area-weighted flux (AWF), and probability of connectivity (PC). ‘High’ indicates maximum magnitude weights (e.g.,1-100) and ’low’ indicates minimum magnitude weights (e.g.,1-5) used to test the sensitivity of analyses to the magnitude of cost values. All correlations were significant (P < 0.05). ..................................................79 Table 4.1. Node weights for graph analysis and cost values for resistance surfaces used to calculate pairwise least cost distance for edge weights. Node weights and cost values are inferred from Carter et al. (2010). Habitat/Non-Habitat values designate land cover positively or negatively correlated with bear presence inferred from Carter et al. 2010. ......................................................................................................................80 Table C1. Black bear genetic diversity at 12 microsatellite loci. No loci showed significant deviations from Hardy-Weinberg equilibrium (Raymond & Rousset 1995). Na, number of alleles per locus; HO, observed heterozygosity; HE, expected heterozygosity; FIS, inbreeding coefficient). ...............................................................93 Table 5.1. Mantel (isolation by distance only; IBD) and partial Mantel correlations (r) between spatial and genetic pairwise distances among individual black bear in the NLP for 2002, 2006, and 2010. Bold indicates competitive models based on causal modeling. ....................................................................................................................................114 Table 5.2. A summary of the landscape variables used to characterize genetic change between 2002, 2006, and 2010 for black bears in the Northern Lower Peninsula of Michigan, USA. ..........................................................................................................................119 ix Table 5.3. Most supported univariate spatial regression models of landscape alteration variables predicting black bear genetic change. Abbreviations are as described in Table 5.1. ....................................................................................................................................124 Table 5.4. Most supported multiple spatial regression models of landscape change characteristics predicting black bear genetic change. ...............................................126 Table D1. Black bear genetic diversity at 12 microsatellite loci. No loci showed significant deviations from Hardy-Weinberg equilibrium (Raymond & Rousset 1995). Na, number of alleles per locus; Ho, observed heterozygosity; HE, expected heterozygosity. ..........................................................................................................133 x LIST OF FIGURES Figure 2.1. Study area in the Northern Lower Peninsula of Michigan showing the location of summer hair (open circles) and fall tissue (filled circles) sample collections for individual black bears sampled during 2003, 2005, and 2009 following removal of duplicate individuals based on multi-locus genotypes. ...............................................12 Figure 2.2. Adaptive kernel density plots showing black bear sample density within Michigan’s Northern Lower Peninsula for the summer and fall of 2003, 2005, and 2009 following removal of duplicate individuals based on multi-locus genotypes. Light shades are areas of high sample density and dark shades are areas of low sample density. ........20 Figure 2.3. Results of K-function analysis characterizing differences in spatial patterns among black bear sample collections made during different seasons (summer and fall) and different years (2003, 2005, and 2009). Solid black lines are the observed differences between K-function estimates at a given distance, and the dashed gray lines indicate the 95% confidence intervals. ..................................................................................... 21 Figure 2.4. Global spatial autocorrelation of summer (squares), fall (triangles) and combined summer/fall (circles) black bear samples for 2003, 2005, and 2009 full (a) and bootstrapped (b) data sets. Correlograms (solid lines) with standard errors. The 95% confidence intervals (dashed lines) about the null hypothesis of no spatial structure were estimated by permutations based on distance classes of 5 km. ..........................23 Figure 2.5. Two-dimensional local spatial autocorrelation analysis of black bear microsatellite data using five nearest neighbors (see Supplemental Materials, Fig. 1 for results of comparable analyses using 3, 7, and 9 nearest neighbors). Circles represent individuals with positive significant local spatial autocorrelation (lr). Size of the circles indicates the strength of the autocorrelation coefficient. .................................25 Figure A1. Two-dimensional local spatial autocorrelation analysis of black bear microsatellite data using (a) three, (b) seven and (c) nine nearest neighbors. Circles represent individuals with positive significant local spatial autocorrelation (lr). Size of the circles indicates the strength of the autocorrelation coefficient. .................................31 Figure 3.1. Study area in the Northern Lower Peninsula of Michigan (NLP) showing locations of black bear harvest samples collected during 2006 (N = 365). Pie graphs were constructed using on individual posterior probabilities of each cluster (K = 2) from Program STRUCTURE. .............................................................................................44 Figure 3.2. Side by side comparison of CCAP and IFMAP land cover maps using black bear classification scheme described in Carter 2007. (A) Land cover maps for the entire study area. (B) Land cover of the same 7 x 10 km area. .............................................48 xi Figure 3.3. Global spatial autocorrelation of black bear samples. Bars represent standard errors. The 95% confidence intervals (dashed lines) about the null hypothesis of no spatial structure were estimated by permutations based on distance classes of 5 km. ...........52 Figure 3.4. Comparison of cumulative resistance maps between pairs of individuals from three competitive isolation-by-resistance models in Table 1:1) land cover only; 2) land cover and rivers; and 3) land cover, rivers, and roads generated from CCAP land cover using program CIRCUITSCAPE. Gradient of colors indicated the probability of black bear movement. Yellow colors indicate high probability; red and pink colors indicate medium probability; and purple and blue representing low probability. ......54 Figure 4.1. Study area in the Northern Lower Peninsula of Michigan (NLP) showing node locations (N = 141) for black bear harvest samples collected during 2002, 2006, and 2010 (N = 569). ...........................................................................................................76 Figure 4.2. Interpolated surfaces of black bear node importance (source strength) based on net flux at the nodes across the Northern Lower Peninsula study area using Inverse Distance Weighting (IDW). Interpolated surfaces of (a) the genetic, (b) Carter et al. 2010 habitat suitability, (c) a habitat/non-habitat, and (d) a bear density model are shown for comparison. Warm colors represent sources and cool colors represent sinks. ...........................................................................................................................85 Figure 4.3. Area-weighted dispersal flux (AWF, a) and probability of connectivity (b) as a function of the number of nodes removed for each four node removal scenarios (random, sources only, sinks only and intermediate only). Error bars represent the standard error (SE) from three independent runs. .......................................................87 Figure C1. Map of the NLP land cover. Land cover data were reclassified into eight land cover classifications from 2006 NOAA Coastal Change Analysis Program Land Cover (CCAP) dataset (resolution = 30m). ...............................................................................94 Figure C2. Map of NLP local black bear harvest density in NLP. Local density was classified into a range of 1-10 (low to high). ..............................................................................95 Figure C3. Map of overlapping Voronoi tessellation point pattern polygons in NLP. Polygons contained at least one sample from 2002, 2006, and 2010. ........................................96 Figure C4. Distribution of the node importance values characterized by inter-individual genetic relatedness, shown are the a) net flux, b) area-weighted flux (AWF) and c) probability of connectivity (PC) metrics. For net flux positive values represent source nodes and negative values represent sink values. ........................................................................97 Figure C5. (a) Distribution of mean relatedness (K) among individuals within the nodes. (b) Mean relatedness at the node according to number of individuals at the node. .........98 Figure 5.1. Study area in the Northern Lower Peninsula of Michigan (NLP) showing areas (n = 141) of consistent sampling for black bear harvested during 2002, 2006, and 2010 (n = 569). .........................................................................................110 xii Figure 5.2. Scatterplot of genetic distance versus geographic distance including least squares regression line (2002, P = 0.01; 2006, P = 0.03; 2010, P = 0.02; Mantel test). Genetic distance was based on pairwise estimates of proportion of shared alleles (Dps) among all individuals. Geographic distances (km) were calculated as the straight-line distance from the center of each sampling location. .................................................121 Figure 5.3. Distribution of (A) number of landscape change patches within a grid cell and (B) the cumulative difference in genetic differentiation within a cell grid for all temporal comparisons (2002-2006, 2006-2010, 2002-2010). For landscape change maps red indicates large number of patches of change and blue indicates low number of patches of change. For genetic change maps red indicates increase in genetic differentiation (GD) from time one (T1) to time two (T2) and blue indicates decrease in genetic differentiation from T1 to T2. ...................................................................................123 Figure D1. Extent and configuration of forested land lost (green) over the sampling period from (A) 2001-2006, (B) 2006-2010, and (C) 2002-2010. ................................................134 Figure D2. Map of overlapping Voronoi tessellation point pattern polygons in NLP. Polygons contained at least one sample from 2002, 2006, and 2010. ......................................135 xiii CHAPTER 1 GENERAL INTRODUCTION With increasing human population numbers and anticipated climatic change, present environments are expected to be altered significantly. Projected human population increases can impact species abundances and distributions as habitat alterations occur due to deforestation, pollution, and other anthropogenic factors (Hoffman & Willi 2008). Land conversion due to human activities places stress on ecosystems to provide requisite requirements for wildlife populations (De Sherbinin et al. 2007). Formerly unaltered continuous habitats are presently fragmented into mosaics of varying habitat qualities creating differential permeabilities within the landscape that can either facilitate or impede movement of organisms between locales (Manel et al. 2003; Spear et al. 2010). Identifying the factors that influence the degree of population connectivity across heterogeneous landscapes provides insight into population dynamics and the evolutionary trajectory of a species (Segelbacher et al. 2010). Understanding demographic and genetic exchange among populations is an important component of developing conservation priorities and management plans for species (Lande 1988; Lande et al. 2003). Therefore, it is important to understand how species use and move through the landscape in order to monitor the consequences of environmental change and predict possible responses under different environmental scenarios. Thus, correctly quantifying the degree of connectivity and gene flow for natural populations is a priority in population conservation and management (Beier & Noss 1998; Bunn et al. 2000; Rabinowitz & Zeller 2010). Landscape genetics seeks to identify how landscape features and environmental factors influence the spatial distribution of genetic variation. For example, landscape genetics assesses how abiotic (e.g., rivers), biotic (e.g. forest structure), or anthropogenic features (e.g. roads) 1 affect gene flow (Anderson et al. 2010; Bolliger et al. 2010; McRae et al. 2005; Spear & Storfer 2010). Unlike traditional theoretical population genetic models (i.e. Wright’s island model and isolation by distance model; Wright 1931, Wright 1943) that assume populations are spatially homogeneous and occupy a uniform landscape (Kimura & Weiss 1964; Wright 1943), a basic tenet of landscape genetics is that natural ecosystems are spatially and temporally heterogeneous. Thus, such heterogeneity may have important ecological implications for organisms inhabiting them (Bolliger et al. 2010; Spear et al. 2010; Storfer et al. 2007). In species that have large space requirements, habitat heterogeneity and connectedness can be important, both ecologically and genetically. For example, wide ranging species such as the black bear (Ursus americanus) prefer large continuous forested areas and can be sensitive disruptions in land cover (Pelton 1992). Habitat fragmentation and pervasiveness of human dominated landscapes have been shown to impede black bear dispersal (Cushman 2006). Historically abundant, black bears inhabited most forested areas of North America. Today, due to habitat alteration, black bears primarily restricted to forested areas with limited human population (Pelton et al., 1982). Within Michigan, the highest density (~85 %) occurs in the UP due to the high quality bear habitat and an abundance of mixed forest and large tracts of forest lands. (Dreher 2004, Carter 2007). In contrast, the Lower Peninsula (LP) is a mosaic of fragmented land and subsequently black bear density is lower. Important factors in habitat selection appear to be availability of food resources, denning sites, cover, and evasion of large human population (Carter et al. 2010). The black bear population in Michigan’s Northern Lower Peninsula (NLP) provides a unique opportunity to apply spatial genetic analytical approaches using multi-locus genotypes in a landscape context. First, by capitalizing on teeth collected annually from all harvested bears, 2 we have acquired a long term data set that spans nine years (2002-2010). Second, the NLP is an isolated population bounded to the south by a matrix of unsuitable urban and agricultural habitat and on all other sides by the Great Lakes. Thus, this population experiences no emigration and immigration that could confound results. Third, the NLP shows considerable landscape heterogeneity, with fragmented areas ranging from high to low quality and has experienced ongoing landscape change (i.e. deforestation and crop conversion). Fourth, NLP bear population is subjected to intensive harvest, where approximately 13-29% of bears are harvested annually (Michigan Department of Natural Resources, Pers. comm.) indicating the potential for rapid population turn over in as little as 4 years. Finally, times series land cover information via Coastal Change Analysis Program (CCAP) for 2001, 2006, and 2010 (covering the entire genetic sampling period) is available for the NLP allowing for investigations into the influence of landscape features on spatial genetic structure across generations. The aim of my dissertation is therefore to utilize multi-locus genotypes, in combination with landscape variables (i.e land cover, roads, and rivers) as tools to address questions pertaining to the influence of landscapes on genetic variation and patterns of genetic differentiation in a solitary, wide-ranging, and elusive species. I examined data from an isolated black bear population in Michigan’s NLP. Genetic sampling of rare and elusive species can be challenging and by necessity, samples are often collected opportunistically from multiple sources that could differ in spatial dispersion and timing of collection. In Chapter 2, I examine the effects of timing and design in collection methodology on estimates of gene flow in NLP black bears. Sampling bias can lead to erroneous interpretations of spatial genetic structure that can subsequently impact conservation efforts and management decisions. I found spatial distributions of individual black bears in 3 Michigan’s NLP sampled during the summer and fall were similar a between seasons and among years. In addition, the timing and method of sampling in different seasons and years did not appear to greatly influence patterns of spatial genetic structure. These results led me to conclude that the overall pattern of spatial genetic structure would be detected regardless of season, year, or collection methodology. In Chapter 3, I quantified levels of spatial genetic structure and the effects of the landscape on connectivity of black bears in NLP. In addition, landscape genetic analyses require generation of resistance surfaces using landscape data layers (e.g. land cover types) that are hypothesized relationships between an organism’s probability of movement through an area and features of the environment. However, there is a degree of uncertainty of how representative a land cover data set portrays reality. Therefore, the relationship between gene flow and landscape features was evaluated using two different land cover data sets to investigate the relative importance of land cover uncertainty on landscape genetic model performance. I found black bears in the NLP exhibit spatial genetic structure and land cover was significantly correlated with genetic distance. However, I found differences in landscape genetic model performance when different land cover data sets were used. My results emphasize the importance of considering land cover classification methodologies and the inherent uncertainty associated with land cover data in landscape genetic studies. In the fourth chapter, I introduce a novel approach that uses individual-based genetic graphs to identify source-sink areas within NLP’s black bear population. Combining population genetics and graph theory offers a compelling framework (Gaggiotti 1996; Murphy et al. 2010; Peery et al. 2008; Pierson et al. 2013) for understanding source sink dynamics and population connectivity. My genetic graphs showed asymmetric connectivity, with areas of high and low net 4 flux (source strength), revealing black bears in the NLP exhibited source-sink dynamics. I also inquired whether source strength was associated with high quality habitat and high local population density two characteristics commonly linked with source areas. I found significant correlations between graph metrics derived from genetic data and ecological models based on habitat quality and local harvest density. My results show genetic graphs provide a viable alternative method for detecting source-sink dynamics in continuously distributed species and a flexible framework for understanding connectivity which could be widely integrated into landscape genetics research and conservation planning at multiple spatial scales. In Chapter 5, I introduce a time series approach in a landscape genetic framework to address how landscape factors affect spatial genetic structure of black bears in the NLP over time and space using multi-locus genetic and land cover data sets sampled in 2002, 2006, and 2010. The promise of time series approaches has been recognized in population genetics (Habel et al. 2014); however, to my knowledge there have no empirical applications in landscape genetics. I found temporal variability in the support of our landscape resistance models which resulted in different conclusions. Had analyses been based solely on samples collected in 2002 we would conclude that land cover has no influence on functional connectivity of black bears in the NLP. Nevertheless, we detected a strong association between black bear SGS and land cover in 2006 and 2010. We found the best model to explain temporal inconsistencies was degree of land cover change. My results demonstrate the importance of conducting long term genetic monitoring for the conservation and management of wildlife in rapidly changing landscapes. 5 REFERENCES 6 REFERENCES Anderson CD, Epperson BK, Fortin MJ, et al. (2010) Considering spatial and temporal scale in landscape-genetic studies of gene flow. Molecular Ecology 19, 3565-3575. Beier P, Noss RF (1998) Do habitat corridors provide connectivity? Conservation Biology 12, 1241-1252. Bolliger J, Le Lay G, Holderegger R (2010) Landscape genetics - how landscapes affect ecological processes. Gaia-Ecological Perspectives for Science and Society 19, 238-240. Bunn AG, Urban DL, Keitt TH (2000) Landscape connectivity: A conservation application of graph theory. Journal of Environmental Management 59, 265-278. Carter NH (2007) Predicting ecological and social suitability of black bear habitat in Michigan's Lower Peninsula, University of Michigan. Carter NH, Brown DG, Etter DR, Visser LG (2010) American black bear habitat selection in northern Lower Peninsula, Michigan, USA, using discrete-choice modeling. Ursus 21, 57-71. Cushman SA (2006) Effects of habitat loss and fragmentation on amphibians: A review and prospectus. Biological Conservation 128, 231-240. De Sherbinin A, Carr D, Cassels S, Jiang L (2007) Population and environment. Annual review of environment and resources 32, 345. Gaggiotti OE (1996) Population genetic models of source-sink metapopulations. Theoretical Population Biology 50, 178-208. Habel JC, Husemann M, Finger A, Danley PD, Zachos FE (2014) The relevance of time series in molecular ecology and conservation biology. Biological Reviews 89, 484-492. Lande R (1988) Genetics and demography in biological conservation. Science(Washington) 241, 1455-1460. Lande R, Engen S, Saether B-E (2003) Stochastic population dynamics in ecology and conservation Oxford University Press Oxford. Manel S, Schwartz MK, Luikart G, Taberlet P (2003) Landscape genetics: combining landscape ecology and population genetics. Trends in Ecology & Evolution 18, 189-197. McRae BH, Beier P, Dewald LE, Huynh LY, Keim P (2005) Habitat barriers limit gene flow and illuminate historical events in a wide-ranging carnivore, the American puma. Molecular Ecology 14, 1965-1977. 7 Murphy MA, Dezzani R, Pilliod DS, Storfer A (2010) Landscape genetics of high mountain frog metapopulations. Molecular Ecology 19, 3634-3649. Peery MZ, Beissinger SR, House RF, et al. (2008) Characterizing source-sink dynamics with genetic parentage assignments. Ecology 89, 2746-2759. Pelton MR (1992) Black Bear (Ursus americanus) In: Wild Mammals of North America: Biology, Management, and Economics (eds. Chapman JA, Feldhamer GA), pp. 504-514. The John Hopkins University Press,USA. Pierson JC, Allendorf FW, Drapeau P, Schwartz MK (2013) Breed locally, disperse globally: fine-scale genetic structure despite landscape-scale panmixia in fire-specialist. Plos One 8, e67248. Rabinowitz A, Zeller KA (2010) A range-wide model of landscape connectivity and conservation for the jaguar, Panthera onca. Biological Conservation 143, 939-945. Segelbacher G, Cushman SA, Epperson BK, et al. (2010) Applications of landscape genetics in conservation biology: concepts and challenges. Conservation Genetics 11, 375-385. Spear SF, Balkenhol N, Fortin MJ, McRae BH, Scribner K (2010) Use of resistance surfaces for landscape genetic studies: considerations for parameterization and analysis. Molecular Ecology 19, 3576-3591. Spear SF, Storfer A (2010) Anthropogenic and natural disturbance lead to differing patterns of gene flow in the Rocky Mountain tailed frog, Ascaphus montanus. Biological Conservation 143, 778-786. Wright S (1931) Evolution in Mendelian populations. Genetics 16, 0097-0159. Wright S (1943) Isolation by distance. Genetics 28, 114-138. 8 CHAPTER 2 EFFECTS OF SPATIAL AND TEMPORAL SAMPLING SCALE ON SAMPLE DISPERSION AND SPATIAL GENETIC STRUCTURE OF A BLACK BEAR POPULATION IN MICHIGAN Draheim, H. M., V. Lopez, D. Etter, S. R. Winterstein, K.T. Scribner ABSTRACT Sampling bias can lead to erroneous interpretations of spatial genetic structure that can subsequently impact conservation efforts and management decisions. Genetic sampling of rare and elusive species can be challenging and by necessity, samples are often collected opportunistically from multiple sources that could differ in spatial dispersion and timing of collection. Here we quantified the effects of timing of sample collection and sampling methods on spatial and temporal variability in sample dispersion and measures of black bear (Ursus americanus) spatial genetic structure. Hair (N = 890) and tissue (N = 1017) samples were collected from Michigan’s Northern Lower Peninsula during the summer using non-invasive (hair snares) and fall harvest methods during of 2003, 2005, and 2009. Point pattern analyses of sample dispersion revealed that sample density did not differ significantly between seasons or among years. Measures of spatial genetic structure (i.e. spatial autocorrelation) revealed significant positive genetic spatial structuring at distances of 0-10 km when samples were analyzed separately by year, season, and when temporally distinct samples were grouped. Local genetic spatial autocorrelation analyses revealed spatial genetic patterns over the entire study area were consistent across seasons and years. Spatial genetic structure is indicative of the extent of dispersal and gene flow which is crucial for developing management plans for harvested 9 species. Collectively, our data show that non-invasive and harvest collection methods similarly capture spatial genetic heterogeneity at the scale and spatial extent appropriate for management. INTRODUCTION Quantifying the degree of connectivity for natural populations is a priority in population conservation and management (Beier & Noss 1998; Bunn et al. 2000; Rabinowitz & Zeller 2010). Genetic sampling of rare and elusive species can be challenging and by necessity, samples are often collected opportunistically from a variety of sources that could differ in timing, effort, and method of collection. To obtain adequate sample numbers, spatial dispersion, and spatial extent, samples may be collected over multiple seasons and/or years (De Barba et al. 2010). For example, samples of many fur bearing species may be collected during direct research efforts and opportunistically from harvests (Dreher et al. 2007; Proctor et al. 2005; Schwartz et al. 2003; Williams & Scribner 2010). Pooling samples collected during different times assumes temporal stability in genetic structure and that distributions of individuals do not vary (Garant et al. 2000). However, this assumption is frequently violated because distributions of individuals typically vary by seasons. Seasonal and inter-annul variation in habitat occupancy may increase the “noise” and mask signatures of spatial genetic structure, thus call into question characterizations of spatial genetic patterns. Measures of fine scale spatial genetic structure can be affected if individuals are outside normal home ranges and breeding areas at the time of collection. Latch and Rhodes (2006) found significant genetic structuring among localized winter flocks of wild turkeys (Meleagris gallopavo) when family groups are spatially aggregated yet no evidence for genetic structure among samples collected in the spring when males move widely while searching for mates. 10 Individuals may be spatially clustered during a particular time of the year like foraging. For example, during the summer several migratory winter elk (Cervus elaphus) herds overlap in Yellowstone National Park to exploit higher-quality forage. In contrast, during breeding seasons geographic overlap is significantly smaller (White et al. 2010). Genetic sampling of American black bear (Ursus americanus) in Michigan’s Northern Lower Peninsula (NLP) has been implemented using two collection methods, field collection of hair during the summer and tissue collected from the fall harvest. If black bears in Michigan’s NLP occupy different habitats seasonally, for example, to areas of high food resources, intraannual variation has potential implications for the interpretation spatial genetic structure. Our objective was to evaluate the effects of timing and heterogeneity in collection methodology on our ability to draw inferences about the extent of gene flow. Spatial and temporal variation of summer hair and fall harvest sample collection methodologies provided an opportunity to evaluate the influence of timing and collection method on sample dispersion and spatial and temporal patterns in measures of black bear spatial genetic structure. METHODS Study Area The study area located in the NLP encompassed approximately 47,739 km2 (Figure 2.1) and contains three Bear Management Units (BMU) (Red Oak, Baldwin, and Gladwin). Approximately 27% of the lands are in public ownership (i.e. state forest and national forest lands) with the majority of public lands located in the northeast and west-central portions of the NLP. The NLP bear population is isolated with the borders of Lake Michigan to the west and north, Lake Huron to the east and north, and intensive human development and agriculture to the 11 Figure 2.1. Study area in the Northern Lower Peninsula of Michigan showing the location of summer hair (open circles) and fall tissue (filled circles) sample collections for individual black bears sampled during 2003, 2005, and 2009 following removal of duplicate individuals based on multi-locus genotypes. 2003 2005 2009 12 south. The NLP is a fragmented mosaic composed of different land cover and land use types including development, agriculture, upland non-forested openings, northern hardwood and mixed hardwood, oak, aspen, pine, forested wetland, non-forested wetland, and open water (Carter et al. 2010). Sampling Sampling was conducted during the summer and fall of 2003, 2005, and 2009. A total of 1907 samples were collected using two methods. Tissue (teeth: N = 1017) samples were obtained from bears brought into hunter check stations during September and October. Registration of all harvested bears is mandatory in Michigan and teeth are collected from > 95% of bears registered annually (D. Etter, MDNR pers. comm). The MDNR issues bear licenses within each BMU annually and the number in each BMU was similar among years. In addition, location within a BMU was influenced by hunter behavior (e.g., establishment of bait location where hunters believe they were most likely to successfully harvest a bear) and site accessibility (Frawley 2001; Frawley 2009). Locations of tissue samples were recorded in township, range, section (e.g. 2.59 km2). Coordinates of harvest locations were georeferenced using the section centroid and converted into UTM coordinates. Hair samples (N = 890) were collected using hair snares (2003 N = 239; 2005 N = 236; 2009 N = 258) during the summer (late May-early July). Hair samples were collected using the protocol outlined in Dreher et al. (2007). Two strands of barbed wire encircled a baited location at uniform strand heights of 20 cm and 50 cm. We collected hair samples from each barb on the top strand (strand at 50 cm). Hair snagged on the 20 cm strand was not collected to minimize the probability of sampling cubs. The location of the snares was similar among years based on 13 criteria established by Dreher et al. (2007). Dreher et al. (2007) found that 403 of 414 harvested bears in 2003 were within a radius of one female bear home-range of a hair-snare and the remainder were within a radius of two female home-ranges. Hair samples were georeferenced to UTM coordinates using a GPS. Individual Identification and Quality Control We extracted DNA from bear hair and teeth using Qiagen DNEasy Tissue Kits following manufacturer protocols (Qiagen Inc., Valencia, CA). We anticipated low qualities of DNA from hair samples, thus we did not quantify DNA concentrations from hair. However, DNA obtained from teeth was quantified using a Nanodrop spectrophotometer (Thermo Scientific, Waltham, MA) and diluted to a 20 ng/μl working concentration. Multilocus genotypes were amplified using polymerase chain reaction (PCR) for five polymorphic microsatellite loci including, G10X, G10L, G10D (Paetkau et al. 1995), UarMU59, and UarMU50 (Taberlet et al. 1997) using methods described in Dreher et al. (2007). We used quality control protocols developed by Paetkau (2003) to assess and minimize genotyping errors. We initially genotyped all samples for all 5 loci and eliminated genotypes that could not be confidently scored by 2 experienced lab personnel. To minimize genotyping error and ensure summer hair samples represented unique individuals and not recaptures we identified duplicate genotypes for all hair and tissue samples within a year using program GENECAP (Wilberg & Dreher 2004). Genotyping errors identified and corrected by reexamining all genotypes that were differed by 1 or 2 loci. To ensure each sample represented one individual, duplicate genotypes from hair samples were assumed to be replicate samples and were not used for further analysis. To assess genotyping error 10% of all samples (N=1805) within each year and season were randomly selected and genotyped a second 14 time for all loci yielding a genotyping error rate of < 2.8%. Probability of identity (PID) and probability of identity for siblings (PSIB) over all loci was calculated using program GenAlEx v. 6.0 (Peakall & Smouse 2006). Point Pattern Analysis of Samples Kernel density estimation is widely used to model sample distribution patterns. We used an adaptive kernel estimator allowing for a varying bandwidth to reduce the variance of areas of low sample size and reduce the bias of areas of high sample size (Silverman 1986). Kernel estimation was weighted to account for sampling effort. We visually compared the spatial distribution and density of individual samples for each sampling period. All maps were produced in spatstat package in R 2.13 (R Development Core Team 2011). To compare the differences between sampling point patterns from different seasons, collection methods, and years we calculated differences between sample K-functions. Ripley’s K function is a measure of clustering or dispersion of observations based on the distance between all pairs of observations, and is used to examine spatial patterns for a range of distances (Dixon 2002). We analyzed differences in mapped spatial distributions between datasets by comparing K-functions for a homogeneous Poisson process using a complete spatial randomness model (CSR) for seasons within a year and between years for summer, fall and a pooled summer/fall sample group. Ripley’s K-function was calculated separately for summer and fall sample datasets for 2003, 2005, and 2009 for distances up to 100 km for each sample dataset. Differences between estimates (e.g. Ksummer – Kfall, K2003 – K2005, K2003 – K2009, and K2005 – K2009) were calculated for sampling group comparisons at all distances. We performed permutations (n = 999) to generated 95% upper and lower confidence intervals to test for significant differences 15 in spatial configurations of samples by comparing the observed difference in K-functions between two sample groups and the difference in K-functions of two randomly drawn data sets (Rowlingson & Diggle 1993). All analyses were performed using splancs package in R 2.13 (R_Development_Core_Team 2011). Genetic Diversity, Hardy-Weinberg, Linkage Disequilibrium Program MICRO-CHECKER (Van Oosterhout et al. 2004), was used to test for presence of null alleles and allelic dropout. We tested for deviations from Hardy-Weinberg using program GENEPOP (Version 3.1d; Raymond and Rousset 1995), and used sequential Bonferroni tests to correct for multiple tests (Rice 1989). We used Bonferroni corrections to test for linkage disequilibrium in program FSTAT 2.93 (Goudet 1995, 2001). We quantified microsatellite genetic diversity in each sampling group using mean number of alleles (A), observed heterozygosity (HO), and expected heterozygosity (HE) over all loci, we calculated the pairwise relatedness coefficient between individuals (rij, following Queller and Goodnight 1989), using program GenAlEx . Differences in number of alleles, expected heterozygosities, and probability of identity between sampling locations were assessed using a Friedman test, which uses loci as the blocking factors, accounting for inter locus variation (Zar 1996). Spatial Genetic Analysis Spatial genetic structure was assessed using a global genetic spatial autocorrelation coefficient (r) calculated using program GenAlEx (Peakall & Smouse 2006). Coefficient r (range, - 1 to + 1) was calculated by comparing pairwise geographical and genetic distance estimates with a null hypothesis of r = 0 when no correlation exists between geographic and genetic distance. 16 Measures of spatial autocorrelation can be influenced by sample number (Epperson et al. 1999). Thus we also implemented a bootstrap resampling procedure to account for discrepancies in sample size between sample groups. One hundred simulated populations of 167, 121, 132 individuals corresponding to the number of individuals in the summer sample for 2003, 2005 and 2009; respectively, were created using a random resampling approach in program POPTOOLS (Hood 2010). Then r was computed using the multi-locus genotypes of all 100 simulated data sets and averaged for 10 distance classes spaced at five km intervals to ensure an adequate number of samples per distance class. The significance of r was determined via 999 permutations to estimate the 2-tailed 95% confidence intervals around zero and the bootstrap procedure was used to provide 95% confidence intervals for the autocorrelation estimated for each distance class. Global spatial autocorrelation estimates are valuable for detecting the presence of spatial genetic structure and distances over which non-random genetic relationships exist. However, these autocorrelation estimates cannot inform researchers of whether spatial patterns are geographically homogeneous or exhibit variation in the relative contribution from different groups of individuals sampled. Local spatial autocorrelation (LSA) analysis is a useful alternative to examine spatial patterns at a fine geographic scale. LSA evaluates subsets of data using nearest neighbors to identify areas occupied by highly related individuals (Anselin 1995). We used methods described in Double et al. (2005) to calculate local autocorrelation coefficient, lr for each individual using 3, 5, 7, and 9 nearest neighbors in program GenAlEx 6 (Peakall & Smouse 2006). Each lr estimate was determined by comparing the estimated value with a null distribution of no spatial structure, using 999 permutations. To visualize LSA we plotted 17 locations of significant groups of related individuals using bubble graphs and superimposed the study area on the graphs. RESULTS Number of bears and genetic diversity We identified 1434 unique genotypes and 473 duplicate samples (hair = 470; tissue = 3; Table 2.1). Estimates of probability of identity (PID and PSIB) were similar across all sampling groups (ranges 2.4 x 10-6 to 6.4 x 10-6, respectively; Friedman Test, P = 0.90, P = 0.80, respectively; Table 2.1). We found no presence of null alleles or allelic dropout, and no loci were found to deviate significantly from Hardy-Weinberg or linkage equilibrium. Expected heterozygosity (range 0.74 to 0.77, Table 2.1) and mean number of alleles (range 8.0 to 9.6, Table 2.1) were not significantly different among sampling groups (Friedman Test, P = 0.58, P = 0.79, respectively; Table 2.1). Point Patterns Analysis Kernel density plots for each sampling group (summer and fall samples for 2003, 2005, and 2009) indicate that areas of highest and lowest sampling density were consistent among sampling periods (Figure 2.2). Kernel density plots reveal two major findings. First, there was a region of high sample density in the east central portion of the study area which was observed for all seasons and years. Second, there was a region of moderately high sample density observed in the south central portion of the study area for most sample groups. Differences in spatial distribution of sampling patterns among years revealed spatial distributions were not significantly different (Figure 2.3). The plots of all K2003 – K2005, K2003 – K2009, and K2005 – K2009, did not fall outside the 95% confidence interval for most distances 18 Table 2.1. Sample size (N), number of resampled individuals (NRES) removed, number of unique individuals (NUNIQ), probability of identity (PID), probability of identity for siblings (PSIB), mean number of alleles (Na), mean observed (HO) and expected heterozygosity (HE) over all loci, and mean relatedness (rxy) for each sample period for the Northern Lower Peninsula, Michigan black bear population. N NRES NIND PID PSIB Na HO HE rxy 302 135 167 2.4 x 10-6 8.4 x 10-3 -6 9.60 0.76 0.77 -0.008 -3 8.20 0.74 0.76 -0.002 2003 Hair Tissue 410 3 407 3.0 x 10 9.0 x 10 Hair 299 178 121 3.8 x 10-6 8.9 x 10-3 8.20 0.76 0.76 -0.008 -6 -3 8.20 0.74 0.76 -0.004 2005 Tissue 248 0 248 4.3 x 10 9.4 x 10 Hair 289 157 132 6.4 x 10-6 1.1 x 10-2 8.00 0.70 0.74 -0.008 Tissue 359 0 359 3.4 x 10-6 9.0 x 10-3 9.00 0.74 0.76 -0.003 2009 19 Figure 2.2. Adaptive kernel density plots showing black bear sample density within Michigan’s Northern Lower Peninsula for the summer and fall of 2003, 2005, and 2009 following removal of duplicate individuals based on multi-locus genotypes. Light shades are areas of high sample density and dark shades are areas of low sample density. 2003 Summer 2003 Fall 2005 Summer 2005 Fall 2009 Summer 2009 Fall High Low 20 Figure 2.3. Results of K-function analysis characterizing differences in spatial patterns among black bear sample collections made during different seasons (summer and fall) and different years (2003, 2005, and 2009). Solid black lines are the observed differences between K-function estimates at a given distance, and the dashed gray lines indicate the 95% confidence intervals. 21 suggesting the dispersion of samples among years was similar (not independent). In contrast, plots of Ksummer – Kfall for all years was greater than zero at short distances (2003 = 10 km; 2005 = 20 km; 2009 = 15 km) indicating summer samples were significantly more clustered up to ~20 km. Spatial Genetic Structure Global spatial autocorrelation analysis indicated a pattern of isolation by distance. However, analyses using complete data sets revealed larger standard errors for the summer samples indicating that sample size could affect significance at larger distance classes (Figure 2.4a). In contrast, after accounting for unequal sample size significant positive genetic autocorrelation was observed for individuals sampled at inter-individual distances of 0-10 km for all sample groups (Figure 2.4b). In general, after bootstrap resampling to account for unequal samples sizes; slightly higher positive correlation values were observed for the summer samples in all years when compared to the fall and pooled samples from the summer and fall. Values from pooling summer and fall samples resulted in intermediate estimates of spatial autocorrelation relative to the estimate from a single sampling period. One dataset did not appear to disproportionately influence the pooled estimate. Three areas (the north, east central and southwest central portion of the research area) of high local autocorrelation were consistently observed in summer and fall in samples from 2003, 2005 and 2009 using five nearest neighbors (P < 0.05; Figure 2.5). The consistency of the pattern was confirmed with similar numbers and distribution of positive autocorrelation values found using 3, 7, and 9 nearest neighbors (Figure A1). 22 Figure 2.4. Global spatial autocorrelation of summer (squares), fall (triangles) and combined summer/fall (circles) black bear samples for 2003, 2005, and 2009 full (a) and bootstrapped (b) data sets. Correlograms (solid lines) with standard errors. The 95% confidence intervals (dashed lines) about the null hypothesis of no spatial structure were estimated by permutations based on distance classes of 5 km. 23 Figure 2.5. Two-dimensional local spatial autocorrelation analysis of black bear microsatellite data using five nearest neighbors (see Supplemental Materials, Figure 1 for results of comparable analyses using 3, 7, and 9 nearest neighbors). Circles represent individuals with positive significant local spatial autocorrelation (lr). Size of the circles indicates the strength of the autocorrelation coefficient. 2003 Summer 6.0E+05 6.0E+05 6.0E+05 2003 Fall 5.5E+05 5.5E+05 NORTHERLY NORTHERLY 5.5E+05 5.0E+05 5.0E+05 5.0E+05 4.5E+05 4.5E+05 4.5E+05 4.0E+05 4.0E+05 4.0E+05 3.5E+05 3.5E+05 3.5E+05 3.0E+05 4.5E+05 5.0E+05 5.5E+05 6.0E+05 6.5E+05 7.0E+05 7.5E+05 5.5E+05 6.0E+05 6.0E+05 2005 Summer 6.5E+05 7.0E+05 7.5E+05 7.0E+05 7.5E+05 2005 Fall 5.5E+05 5.5E+05 5.0E+05 5.0E+05 5.0E+05 4.5E+05 4.5E+05 4.5E+05 4.0E+05 4.0E+05 4.0E+05 3.5E+05 3.5E+05 3.5E+05 3.0E+05 3.0E+05 3.0E+05 4.5E+05 5.0E+05 5.5E+05 6.0E+05 6.5E+05 7.0E+05 4.5E+05 7.5E+05 5.0E+05 5.5E+05 6.0E+05 6.5E+05 6.0E+05 6.0E+05 2009 Summer 2009 Fall NORTHERLY NORTHERLY 5.5E+05 5.5E+05 5.0E+05 5.0E+05 5.0E+05 4.5E+05 4.5E+05 4.5E+05 4.0E+05 4.0E+05 4.0E+05 3.5E+05 3.5E+05 3.5E+05 3.0E+05 4.5E+05 6.0E+05 EASTERLY EASTERLY 5.5E+05 6.0E+05 NORTHERLY NORTHERLY NORTHERLY 5.5E+05 5.0E+05 EASTERLY EASTERLY 6.0E+05 3.0E+05 3.0E+05 4.5E+05 5.0E+05 5.5E+05 6.0E+05 6.5E+05 7.0E+05 7.5E+05 3.0E+05 3.0E+05 4.5E+05 EASTERLY 5.0E+05 5.5E+05 6.0E+05 6.5E+05 EASTERLY 24 7.0E+05 7.5E+05 DISCUSSION Spatial distributions of individual black bears sampled during the summer and fall were similar and analogous spatial genetic structure exists between seasons and among years for black bears in Michigan’s NLP. Sample distribution was concordant across seasons, years and collection methods. Kernel density plots show areas of high sampling density to be consistent among sampling groups (Figure 2.2). Landscape in Michigan’s NLP is a highly fragmented and black bears inhabit large contiguous tracts of forested land (Manville 1987). Areas of high sample density were located in Montmorency, Oscoda, Alpena and Alcona counties on private land and forested public lands (Ausable, Pere Marquette, Mackinaw States Forest and Huron National Forest lands). Point pattern analyses using differences in Ripley’s K-function revealed that summer samples were more clustered at short distances (< 20 km), indicating some differences in the spatial distributions (Figure 2.3). Because individual hair snares capture multiple bears, locations are likely to have multiple observations with the same coordinates, thus inflating the degree of clustering. In contrast, greater variation in distances among fall samples was expected as more than one bear will not generally be harvested in the same location. Estimated bear density in the NLP was < 0.24 bear/km2 in any year sampled (MDNR, unpublished data). The timing and method of sampling in different seasons and years did not appear to greatly influence patterns of spatial genetic structure. We observed the same overall pattern of significant positive spatial genetic correlation for distances up to 10 km in all sampling periods. Positive spatial autocorrelation is expected to accrue due to restricted dispersal and beyond these distance classes individuals are no more related than expected by chance (Epperson 2005). Positive spatial genetic autocorrelation over short distances is likely attributed to female natal philopatry commonly exhibited in black bears where female offspring establish home ranges 25 adjacent to the mothers whereas male offspring disperse from the natal area (Costello 2010; Rogers 1977; Rogers 1987; Schwartz & Franzmann 1992). The isolation by distance pattern is consistent with those reported in bears and other wide ranging carnivores (Brown et al. 2009; Paetkau et al. 1997; Rueness et al. 2003b). Radio-tracking studies have demonstrated adult male bears are more active during the summer breeding season while seeking out perspective mates (Amstrup & Beecham 1976; Garshelis & Pelton 1981). In addition, sub-adult males move large distances to find and establish new home ranges away from natal territories (Lee & Vaughan 2003). These bear movement patterns suggest genetic structure should be comparatively weaker during the summer than in the fall which is inconsistent with our results. We observed slightly higher spatial genetic autocorrelation during summer in the shortest distance class relative to fall and pooled samples. Greater spatial autocorrelation during the summer could be attributed to capturing related individuals (females with yearlings) at the same hair snag, thereby inflating autocorrelation. However, the proportion of individuals that could be related as parent/offspring (rij > 0.5) was similar among hair and tissue samples within the first distance class (hair: range = 0.07 -0.08; tissue: range = 0.06-0.08; Table 2.1). Lower autocorrelation values during the fall may be attributed to misreporting of harvest locations by hunters. However, this is unlikely due to our large sample sizes for all tissue samples (N = 248-407) and after performing bootstrap analyses (N = 121-177). We see the same spatial genetic pattern for either data set (Figure 2.4). Alternatively, a possible explanation for the lower autocorrelation values observed during the fall involves sampling individuals during seasonal movements associated with food resources, which has been previously observed in black bears (Noyce & Garshelis 2011). 26 Low statistical power may limit our ability to precisely measure fine scale spatial structure (Kelly et al. 2010). Limited DNA quantity of the summer hair samples restricted our analyses to five microsatellite loci. However, studies examining the effects of sample size on spatial genetic structure suggest the product of K (total number of alleles) and n (sample size) should be greater than 2,000 to obtain sufficient statistical power (Epperson 2005). The total number of alleles in our analysis was K = 50. Therefore, the product Kn = 6,050 for a simulated population of n = 121 (smallest sample size employed in our analysis) indicates our multilocus data set is adequate to detect and compare annual and seasonal patterns in spatial genetic structure in Michigan black bears. Local genetic spatial autocorrelation analyses revealed spatial genetic structure is geographically localized and consistent across seasons and years. Plots superimposed over the study area revealed three areas characterized by a larger number of clusters of genetically similar individuals. Essentially, the entire NLP is accessible by motorized vehicles and bear density is cited as the primary reason hunters choose a harvest site (Frawley 2001; Frawley 2009). Thus, if we assume sampling density is a surrogate of black bear density (Moore et al. 2014), sampling areas with high genetic relatedness occur in the areas of high black bear density (Figs. 2 and 5). Spatial heterogeneity in degree of spatial genetic clustering of related individuals likely reflects source lineages and sink habitats where immigration reduces individual relatedness in sink environments (Dias 1996). Spatial genetic structure suggest restricted dispersal in NLP black bears; however, limited statistical power precluded further exploration of spatial genetic structure using more data intensive approaches (e.g. clustering assignment; Prichard et al. 2000; Guillot et al. 2005). It is uncertain what is driving genetic heterogeneity in the NLP. The isolation by distance pattern we 27 observed is common in many continuously distributed species (Kopatz et al. 2012; Rueness et al. 2003a; Wasserman et al. 2013). However, an increasing number of landscape genetic studies have found landscape features can impede or facilitate gene flow (Beier & Noss 1998; Coulon et al. 2004; Funk et al. 2005; McRae & Beier 2007; Murphy et al. 2010). MANAGEMENT IMPLICATIONS Few empirical population genetic or landscape genetic studies have explored the impact of sample design (Anderson et al. 2010), specifically how sampling might affect measures of spatial genetic structure. Genetic sampling of rare and elusive species may preclude sampling at adequate spatial and temporal scales, and studies often rely on data collected over extended periods of time. Results of this study show significant positive spatial autocorrelation was concordant across seasons and years indicating the overall pattern of spatial genetic structure would be detected regardless of season, year, or collection methodology. However, our findings, while supported in Michigan black bears, should be evaluated empirically for other species, habitats and locales. Combining samples from different temporal scales that capture seasonal movements or overlapping generations could result in temporal instability in allele frequencies (Ryman 1997). Therefore, conservation management research should conduct analyses as described here before samples are aggregated indiscriminately. ACKNOWLEDGMENTS We thank the many MDNR staff for collecting summer hair samples and Michigan bear hunters for providing DNA from harvested bears. We are grateful to K. Filcek, B. Dreher, and S. Libants for assisting in data generation. We also thank J. Messina, B. Epperson, K. Holecamp, J. 28 Moore, J. Kanefsky, J. McGuire, for insightful discussion and helpful comments on the analysis and manuscript. Support for this project was provided by the Michigan Department of Natural Resources through the Wildlife and Sportfish Restoration Program F11AF00640 and by the Department of Fisheries and Wildlife at Michigan State University. 29 APPENDIX 30 APPENDIX Figure A1. Two-dimensional local spatial autocorrelation analysis of black bear microsatellite data using (a) three, (b) seven and (c) nine nearest neighbors. Circles represent individuals with positive significant local spatial autocorrelation (lr). Size of the circles indicates the strength of the autocorrelation coefficient. 2003 Summer 2003 Fall a .. b c 31 Figure A1 (cont’d) 2005 Summer 2005 Fall a .. b c 32 Figure A1 (cont’d) 2009 Summer 2009 Fall a .. b .. c 33 REFERENCES 34 REFERENCES Amstrup SC, Beecham J (1976) Activity patterns of radio-collared black bears in idaho. The Journal of Wildlife Management 40, 340-348. Anderson CD, Epperson BK, Fortin MJ, et al. (2010) Considering spatial and temporal scale in landscape-genetic studies of gene flow. Molecular Ecology 19, 3565-3575. Anselin L (1995) Local indicators of spatial association - Lisa. Geographical Analysis 27, 93115. Beier P, Noss RF (1998) Do habitat corridors provide connectivity? Conservation Biology 12, 1241-1252. Brown SK, Hull JM, Updike DR, Fain SR, Ernest HB (2009) Black bear population genetics in california: signatures of population structure, competitive release, and historical translocation. Journal of Mammalogy 90, 1066-1074. Bunn AG, Urban DL, Keitt TH (2000) Landscape connectivity: A conservation application of graph theory. Journal of Environmental Management 59, 265-278. Carter NH, Brown DG, Etter DR, Visser LG (2010) American black bear habitat selection in northern Lower Peninsula, Michigan, USA, using discrete-choice modeling. Ursus 21, 57-71. Costello CM (2010) Estimates of Dispersal and Home-Range Fidelity in American Black Bears. Journal of Mammalogy 91, 116-121. Coulon A, Cosson JF, Angibault JM, et al. (2004) Landscape connectivity influences gene flow in a roe deer population inhabiting a fragmented landscape: an individual-based approach. Molecular Ecology 13, 2841-2850. De Barba M, Waits LP, Genovesi P, et al. (2010) Comparing opportunistic and systematic sampling methods for non-invasive genetic monitoring of a small translocated brown bear population. Journal of Applied Ecology 47, 172-181. Dias PC (1996) Sources and sinks in population biology. Trends in Ecology & Evolution 11, 326-330. Dixon PM (2002) Ripley's K function. In: Encyclopedia of Environmetrics, 2 (eds. EIShaarawi AH, Piegorsch WW), pp. 1796-1803. John Wiley & Sons, Ltd., Chichester. Double MC, Peakall R, Beck NR, Cockburn A (2005) Dispersal, philopatry, and infidelity: Dissecting local genetic structure in superb fairy-wrens (Malurus cyaneus). Evolution 59, 625-635. 35 Dreher BP, Winterstein SR, Scribner KT, et al. (2007) Noninvasive estimation of black bear abundance incorporating genotyping errors and harvested bear. Journal of Wildlife Management 71, 2684-2693. Epperson BK (2005) Estimating dispersal from short distance spatial autocorrelation. Heredity 95, 7-15. Epperson BK, Huang Z, Li TQ (1999) Measures of spatial structure in samples of genotypes for multiallelic loci. Genetical Research 73, 251-261. Frawley BJ (2001) 2000 Michigan black bear hunter survey, p. 14. MDNR, Wildl. Div. Report. Frawley BJ (2009) 2008 Michigan black bear hunter survey, p. 24. MDNR, Wildl. Div. Funk WC, Blouin MS, Corn PS, et al. (2005) Population structure of Columbia spotted frogs (Rana luteiventris) is strongly affected by the landscape. Molecular Ecology 14, 483-496. Garant D, Dodson JJ, Bernatchez L (2000) Ecological determinants and temporal stability of the within-river population structure in Atlantic salmon (Salmo salar L.). Molecular Ecology 9, 615-628. Garshelis DL, Pelton MR (1981) Movements of black bears in the Great Smoky Mountains National Park. The Journal of Wildlife Management 45, 912-925. Goudet J (1995) FSTAT (Version 1.2): A computer program to calculate F-statistics. Journal of Heredity 86, 485-486. Goudet J (2001) FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3). Guillot G, Estoup A, Mortier F, Cosson JF (2005) A spatial statistical model for landscape genetics. Genetics 170, 1261-1280. Hood GM (2010) PopTools version 3.2.5 Available on the internet. Kelly RP, Oliver TA, Sivasundar A, Palumbi SR (2010) A method for detecting population genetic structure in diverse, high gene-flow species. Journal of Heredity 101, 423-436. Kopatz A, Eiken HG, Hagen SB, et al. (2012) Connectivity and population subdivision at the fringe of a large brown bear (Ursus arctos) population in North Western Europe. Conservation Genetics 13, 681-692. Latch EK, Rhodes E (2006) Evidence for bias in estimates of local genetic structure due to sampling scheme. Animal Conservation 9, 308-315. Lee DJ, Vaughan MR (2003) Dispersal Movements by Subadult American Black Bears in Virginia. Ursus 14, 162-170. Manville AM (1987) Den selection and use by black bears in Michigan’s Northern Lower 36 Peninsula. International Conference of Bear Research and Management 7, 317-322. McRae BH, Beier P (2007) Circuit theory predicts gene flow in plant and animal populations. Proceedings of the National Academy of Sciences of the United States of America 104, 19885-19890. Moore JA, Draheim HM, Etter D, Winterstein S, Scribner KT (2014) Application of large-scale parentage analysis for investigating natal dispersal in highly vagile vertebrates: a case study of american black bears (Ursus americanus). Plos One 9, e91168. Murphy MA, Dezzani R, Pilliod DS, Storfer A (2010) Landscape genetics of high mountain frog metapopulations. Molecular Ecology 19, 3634-3649. Noyce K, Garshelis D (2011) Seasonal migrations of black bears (Ursus americanus): causes and consequences. Behavioral Ecology and Sociobiology 65, 823-835. Paetkau D (2003) An empirical exploration of data quality in DNA-based population inventories. Molecular Ecology 12, 1375-1387. Paetkau D, Calvert W, Stirling I, Strobeck C (1995) Microsatellite analysis of populationstructure in Canadian polar bears. Molecular Ecology 4, 347-354. Paetkau D, Waits LP, Clarkson PL, Craighead L, Strobeck C (1997) An empirical evaluation of genetic distance statistics using microsatellite data from bear (Ursidae) populations. Genetics 147, 1943-1957. Peakall R, Smouse PE (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes 6, 288-295. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155, 945-959. Proctor MF, McLellan BN, Strobeck C, Barclay RMR (2005) Genetic analysis reveals demographic fragmentation of grizzly bears yielding vulnerably small populations. Proceedings of the Royal Society B-Biological Sciences 272, 2409-2416. Queller DC, Goodnight KF (1989) Estimating relatedness using genetic markers. Evolution 43, 258-275. R_Development_Core_Team (2011) R: A Language and Enviroment for Statistical Computing (ed. Team RDC). R Foundation for Statistical Computing, Vienna, Austria. Rabinowitz A, Zeller KA (2010) A range-wide model of landscape connectivity and conservation for the jaguar, Panthera onca. Biological Conservation 143, 939-945. Raymond M, Rousset F (1995) GENEPOP (VERSION-1.2) - Population-genetics software for exact tests and ecumenicism. Journal of Heredity 86, 248-249. 37 Rice WR (1989) Analyzing tables of statistical tests. Evolution 43, 223-225. Rogers LL (1977) Social relationships, movements, and population dynamics of black bears in northeastern Minnesota, University of Minnesota. Rogers LL (1987) Effects of food-supply and kinship on social-behavior, movements, and population-growth of black bears in northeastern Minnesota. Wildlife Monographs, 1-72. Rowlingson BS, Diggle PJ (1993) SPLANCS: spatial point pattern analysis code in S-Plus. Comput. Geosci. 19, 627-655. Rueness EK, Jorde PE, Hellborg L, et al. (2003a) Cryptic population structure in a large, mobile mammalian predator: the Scandinavian lynx. Molecular Ecology 12, 2623-2633. Rueness EK, Stenseth NC, O'Donoghue M, et al. (2003b) Ecological and genetic spatial structuring in the Canadian lynx. Nature 425, 69-72. Ryman N (1997) Minimizing adverse effects of fish culture: understanding the genetics of populations with overlapping generations. ICES J. Mar. Sci. 54, 1149-1159. Schwartz CC, Franzmann AW (1992) Dispersal and survival of subadult black bears from the Kenai Peninsula, Alaska. Journal of Wildlife Management 56, 426-431. Schwartz MK, Mills LS, Ortega Y, Ruggiero LF, Allendorf FW (2003) Landscape location affects genetic variation of Canada lynx (Lynx canadensis). Molecular Ecology 12, 18071816. Silverman BW (1986) Density estimation for statistics and data Analyses. In: Monographs on Statistics and Applied Probability. Chapman and Hall, London. Taberlet P, Camarra JJ, Griffin S, et al. (1997) Noninvasive genetic tracking of the endangered Pyrenean brown bear population. Molecular Ecology 6, 869-876. Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes 4, 535-538. Wasserman TN, Cushman SA, Littell JS, Shirk AJ, Landguth EL (2013) Population connectivity and genetic diversity of American marten (Martes americana) in the United States northern Rocky Mountains in a climate change context. Conservation Genetics 14, 529541. White PJ, Proffitt KM, Mech LD, et al. (2010) Migration of northern Yellowstone elk: implications of spatial structuring. Journal of Mammalogy 91, 827-837. Wilberg MJ, Dreher BP (2004) GENECAP: a program for analysis of multilocus genotype data for non-invasive sampling and capture-recapture population estimation. Molecular Ecology Notes 4, 783-785. 38 Williams BW, Scribner KT (2010) Effects of multiple founder populations on spatial genetic structure of reintroduced American martens. Molecular Ecology 19, 227-240. Zar JH (1996) Biostatistical analysis, 3rd edn. Prentice-Hall, Upper Saddle River, New Jersey, USA. 39 CHAPTER 3 SPATIAL GENETIC STRUCTURE AND LANDSCAPE CONNECTIVITY IN BLACK BEARS: INVESTIGATING THE SIGNIFICANCE OF DIFFERENT LAND COVER DATA SETS IN LANDSCAPE GENETICS ANALYSES. Draheim, H. M., J. A. Moore, S. R. Winterstein, K. T. Scribner ABSTRACT Landscape genetic analyses allow detection of fine scale spatial genetic structure (SGS) and quantification of effects of landscape features on gene flow and connectivity. Typically, analyses require generation of resistance surfaces that are hypothesized relationships between an organism’s movement and features of the environment. How accurately hypotheses predict associations is determined in large part on 1) the landscape features used, 2) the resistance values assigned to features, and 3) how accurately resistance surfaces represent actual landscape permeability. Our objectives were to perform individual-based analyses to examine fine-scale spatial genetic patterns and the influence of landscape features on functional connectivity (i.e. gene flow) of black bears (Ursus americanus) in the Northern Lower Peninsula (NLP) of Michigan, USA. The relationship between gene flow and landscape features was evaluated using two different land cover data sets to investigate the relative importance of land cover classification and accuracy for landscape resistance model performance. We genotyped 365 individuals (from 2006) at 12 microsatellite loci, and used spatial genetic autocorrelation and a Bayesian approach to infer SGS. Black bears in Michigan’s NLP are not genetically homogenous and land cover was significantly correlated with genetic distance. However, despite using the same ecological model and biological assumptions to generate resistance surfaces, we observed 40 differences in model performance when different land cover data sets were used. Results indicate the importance of considering land cover uncertainty in landscape genetic studies. INTRODUCTION Understanding the mechanisms that dictate habitat occupancy and form movement patterns is a cornerstone to development of successful conservation, restoration, and management strategies. Natural ecosystems are spatially and temporally heterogeneous, such heterogeneity has important ecological implications for the organisms inhabiting them and dispersing through them (Bolliger et al. 2010; Spear et al. 2010; Storfer et al. 2007). As natural habitats are altered and fragmented, dispersal and colonization potential decreases (Kool et al. 2012), and patterns of gene flow and movement are likely to be disrupted, leading to increased spatial genetic structure, reduced abundance, and a decreased likelihood of population persistence (Frankham et al. 2002). For example, species that have large space requirements may be particularly susceptible to habitat alternation; therefore habitat and connectedness can be important (Opdam & Wascher 2004). Identifying the factors that influence the degree of population connectivity across heterogeneous landscapes provides insight into population dynamics and the evolutionary trajectory of a species (Segelbacher et al. 2010). For example, wide ranging species such as the black bear (Ursus americanus) can be sensitive to land cover (Pelton 1992). Interruptions in forest cover and prevalence of human dominated landscapes have been shown to impede black bear dispersal (Cushman 2006; Short Bull et al. 2011). Understanding the synergistic effects of land-use/land-cover on species movements and abundance will enable managers to target regions or habitat types that are important to maintain connectivity across anthropogenically-altered habitats and ultimately assess the impacts of future change. Landscape genetics seeks to identify 41 how landscape elements and environmental factors influence the spatial distribution of genetic variation. Landscape geneticists face significant challenges to understand how animals move through landscapes. Landscape genetic analyses require generation of resistance surfaces that are hypothesized relationships between an organism’s movement and features of the environment. How accurately hypotheses predict associations is determined in large part on the 1) landscape features used, 2) the resistance values assigned to features, and 3) how accurately resistance surfaces represent actual landscape permeability. For example, land cover is an indicator of ecological condition (e.g. available resources) and a common predictor variable. Land cover resistance surfaces are commonly constructed using freely-available land cover data layers (i.e. National Land Cover Database (NLCD), Global Land Cover (GLC) database and Moderate Resolution Imaging Spectroradiometer (MODIS) land cover (Bartholomé & Belward 2005; Friedl et al. 2002; Homer et al. 2007). However, these data sets differ in the degree of uncertainty. How representative a land cover data set portrays reality may result in markedly different resistance surfaces (Liu et al. 2004) and ability to perform hypothesis testing. Uncertainty in land cover data is widely recognized in the remote sensing community (Foody 2002), but there has been little discussion of potential implications within landscape genetic literature. In light of this uncertainty practitioners need to understand the consequences of using different land cover datasets of varying uncertainty and how it may affect landscape genetic analyses. We investigated the effects of the landscape on functional connectivity of black bears in Michigan’s Northern Lower Peninsula. Our objectives were to quantify fine-scale spatial genetic structure (SGS) and the influence of landscape features (roads, river, and land cover) on 42 connectivity of black bears in the NLP. We first used Bayesian clustering methods to identify population subdivision and spatial autocorrelation to examine patterns of SGS. We then performed individual-based landscape genetic analyses to examine whether landscape permeability more accurately explained genetic distance than Euclidean distance. Inconsistencies among resistance surfaces generated by different land cover data sets may ultimately influence landscape genetic model performance. Therefore, is important to understand the consequences of choosing one data layer over another. Thus, we compared landscape genetic models using the two land cover data sets. METHODS Sampling, DNA extraction, PCR amplification Tissue samples (teeth) (N= 365) were collected from harvested bears registered at hunter check stations during the fall harvest (September-October) of 2006 in Michigan’s Northern Lower Peninsula (NLP, 47,739 km2; Figure 3.1). Bear harvest locations were recorded as township, range, section (e.g. 2.59 km2). Coordinates of harvest locations were georeferenced using the section centroid and converted into UTM coordinates. We extracted DNA from bear teeth using Qiagen DNEasy Tissue Kits (Qiagen Inc., Valencia, CA) following manufacturer protocols. DNA was quantified using a Nanodrop spectrophotometer (Thermo Scientific, Waltham, MA) and diluted to a 20 ng/μl working concentration. We used PCR to amplify 12 variable microsatellite loci including: G10X, G10L, G10D, G10B, G10M (PCR annealing temperature TA = 58°C; Paetkau et al. 1995) UarMU59, UarMU50 (TA = 58°C, Taberlet et al. 1997), ABB1, ABB4 (TA = 54°C; Wu et al. 2010), UT29, UT35, and UT38 (TA = 54°C; Shih et al. 2009). DNA was amplified using a PCR profile with 43 Figure 3.1. Study area in the Northern Lower Peninsula of Michigan (NLP) showing locations of black bear harvest samples collected during 2006 (N = 365). Pie graphs were constructed using on individual posterior probabilities of each cluster (K = 2) from Program STRUCTURE. 44 the following steps: initial denaturation 94°C for four minutes, followed by 25-42 cycles of 30s at 94°C, 30s at locus specific annealing temperature, 1 min at 72°C, then an additional 10 min extension step at 72°C. Ten microliter reactions were prepared using 40 ng of DNA in 10 mM Tris-HCl, 50 mM KCl, 2.0 mM MgCl2, 0.2 mM of each dNTP, 1 pmol of each primer, and 0.5 units of Taq polymerase. Amplified products were sized on 6.5 % denaturing acrylamide gels for electrophoresis and visualized on a LI-COR 4200 Global IR2 System (LI-COR Inc., Lincoln, NE). All alleles were scored independently by two experienced lab personnel using SAGA genotyping software (LI-COR Inc., Lincoln, NE). To assess genotyping error 10% of samples were randomly selected and genotyped twice to yield a genotyping error rate of < 2%. Population genetic analysis Program MICRO-CHECKER (Van Oosterhout et al. 2004), was used to test for the presence of null alleles and allelic dropout. We tested for deviations from Hardy-Weinberg equilibrium using program GENEPOP (Version 3.1d; Raymond and Rousset 1995), and used sequential Bonferroni tests to correct for multiple tests (Rice 1989). We used Bonferroni corrections (Goudet 1995, 2001) to test for linkage disequilibrium using program FSTAT 2.93. We quantified microsatellite genetic diversity using mean number of alleles (A), observed heterozygosity (HO), and expected heterozygosity (HE) over all loci, using program GenAlEx (v. 6.0, Peakall and Smouse 2006,). Spatial genetic structure was assessed using a global genetic spatial autocorrelation coefficient (r) calculated using program GenAlEx (Peakall and Smouse 2006, v. 6.0). Coefficient r (range, - 1 to + 1) was calculated by comparing pairwise geographical inter-individual and genetic distance. We used program STUCTURE v. 2.3.4 (Pritchard et al. 2000) to evaluate population genetic structure. We estimated the number of genetic clusters (K) based on Hardy- 45 Weinberg expectations using no geographic location information and the posterior probability of each individual belonging to each cluster. We performed 10 independent runs of K = 1-10 using simulations of 2 x 106 iterations after a burn-in period of 5 x 105 Markov chain Monte Carlo (MCMC) iterations. The most likely number of clusters was determined by the log likelihood of K and the posterior probability of K (P(K|X)) as determined by the method described in Pritchard et al. (2000) and estimates of delta K (Dk) (Evanno et al. 2005) using program STRUCTURE HARVESTER v. 06.93 (Earl 2012). To visualize the spatial distribution of the genetic clusters we used ArcGIS 10.1 to plot posterior probability of cluster membership for each individual to each cluster. To measure genetic differentiation among individuals, we calculated the ratio proportion of shared alleles (Dps; Bowcock et al. 1994) for each pairwise combination of individuals using GenAlEx ver. 6.0 (Peakall & Smouse 2006). Landscape Genetic Analysis We used two land use/land cover digital coverage maps derived from Landstat TM imagery; 1) Michigan Resource Information Systems (MIRIS) Integrated Forest Monitoring and Assessment Prescription Project (IFMAP) land cover data and 2) NOAA’s Coastal Change Analysis Program (CCAP) Land Cover dataset. The IFMAP data set was categorized into 34 land cover classes (Table B1). Reported accuracy assessment (based on 2817 reference points) for IFMAP was 77% for overall accuracy and 68% for division among major forest types (i.e. mixed, deciduous, coniferous). Michigan Department of Natural Resources (MDNR) defined a stand as coniferous, deciduous, or mixed using a 60% purity rule. That is, a forest stand is considered of mixed composition unless 60% of the forest stand is either dominated by coniferous or deciduous trees. In contrast, CCAP land cover data was classified in to 25 land cover classes (Table B1) with an 46 overall accuracy of 87% and 83% for major forest types based on 900 reference points. NOAA used a much higher purity requirement (75%) for defining coniferous and deciduous forest types resulting in larger of proportion of forest stands to be classified as mixed forest (Figure 3.2). In order to understand how genetic differentiation was influenced by landscape features, resistance surfaces were generated using Spatial Analyst in ARCGIS 10.1 by weighing each land cover class according to positive or negative associations with black bear presence based on habitat suitability estimates derived by Carter et al. (2010). Landscape resistance surfaces were based on three factors that have previously been reported to influence habitat selection by black bears in the NLP (Carter et al. 2010). These included roads, rivers, and land cover. The first set of models included resistance surfaces using a fine-scale land cover classification from the IFMAP data set (resolution = 150 m). We reclassified the IFMAP land cover data according to bear habitat suitability (cost = 1-100, least to most resistant to bear movement): northern hardwood mixed forest (NHM, 1), aspen (A, 1), forested wetland (FWL, 1), pine (P, 10), oak (O, 20), non-forested upland (NFU, 20), agriculture (AG, 50), non-forested wetland (NFW, 100), and developed (DEV, 100). Our second set of land cover models included resistance surfaces based on the CCAP Land Cover dataset (resolution = 150 m) which we reclassified based on bear habitat suitability (1-100, least to most resistant to bear movement): mixed deciduous forest (MF, 1), forested wetland (FW, 1), evergreen forest (EF, 10), non-forested upland (NFU, 20), agriculture (AG, 50), non-forested wetland (NFW, 100), and developed (DEV, 100). We then modeled resistance as a function of land cover according to two alternative hypotheses using a coarser classification of land cover: 1) bear habitat (significantly positively correlated with bear presence; MF, FW or non-habitat (either positively associated but not significant or significantly negatively correlated with bear presence; EF, NFU, AG, NFW, DEV; 100), and 2) forested areas 47 Figure 3.2. Side by side comparison of CCAP and IFMAP land cover maps using black bear classification scheme described in Carter 2007. (A) Land cover maps for the entire study area. (B) Land cover of the same 7 x 10 km area. IFMAP CCAP A B 48 (MF, FW, EF; 1) or non-forested areas (NFU, AG, NFW, DEV; 100). Major rivers (10) and all roads (I-75, 100; state roads, 50; all other roads, 10; weight based on traffic patterns) were included as predictor variables as potential physical barriers to dispersal. We proposed 17 alternative landscape hypotheses (Table 3.1). Using raster algebra in Spatial Analyst ARCGIS 10.1 we combined each land cover classification with roads, rivers, and roads + rivers. To ensure resistance distances are not unduly influenced upon the relative cost values assigned to land cover types (Rayfield et al. 2010), we performed a sensitivity analysis on the range of relative cost values between resistance surface land cover types by comparing resistance distances between rasters weighted using large differences (i.e. maximum weights = 1, 10, 20, 50, 100) and small differences (minimum weights = 0.1, 1, 2, 5, 10) in the relative resistance values for land cover types. Pairwise resistance distances among individuals were calculated using program CIRCUITSCAPE 3.5 (McRae 2006). CIRCUITSCAPE incorporates circuit theory that evaluates the total landscape resistance between individuals via multiple potential paths of resistance (McRae & Beier 2007). We used Mantel (Mantel 1967) and partial Mantel tests (Smouse et al. 1986) to quantify the associations between genetic and geographic/resistance distances. Legendre and Fortin (2010) noted that Mantel and partial Mantel testing is an appropriate framework when hypotheses are explicitly defined in terms of distances, and are appropriate for individual-based analyses. We tested for the effect of geographic distance alone using a resistance distance in CIRCUITSCAPE with a raster surface value of 1 bounded the same as all other surfaces then partial out the pairwise distance only matrix to assess the influence of landscape features. The significance of the matrix correlation was evaluated by permutation test, based on 10,000 randomizations. The top candidate model corresponds to the largest partial 49 Table 3.1. Mantel (r, IBD) and partial Mantel correlations (r, all resistance models) between geographic and genetic pairwise distances among 365 individual black bear in the NLP. Bold indicates significant correlations. Mantel and Partial Mantel test r P IBD I75 only State Roads only Roads Rivers Roads + Rivers CCAP HSI land cover only HSI + Roads HSI + Rivers HSI + Roads + River Hab/Non-Hab land cover only Hab/Non-Hab + Roads Hab/Non-Hab + River Hab/Non-Hab + Roads + River Forest/Non Forest Forest/Non Forest + Roads Forest/Non Forest + River Forest/Non Forest + Roads + River IFMAP HSI land cover only HSI + Roads HSI + Rivers HSI + Roads + River Hab/Non-Hab land cover only Hab/Non-Hab + Roads Hab/Non-Hab + River Hab/Non-Hab + Roads + River Forest/Non Forest Forest/Non Forest + Roads Forest/Non Forest + River Forest/Non Forest + Roads + River 0.052 0.008 0.019 0.025 0.011 0.029 0.004 0.051 0.480 0.157 0.334 0.131 0.044 0.041 0.053 0.056 0.008 0.033 0.015 0.041 0.015 0.037 0.036 0.041 0.050 0.081 0.039 0.021 0.334 0.082 0.226 0.075 0.130 0.091 0.092 0.076 0.017 0.020 0.023 0.012 -0.011 -0.013 -0.012 -0.014 0.003 -0.009 -0.009 -0.009 0.224 0.202 0.166 0.297 0.679 0.712 0.695 0.712 0.402 0.676 0.673 0.676 CCAP = Costal Change Analysis Program land cover data set. IFMAP = Integrated Forest Monitoring and Assessment Prescription Project land cover data set. HIS = Habitat Suitability Index, designating models in which land cover is reclassified using Carter et al. 2010 ecological model. Hab = Habitat 50 Mantel r values. Potential corridors of for black bear gene flow were generated using cumulative current map option in CIRCUITSCAPE. Next, to eliminate the potential for Type I error and spurious correlations we compared a set of models against one another rather than the null model of isolation by distance (IBD). We used an additional robust causal modeling frame work proposed by Cushman et al. (2013) that is based on the relative support (RS) of each candidate model. The relative support (RS) of each model is the difference in the partial Mantel correlation after partialling out the resistance distance of the other model such that: RS1|2 = (GD ~ RD1|RD2) – (GD ~ RD2|RD1) where, GD is the genetic distance matrix and RD is the resistance distance matrix from model. The model with significant and positive RS in all comparisons is considered the best candidate model (Cushman et al. 2013). RESULTS Genetic Diversity and Spatial Genetic Structure We found no evidence for null alleles or allelic dropout, and no loci were found to deviate significantly from Hardy-Weinberg or linkage equilibrium, so all 12 loci were retained for further analyses. For all loci expected heterozygosity ranged from 0.71 to 0.91 and number of alleles per locus ranged from 6 to 26. The NLP black bear population is not genetically homogenous. Global spatial autocorrelation analysis revealed that genetically similar individuals were not randomly distributed. Correlograms indicated a pattern of isolation by distance (Figure 3.3). Significant positive genetic autocorrelation was observed for individuals sampled at inter-individual 51 Figure 3.3. Global spatial autocorrelation of black bear samples. Bars represent standard errors. The 95% confidence intervals (dashed lines) about the null hypothesis of no spatial structure were estimated by permutations based on distance classes of 5 km. 52 distances of 0-30 km (Figure 3.3). Our STRUCTURE results indicated the highest average loglikelihood value (-15714.50) was observed for K = 4. However, using Delta K as suggested by Evanno et al (2005) the greatest support was for K = 2 which suggests the NLP consist of two genetic groupings. Distribution of individual posterior probabilities indicated a west to east gradient (Figure 3.1). Landscape Genetic Analysis The best models based on Mantel correlations included land cover weighted according to the habitat suitability model from Carter et al. 2010 (LCHSI) for both CCAP and IFMAP land cover data sets. However, after partialling out Euclidean distance only the CCAP models were significant (Table 3.1). The model that included LCHSI, rivers, and roads was the best supported model (r = 0.56, P = 0.021; Tabel 3.1). Equally well supported were two reduced models that included LCHSI only (r = 0.044, P = 0.05) and LCHSI with rivers (r = 0.053, P = 0.039). Similarly, the best multivariate model based on RS included LCHSI, rivers, and roads. However, the model with the highest RS was not significant after partialling out each of the 2 other alternative models. Therefore, we were unable to conclude which single model is the top model. Indeed, the cumulative resistance “heat” maps between pairs of individuals from the three competitive isolation-by-resistance models appear similar (Figure 3.4). It should be noted that all competitive models included land cover confirming that land cover is extremely important for black bear connectivity in the NLP. None of the models that classified land cover using habitat/non-habitat or forest/non-forest classes were significant after partialling out Euclidean distance. 53 Figure 3.4. Comparison of cumulative resistance maps between pairs of individuals from three competitive isolation-by-resistance models in Table 3.1: 1) land cover only; 2) land cover and rivers; and 3) land cover, rivers, and roads generated from CCAP land cover using program CIRCUITSCAPE. Gradient of colors indicated the probability of black bear movement. Yellow colors indicate high probability; red and pink colors indicate medium probability; and purple and blue representing low probability. A B C High Low Black bear spatial genetic structure 54 DISCUSSION We evaluated fine-scale genetic structure within a continuously distributed black bear population inhabiting a heterogeneous landscape. Black bears in the NLP exhibited significant positive spatial genetic autocorrelation for distances up to 30 km, such spatial structuring may be interpreted under the concept of IBD (Coulon et al. 2004; Wright 1943). Positive spatial autocorrelation is expected to accrue due to restricted dispersal and beyond these distance classes individuals are no more related than expected by chance (Epperson 2005). Positive spatial genetic autocorrelation over short distances is likely attributed to male-biased dispersal and female natal philopatry, commonly exhibited in black bears, where female offspring establish home ranges adjacent to the mothers whereas male offspring disperse from the natal area (Costello 2010; Moore et al. 2014; Rogers 1977; Rogers 1987; Schwartz & Franzmann 1992). Indeed, black bear dispersal in Michigan’s NLP is strongly male biased. IBD is consistent with studies on bears and other wide ranging carnivores (Brown et al. 2009; Paetkau et al. 1997; Rueness et al. 2003b). Our STRUCTURE results indicate that two genetically distinct groups exist in the NLP, defined as a western and an eastern genetic cluster. The distinction between the two genetic clusters was modest. When the individual cluster membership probabilities are plotted over the study area we observed a gradual west-east pattern with a number of individuals with mixed ancestry in the middle of the study area and a number of “western” individuals scattered in the eastern portion and vice versa (Figure 3.1). It is uncertain what is driving genetic heterogeneity in the NLP. Running proximal to the primary shift from western to eastern individuals is Interstate-75 indicating major roads could be a barrier to dispersal. However, we found I-75 and/or state roads alone did not significantly explain genetic distance (Table 3.1). Other studies 55 have shown major roads to impede dispersal in both black bears (Lee & Vaughan 2003; Thompson et al. 2005) and brown bears (Proctor et al. 2005). Effect of landscape features on black bear connectivity In our study area isolation by landscape resistance was more strongly supported than IBD. Our results using univariate landscape models suggest that land cover was the landscape factor most strongly correlated with genetic distance for black bears in NLP followed by roads and rivers, with the latter two variables not significantly influencing genetic distance after partialling out geographic distance. The importance of land cover in our models reflects habitats required to meet nutritional needs, and land cover largely reflects food abundance and available refugia (Amstrup & Beecham 1976; Davis et al. 2006; Noyce & Garshelis 2011; Rogers 1987). For example, understory communities of deciduous and mixed forest provides essential food resources for black bear diets such as aspen catkins and leaves, wild berries and hard mast (e.g. acorns). However, land cover was only identified as significantly associated with genetic differentiation when classified using the habitat suitability model from Carter et al. (2010) (Table 3.1). That is, when land cover was classified using two broad definitions:1) habitat/non-habitat or 2) forest/non-forest land cover was not significant after partialling out distance. Our results are inconsistent with previous studies that have found genetic structure related to land cover when classified broadly as forest/non- forest (Cushman et al. 2006; Short Bull et al. 2011). One possible explanation for this disparity is the difference in spatial distribution of food availability between study areas (Rocky Mountains vs. NLP). Unlike conifer dominated forest types that occur in the Rocky Mountains of the United States, Michigan’s NLP forests are a heterogeneous 56 mix of deciduous, coniferous, and mixed stands with different understory communities and thus food availabilities. Comparisons between the competing landscape models parameterized using multiple variables identified significant partial Mantel correlation after partialling out distance for three equally supported models: 1) land cover only; 2) land cover and rivers; and land cover, rivers, and roads. Visualization of potential dispersal pathways based on these three models revealed similar patterns of high, medium, and low probability of bear movement (Figure 3.4). Heat maps show two clear areas of high flow; a large patch located in the east central region and a smaller patch located in the south central region (near Houghton and Higgins Lake) of the study area. Connecting these two regions are many areas of medium flow suggesting dispersal within the NLP can occur via multiple pathways. Corridors, if adequate in arrangement and number, can in theory offset the negative consequences of population fragmentation and are necessary to maintain occupancy in a population that exhibits source-sink dynamics (Dixo et al. 2009). Black bears exhibit asymmetric dispersal among areas within the NLP suggesting loss of movement corridors could have wide ranging effects (Chapter 4). For example, long-term persistence of black bears in sink areas depends on the rate of extinction and the rate of movement between areas (Fahrig & Merriam 1994). Thus, loss of corridors by habitat alternation in the NLP can decrease local bear occupancy. Land cover dataset comparison We observed differences in best fit model when landscape resistances were weighted using different land cover data sets. In contrast to the CCAP models all landscape resistance models using IFMAP land cover data failed to improve model fit relative to the null model of isolation 57 by distance. These results are surprising because the spatial dispersion of land cover types appear very similar qualitatively (Figure 3.2A). However, on closer inspection, when the data sets are classified under the same land cover categories approximately 12% of raster cells differed (Figure 3.2B). Differences may be attributed to a number of factors. First, despite the potential of high-resolution satellite imagery to pick up fine-scale differences in land cover there is uncertainty. Often researchers take GIS land cover classification for granted and assume that classifications represent an accurate representation of physical land cover attributes when digitized land cover maps have inherent limitations. Uncertainty or error in remote sensing is primarily associated with land cover class misidentification (Foody 2002; Liu et al. 2004). Land cover classification can be subjective (Thomlinson et al. 1999) and can have a propensity for high error rates, despite recommended guidelines to ensure accuracy of classifications (e.g. overall accuracy of 85% with no class less than 70% accurate) (Congalton & Green 2008). For example, the same satellite images classified by two independent researchers but applying the same definition of land cover types can produce different land cover maps. The IFMAP and CCAP data sets differed considerably in their overall misclassifications (e.g. accuracy assessments; IFMAP > 77 %, CCAP > 85%), and propensity for misclassification increased when unraveling closely-related land cover types (e.g. accuracy assessment among forest types; IFMAP > 67 %, CCAP > 83%; Beier et al 2008; MDNR 2004; NOAA 2006) . Second, the data layers we used in our evaluation use remarkably different thresholds for stand dominance (IFMAP = 60%, CCAP = 75%) for classifying coniferous/deciduous forest cover types (MDNR 2004). For the IFMAP data, the MDNR lowered the thresholds to 60% because much of Michigan’s forests are mixed and they wanted to emphasis species predominance in cover type for forestry practices. Thereby reducing the number of forest cover types that fall within mixed 58 forest classes (MDNR 2004). Finally, the IFMAP data set was published 5 years before the CCAP data. Within the NLP there has been landscape modification from ongoing anthropogenic activities (primarily deforestation). However, we are unsure as to whether the landscape change is contributing to our results. Further investigation using time series data is need to address this question. While the disparity observed in the land cover differences may not seem extensive, errors can be amplified when parameterizing our models to reflect landscape resistance. For instance, our land cover data sets were generated using disparate forest stand purity requirements (IFMAP = 60%, CCAP = 75%). As a result, mixed forests are more likely to be misclassified as either deciduous or coniferous forest using IFMAP classification rule. The Carter et al. 2010 habitat suitability model we used in our analyses found black bears in the NLP are positively associated with aspen, forested wetlands, and northern hardwood/mixed forest (females only) therefore were assigned a low resistance value. In contrast, no significant association was observed for coniferous forest thus assigned a higher cost value. The intrinsic problem in this scenario would be if areas classified as coniferous but are in reality mixed forest (i.e. exhibit low resistance to bear movement) are then assigned an erroneous higher cost value resulting in spurious resistance surface. Had analyses been based on IFMAP land cover resistances surfaces we would conclude that land cover has no influence on functional connectivity of black bears in the NLP. Nevertheless, we detected a strong association between black bear SGS and land cover using the CCAP data set. Our results suggest minor differences in land cover classification schemes and data quality can lead to markedly different depictions of landscape resistance and can impact our conclusions. 59 CONCLUSIONS Landscape configuration and composition strongly influence population dynamics (Allendorf & Luikart 2012). Thus, identifying factors that influence the degree of population connectivity across heterogeneous landscapes can provide insight into population dynamics and the evolutionary trajectory of a species (Manel et al. 2010; Manel et al. 2003). Our results indicate black bear dispersal in Michigan’s NLP is impacted by land cover. Improved model performance when land cover data are categorized using a more refine classification scheme (more classes) suggests despite being habitat generalist forest/non-forest classifications are too broad for black bears in Michigan. However, our ability to relate genetic distance to landscape elements is largely dependent on hypothesized resistance surfaces. While there is no single way of generating resistance surfaces we strongly suggest that, when planning landscape genetic studies, researchers carefully evaluate which landscape variables they use to meet their objectives and important biological assumptions (Reding et al. 2013). However, we obtained spurious correlations using two different land cover data sets; weighted using the same ecological model, but differed greatly in their classification schemes and how accurately landscapes were represented (i.e. uncertainty). Our analyses suggest different land cover data sets can result in disparate or contradictory conclusions. With increasing availability of high-resolution satellite imagery there are a number of land cover data layers from which to choose. However, agencies differ in the way in which land cover data is classified; thus there are inherent inconsistencies among them. Our results emphasize that landscape genetic researchers should not only carefully consider which landscape variables to use when constructing resistances surfaces, but also, the accuracy and classification schemes used to generate the underlying datasets. 60 ACKNOWLEDGMENTS We are grateful to J. Fierke, K. Filcek, B. Dreher, and S. Libants for assisting in data generation. We also thank J. Messina, B. Epperson, K. HoleKamp, J. Kanefsky, J. McGuire, Robert Goodwin for insightful discussion and helpful comments on the analysis and manuscript. All samples used in analyses were provided by Michigan bear hunters and collected by the many MDNR staff at registration stations. Support for this project was provided by the Michigan Department of Natural Resources through the Wildlife and Sportfish Restoration Program F11AF00640 and by the Department of Fisheries and Wildlife at Michigan State University. 61 APPENDIX 62 APPENDIX Table B1. Side by side comparisons of land cover types according to IFMAP and CCAP classification schemes IFMAP and CCAP data the class number and descriptions pulled from original data sets. Class numbers and descriptions for the black bear ecological model were originated from Carter 2007. IFMAP Class 110 123 121 122 2111 2112 2113 222 310 320 330 350 411 412 413 414 419 421 423 429 Description Low Intensity Urban High Intensity Urban Airports Roads / Paved Non-vegetated Agriculture Row Crops Forage / Pasture Orchards / Vineyards / Nursery Herbaceous Openland Upland Shrub Low-Density Trees Parks / Golf Courses Northern Hardwoods Oak Type Aspen Type Other Upland Hardwoods Mixed Upland Deciduous Pines Other Upland Conifers Mixed Upland Conifers CCAP Class 3 2 2 2 20 6 6 6 8 12 12 8 9 9 9 9 9 10 10 10 Description Low-Intensity Developed High-Intensity Developed High-Intensity Developed High-Intensity Developed Bare Land Cultivated Land Cultivated Land Cultivated Land Grassland Scrub/Shrub Scrub/Shrub Grassland Deciduous Forest Deciduous Forest Deciduous Forest Deciduous Forest Deciduous Forest Evergreen Forest Evergreen Forest Evergreen Forest 63 Black Bear Ecological Model Class 1 1 1 1 2 2 2 2 3 3 3 1 4 5 6 4 4 7 7 7 Description Human Development Human Development Human Development Human Development Agriculture Agriculture Agriculture Agriculture Upland Nonforested Upland Nonforested Upland Nonforested Upland Nonforested Northern Hardwood and Mixed Forest Oak Forest Aspen Forest Northern Hardwood and Mixed Forest Northern Hardwood and Mixed Forest Pine Pine Pine Table B1 (cont’d) 431 500 611 612 613 621 622 623 629 710 720 730 790 Upland Mixed Forest Water Lowland Deciduous Forest Lowland Coniferous Forest Lowland Mixed Forest Floating Aquatic Lowland Shrub Emergent Wetland Mixed Non-Forest Wetland Sand / Soil Exposed Rock Mud Flats Other Bare /Sparsely Vegetated 11 21 13 13 13 22 14 15 15 20 20 19 20 Mixed Forest Water Palustrine Forested Wetland Palustrine Forested Wetland Palustrine Forested Wetland Palustrine Aquatic Bed Palustrine Shrub/Scrub Wetland Palustrine Emergent Wetland Palustrine Emergent Wetland Bare Land Bare Land Unconsolidated Shore Bare Land CCAP = Costal Change Analysis Program land cover data set. IFMAP = Integrated Forest Monitoring and Assessment Prescription Project land cover data set. 64 4 10 8 8 8 9 9 9 9 2 2 2 2 Northern Hardwood and Mixed Forest Water Forested Wetland Forested Wetland Forested Wetland Nonforested Wetland Nonforested Wetland Nonforested Wetland Nonforested Wetland Agriculture Agriculture Agriculture Agriculture REFERENCES 65 REFERENCES Allendorf FW, Luikart G (2012) Conservation and the Genetics of Populations, 2nd edn. WileyBlackwel. Amstrup SC, Beecham J (1976) Activity patterns of radio-collared black bears in idaho. The Journal of Wildlife Management 40, 340-348. Bartholomé E, Belward A (2005) GLC2000: a new approach to global land cover mapping from Earth observation data. International Journal of Remote Sensing 26, 1959-1977. Beier P, Majka DR, Spencer WD (2008) Forks in the road: choices in procedures for designing wildland linkages. Conservation Biology 22, 836-851. Bolliger J, Le Lay G, Holderegger R (2010) Landscape genetics - how landscapes affect ecological processes. Gaia-Ecological Perspectives for Science and Society 19, 238-240. Bowcock AM, Ruizlinares A, Tomfohrde J, et al. (1994) High resolution of human evolutionary threes with polymorphic microsatellites. Nature 368, 455-457. Brown SK, Hull JM, Updike DR, Fain SR, Ernest HB (2009) Black bear population genetics in california: signatures of population structure, competitive release, and historical translocation. Journal of Mammalogy 90, 1066-1074. Carter NH (2007) Predicting ecological and social suitability of black bear habitat in Michigan's Lower Peninsula, University of Michigan. Carter NH, Brown DG, Etter DR, Visser LG (2010) American black bear habitat selection in northern Lower Peninsula, Michigan, USA, using discrete-choice modeling. Ursus 21, 57-71. Congalton RG, Green K (2008) Assessing the accuracy of remotely sensed data: principles and practices CRC press. Costello CM (2010) Estimates of Dispersal and Home-Range Fidelity in American Black Bears. Journal of Mammalogy 91, 116-121. Coulon A, Cosson JF, Angibault JM, et al. (2004) Landscape connectivity influences gene flow in a roe deer population inhabiting a fragmented landscape: an individual-based approach. Molecular Ecology 13, 2841-2850. Cushman SA (2006) Effects of habitat loss and fragmentation on amphibians: A review and prospectus. Biological Conservation 128, 231-240. Cushman SA, McKelvey KS, Hayden J, Schwartz MK (2006) Gene flow in complex landscapes: Testing multiple hypotheses with causal modeling. American Naturalist 168, 486-499. Cushman SA, Wasserman TN, Landguth EL, Shirk AJ (2013) Re-evaluating causal modeling with Mantel tests in landscape genetics. Diversity 5, 51-72. 66 Davis H, Weir RD, Hamilton AN, Deal JA (2006) Influence of phenology on site selection by female American black bears in coastal British Columbia. Ursus 17, 41-51. Dixo M, Metzger JP, Morgante JS, Zamudio KR (2009) Habitat fragmentation reduces genetic diversity and connectivity among toad populations in the Brazilian Atlantic Coastal Forest. Biological Conservation 142, 1560-1569. Earl DA (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genetics Resources 4, 359-361. Epperson BK (2005) Estimating dispersal from short distance spatial autocorrelation. Heredity 95, 7-15. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software structure: a simulation study. Molecular Ecology 14, 2611-2620. Fahrig L, Merriam G (1994) Conservation of fragmented populations. Conservation Biology 8, 50-59. Foody GM (2002) Status of land cover classification accuracy assessment. Remote Sensing of Environment 80, 185-201. Frankham R, Ballou JD, Briscoe DA (2002) Introduction to conservation genetics Cambridge University Press. Friedl MA, McIver DK, Hodges JC, et al. (2002) Global land cover mapping from MODIS: algorithms and early results. Remote Sensing of Environment 83, 287-302. Goudet J (1995) FSTAT (Version 1.2): A computer program to calculate F-statistics. Journal of Heredity 86, 485-486. Goudet J (2001) FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3). Homer C, Dewitz J, Fry J, et al. (2007) Completion of the 2001 National Land Cover Database for the Counterminous United States. Photogrammetric Engineering and Remote Sensing 73, 337. Kool J, Moilanen A, Treml E (2012) Population connectivity: recent advances and new perspectives. Landscape Ecology, 1-21. Lee DJ, Vaughan MR (2003) Dispersal Movements by Subadult American Black Bears in Virginia. Ursus 14, 162-170. Legendre P, Fortin M-J (2010) Comparison of the Mantel test and alternative approaches for detecting complex multivariate relationships in the spatial analysis of genetic data. Molecular Ecology Resources 10, 831-844. 67 Liu W, Gopal S, Woodcock CE (2004) Uncertainty and confidence in land cover classification using a hybrid classifier approach. Photogrammetric Engineering and Remote Sensing 70, 963-972. Manel S, Joost S, Epperson BK, et al. (2010) Perspectives on the use of landscape genetics to detect genetic adaptive variation in the field. Molecular Ecology 19, 3760-3772. Manel S, Schwartz MK, Luikart G, Taberlet P (2003) Landscape genetics: combining landscape ecology and population genetics. Trends in Ecology & Evolution 18, 189-197. Mantel N (1967) Detection of disease clustering and a generalized regression approach. Cancer Research 27, 209-&. McRae BH (2006) Isolation by resistance. Evolution 60, 1551-1561. McRae BH, Beier P (2007) Circuit theory predicts gene flow in plant and animal populations. Proceedings of the National Academy of Sciences of the United States of America 104, 19885-19890. MDNR (2004) Review of remote sensing technology used in the IFMAP project, pp. 1-49. Michigan Department of Natural Resources. Moore JA, Draheim HM, Etter D, Winterstein S, Scribner KT (2014) Application of large-scale parentage analysis for investigating natal dispersal in highly vagile vertebrates: a case study of american black bears (Ursus americanus). Plos One 9, e91168. NOAA (2014) Western Great Lakes 2010 Coastal Change Analysis Program accuracy assessment, pp. 1-11, Charleston, South Carolina. Noyce K, Garshelis D (2011) Seasonal migrations of black bears (Ursus americanus): causes and consequences. Behavioral Ecology and Sociobiology 65, 823-835. Opdam P, Wascher D (2004) Climate change meets habitat fragmentation: linking landscape and biogeographical scale levels in research and conservation. Biological Conservation 117, 285-297. Paetkau D, Calvert W, Stirling I, Strobeck C (1995) Microsatellite analysis of populationstructure in Canadian polar bears. Molecular Ecology 4, 347-354. Paetkau D, Waits LP, Clarkson PL, Craighead L, Strobeck C (1997) An empirical evaluation of genetic distance statistics using microsatellite data from bear (Ursidae) populations. Genetics 147, 1943-1957. Peakall R, Smouse PE (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes 6, 288-295. Pelton MR (1992) Black Bear (Ursus americanus) In: Wild Mammals of North America: Biology, Management, and Economics (eds. Chapman JA, Feldhamer GA), pp. 504-514. The John Hopkins University Press,USA. 68 Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155, 945-959. Proctor MF, McLellan BN, Strobeck C, Barclay RMR (2005) Genetic analysis reveals demographic fragmentation of grizzly bears yielding vulnerably small populations. Proceedings of the Royal Society B-Biological Sciences 272, 2409-2416. Rayfield B, Fortin M-J, Fall A (2010) The sensitivity of least-cost habitat graphs to relative cost surface values. Landscape Ecology 25, 519-532. Raymond M, Rousset F (1995) GENEPOP (VERSION-1.2) - Population-genetics software for exact tests and ecumenicism. Journal of Heredity 86, 248-249. Reding D, Cushman S, Gosselink T, Clark W (2013) Linking movement behavior and fine-scale genetic structure to model landscape connectivity for bobcats (Lynx rufus). Landscape Ecology 28, 471-486. Rice WR (1989) Analyzing tables of statistical tests. Evolution 43, 223-225. Rogers LL (1977) Social relationships, movements, and population dynamics of black bears in northeastern Minnesota, University of Minnesota. Rogers LL (1987) Effects of food-supply and kinship on social-behavior, movements, and population-growth of black bears in northeastern Minnesota. Wildlife Monographs, 1-72. Rueness EK, Stenseth NC, O'Donoghue M, et al. (2003) Ecological and genetic spatial structuring in the Canadian lynx. Nature 425, 69-72. Schwartz CC, Franzmann AW (1992) Dispersal and survival of subadult black bears from the Kenai Peninsula, Alaska. Journal of Wildlife Management 56, 426-431. Segelbacher G, Cushman SA, Epperson BK, et al. (2010) Applications of landscape genetics in conservation biology: concepts and challenges. Conservation Genetics 11, 375-385. Shih CC, Huang CC, Li SH, Hwang MH, Lee LL (2009) Ten novel tetranucleotide microsatellite DNA markers from Asiatic black bear, Ursus thibetanus. Conservation Genetics 10, 1845-1847. Short Bull R, Cushman S, Mace R, et al. (2011) Why replication is important in landscape genetics: American black bear in the Rocky Mountains. Molecular Ecology 20, 10921107. Smouse PE, Long JC, Sokal RR (1986) Multiple regression and correlation extensions of the Mantel test of matrix correspondence. Systematic Zoology, 627-632. Spear SF, Balkenhol N, Fortin MJ, McRae BH, Scribner K (2010) Use of resistance surfaces for landscape genetic studies: considerations for parameterization and analysis. Molecular Ecology 19, 3576-3591. Storfer A, Murphy MA, Evans JS, et al. (2007) Putting the 'landscape' in landscape genetics. Heredity 98, 128-142. 69 Taberlet P, Camarra JJ, Griffin S, et al. (1997) Noninvasive genetic tracking of the endangered Pyrenean brown bear population. Molecular Ecology 6, 869-876. Thomlinson JR, Bolstad PV, Cohen WB (1999) Coordinating methodologies for scaling landcover classifications from site-specific to global: Steps toward validating global map products. Remote Sensing of Environment 70, 16-28. Thompson LM, van Manen FT, King TL (2005) Geostatistical analysis of allele presence patterns among American black bears in eastern North Carolina. Ursus 16, 59-69. Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes 4, 535-538. Wright S (1943) Isolation by distance. Genetics 28, 114-138. Wu H, Zhang SN, Wei FW (2010) Twelve novel polymorphic microsatellite loci developed from the Asiatic black bear (Ursus thibetanus). Conservation Genetics 11, 1215-1217. 70 CHAPTER 4 DETECTING BLACK BEAR SOURCE-SINK DYNAMICS USING INDIVIDUAL-BASED GENETIC GRAPHS Draheim, H. M., J. A. Moore, S. R. Winterstein, D. Etter, K. T. Scribner ABSTRACT Source-sink dynamics affect population connectivity, spatial genetic structure, and population viability for many species. We introduce a novel approach that uses genetic graphs to identify source-sink areas within a continuously-distributed population of black bears (Ursus americanus) in Michigan, USA. Black bear harvest samples (N=569, from 2002, 2006, and 2010) were genotyped at 12 microsatellite loci and locations were compared across years to identify areas of consistent occupancy over time. We compared graph metrics estimated for 10 ecological models, with comparable metrics from a genetic model, to identify ecological factors that are associated with sources and sinks. Genetic graphs showed asymmetric connectivity, with areas of high and low net flux (source strength), revealing black bears in the NLP exhibited source-sink dynamics. Source strength was associated with black bear local harvest density (a proxy for bear density). Simulated removal of nodes revealed that nodes removed in source areas substantially reduce connectivity. Findings demonstrate that through the use of individual-based genetic graphs it is possible to characterize source-sink dynamics in a continuously distributed species in the absence of discrete habitat patches. In addition, because black bears are a harvested game species, findings stress the importance of considering source-sink dynamics in harvest management. Keywords: source-sink dynamics, graph theory, genetic relatedness, black bear, connectivity 71 INTRODUCTION When habitat quality varies across heterogeneous landscapes, occupancy and relative abundance can concomitantly vary (Knowlton & Graham 2010). As population densities increase in higher quality habitats, intraspecific competition for resources can force individuals to disperse to lower density and lower quality habitats (Clobert et al. 2009). This phenomenon is commonly referred to as source-sink dynamics (Pulliam 1988). Source areas are characterized by higher quality habitat and a surplus of individuals that emigrate to other areas. In contrast, sink areas are characterized by lower quality habitat, have negative population growth rates and require an influx of dispersing individuals from source areas to maintain population persistence (Boughton 1999; Minor & Urban 2008). Source-sink dynamics fundamentally affect population connectivity and spatial genetic structure of species, and have evolutionary and ecological implications that could influence population viability, species persistence and evolutionary potential (Dias 1996; Gaggiotti 1996; Pulliam 1988). Furthermore, source-sink dynamics may influence the evolutionary trajectory of species by disrupting local adaptation to low quality (sink) habitats as directional selection is mitigated by an influx of individuals from high quality (source) habitats (Anderson & Geber 2010; Kawecki 1995; Kawecki 2008). Genetic variation may be lost as alleles present in the source populations, which have higher mean fitness across all populations, tend to be fixed (Lenormand 2002). Thus, detecting sources and sinks is critical for understanding microevolutionary processes that affect populations and for managing wildlife and conserving at-risk species. Traditionally, source-sink dynamics have been detected using direct methods including field observations or capture-mark-recapture studies (Peery et al. 2008). An alternative method is to apply graph theory to estimate asymmetrical connectivity (Dale & Fortin 2010). Graphs are 72 simple visual representations of real systems (Peterman et al. 2013) that are composed of nodes (representing biological units like individuals, populations, or habitat patches) that are characterized by attributes like size or quality. Nodes are connected by edges that represent biological phenomena like dispersal. Edges can be symmetric, with uniform weighting between nodes. Alternatively, asymmetric edges are directionally weighted where, for instance, dispersal from patch i to patch j is larger than from patch j to patch i, a pattern commonly found in sourcesink dynamics (Rayfield et al. 2011). Combining population genetics and graph theory offers a compelling framework (Gaggiotti 1996; Murphy et al. 2010; Peery et al. 2008; Pierson et al. 2013) for understanding source sink dynamics and population connectivity. However, for species that are continuously rather than patchily distributed, using genetic data to estimate gene flow has been hampered by an inability to delineate populations. Use of individual-based measures of pairwise genetic relatedness allows for investigations of fine-scale genetic structure and dispersal patterns of continuously distributed species without a priori delineation of populations (Cushman et al. 2006; Rousset 2000). Pairwise relatedness represents the degree of shared ancestry, or the probability that two alleles at a locus, one taken at random from each of two individuals, are identical-by-descent (Queller & Goodnight 1989). When applied in a graph context, nodes can represent individuals or populations connected by edges representing shared ancestry (Dale & Fortin 2010; Garroway et al. 2011; Murphy et al. 2010). Continuously distributed species commonly exhibit a pattern of isolationby-distance (IBD), where individuals closer in geographic proximity are more genetically similar (Wright 1943). Alternatively, if habitats are highly heterogeneous, habitat quality can more strongly impact genetic structure than distance alone (Manel et al. 2003). Differential occupancy 73 in heterogeneous landscapes will ultimately affect patterns of relatedness. For example, within source areas with high recruitment and little immigration, average pairwise relatedness will be higher than in sinks due to accumulation of multi-generational family groups. In contrast, average pairwise relatedness will be lower in sink areas that are characterized by high rates of immigration of individuals with mixed ancestry (Peery et al. 2008). Graph theory also provides metrics that quantify the contributions of individual nodes and edges to overall graph connectivity. When applied to management questions, such metrics can provide a novel and powerful means for managers and conservationists to predict the impact of future policy decisions or land use changes. For example, wildlife managers can predict how increasing regional harvest quotas will impact global population connectivity. A need for landscape genetic methods that move beyond assessing the influences of landscape features on gene flow to predicting the impacts of management decisions and landscape change has recently been recognized (Keller et al. 2014; Manel & Holderegger 2013; van Strien et al. 2014). We use an individual-based graph approach to investigate source-sink dynamics in a population of American black bears (Ursus americanus) in the Northern Lower Peninsula (NLP) of Michigan, USA. Black bears have home ranges that can overlap depending on resource availability (Carter et al. 2010; Costello 2010; Horner & Powell 1990; Mitchell & Powell 2007) and habitat quality (Figure C1). Bears in the NLP are geographically isolated as the population is bounded to the south by a matrix of unsuitable urban and agricultural habitat and on all other sides by the Great Lakes, thus the population experiences no emigration and immigration that could confound results. We applied graph theory to a dataset spanning a nine year period (20022010) to identify specific regions in the NLP that function as source and sink areas. Our objectives were to examine the importance of landscape and ecological features on directional 74 dispersal, to assess the relative influence of specific nodes on overall connectivity, and to test the hypothesis that source areas are associated with high quality habitat and high local population density. METHODS Study Area Michigan’s Northern Lower Peninsula (NLP) is a 47,739 km2 (Figure 4.1) fragmented mosaic of variable land cover and land use types including development, agriculture, upland non-forested openings, northern hardwood and mixed hardwood, oak, aspen, pine, forested wetland, and nonforested wetland (Carter et al. 2010). During the fall harvest (September-October) of 2002 (N=263), 2006 (N=385), and 2010 (N = 336), teeth were collected from harvested black bears that were registered at hunter check stations. Locations of bear harvest samples were recorded as township, range, and sections, which were then georeferenced to the section centroid and converted into UTM coordinates. Laboratory Analysis We extracted DNA from bear teeth using Qiagen DNEasy Tissue Kits following manufacturer protocols (Qiagen Inc., Valencia, CA). DNA was quantified using a Nanodrop spectrophotometer (Thermo Scientific, Waltham, MA) and diluted to a 20 ng/μl working concentration. We amplified twelve microsatellite loci using polymerase chain reaction (PCR), including; G10X, G10L, G10D, G10M, G10B (Paetkau et al. 1995), UarMU50, and UarMU59 (Taberlet et al. 1997), UT29, UT35, UT38 (Shih et al. 2009), ABB1, and ABB4 (Wu et al. 2010) using methods described in Moore et al. (Moore et al. 2014). 10% of samples were randomly 75 Figure 4.1. Study area in the Northern Lower Peninsula of Michigan (NLP) showing node locations (N = 141) for black bear harvest samples collected during 2002, 2006, and 2010 (N = 569). 76 selected and genotyped twice to provide a genotyping error rate of <2%. Using MICROCHECKER (Van Oosterhout et al. 2004) and GenAlEx ver. 6.0, (Peakall & Smouse 2006) , no loci were found to deviate significantly from Hardy-Weinberg and linkage equilibrium, so all twelve were retained for further analyses (Table C1). Defining the Nodes Graphs are composed of nodes connected by edges, and both of these can be assigned weights. When applied in a landscape context, nodes are traditionally defined as the spatial centroid of a patch or population, and are generally weighted by attributes representing measures of quality (e.g. patch size or habitat suitability). Since it is not feasible to assign node weights representing local landscape or ecological characteristics to a location representing a single individual, we defined nodes as areas that are consistently occupied in both space and time. Our sample includes individuals collected from 2002, 2006, and 2010 representing three generations. We used the harvest locations from each year to generate polygon layers using Voronoi tessellations. Briefly, Voronoi tessellations are polygons created around each point such that the borders of each polygon are equidistant from their two nearest points. Annual polygon layers were then overlaid and areas of consistent occupancy were identified where polygons from 2002, 2006, and 2010 overlapped and contained at least one bear from each of those years. We used STAMP (spatial–temporal analysis of moving polygons (Robertson et al. 2007) in ArcGIS 10.0 to identify overlap areas. Nodes were thus defined as areas of consistent polygon overlap (i.e. occupancy) across the three time points. Node locations were defined as the centroids of overlapping polygons that contained at least one bear from 2002, 2006, and 2010. 77 Graph construction We generated a genetic graph and a set of ecological graphs to explore genetic, habitat, and demographic connectivity of the NLP black bear population: 1. Genetic Modeling-- We calculated pairwise relatedness among individuals using a maximum likelihood estimator of relatedness (K; 38), in the program ML-Relate (Kalinowski et al. 2006). We constructed a genetic ‘saturated graph’ in which each node was connected to every other node. We assigned node weights as the mean relatedness among individuals within each node, and edge weights as the mean pairwise relatedness between individuals at node i and node j. We tested for effects of sample size on mean relatedness at the nodes (node weight) and found no significant association (simple liner regression; R2 = 0.004, P = 0.50). 2. Ecological Modeling--To test associations between factors that affect node importance (strength) and overall graph connectivity we created ecological graphs. Edges were calculated using Euclidean or least-cost distance (based on habitat permeability) from node i to node j. We assigned node weights based on habitat suitability or local harvest density. We performed Spearman’s rank correlation analyses between the genetic model and our set of ecological candidate models (Table 4.1) in R (R Development Core Team 2012). a. Weighting nodes--To obtain habitat suitability weights for each node, we reclassified the 2006 NOAA Coastal Change Analysis Program (C-CAP) Land Cover dataset (resolution = 30m) into eight land cover classifications and ranked them based on bear habitat suitability (1-100, least to most suitable; Table 4.2) including: mixed deciduous forest (MF, 100), forested wetland (FW, 100), evergreen forest (EF, 90), non-forested upland (NFU, 80), agriculture (AG, 50), non-forested wetland (NFW, 1), developed (DEV, 1), and open water (NA) using an ecological model developed independently for Michigan’s NLP black bears 78 Table 4.1. Spearman's rank correlations between genetic and ecological measures of node importance (source strength). The three graph metrics compared are net flux (Fout-Fin), area-weighted flux (AWF), and probability of connectivity (PC). ‘High’ indicates maximum magnitude weights (e.g.,1-100) and ’low’ indicates minimum magnitude weights (e.g.,1-5) used to test the sensitivity of analyses to the magnitude of cost values. All correlations were significant (P < 0.05). Ecological Model Node weight Edge Weight Habitat Suitability (high) Habitat Suitability (low) Habitat Suitability (high) Habitat Suitability (low) Habitat Suitability (high) Graph Metrics Fout-Fin AWF PC Habitat Suitability (high) Habitat Suitability (low) road (high) road (low) distance only 0.1865 0.1947 0.1836 0.1809 0.1799 0.2040 0.2160 0.2270 0.2150 0.2270 0.2094 0.2050 0.2160 0.2040 0.1980 Habitat/Non-habitat (high) Habitat/Non-habitat (low) Habitat/Non-habitat (high) Habitat/Non-habitat (low) Habitat/Non-habitat (high) Habitat/Non-habitat (high) Habitat/Non-habitat (low) road (high) road (low) distance only 0.1712 0.1690 0.1681 0.1581 0.1490 0.2030 0.2100 0.2120 0.2040 0.2100 0.1985 0.2081 0.1970 0.2030 0.1970 Density Density Density Density Density Density Density Habitat Suitability (high) Habitat Suitability (low) Habitat/Non-habitat (high) Habitat/Non-habitat (low) road (high) road (low) distance only 0.2021 0.2032 0.2020 0.2020 0.2034 0.2020 0.2021 0.2630 0.2780 0.2740 0.2820 0.2800 0.2830 0.2820 0.2590 0.2750 0.2700 0.2770 0.2760 0.2790 0.2780 79 Table 4.2. Node weights for graph analysis and cost values for resistance surfaces used to calculate pairwise least cost distance for edge weights. Node weights and cost values are inferred from Carter et al. (2010). Habitat/Non-Habitat values designate land cover positively or negatively correlated with bear presence inferred from Carter et al. (2010). Graph Attribute Node Weight Max Min Cost Values Max Min Habitat Suitability Development Agriculture Non-Forest Upland Mixed Forest Conifer Forest Forested Wetland Non-Forested Wetland Water Open water (lakes) Non-state roads State roads 1 50 80 100 90 100 1 1 N/A 90 1 1 2 3 5 4 5 1 1 N/A 4 1 100 50 20 1 10 1 100 100 N/A 10 100 5 4 3 1 2 1 5 5 N/A 2 5 Habitat/Non-Habitat Habitat Non-Habitat 100 1 2 1 1 100 1 3 N/A N/A N/A N/A N/A N/A 1 10 100 1 2 3 Roads No road Non-state roads State roads 80 (Carter et al. 2010). In addition, we reclassified land cover data into a coarser habitat suitability classification as either bear habitat (MF, FW; 100) or non-habitat (EF, NFU, AG, NFW, DEV; 1). We used density of harvest locations as a surrogate for population density. To estimate local harvest density we used annual harvest locations from 2002-2010 and generated kernel density grids in ArcGIS 10.0 and reclassified them into a range of 1-10 (low to high as previously defined (Moore et al. 2014)). We then created a harvest grid by calculating the median values over the 9 annual harvest grids (Figure C2). We extracted a local harvest density and land cover value based on majority of grid cells within the nodes. b. Weighting edges-Edges for ecological graphs were weighted by either Euclidean distance or least cost distance (one path of least resistance through a resistance surface, reflecting permeability through habitat types or boundaries) in ARCGIS 10.0. Least cost distance was calculated using three resistance surfaces to model habitat permeability (Table 4.1). The 2006 CCAP data set was reclassified into eight cover types reflecting the ecological model (Carter et al. 2010) and two cover types for bear habitat vs. non-habitat (Table 4.1). In addition, we obtained and reclassified a road only resistance surface from the MI Geographic Data Library (Center for Geographic Information, Table 4.1). Pairwise least cost distances were calculated using the landscape genetic toolbox in ArcGIS (Etherington 2011). Least cost analysis and node weights may be sensitive to relative cost values assigned to land cover types (Rayfield et al. 2010) and could potentially influence our results. Thus, we performed a sensitivity analysis on the range of relative cost values between resistance surface land cover types by comparing least cost distances between rasters weighted using large differences (i.e. maximum weights = 1, 10, 20, 50, 100) and 81 small differences (minimum weights = 1, 2, 3, 4, 5) in the relative cost values for land cover types (Table 4.2). Node importance and source-sink analysis Graph metrics were calculated to define node importance to identify source areas. We calcuated dispersal probability (pij) for each pair of nodes: Eq 1 pij = exp(k x dij) where k is a decay coefficient calculated from a tail dispersal distance (an arbitrarily selected point on the flat tail of a negative-exponential dispersal-distance function) and dij is the distance from node i to j (Urban & Keitt 2001). Area-weighted flux (AWF) was estimated using the previous metric (dispersal probability) multiplied by the quality area (i.e. genetic graph, qi = mean relatedness; ecological graph, qi = density or habitat quality) of the nodes; Eq 2 AWF = qi x pij (Minor & Urban 2007) Probability of connectivity (PC) is a graph metric that combines the attributes of the nodes with the maximum product probability of all the possible paths between every pair of nodes; Eq 3 where AL is the total quality measure of the entire graph and ai and aj are the qualities of node i and j, respectively (Saura & Pascual-Hortal 2007). AWF and PC were calculated using program CONEFOR 2.6 (Saura & Torne 2009). 82 Influx (flux entering a node) and outflux (flux exiting a node) were calculated to account for asymmetric dispersal as there are two distinct edges for each pair of nodes. Edges were based on distances so that pij = pji. However, influx and outflux at a particular node may differ, as flux also incorporates quality area qi of the donor node i (i.e., fij = qi x pij vs. fji = qj x pij ). Total flux out of node i and total flux into node i was calculated as: Eq 4 Fout i = Σ fji Eq 5 Fin i = Σ fji Source nodes were characterized by high total outflux and sink nodes were characterized by high influx and low outflux. To determine the magnitude and direction of the net flux associated with each node, we calculated the flux difference (Fout i - Fin i) (Minor & Urban 2007). To visualize the spatial configuration of source and sink nodes we generated interpolated surfaces of node importance (source strength) using Inverse Distance Weighting (IWD). Node removal--Node removal is a means to examine the relative importance of a node for maintaining connectivity in the landscape. For our study, node removal is a way to predict the effects of different harvest or landscape change scenarios. We used node removal to quantify the effects of different removal scenarios on changes in connectivity over the entire landscape (NLP). Ten nodes were removed randomly, in a stepwise manner (one node per step). After each step of the removal process the graph was reanalyzed to determine the importance of the node to overall graph connectivity by comparing graph global AWF and PC values. Next, we sorted the nodes according to net flux and randomly removed 10 nodes from within the highest and lowest 10% net flux values to evaluate the impact of removing important sources or extreme sink areas, respectively. In addition, we removed 10 nodes with intermediate net flux values (10% of net 83 flux near zero). We performed three replicates of the four node removal scenarios (source nodes, sink nodes, intermediate nodes (~ 0 net flux) and random nodes). RESULTS Based on our saturated genetic graph using 141 nodes (areas of consistent occupancy) Figure 4.2 presents the juxtaposition of inferred source and sink nodes over the study area. Of the 141 nodes in the NLP, 62 are considered source nodes (positive net flux) and 16 nodes had net flux > 0.7. In contrast, 79 nodes were classified as sinks (negative net flux). We estimated two additional graph metrics that measure node importance (i.e., source strength); 1) area-weighted flux (AWF) and 2) probability of connectivity (PC). AWF and PC were highly correlated with net flux (Spearman’s rank correlations ρ > 0.90). We constructed additional saturated graphs based on the same 141 nodes using purely ecological data (only habitat or local harvest density). Correlation tests of node metrics derived from genetic and ecological graphs revealed four major findings. First, correlations between graph metrics (net flux, AWF, and PC) and habitat/local harvest abundance graph metrics were all significant (ρ = 0.149 – 0.279, P < 0.05; Table 4.1). Second, local bear harvest density was most highly correlated with genetic graph metrics (net flux, mean ρ = 0.2024; AWF, mean ρ = 0.2812; PC, mean ρ = 0.2730, P < 0.05). Third, ecological models where nodes were weighted by habitat suitability revealed similar correlation values irrespective of how finely (habitat suitability model; net flux, mean ρ = 0.1842; AWF, mean ρ = 0.2182; PC, mean ρ = 0.2045) or coarsely (habitat vs. non-habitat; net flux, mean ρ = 0.1620; AWF, mean ρ = 0.2100; PC, mean ρ = 0.1993) the land cover was defined (Table 4.1). Lastly, because least cost distances can be sensitive to the magnitude of cost values, we performed a sensitivity analysis of cost values 84 Figure 4.2. Interpolated surfaces of black bear node importance (source strength) based on net flux at the nodes across the Northern Lower Peninsula study area using Inverse Distance Weighting (IDW). Interpolated surfaces of (a) the genetic, (b) Carter et al. 2010 habitat suitability, (c) a habitat/non-habitat, and (d) a bear density model are shown for comparison. Warm colors represent sources and cool colors represent sinks. a b c d 85 assigned to land cover and found graph metrics were not sensitive to the magnitude of land cover cost values (Table 4.1). Node removal To investigate the effects of node removal that reflects potential harvest scenarios, we removed nodes and examined their impact on overall population connectivity. Figure 4.3 illustrates four scenarios of node removal including removal of random nodes, source nodes, sink nodes, and intermediate nodes (net flux ≈ 0). Overall, global graph connectivity measures AWF and PC declined as successive nodes were systematically removed (Figure 4.3). When nodes were removed at random with respect to net flux, overall connectivity of the graph declined. When intermediate nodes were removed (i.e., those that exhibited no net flux difference) the decline in connectivity did not differ from random. However, when high net flux nodes, classified as important sources, were removed, connectivity declined rapidly. In contrast, when low net flux nodes, classified as sinks, were removed, connectivity was virtually unchanged (Figure 4.3). DISCUSSION By combining genetic and ecological graphs we have identified source and sink areas using individual-based genetic data in a landscape genetic framework. Our methodological challenge was to find a way to apply genetic graphs to a continuously-distributed population exhibiting clinal isolation-by-distance. If populations are patchily distributed and/or strongly genetically structured, allowing for a priori delineation of populations or patches, genetic weights can easily be assigned by calculating the mean genetic similarity (e.g. relatedness) over 86 Figure 4.3. Area-weighted dispersal flux (AWF, a) and probability of connectivity (b) as a function of the number of nodes removed for each four node removal scenarios (random, sources only, sinks only and intermediate only). Error bars represent the standard error (SE) from three independent runs. a b 87 the entire population or patch (Clobert et al. 2009; Murphy et al. 2010). However, for continuously distributed species that cannot be unambiguously placed into defined populations or patches, such as black bears, an alternative, individual-based method of weighting nodes is necessary. We overcame this methodological obstacle by using point pattern and spatial analyses to identify areas of consistent spatial occupancy over time and assigned node weights based on average pairwise relatedness of individuals within the nodes. Genetic graphs can be applied to detect source-sink dynamics in other systems where populations are not patchily distributed or strongly genetically structured into discrete populations. For example, our method may find utility in wide ranging and elusive species where source-sink dynamics have gone undetected (e.g. fur bearing species, carnivores, marine mammals; 47-49). Undetected source-sink dynamics have potentially negative consequences, particularly for harvested species. For example, if the management goal is to reduce abundance in a specific sink area, through harvest, dispersal from surrounding source areas will confound efforts. Alternatively a source area with abundant resources and subsequent high reproduction may become an “attractive sink” if harvest pressure is high (Robinson et al. 2008). Our analyses revealed patterns of asymmetric connectivity in NLP black bears based on areas of high and low net flux. Findings are consistent with previous local spatial autocorrelation analyses. Chapter 3 shows evidence of fine-scale spatial genetic structure in NLP black bears but no strong population structuring across the region. We identified 62 source areas including 16 nodes we consider to be important source areas (net flux > 0.7) in the NLP. Interpolated surfaces of net flux values (Figure 4.2) revealed important source areas were patchily distributed with clusters of adjacent high net flux nodes occurring on private and forested public lands in the northwestern, northern and central portions of the study area. Net flux estimates along the study 88 area boundary should be interpreted with caution as interpolations may be prone to edge effects (Okabe et al. 2000). In addition, this map is interpolated from point data (N = 141 nodes), and observed patterns are dependent on the spatial dispersion of nodes. Therefore, interpolated maps should be interpreted as displaying regional patterns and not the exact locations of source and sink nodes. We found significant correlations between graph metrics derived from genetic data and ecological models based on habitat quality and local harvest density. Graphs with nodes connected by edges that incorporated measures of habitat permeability (least-cost distance) had slightly higher correlations with genetic graph metrics characterizing source strength than Euclidean distance, suggesting black bear connectivity is influenced by landscape features. Local harvest density at the nodes better explained source strength than did habitat quality. Source-sink dynamics, as defined by Pullman (1988), result from density-dependent dispersal (Amarasekare 2004; Pulliam 1988). High densities often lead to density-dependent dispersal (Amarasekare 2004) and are probably indicative of productive areas with high resource availability. Thus, we predicted source areas identified using genetic data would be characterized by areas of high local harvest abundance which is a surrogate for population density (Cattadori et al. 2003; Kitson 2004; Royama 1992). Results are consistent with previously observed dispersal patterns for black bears (Moore et al. 2014; Roy et al. 2012). Node-removal analysis allowed assessment of relative contributions of individual source and sink nodes to connectivity across the study site (Minor & Urban 2007; Urban et al. 2009). Overall, genetic connectivity decreased as nodes were removed; however, the impact varied dramatically based on net flux. Graph connectivity is sensitive to loss of nodes characterized by high net flux (sources), and connectivity declined rapidly when source nodes (net flux > 0.7) 89 were removed compared to removal of sink nodes or nodes with intermediate net flux. The importance of maintaining genetic connectivity cannot be overstated. Decreased genetic connectivity can result in reduced genetic diversity via genetic drift and inbreeding, and consequently can lead to higher extinction rates and loss of evolutionary potential (Allendorf et al. 2012; Frankham et al. 2002). Node removal results have management implications as data indicate loss of source areas will impact genetic connectivity across the entire region. Node removal provides a practical methodology for assessing the potential impacts of harvest scenarios or landscape change on population connectivity, and therefore has strong potential for management and conservation applications. CONCLUSIONS Graphs are simplified portrayals of natural systems that we have shown can be powerful tools for assessing and detecting source-sink dynamics (Bunn et al. 2000; Minor & Urban 2007; Urban & Keitt 2001; Urban et al. 2009). Our method does not estimate migration rates. Rather, we derived graph metrics based on pairwise relatedness within and among nodes to evaluate plausible spatially-explicit hypotheses about degree and directionally of connectivity among nodes (source-sink dynamics) and the ecological and environmental factors that affect connectivity. Rigorous identification of source and sink areas traditionally required a great deal of detailed information on population survival rates, fecundity, immigration and emigration (Dias 1996; Lowe & Allendorf 2010; Pulliam 1988). However, such detailed information is difficult to acquire for an elusive, highly mobile animal such as black bears. Dispersal studies that use traditional methods (i.e., capture-mark-recapture or radio-telemetry), are expensive and labor intensive (Solberg et al. 2006) and may lead to underestimates due to small sample size and 90 limited geographic scope. Genetic graphs provide a viable alternative method for detecting source-sink dynamics in continuously distributed species. Because genetic graphs provide a flexible framework for understanding connectivity, they could be widely integrated into landscape genetics research and conservation planning at multiple spatial scales. While genetic graphs may be broadly applicable across a range of taxa, the graph models need to be parameterized based on the specific life history characteristics of the species or population of interest. ACKNOWLEDGMENTS We are grateful to J. Fierke, K. Filcek, B. Dreher, and S. Libants for assisting in data generation. We also thank M. Fortin, J. Messina, B. Epperson, K. Holekamp, J. Kanefsky, J. McGuire, for insightful discussion and helpful comments on the analysis and manuscript. All samples used in analyses were provided by Michigan bear hunters and collected by the many MDNR staff at registration stations. Support for this project was provided by the Michigan Department of Natural Resources through the Wildlife and Sportfish Restoration Program F11AF00640 and by the Department of Fisheries and Wildlife at Michigan State University. 91 APPENDIX 92 APPENDIX Table C1. Black bear genetic diversity at 12 microsatellite loci. No loci showed significant deviations from Hardy-Weinberg equilibrium (Raymond & Rousset 1995). Na, number of alleles per locus; HO, observed heterozygosity; HE, expected heterozygosity; FIS, inbreeding coefficient). Locus X M D L B Uar50 Uar59 ABB1 ABB4 UT35 UT38 UT29 Na 11 10 9 13 7 9 7 9 8 6 26 12 HO HE 0.697 0.756 0.766 0.862 0.696 0.712 0.695 0.712 0.815 0.764 0.921 0.877 0.684 0.814 0.759 0.861 0.766 0.796 0.702 0.732 0.799 0.781 0.907 0.896 FIS -0.018 0.071 -0.010 -0.001 0.091 0.105 0.010 0.028 -0.019 0.021 -0.016 0.020 93 Figure C1. Map of the NLP land cover. Land cover data were reclassified into eight land cover classifications from 2006 NOAA Coastal Change Analysis Program Land Cover (C-CAP) dataset (resolution = 30m). 94 Figure C2. Map of NLP local black bear harvest density in NLP. Local density was classified into a range of 1-10 (low to high). 95 Figure C3. Map of overlapping Voronoi tessellation point pattern polygons in NLP. Polygons contained at least one sample from 2002, 2006, and 2010. 96 Figure C4. Distribution of the node importance values characterized by inter-individual genetic relatedness, shown are the a) net flux, b) area-weighted flux (AWF) and c) probability of connectivity (PC) metrics. For net flux positive values represent source nodes and negative values represent sink values. b a c 97 Figure C5. (a) Distribution of mean relatedness (K) among individuals within the nodes. (b) Mean relatedness at the node according to number of individuals at the node. a b 98 REFERENCES 99 REFERENCES Allendorf FW, Luikart GH, Aitken SN (2012) Conservation and the genetics of populations Wiley. com. Amarasekare P (2004) The role of density-dependent dispersal in source–sink dynamics. Journal of Theoretical Biology 226, 159-168. Anderson JT, Geber MA (2010) Demographic source-sink dynamics restrict local adaptation in Elliott’s blueberry (Vaccinium elliottii). Evolution 64, 370-384. Boughton DA (1999) Empirical evidence for complex source-sink dynamics with alternative states in a butterfly metapopulation. Ecology 80, 2727-2739. Bunn AG, Urban DL, Keitt TH (2000) Landscape connectivity: A conservation application of graph theory. Journal of Environmental Management 59, 265-278. Carter NH, Brown DG, Etter DR, Visser LG (2010) American black bear habitat selection in northern Lower Peninsula, Michigan, USA, using discrete-choice modeling. Ursus 21, 57-71. Cattadori IM, Haydon DT, Thirgood SJ, Hudson PJ (2003) Are indirect measures of abundance a useful index of population density? The case of red grouse harvesting. Oikos 100, 439446. Clobert J, Le Galliard JF, Cote J, Meylan S, Massot M (2009) Informed dispersal, heterogeneity in animal dispersal syndromes and the dynamics of spatially structured populations. Ecology Letters 12, 197-209. Costello CM (2010) Estimates of dispersal and home-range fidelity in American black bears. Journal of Mammalogy 91, 116-121. Cushman SA, McKelvey KS, Hayden J, Schwartz MK (2006) Gene flow in complex landscapes: Testing multiple hypotheses with causal modeling. American Naturalist 168, 486-499. Dale MRT, Fortin M-J (2010) From graphs to spatial graphs. Annual Review of Ecology, Evolution, and Systematics 41, 21-38. Dias PC (1996) Sources and sinks in population biology. Trends in Ecology & Evolution 11, 326-330. Draheim HM (2015) Landscape genetics of black bears (Ursus Americanus) in the Northern Lower Peninsula (NLP) of Michigan, USA Dissertation, Michigan State University. 100 Etherington TR (2011) Python based GIS tools for landscape genetics: visualising genetic relatedness and measuring landscape connectivity. Methods in Ecology and Evolution 2, 52-55. Frankham R, Ballou JD, Briscoe DA (2002) Introduction to conservation genetics Cambridge University Press. Gaggiotti OE (1996) Population genetic models of source-sink metapopulations. Theoretical Population Biology 50, 178-208. Garroway CJ, Bowman J, Wilson PJ (2011) Using a genetic network to parameterize a landscape resistance surface for fishers, Martes pennanti. Molecular Ecology 20, 3978-3988. Hoffman JI, Matson CW, Amos W, Loughlin TR, Bickham JW (2006) Deep genetic subdivision within a continuously distributed and highly vagile marine mammal, the Steller's sea lion (Eumetopias jubatus). Molecular Ecology 15, 2821-2832. Horner MA, Powell RA (1990) Internal structure of home ranges of black bears and analyses of home-range overlap. Journal of Mammalogy 71, 402-410. Kalinowski ST, Wagner AP, Taper ML (2006) ml‐relate: a computer program for maximum likelihood estimation of relatedness and relationship. Molecular Ecology Notes 6, 576579. Kawecki T (1995) Demography of source—sink populations and the evolution of ecological niches. Evolutionary Ecology 9, 38-44. Kawecki TJ (2008) Adaptation to marginal habitats. Annual Review of Ecology, Evolution, and Systematics 39, 321-342. Keller D, Holderegger R, van Strien MJ, Bolliger J (2014) How to make landscape genetics beneficial for conservation management? Conservation Genetics, 1-10. Kitson JC (2004) Harvest rate of sooty shearwaters (Puffinus griseus) by Rakiura Māori: a potential tool to monitor population trends? Wildlife Research 31, 319-325. Knowlton JL, Graham CH (2010) Using behavioral landscape ecology to predict species' responses to land-use and climate change. Biological Conservation 143, 1342-1354. Lenormand T (2002) Gene flow and the limits to natural selection. Trends in Ecology & Evolution 17, 183-189. Lowe WH, Allendorf FW (2010) What can genetics tell us about population connectivity? (vol 19, pg 3038, 2010). Molecular Ecology 19, 5320-5320. Manel S, Holderegger R (2013) Ten years of landscape genetics. Trends in Ecology & Evolution 28, 614-621. 101 Manel S, Schwartz MK, Luikart G, Taberlet P (2003) Landscape genetics: combining landscape ecology and population genetics. Trends in Ecology & Evolution 18, 189-197. Minor ES, Urban DL (2007) Graph theory as a proxy for spatially explicit population models in conservation planning. Ecological Applications 17, 1771-1782. Minor ES, Urban DL (2008) A graph-theory frarmework for evaluating landscape connectivity and conservation planning. Conservation Biology 22, 297-307. Mitchell MS, Powell RA (2007) Optimal use of resources structures home ranges and spatial distribution of black bears. Animal Behaviour 74, 219-230. Moore JA, Draheim HM, Etter D, Winterstein S, Scribner KT (2014) Application of large-scale parentage analysis for investigating natal dispersal in highly vagile vertebrates: a case study of american black bears (Ursus americanus). Plos One 9, e91168. Murphy MA, Dezzani R, Pilliod DS, Storfer A (2010) Landscape genetics of high mountain frog metapopulations. Molecular Ecology 19, 3634-3649. Okabe A, Boots B, Sugihara K, SN C (2000) Spatial tessellations: concepts and applications of voronoi diagrams., 2nd ed. edn. Wiley, Chichester. Paetkau D, Calvert W, Stirling I, Strobeck C (1995) Microsatellite analysis of populationstructure in Canadian polar bears. Molecular Ecology 4, 347-354. Peakall R, Smouse PE (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes 6, 288-295. Peery MZ, Beissinger SR, House RF, et al. (2008) Characterizing source-sink dynamics with genetic parentage assignments. Ecology 89, 2746-2759. Peterman W, Rittenhouse TG, Earl J, Semlitsch R (2013) Demographic network and multiseason occupancy modeling of Rana sylvatica reveal spatial and temporal patterns of population connectivity and persistence. Landscape Ecology, 1-13. Pierson JC, Allendorf FW, Drapeau P, Schwartz MK (2013) Breed locally, disperse globally: fine-scale genetic structure despite landscape-scale panmixia in fire-specialist. Plos One 8, e67248. Pulliam HR (1988) Sources, sinks, and population regulation. American Naturalist 132, 652-661. Queller DC, Goodnight KF (1989) Estimating relatedness using genetic markers. Evolution 43, 258-275. Rayfield B, Fortin M-J, Fall A (2010) The sensitivity of least-cost habitat graphs to relative cost surface values. Landscape Ecology 25, 519-532. 102 Rayfield B, Fortin M-J, Fall A (2011) Connectivity for conservation: a framework to classify network measures. Ecology 92, 847-858. Raymond M, Rousset F (1995) GENEPOP (VERSION-1.2) - Population-genetics software for exact tests and ecumenicism. Journal of Heredity 86, 248-249. Robertson C, Nelson T, Boots B, Wulder M (2007) STAMP: spatial–temporal analysis of moving polygons. Journal of Geographical Systems 9, 207-227. Robinson HS, Wielgus RB, Cooley HS, Cooley SW (2008) Sink populations in carivore managment: cougar demography and immigration in a hunted population. Ecological Applications 18, 1028-1037. Rousset (2000) Genetic differentiation between individuals. Journal of Evolutionary Biology 13, 58-62. Roy J, Yannic G, Côté SD, Bernatchez L (2012) Negative density-dependent dispersal in the American black bear (Ursus americanus) revealed by noninvasive sampling and genotyping. Ecology and Evolution 2, 525-537. Royama T (1992) Analytical Population Dynamics. Chapman and Hall, London. Saura S, Pascual-Hortal L (2007) A new habitat availability index to integrate connectivity in landscape conservation planning: Comparison with existing indices and application to a case study. Landscape and Urban Planning 83, 91-103. Saura S, Torne J (2009) Conefor Sensinode 2.2: A software package for quantifying the importance of habitat patches for landscape connectivity. Environmental Modelling & Software 24, 135-139. Shih CC, Huang CC, Li SH, Hwang MH, Lee LL (2009) Ten novel tetranucleotide microsatellite DNA markers from Asiatic black bear, Ursus thibetanus. Conservation Genetics 10, 1845-1847. Solberg KH, Bellemain E, Drageset OM, Taberlet P, Swenson JE (2006) An evaluation of field and non-invasive genetic methods to estimate brown bear (Ursus arctos) population size. Biological Conservation 128, 158-168. Taberlet P, Camarra JJ, Griffin S, et al. (1997) Noninvasive genetic tracking of the endangered Pyrenean brown bear population. Molecular Ecology 6, 869-876. Tammeleht E, Remm J, Korsten M, et al. (2010) Genetic structure in large, continuous mammal populations: the example of brown bears in northwestern Eurasia. Molecular Ecology 19, 5359-5370. Urban D, Keitt T (2001) Landscape connectivity: A graph-theoretic perspective. Ecology 82, 1205-1218. 103 Urban DL, Minor ES, Treml EA, Schick RS (2009) Graph models of habitat mosaics. Ecology Letters 12, 260-273. Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes 4, 535-538. van Strien MJ, Keller D, Holderegger R, et al. (2014) Landscape genetics as a tool for conservation planning: predicting the effects of landscape change on gene flow. Ecological Applications 24, 327-339. Wagner AP, Creel S, Kalinowski ST (2006) Estimating relatedness and relationships using microsatellite loci with null alleles. Heredity 97, 336-345. Wright S (1943) Isolation by distance. Genetics 28, 114-138. Wu H, Zhang SN, Wei FW (2010) Twelve novel polymorphic microsatellite loci developed from the Asiatic black bear (Ursus thibetanus). Conservation Genetics 11, 1215-1217. 104 CHAPTER 5 BEYOND THE SNAPSHOT: LANDSCAPE GENETIC ANALYSIS OF TIME SERIES DATA REVEALS RESPONSES OF AMERICAN BLACK BEARS TO LANDSCAPE CHANGE Draheim, H. M., J. A. Moore, M-J., Fortin, K. T. Scribner ABSTRACT Landscape genetic studies typically focus on the spatial pattern and evolutionary process occurring at a single point in time. Although landscape change is widely recognized as a strong driver of microevolutionary processes, few landscape genetic studies have directly evaluated the change in spatial genetic structure (SGS) across generations and the concurrent change in landscape pattern. We introduce time series approach in a landscape genetic framework for a large black bear (Ursus americanus) population in Michigan, USA. We genotyped black bear samples (n = 569) harvested during three years (2002, 2006, and 2010), representing three consecutive generations. We identified areas that were consistently occupied and quantified temporal variation in SGS during this time period. We examined associations between gene flow and landscape change over time, and evaluated effects of ongoing landscape alteration (e.g. deforestation or crop conversion) using spatial autoregressive modeling and model selection to test alternative hypotheses about changes in SGS. Land cover was significantly correlated with gene flow; however, comparisons among sampling years revealed temporal variability in the predictive power and performance of landscape resistance models which resulted in different conclusions. We found the best model to explain temporal inconsistencies was land cover 105 change. Our results demonstrate the importance of conducting long term genetic monitoring for the conservation and management of wildlife inhabiting rapidly changing landscapes. INTRODUCTION With increasing pressure associated with expanding human populations and climate change, current environments have significantly changed and are expected to continue to change (Millennium Ecosysterm Assessment 2005). Ecosystems will undergo significant landscape alterations due to anthropogenic factors such as urbanization, land conversion, and energy production (Hoffmann & Willi 2008). Understanding how landscape alteration will influence species distributions and connectivity is the cornerstone to development of successful conservation, restoration, and management strategies. As natural habitats are altered and fragmented, dispersal and the potential for colonization decreases (Kool et al. 2012). Thus, gene flow and dispersal are likely to be disrupted, leading to increased genetic differentiation, reduced genetic variation and population abundance, and a decreased likelihood of population persistence (Frankham et al. 2002). As a result, landscape change will alter patterns of fine-scale spatial genetic structure (SGS, Epperson 2003). Successful conservation strategies require an understanding of the biotic and abiotic forces that structure populations at multiple spatio-temporal scales (Anderson et al. 2010). Processes affecting the strength and distribution of SGS, such as genetic drift or migration, are both spatially and temporally stochastic (Epperson 1993). Therefore, SGS should be quantified over both space and time. However, traditional population genetic sampling involves a single “snapshot” of genetic data collected at one point in time. Without modeling or accounting for temporal stochasticity, a single snapshot provides an incomplete understanding of natural or 106 anthropogenic processes influencing population genetic structure (Crispo & Chapman 2009). Thus time series analyses based on long term data sets are inherently valuable (Schwartz et al. 2007). Time series data are increasingly applied to understand biological processes, most notably in the field of ecology and population genetics (Lindenmayer et al. 2012; Schwartz et al. 2007). A number of empirical studies have used temporal genetic data to contrast historical and recent genetic diversity using contemporary samples and natural history collections (Wandeler et al. 2007). However, the use of time series analyses as a tool for quantifying temporal changes in population genetic structure on the same population(s) sampled at regular intervals is generally lacking despite the potential value for conservation and management. For example, time series data can be used to examine demographic histories of populations, test for loss of genetic diversity or documenting persistence of spatial genetic structure (Dowling et al. 2014; Waples 1989). While the promise of time series approaches has been recognized in population genetics (Habel et al. 2014) to our knowledge there have been no empirical applications in landscape genetics. Landscape genetics addresses how landscape features influence population genetic structure (Manel et al. 2003). Unlike traditional population genetic models (i.e. Wright’s island model and isolation by distance model; Wright 1943) that assume populations are spatially homogeneous or occupy a uniform landscape (Kimura & Weiss 1964; Wright 1943), a basic tenet of landscape genetics is that natural ecosystems are spatially and temporally heterogeneous. Thus, such heterogeneity may have important ecological implications (Bolliger et al. 2010; Spear et al. 2010). By quantifying concurrent changes in genetic and land cover data over time, we can identify associations that allow for greater insights into underlying processes (Anderson et al. 107 2010). In addition, if the extent and dispersion of samples are consistent among temporal periods, time series data can be used to evaluate the strengths and weaknesses of current landscape genetic analytical approaches to infer patterns beyond a single point in time (Short Bull et al. 2011). Information will allow researchers to evaluate cumulative the impacts of ongoing anthropogenic disturbance and better predict how populations will respond to future landscape change. In this study, we used a novel application of time series analyses to detect changes in spatial genetic structure and relate those changes to underlying landscape change for a population of American black bears (Ursus americanus) in the Northern Lower Peninsula (NLP) of Michigan, USA. The NLP black bears provided a unique opportunity to apply time series analyses in a landscape context. First, by capitalizing on harvest samples, we have acquired a large, long-term data set that spans three consecutive generations of bears. Second, because they occur on a peninsula and are bounded to the south by unsuitable habitat, the NLP black bears are a closed population unaffected by emigration and immigration. Third, Michigan’s NLP is a heterogeneous landscape that has experienced past and current landscape change (i.e. deforestation and agricultural conversion). Fourth, the NLP bear population is subjected to intensive annual harvest (13-29% of the population is harvested annually, Michigan Department of Natural Resources) indicating the potential for complete population turn over in as little as four years. Finally, remotely-sensed land cover data (Coastal Change Analysis Program (CCAP; http://www.csc.noaa.gov/digitalcoast) are publicly available every five years (2001, 2006, 2011) for the entirety of the NLP during the time period over which our genetic samples were collected. 108 We used three temporal genetic data sets with corresponding time series land cover data to address three objectives: (1) Quantify temporal variability in spatial genetic structure between each of three time periods; (2) Quantify the relationship between landscape features and spatial genetic structure for black bears in the NLP within each time period and compare results inferred from each time period to assess consistency of landscape resistance model selection; (3) Quantify associations between landscape change and changes in spatial genetic structure. METHODS Study Area and Sample Collection Our study area covers the northern two-thirds of the lower peninsula of Michigan (47,739 km2) (Figure 5.1). The NLP is a fragmented mosaic composed of different land cover and land use types including development, agriculture, upland non-forested openings, northern hardwood and mixed hardwood, oak, aspen, pine, forested wetland, and non-forested wetland. The NLP experiences ongoing landscape alteration primarily due to forestry practices and crop conversion resulting in many small fragmented patches of land cover change (Figure D1). Sampling occurred during fall harvest (September-October) of 2002 (n = 263), 2006 (n = 385), and 2010 (n = 336). Annually, Michigan’s Department of Natural Resources (MDNR) requires all harvested bears to be registered at hunter check stations. During registration,a pre-molar tooth is extracted for aging and DNA extraction, and hunters report the bear’s harvest location and sex. Locations of bear harvest samples were recorded as township, range, section, which were then georeferenced to the section centroid and converted into UTM coordinates. 109 Figure 5.1. Study area in the Northern Lower Peninsula of Michigan (NLP) showing areas (n = 141) of consistent sampling for black bear harvested during 2002, 2006, and 2010 (n = 569). 110 Microsatellite Genotyping We extracted DNA from bear teeth using Qiagen DNEasy Tissue Kits (Qiagen Inc., Valencia, CA) following manufacturer protocols. We quantified DNA using a Nanodrop spectrophotometer (Thermo Scientific, Waltham, MA) and diluted samples to a 20 ng/μl working concentration. We used PCR to amplify 12 variable microsatellite loci including: G10X, G10L, G10D, G10B, G10M (PCR annealing temperature TA = 58°C; Paetkau et al. 1995) UarMU59, UarMU50 (TA = 58°C, Taberlet et al. 1997), ABB1, ABB4 (TA = 54°C; Wu et al. 2010), UT29, UT35, and UT38 (TA = 54°C; Shih et al. 2009). We amplified DNA according to conditions outlined in Moore et al. 2014. Amplified products were sized on 6.5 % denaturing acrylamide gels for electrophoresis and visualized on a LI-COR 4200 Global IR2 System (LI-COR Inc., Lincoln, NE). All alleles were scored independently by two experienced lab personnel using SAGA genotyping software (LI-COR Inc., Lincoln, NE). To assess genotyping error 10% of samples were randomly selected and genotyped twice to yield a genotyping error rate of < 2%. Defining the Temporal Sampling Locations Harvest samples are opportunistically collected and may vary in spatial dispersion and extent among sampling years which can lead to sampling error. However, the spatial and temporal nature of our data allowed us to account for spatial heterogeneity among years. We defined areas that were consistently occupied in both space and time (sampling from 2002, 2006, and 2010 representing approximately three continuous generations). For each year, these areas were delineated using Voronoi polygons based on the black bear GPS harvest locations. We overlaid the three Voronoi polygon layers from 2002, 2006, and 2010 (Figure D2) and used STAMP 111 (spatial–temporal analysis of moving polygons) in ArcGIS 10.1 to identify and define areas where the three polygon layers overlapped and contained at least one bear from each year. We used program MICRO-CHECKER (Van Oosterhout et al. 2004) to test for the presence of null alleles and allelic dropout. We tested for deviations from Hardy-Weinberg equilibrium using program GENEPOP (Version 3.1d; Raymond and Rousset 1995), and used sequential Bonferroni tests to correct for multiple comparisons (Rice 1989). We used Bonferroni corrections (Goudet 1995, 2001) to test for linkage disequilibrium in program FSTAT 2.93. We quantified microsatellite genetic diversity using mean number of alleles (A), observed heterozygosity (HO), and expected heterozygosity (HE) over all loci, using program GenAlEx (v. 6.0, Peakall and Smouse 2006). To measure genetic differentiation (GD), we calculated the ratio proportion of shared alleles (Dps; Bowcock et al. 1994) for each pairwise combination of individuals using GenAlEx v. 6.0 (Peakall & Smouse 2006). Dps is a commonly used genetic distance estimate in landscape genetic studies and has been shown to accurately reflect population genetic structure and connectivity at small spatial scales (Murphy et al. 2010). Landscape Genetic Analysis We characterized land cover for our sampling periods using three NOAA Coastal Change Analysis Program (CCAP) Land Cover digital coverage maps, derived from Landstat TM imagery for 2001, 2006, and 2011. In order to understand how functional connectivity was influenced by landscape features, we generated resistance surfaces for each temporal land cover data set using Spatial Analyst in ARCGIS 10.1 by weighting each land cover class according to positive or negative associations with black bear presence based on habitat suitability indices (HIS) developed independently from telemetry data by Carter et al. (2010). Landscape resistance 112 surfaces were based on three features (roads, rivers and land cover) that have previously been reported to influence habitat selection by black bears in the NLP (Carter et al. 2010). The models included resistance surfaces using a land cover classification from CCAP Land Cover datasets (resolution = 150 m, 25 classes) which we reclassified based on bear habitat suitability (Carter et al. 2010) (HSI; cost = 1 - 100, least to most resistant to bear movement): mixed deciduous forest (MF = 1), forested wetland (FW = 1), evergreen forest (EF = 10), non-forested upland (NFU = 20), agriculture (AG = 50), non-forested wetland (NFW = 100), and developed (DEV = 100). We also included major rivers (10) and roads weighted based on traffic patterns (Interstate-75 = 100; state roads = 50; all other roads = 10), as they represent putative physical barriers to dispersal. In addition, we created distance only landscape with a raster surface value of 1 bounded by the same spatial extent as all other resistance surfaces to reflect isolation by distance (IBD). Based on those predictors, we proposed 11 alternative landscape hypotheses (Table 5.1). Using the sum function of raster algebra in Spatial Analyst ARCGIS 10.1 we combined land cover with roads, rivers, and roads + rivers. To ensure resistance distances are not dependent upon the relative weight assigned to land cover types (Rayfield et al. 2010), we performed a sensitivity analysis on the range of relative weight between resistance surface land cover types by comparing resistance distances between raster weighted layers using rescaled values according to a suite of weights having either large differences (i.e. maximum weights = 1, 10, 20, 50, 100) or small differences (minimum weights = 0.1, 1, 2, 5, 10) in the relative cost values for land cover types. We calculated pairwise resistance distances using CIRCUITSCAPE 3.5 (McRae 2006). CIRCUITSCAPE incorporates circuit theory to quantify the total landscape resistance between individuals via multiple potential paths of least resistance (McRae & Beier 2007). We used Mantel (Mantel 1967) and partial Mantel tests (Smouse et al. 1986) to correlate genetic and 113 Table 5.1. Mantel (isolation by distance only; IBD) and partial Mantel correlations (r) between spatial and genetic pairwise distances among individual black bear in the NLP for 2002, 2006, and 2010. Bold indicates competitive models based on causal modeling. 2002 Mantel and Partial Mantel test Isolation by Distance r 2006 P r 2010 P r P 0.123 0.013 0.101 0.034 0.110 0.020 -0.060 -0.089 -0.043 -0.051 -0.049 0.990 0.986 0.806 0.858 0.858 -0.027 -0.013 -0.039 -0.009 -0.032 0.672 0.566 0.184 0.509 0.752 -0.033 -0.031 -0.037 0.106 -0.004 0.714 0.738 0.801 0.037 0.440 0.008 0.008 0.019 0.011 0.002 0.026 0.390 0.361 0.219 0.257 0.316 0.200 0.214 0.143 0.192 0.159 0.159 0.187 0.005 0.011 0.005 0.009 0.009 0.005 0.097 0.082 0.101 0.081 0.088 0.104 0.031 0.041 0.027 0.035 0.041 0.026 Resistance models State Roads (STR) Interstate 75 (I75) All Roads Rivers Roads + Rivers Land cover (LC) LC cover only LC + Roads LC + Rivers LC + Roads + Rivers LC + STR LC+ I75 STR = State Roads. I75 = Interstate-75. LC = Land cover 114 geographic/resistance distances. Legendre and Fortin (2010) noted that Mantel and partial Mantel testing is an appropriate framework when hypotheses are explicitly defined in terms of distances, and are appropriate for individual-based analyses. The significance of the matrix correlations was evaluated by permutation tests, based on 10,000 randomizations. The top candidate model corresponded to the largest partial Mantel r values. Next, to eliminate the potential of high Type I error and spurious correlations, we compared a set of competitive models (models with correlation values close to the top model) against each other rather than the null model of IBD. We used an additional causal modeling frame work proposed by Cushman et al. (2013) that is based on the relative support (RS) of each candidate model. The relative support of each model is the difference in the partial Mantel correlation after partialling out the resistance distance of the other model. The model with significant and positive RS values for all comparisons is considered the best candidate model (Cushman et al. 2006; Cushman et al. 2013). Genetic and Landscape Change Analysis Defining genetic change --To assess the relationship between genetic and landscape change over time, we developed a novel approach to analyze data using a combination of GIS and spatial autoregressive modeling. First, we created spatial genetic data layers using the Genetic Landscapes Toolbox (Vandergast et al. 2011) in ARCGIS 10.1. We used Delaunay triangulation to determine neighboring individuals then assigned genetic distances to midpoints of the lines halfway between two individual locations. We interpolated a surface of the genetic distances using Inverse Distance Weighting (IDW; power = 2, variable search radius with 12 points) to form a raster layer. To avoid extrapolating beyond the original collection locations, genetic surfaces were clipped to the spatial extent of collection locations (Vandergast et al. 2013). We 115 calculated the difference in genetic differentiation among two temporal periods using math algebra function in Spatial Analyst ARCGIS 10.1 to obtain a genetic change surface (i.e. genetic differentiation at time 1 – genetic differentiation at time 2). We generated three genetic change surfaces by subtracting the earlier genetic differentiation values from the later period (GD2002GD2006, GD2006-GD2010, GD2002-GD2010). Positive values indicate genetic differentiation was greater during an earlier sampling period, whereas negative values indicate genetic differentiation is greater during the later sampling period. We then overlaid a sampling grid of 2.6 km2 cells and quantified our response variable “cumulative genetic change” (GC) by calculating the cumulative difference of genetic differentiation that fell within a cell. We chose sampling grid cells of 2.6 km2 (size of a section) because they represent sampling uncertainty in the genetic sampling and are larger enough to capture landscape change patterns yet yield sufficiently large sample sizes for modeling. Defining landscape change --To quantify the land cover change among time periods we first reclassified NOAA’s Coastal Change Analysis Program (CCAP) Land Cover change coverage maps (resolution 30 m); 2001-2006; 2006-2011; 2001-2011, to reflect change in black bear habitat suitability among temporal periods. First, we assigned raster cells with no land cover change a value of 0, then weighted areas of change to reflect changes in permeability to black bear movement by calculating the difference in cost from time one to time two. The values can either be positive or negative depending on the type of land cover change (negative = decreased cost or gain of habitat; positive = increased cost or loss of habitat; range, min = -99, max = 99). For example, if a raster cell changed from deciduous forest (low cost to bear movement = 1) to non-forest (high cost to movement = 20) the cell was assigned a landscape change value of 19. Similarly, if a raster cell changed from non-forest to deciduous forest it was assigned a landscape 116 change value of -19. We then overlaid the same sampling grid of 2.6 km2 cells that quantified genetic change (response variable). For the land cover change layers, we used 13 landscape metrics to quantify the configuration and composition of landscape change within each 2.6 km2 sampling cell (definitions of variables and data sources are listed in Table 5.2): (i) number of patches (i.e., contiguous grid cells) of landscape change (NP); (ii) proportion of area with landscape change (PLAND); (iii) aggregation index (AI), which measures the extent to which land cover change is aggregated or dispersed; (iv) patch cohesion index (COH), which quantifies connectivity of land cover change; (v) nearest neighbor distance (NNDIST), which is the shortest distance from one landscape change patch to the next (Table 5.2) in program FRAGSTATS (McGarigal et al. 2002). In addition, using ARCGIS 10.1 we calculated the cumulative difference of landscape change within each 2.6 km2 sampling cell (DEG) (this measure captures overall degree of change) and between neighboring sampling cells (NDEG). Because deforestation is the predominate type of landscape change we also described landscape change in terms of proportion of area with forest loss (FL). Also, because bear presence is not significantly correlated with coniferous stands (Carter et al. 2010) we described landscape change in terms of proportion of deciduous forest and mixed forest loss (DMFL). However, genetic structure can have a substantial time lag in its response to changes in gene flow resulting from a barrier. Thus, we also calculated the distance to landscape features that may act as a putative barrier to functional connectivity: (i) distance to major human population center (i.e. city or town; HPDIST); (ii) distance to major road (RDDIST); (iii) distance to major river (RVDIST). Spatial Modeling Analyses— We compared multiple spatial autoregressive lag models to determine whether landscape change best predicted genetic change. These models are linear 117 regressions that account for spatial non-independence in the response variable (genetic change) layer within an additional spatial lag term as an explanatory variable (Ward & Gleditsch 2008). Models without significant regression coefficients were discarded. Because the inherent relationships between landscape characteristics NP, PLAND, DEG, FL, and DMFL (Table 5.2), univariate models with these explanatory variables were compared using Akaike's information criterion (AIC) (Burnham & Anderson 2002; Johnson & Omland 2004), where Akaike weights represent the probability that a model is the best of the set (Ward & Gleditsch 2008). The best landscape alteration variable of the univariate regression models was then used in multiple regression analyses. All tests were performed in spdep package in R v. 3.0.1 (R Development Core Team 2013). RESULTS Spatial Genetic Structure We identified 141 areas of consistent occupancy across the three time periods, and thus retained 569 individuals within those areas for subsequent landscape genetic analyses (2002, n = 204; 2006, n = 199; 2010, n= 166). We did not find evidence for null alleles or allelic dropout and no loci were found to deviate significantly from Hardy-Weinberg or linkage equilibrium, so all 12 loci were retained for further analyses. For all loci, expected heterozygosity ranged from 0.67 to 0.94, number of alleles per locus ranged from 6 to 28 and genetic diversity was similar among years (Table D1). Inter-individual genetic distances were correlated with Euclidean geographic distances for all sampling years (2002, r = 0.123, P = 0.01, 2006, r = 0.101, P = 0.03; 2010, r = 0.106, P = 0.02; Figure 5.2). 118 Table 5.2. A summary of the landscape variables used to characterize genetic change between 2002, 2006, and 2010 for black bears in the Northern Lower Peninsula of Michigan, USA. Variable Habitat Alteration Quantity of Change Percent Area of Change Degree of Change Deforestation Deforestation (deciduous or mixed) Neighborhood Nearest Neighborhood Change Nearest Neighbor Distance Deforestation (all forest) Deforestation (deciduous and mixed) Abbreviation Description Program NP PLAND DEG FL Total number of landscape change patches Percent of area occupied by landscape change Cumulative difference in landscape change Percent of area occupied by deforestation FRAGSTATS FRAGSTATS ArcGIS FRAGSTATS DMFL Percent of area occupied by loss of deciduous or mixed forest FRAGSTATS Cumulative difference in landscape change in neighboring sampling cells Distance to nearest landscape change patch FRAGSTATS Percent of area occupied by deforestation FRAGSTATS Percent of area occupied by habitat loss FRAGSTATS Extent to which land cover change is aggregated Connectivity of landscape change FRAGSTATS FRAGSTATS NNDEG NNDIST FL DMFL Heterogeneity Aggregation Cohesion AI COH Distance to Barriers 119 Table 5.2 (cont’d) Human Population (city) Roads Rivers HPDIST RDDIST RVDIST Distance to nearest human population center Distance to nearest major road Distance to major river 120 ArcGIS ArcGIS ArcGIS Figure 5.2. Scatterplot of genetic distance versus geographic distance including least squares regression line (2002, P = 0.01; 2006, P = 0.03; 2010, P = 0.02; Mantel test). Genetic distance was based on pairwise estimates of proportion of shared alleles (dps) among all individuals. Geographic distances (km) were calculated as the straight-line distance from the center of each sampling location. 2002 2006 2010 121 Landscape Genetic Analyses Landscape genetic analyses revealed discrepancies in landscape resistance model performance among years. Only 2006 and 2010 had statistically significant landscape resistance models (P<0.05) after partialling out the effects of geographic distance (Table 5.1). The best supported models differed among years. For 2006, the best supported model included land cover only (LC) based on the highest partial Mantel correlation (r = 0.214, P = 0.005; Table 5.1) and relative support (RS). In contrast, the resistance model with the highest correlation for 2010 included land cover and interstate-75 (I-75) as a barrier (partial Mantel, r = 0.104, P = 0.026; Table 5.1); however, causal modeling revealed this model was not significant after partialling out three competing models. Therefore, we are unable to conclude which single model is the best supported model for this year. Equally well supported were three alternative models: (1) LC only (r = 0.097, P = 0.031) (2); LC + rivers (r = 0.101, P = 0.027); and (3) LC + state roads (STR) (r = 0.088, P = 0.041; Table 5.1). All competitive models included land cover, confirming that land cover is important for black bear connectivity in the NLP. In addition, landscape genetic analyses were not sensitive to magnitude of land cover cost values. Genetic Change and Landscape Change “Landscape alteration” univariate modeling predicting genetic change revealed the number of landscape change patches (NP) within a sampling unit was the most supported model for all temporal comparisons (2002-2006, wAIC = 0.48; 2006-2010, wAIC = 0.65; 2002-2010, wAIC = 0.29; Table 5.3; Figure 5.3). Therefore, NP was used in multiple regression analyses. Multivariate models with significant coefficients are described in Table 5.4. All temporal 122 Figure 5.3. Distribution of (A) number of landscape change patches within a grid cell and (B) the cumulative difference in genetic differentiation within a cell grid for all temporal comparisons (2002-2006, 2006-2010, 2002-2010). For landscape change maps red indicates large number of patches of change and blue indicates low number of patches of change. For genetic change maps red indicates increase in genetic differentiation (GD) from time one (T1) to time two (T2) and blue indicates decrease in genetic differentiation from T1 to T2. 2006-2010 2002-2010 Genetic change 2002-2006 A 0 No Change Landscape change -1 T1GD < T2GD 1 T1GD < T2GD B 0 81 123 Table 5.3. Most supported univariate spatial regression models of landscape alteration variables predicting black bear genetic change. Abbreviations are as described in Table 1. 2002-2006 Model (GC ∼) βNP βPLAND βDEG βHL βFL Intercept 0.353 0.212 0.194 0.255 0.213 Coeff. Z-value 0.011 -0.045 -0.112 -0.398 -0.052 2.634 1.125 -0.734 -2.595 -1.266 Coeff. Z-value Ρ AIC 0.955 0.955 0.955 0.955 0.955 17442 17448 17448 17443 17447 Ρ AIC 0.872 2.131 0.871 0.871 0.872 18667 18670 18674 18670 18672 ρ AIC 0.959 0.959 0.959 0.959 0.959 16535 16536 16537 16537 16537 ΔAIC wAIC 0 6 6 1 5 0.48 0.03 0.19 0.03 0.43 ΔAIC wAIC 0 3 7 3 5 0.65 0.12 0.02 0.06 0.16 ΔAIC wAIC 2006-2010 Model (GC ∼) βNP βPLAND βDEG βHL βFL Intercept -0.613 -0.470 -2.465 -2.467 -0.412 0.010 0.092 0.188 0.273 0.091 2.806 2.132 0.792 2.249 1.789 2002-2010 Model (GC ∼) βNP βPLAND βDEG βHL βFL Intercept 0.007 0.128 0.111 0.103 0.131 Coeff. Z-value 0.011 -0.025 -0.098 -0.039 -0.032 0.539 -1.185 -0.883 -0.776 -1.354 0 1 2 2 2 0.29 0.23 0.17 0.16 0.13 Coefficients (Coeff) and their corresponding Z-values refer to genetic change layer coefficients, whereas ρ is the spatial lag coefficient. All values of ρ were significant (P < 0.01). AIC, ΔAIC, and weighted (w)AIC values are reported . The best model is in boldface type. 124 comparisons that included 2010, genetic change was better predicted by inclusion of landscape change heterogeneity (Table 5.4). The best model for 2006-2010 included NP and extent of heterogeneity calculated as AI (wAIC = 0.95). Similarly, the most probable model in the 20022010 comparison included NP, AI, and an additional heterogeneity variable COH (2002-2010, wAIC = 1.0 Table 5.4). DISCUSSION As the field of landscape genetics advances beyond characterizing static landscape patterngenetic process relationships, incorporating landscape change and quantifying its effects on spatial genetic structure over time is vital to improve understanding of anthropogenic impacts and functional connectivity (Bolliger et al. 2010; Segelbacher et al. 2010; Storfer et al. 2010). Our study illustrates the advantages of time-series genetic and landscape data when analyzed using established landscape genetics methods. We have taken advantage of a long term genetic data set and land cover imagery covering three generations of black bears to show inconsistencies in landscape resistance models, using different “snapshots”. Our method showed that black bear spatial genetic structure is strongly influenced by landscape change, which has broad implications for conservation and management. We found isolation by landscape resistance was more strongly supported than IBD. Land cover was the landscape factor most strongly correlated with genetic distance in black bears in the NLP after partialling out geographic distance. In contrast, roads and rivers alone did not appear to significantly influence genetic distance. However, the univariate model which includes rivers only was significantly correlated with genetic differentiation in 2010. Our results are consistent with a previous study that performed similar analyses on black bears sampled during 125 Table 5.4. Most supported multiple spatial regression models of landscape change characteristics predicting black bear genetic change.. 2002-2006 Model (GC ∼) βNP ρ 0.955 AIC ΔAIC wAIC ΔAIC wAIC 0 4 1 0 ΔAIC wAIC 17442 2006-2010 Model (GC ∼) βNP +βAI βNP ρ 0.93415 0.959 AIC 16533 16537 2006-2010 Model (GC ∼) βNP+ βAI + βCOH βNP +βAI βNP βAI + βCOH βAI ρ AIC 0.900 16683 0.951 0.872 0.945 0.945 16689 16886 16899 17001 0 6 203 0.95 0.05 0 216 318 0 0 ρ is the spatial lag coefficient. All values of ρ were significant (P < 0.01). AIC, ΔAIC, and weighted (w)AIC values are reported . The best model is in boldface type. 126 2006 but used a much larger sample size (N = 365) and spatial distribution (Chapter 3). The influence of land cover is not surprising as land cover is fundamentally tied to resources such as food availability and forest cover for security and resting (Amstrup & Beecham 1976; Davis et al. 2006; Noyce & Garshelis 2011; Rogers 1977) We found temporal inconsistencies in support of our landscape resistance models which resulted in different conclusions. Had analyses been based solely on samples collected in 2002 we would conclude that land cover has no influence on functional connectivity of black bears in the NLP. Nevertheless, we detected a strong association between black bear SGS and land cover in 2006 and 2010. Results suggest variations in habitat quality and subsequent different habitat permeabilities either facilitate or impede black bear gene flow. We propose the cumulative effects of landscape alteration over increasing lengths of time likely explains the differences in relative effects of landscape features on spatial genetic structure in NLP black bears. Today the NLP is a patchy landscape with areas of deciduous, mixed and coniferous forest fragmented by areas of intensive agriculture or human development. If current patterns of human land-use continue, the remaining forested areas could undergo further modification in the near future (Public Service Consultants 2001). Based on our landscape change raster layers approximately 1.3% of black bear habitat is being lost every 5 years (20022006 = 0.9%, 2006-2010 = 1.6%; Figure D1). On the whole, this habitat loss may not seem extensive; however, the effects of habitat loss can be amplified by concurrent fragmentation (Fahrig 1997) resulting in an increase of dispersed, complex, habitat patches. We found significant effects of landscape change on measures of change in genetic differentiation (GC; Tables 5.3 and 5.4). The magnitude of landscape change within a sampling grid cell was predictive of the genetic change. However, the amount of variation explained in models that 127 included variables associated with habitat heterogeneity were the best predictors of genetic change for temporal comparisons that included 2010. This discrepancy among years reflects an increasing amount of habitat fragmentation due to deforestation over the nine year sampling period. Therefore, lack of an association between land cover and spatial genetic structure in bears sampled during 2002 could be due to the degree of land cover heterogeneity present in 2002 which may have been insufficient to cause a response. Our result is consistent with a recent study in the Rocky Mountains of the United States that evaluated the influence of forest cover on black bear gene flow based on twelve landscapes that varied in degree of fragmentation. Short Bull et al. (2011) found higher correlations between land cover and black bear gene flow in landscapes where forest cover was highly fragmented compared to landscapes of contiguous forest. Our temporal analyses may also have been influenced by sampling error due to the use of harvest data. However, this is unlikely as our temporal sampling approach and use of consistently occupied areas controlled for spatial heterogeneity among years. In addition, Chapter 1 results shows spatial dispersion of harvested bears sampled in the NLP is concordant across years. Also, our samples sizes for this study are suitable (N = 166-204), thus bias due to small sample size is minor. Black bears in the NLP may be experiencing rapid rates of population turn over and subsequent local population fluctuations due to intensive annual harvest (~13-29% of the population harvested annually between 2002-2010, MDNR unpublished data). However, we are unsure as to how genetic drift is contributing to our results. Further investigations using local harvest abundance as a proxy for bear density (e.g. Moore et al. 2014) are needed to address this question. While our method may be broadly applicable across a range of taxa, the approach needs to be parameterized based on specific life history characteristics of the species or population of 128 interest. Our findings, while supported in Michigan black bears in the NLP, should be evaluated empirically for other species and locales. For example, we would predict similar results based on the degree of concordance between our results and previous landscape genetic studies of black bears in the Rocky Mountains (Short Bull et al. 2011). However, there are a number of possible factors that could confound results in another bear population. For example, the NLP is an isolated population thus we could assume changes in spatial genetic structure were not due to immigration or emigration. Also, it is important that the rate of landscape change be sufficient to influence spatial genetic structure over relatively short time intervals. To the best of our knowledge our study provides the first time series landscape genetic analyses performed across multiple generations of the same population. Here we have shown the importance of time series data to test for consistencies in landscape-genetic relationships over time. Our ability to relate gene flow to landscape features is largely dependent on the temporal scale of sampling. Our finding that genetic change is associated with landscape change highlights the synergistic effects of habitat loss and fragmentation on black bear gene flow. Our data enables managers to target regions or habitat types that are important for maintaining connectivity across anthropogenically-altered habitats and assess impacts of future change. Our analysis lays the groundwork for using landscape genetic analyses to predict the genetic implications of impending land use/cover change. For example, our results can be used to parametrize simulations to model different landscape genetic conditions using multiple resistance grids that mimic possible landscape change scenarios to predict changes in genetic connectivity or estimate the time required for the spatial patterns of genetic structure to equilibrate (Epperson et al. 2010; Landguth & Cushman 2010). Also, as a harvested species we would be able to simulate different rates of bear mortality using different harvest forecasts coupled with different 129 landscape configurations that will mimic varying anthropogenic pressures and evaluate their effects on spatial genetic structure. ACKNOWLEDGMENTS We are grateful to J. Fierke, K. Filcek, B. Dreher, and S. Libants for assisting in data generation. We also thank J. Messina, B. Epperson, K. Holekamp, J. Kanefsky, J. McGuire, for insightful discussion and helpful comments on the analysis and manuscript. All samples used in analyses were provided by Michigan bear hunters and collected by the many MDNR staff at registration stations. Support for this project was provided by the Michigan Department of Natural Resources through the Wildlife and Sportfish Restoration Program F11AF00640 and by the Department of Fisheries and Wildlife at Michigan State University. 130 APPENDIX 131 APPENDIX Table D1. Black bear genetic diversity at 12 microsatellite loci. No loci showed significant deviations from Hardy-Weinberg equilibrium (Raymond & Rousset 1995). Na, number of alleles per locus; Ho, observed hererozygosity; HE, expected heterozygosity. Locus X M D L B Uar50 Uar59 ABB1 ABB4 UT35 UT38 UT29 Locus X M D L B Uar50 Uar59 ABB1 ABB4 UT35 UT38 UT29 2002 (n = 204) Na HO 8 9 9 10 6 8 6 8 7 6 21 12 0.658 0.744 0.776 0.892 0.709 0.756 0.688 0.710 0.789 0.760 0.951 0.859 2006 (n = 199) Na HO 15 10 11 10 8 10 9 9 8 7 28 12 0.668 0.793 0.790 0.818 0.715 0.656 0.667 0.689 0.823 0.773 0.888 0.850 HE 0.617 0.757 0.752 0.856 0.737 0.783 0.707 0.711 0.795 0.771 0.908 0.844 HE 0.702 0.787 0.762 0.860 0.727 0.787 0.677 0.737 0.798 0.784 0.898 0.816 132 Table D1 (cont’d) Locus X M D L B Uar50 Uar59 ABB1 ABB4 UT35 UT38 UT29 2010 (n = 166) Na HO 9 10 11 12 6 9 7 10 10 7 22 11 0.776 0.759 0.735 0.878 0.680 0.747 0.721 0.738 0.822 0.753 0.937 0.911 HE 0.721 0.765 0.754 0.862 0.716 0.816 0.691 0.742 0.790 0.776 0.903 0.941 133 Figure D1. Extent and configuration of forested land lost (green) over the sampling period from (A) 2001-2006, (B) 2006-2010, and (C) 2002-2010. A B C 134 Figure D2. Map of overlapping Voronoi tessellation point pattern polygons in NLP. Polygons contained at least one sample from 2002, 2006, and 2010. 135 REFERENCES 136 REFERENCES Allendorf FW, England PR, Luikart G, Ritchie PA, Ryman N (2008) Genetic effects of harvest on wild animal populations. Trends in Ecology & Evolution 23, 327-337. Allendorf FW, Luikart G (2012) Conservation and the Genetics of Populations, 2nd edn. WileyBlackwel. Allendorf FW, Luikart GH, Aitken SN (2012) Conservation and the genetics of populations Wiley. com. Amarasekare P (2004) The role of density-dependent dispersal in source–sink dynamics. Journal of Theoretical Biology 226, 159-168. Amstrup SC, Beecham J (1976) Activity patterns of radio-collared black bears in idaho. The Journal of Wildlife Management 40, 340-348. Anderson CD, Epperson BK, Fortin MJ, et al. (2010) Considering spatial and temporal scale in landscape-genetic studies of gene flow. Molecular Ecology 19, 3565-3575. Anderson JT, Geber MA (2010) Demographic source-sink dynamics restrict local adaptation in Elliott’s blueberry (Vaccinium elliottii). Evolution 64, 370-384. Bartholomé E, Belward A (2005) GLC2000: a new approach to global land cover mapping from Earth observation data. International Journal of Remote Sensing 26, 1959-1977. Beier P, Majka DR, Spencer WD (2008) Forks in the road: choices in procedures for designing wildland linkages. Conservation Biology 22, 836-851. Beier P, Noss RF (1998) Do habitat corridors provide connectivity? Conservation Biology 12, 1241-1252. Bolliger J, Le Lay G, Holderegger R (2010) Landscape genetics - how landscapes affect ecological processes. Gaia-Ecological Perspectives for Science and Society 19, 238-240. Boughton DA (1999) Empirical evidence for complex source-sink dynamics with alternative states in a butterfly metapopulation. Ecology 80, 2727-2739. Bowcock AM, Ruiz-Linares A, Tomfohrde J, et al. (1994a) High resolution of human evolutionary trees with polymorphic microsatellites. Nature 368, 455-457. Bowcock AM, Ruizlinares A, Tomfohrde J, et al. (1994b) High resolution of human evolutionary threes with polymorphic microsatellites. Nature 368, 455-457. Brown SK, Hull JM, Updike DR, Fain SR, Ernest HB (2009) Black bear population genetics in california: signatures of population structure, competitive release, and historical translocation. Journal of Mammalogy 90, 1066-1074. 137 Bunn AG, Urban DL, Keitt TH (2000) Landscape connectivity: A conservation application of graph theory. Journal of Environmental Management 59, 265-278. Burnham KP, Anderson DR (2002) Model selection and multimodel inference: a practical information-theoretic approach Springer. Carter NH (2007) Predicting ecological and social suitability of black bear habitat in Michigan's Lower Peninsula, University of Michigan. Carter NH, Brown DG, Etter DR, Visser LG (2010) American black bear habitat selection in northern Lower Peninsula, Michigan, USA, using discrete-choice modeling. Ursus 21, 57-71. Cattadori IM, Haydon DT, Thirgood SJ, Hudson PJ (2003) Are indirect measures of abundance a useful index of population density? The case of red grouse harvesting. Oikos 100, 439446. Clobert J, Le Galliard JF, Cote J, Meylan S, Massot M (2009) Informed dispersal, heterogeneity in animal dispersal syndromes and the dynamics of spatially structured populations. Ecology Letters 12, 197-209. Congalton RG, Green K (2008) Assessing the accuracy of remotely sensed data: principles and practices CRC press. Costello CM (2010) Estimates of Dispersal and Home-Range Fidelity in American Black Bears. Journal of Mammalogy 91, 116-121. Coulon A, Cosson JF, Angibault JM, et al. (2004) Landscape connectivity influences gene flow in a roe deer population inhabiting a fragmented landscape: an individual-based approach. Molecular Ecology 13, 2841-2850. Crispo E, Chapman LJ (2009) Temporal variation in population genetic structure of a riverine African cichlid fish. Journal of Heredity, esp078. Cushman SA (2006) Effects of habitat loss and fragmentation on amphibians: A review and prospectus. Biological Conservation 128, 231-240. Cushman SA, McKelvey KS, Hayden J, Schwartz MK (2006) Gene flow in complex landscapes: Testing multiple hypotheses with causal modeling. American Naturalist 168, 486-499. Cushman SA, Wasserman TN, Landguth EL, Shirk AJ (2013) Re-evaluating causal modeling with Mantel tests in landscape genetics. Diversity 5, 51-72. Dale MRT, Fortin M-J (2010) From graphs to spatial graphs. Annual Review of Ecology, Evolution, and Systematics 41, 21-38. Davis H, Weir RD, Hamilton AN, Deal JA (2006) Influence of phenology on site selection by female American black bears in coastal British Columbia. Ursus 17, 41-51. De Sherbinin A, Carr D, Cassels S, Jiang L (2007) Population and environment. Annual review of environment and resources 32, 345. 138 Dias PC (1996) Sources and sinks in population biology. Trends in Ecology & Evolution 11, 326-330. Dixo M, Metzger JP, Morgante JS, Zamudio KR (2009) Habitat fragmentation reduces genetic diversity and connectivity among toad populations in the Brazilian Atlantic Coastal Forest. Biological Conservation 142, 1560-1569. Dowling TE, Turner TF, Carson EW, et al. (2014) Time-series analysis reveals genetic responses to intensive management of razorback sucker (Xyrauchen texanus). Evolutionary Applications 7, 339-354. Draheim HM, Lopez V, Winterstein SR, Etter DR, Scribner KT (in prep) Effects of spatial and temporal sampling scale on sample dispersion and spatial genetic structure. Earl DA (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conservation Genetics Resources 4, 359-361. Epperson BK (1993) Spatial and Space-Time Correlations in Systems of Subpopulations with Genetic Drift and Migration. Genetics 133, 711-727. Epperson BK (2005) Estimating dispersal from short distance spatial autocorrelation. Heredity 95, 7-15. Epperson BK, McRae BH, Scribner K, et al. (2010) Utility of computer simulations in landscape genetics. Molecular Ecology 19, 3549-3564. Etherington TR (2011) Python based GIS tools for landscape genetics: visualising genetic relatedness and measuring landscape connectivity. Methods in Ecology and Evolution 2, 52-55. Evanno G, Regnaut S, Goudet J (2005) Detecting the number of clusters of individuals using the software structure: a simulation study. Molecular Ecology 14, 2611-2620. Fahrig L (1997) Relative effects of habitat loss and fragmentation on population extinction. The Journal of Wildlife Management 61, 603-610. Fahrig L, Merriam G (1994) Conservation of fragmented populations. Conservation Biology 8, 50-59. Foody GM (2002) Status of land cover classification accuracy assessment. Remote Sensing of Environment 80, 185-201. Frankham R, Ballou JD, Briscoe DA (2002) Introduction to conservation genetics Cambridge University Press. Frawley BJ (2001) 2000 Michigan black bear hunter survey, p. 14. MDNR, Wildl. Div. Report. Frawley BJ (2009) 2008 Michigan black bear hunter survey, p. 24. MDNR, Wildl. Div. 139 Friedl MA, McIver DK, Hodges JC, et al. (2002) Global land cover mapping from MODIS: algorithms and early results. Remote Sensing of Environment 83, 287-302. Funk WC, Blouin MS, Corn PS, et al. (2005) Population structure of Columbia spotted frogs (Rana luteiventris) is strongly affected by the landscape. Molecular Ecology 14, 483-496. Gaggiotti OE (1996) Population genetic models of source-sink metapopulations. Theoretical Population Biology 50, 178-208. Garroway CJ, Bowman J, Wilson PJ (2011) Using a genetic network to parameterize a landscape resistance surface for fishers, Martes pennanti. Molecular Ecology 20, 3978-3988. Garshelis DL, Pelton MR (1981) Movements of black bears in the Great Smoky Mountains National Park. The Journal of Wildlife Management 45, 912-925. Goudet J (1995) FSTAT (Version 1.2): A computer program to calculate F-statistics. Journal of Heredity 86, 485-486. Goudet J (2001) FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3). Guillot G, Estoup A, Mortier F, Cosson JF (2005) A spatial statistical model for landscape genetics. Genetics 170, 1261-1280. Habel JC, Husemann M, Finger A, Danley PD, Zachos FE (2014) The relevance of time series in molecular ecology and conservation biology. Biological Reviews 89, 484-492. Hoffman JI, Matson CW, Amos W, Loughlin TR, Bickham JW (2006) Deep genetic subdivision within a continuously distributed and highly vagile marine mammal, the Steller's sea lion (Eumetopias jubatus). Molecular Ecology 15, 2821-2832. Hoffmann AA, Willi Y (2008) Detecting genetic responses to environmental change. Nat Rev Genet 9, 421-432. Homer C, Dewitz J, Fry J, et al. (2007) Completion of the 2001 National Land Cover Database for the Counterminous United States. Photogrammetric Engineering and Remote Sensing 73, 337. Horner MA, Powell RA (1990) Internal structure of home ranges of black bears and analyses of home-range overlap. Journal of Mammalogy 71, 402-410. Johnson JB, Omland KS (2004) Model selection in ecology and evolution. Trends in Ecology & Evolution 19, 101-108. Kalinowski ST, Wagner AP, Taper ML (2006) ml‐relate: a computer program for maximum likelihood estimation of relatedness and relationship. Molecular Ecology Notes 6, 576579. Kawecki T (1995) Demography of source—sink populations and the evolution of ecological niches. Evolutionary Ecology 9, 38-44. 140 Kawecki TJ (2008) Adaptation to marginal habitats. Annual Review of Ecology, Evolution, and Systematics 39, 321-342. Keller D, Holderegger R, van Strien MJ, Bolliger J (2014) How to make landscape genetics beneficial for conservation management? Conservation Genetics, 1-10. Kelly RP, Oliver TA, Sivasundar A, Palumbi SR (2010) A method for detecting population genetic structure in diverse, high gene-flow species. Journal of Heredity 101, 423-436. Kimura M, Weiss GH (1964) The stepping stone model of population structure and the decrease of genetic correlation with distance. Genetics 49, 561. Kitson JC (2004) Harvest rate of sooty shearwaters (Puffinus griseus) by Rakiura Māori: a potential tool to monitor population trends? Wildlife Research 31, 319-325. Knowlton JL, Graham CH (2010) Using behavioral landscape ecology to predict species' responses to land-use and climate change. Biological Conservation 143, 1342-1354. Kool J, Moilanen A, Treml E (2012) Population connectivity: recent advances and new perspectives. Landscape Ecology, 1-21. Kopatz A, Eiken HG, Hagen SB, et al. (2012) Connectivity and population subdivision at the fringe of a large brown bear (Ursus arctos) population in North Western Europe. Conservation Genetics 13, 681-692. Lande R (1988) Genetics and demography in biological conservation. Science(Washington) 241, 1455-1460. Lande R, Engen S, Saether B-E (2003) Stochastic population dynamics in ecology and conservation Oxford University Press Oxford. Landguth EL, Cushman SA (2010) CDPOP: A spatially explicit cost distance population genetics program. Molecular Ecology Resources 10, 156-161. Lee DJ, Vaughan MR (2003) Dispersal Movements by Subadult American Black Bears in Virginia. Ursus 14, 162-170. Legendre P, Fortin M-J (2010) Comparison of the Mantel test and alternative approaches for detecting complex multivariate relationships in the spatial analysis of genetic data. Molecular Ecology Resources 10, 831-844. Lenormand T (2002) Gene flow and the limits to natural selection. Trends in Ecology & Evolution 17, 183-189. Lindenmayer DB, Likens GE, Andersen A, et al. (2012) Value of long‐term ecological studies. Austral Ecology 37, 745-757. Liu W, Gopal S, Woodcock CE (2004) Uncertainty and confidence in land cover classification using a hybrid classifier approach. Photogrammetric Engineering and Remote Sensing 70, 963-972. 141 Lowe WH, Allendorf FW (2010) What can genetics tell us about population connectivity? (vol 19, pg 3038, 2010). Molecular Ecology 19, 5320-5320. Manel S, Holderegger R (2013) Ten years of landscape genetics. Trends in Ecology & Evolution 28, 614-621. Manel S, Joost S, Epperson BK, et al. (2010) Perspectives on the use of landscape genetics to detect genetic adaptive variation in the field. Molecular Ecology 19, 3760-3772. Manel S, Schwartz MK, Luikart G, Taberlet P (2003) Landscape genetics: combining landscape ecology and population genetics. Trends in Ecology & Evolution 18, 189-197. Mantel N (1967) Detection of disease clustering and a generalized regression approach. Cancer Research 27, 209-&. Manville AM (1987) Den selection and use by black bears in Michigan’s Northern Lower Peninsula. International Conference of Bear Research and Management 7, 317-322. McGarigal K, Cushman SA, Neel MC, Ene E (2002) FRAGSTATS: spatial pattern analysis program for categorical maps. McRae BH (2006) Isolation by resistance. Evolution 60, 1551-1561. McRae BH, Beier P (2007) Circuit theory predicts gene flow in plant and animal populations. Proceedings of the National Academy of Sciences of the United States of America 104, 19885-19890. McRae BH, Beier P, Dewald LE, Huynh LY, Keim P (2005) Habitat barriers limit gene flow and illuminate historical events in a wide-ranging carnivore, the American puma. Molecular Ecology 14, 1965-1977. MDNR (2004) Review of remote sensing technology used in the IFMAP project, pp. 1-49. Michigan Department of Natural Resources. Millennium Ecosysterm Assessment (2005) Ecosystems and human well-being: desertification synthesis. In: Millennium Ecosysterm Assessment. World Resources Institute, Washington DC. Minor ES, Urban DL (2007) Graph theory as a proxy for spatially explicit population models in conservation planning. Ecological Applications 17, 1771-1782. Minor ES, Urban DL (2008) A graph-theory frarmework for evaluating landscape connectivity and conservation planning. Conservation Biology 22, 297-307. Mitchell MS, Powell RA (2007) Optimal use of resources structures home ranges and spatial distribution of black bears. Animal Behaviour 74, 219-230. Moore JA, Draheim HM, Etter D, Winterstein S, Scribner KT (2014) Application of large-scale parentage analysis for investigating natal dispersal in highly vagile vertebrates: a case study of american black bears (Ursus americanus). Plos One 9, e91168. 142 Murphy MA, Dezzani R, Pilliod DS, Storfer A (2010) Landscape genetics of high mountain frog metapopulations. Molecular Ecology 19, 3634-3649. NOAA (2014) Western Great Lakes 2010 Coastal Change Analysis Program accuracy assessment, pp. 1-11, Charleston, South Carolina. Noyce K, Garshelis D (2011) Seasonal migrations of black bears (Ursus americanus): causes and consequences. Behavioral Ecology and Sociobiology 65, 823-835. Okabe A, Boots B, Sugihara K, SN C (2000) Spatial tessellations: concepts and applications of voronoi diagrams., 2nd ed. edn. Wiley, Chichester. Opdam P, Wascher D (2004) Climate change meets habitat fragmentation: linking landscape and biogeographical scale levels in research and conservation. Biological Conservation 117, 285-297. Paetkau D, Calvert W, Stirling I, Strobeck C (1995) Microsatellite analysis of populationstructure in Canadian polar bears. Molecular Ecology 4, 347-354. Paetkau D, Waits LP, Clarkson PL, Craighead L, Strobeck C (1997) An empirical evaluation of genetic distance statistics using microsatellite data from bear (Ursidae) populations. Genetics 147, 1943-1957. Peakall R, Smouse PE (2006) GENALEX 6: genetic analysis in Excel. Population genetic software for teaching and research. Molecular Ecology Notes 6, 288-295. Peery MZ, Beissinger SR, House RF, et al. (2008) Characterizing source-sink dynamics with genetic parentage assignments. Ecology 89, 2746-2759. Pelton MR (1992) Black Bear (Ursus americanus) In: Wild Mammals of North America: Biology, Management, and Economics (eds. Chapman JA, Feldhamer GA), pp. 504-514. The John Hopkins University Press,USA. Peterman W, Rittenhouse TG, Earl J, Semlitsch R (2013) Demographic network and multiseason occupancy modeling of Rana sylvatica reveal spatial and temporal patterns of population connectivity and persistence. Landscape Ecology, 1-13. Pierson JC, Allendorf FW, Drapeau P, Schwartz MK (2013) Breed locally, disperse globally: fine-scale genetic structure despite landscape-scale panmixia in fire-specialist. Plos One 8, e67248. Pritchard JK, Stephens M, Donnelly P (2000) Inference of population structure using multilocus genotype data. Genetics 155, 945-959. Proctor MF, McLellan BN, Strobeck C, Barclay RMR (2005) Genetic analysis reveals demographic fragmentation of grizzly bears yielding vulnerably small populations. Proceedings of the Royal Society B-Biological Sciences 272, 2409-2416. 143 Public Service Consultants (2001) Michigan land resource project. Prepared for the Frey Foundation and WK Kellogg Foundation on behalf of the Michigan Economic and Environmental Roundtable (ed. Consultants PS), Lansing, Michigan, USA. Pulliam HR (1988) Sources, sinks, and population regulation. American Naturalist 132, 652-661. Queller DC, Goodnight KF (1989) Estimating relatedness using genetic markers. Evolution 43, 258-275. Rabinowitz A, Zeller KA (2010) A range-wide model of landscape connectivity and conservation for the jaguar, Panthera onca. Biological Conservation 143, 939-945. Rayfield B, Fortin M-J, Fall A (2010) The sensitivity of least-cost habitat graphs to relative cost surface values. Landscape Ecology 25, 519-532. Rayfield B, Fortin M-J, Fall A (2011) Connectivity for conservation: a framework to classify network measures. Ecology 92, 847-858. Raymond M, Rousset F (1995) GENEPOP (VERSION-1.2) - Population-genetics software for exact tests and ecumenicism. Journal of Heredity 86, 248-249. Reding D, Cushman S, Gosselink T, Clark W (2013) Linking movement behavior and fine-scale genetic structure to model landscape connectivity for bobcats (Lynx rufus). Landscape Ecology 28, 471-486. Rice WR (1989) Analyzing tables of statistical tests. Evolution 43, 223-225. Robertson C, Nelson T, Boots B, Wulder M (2007) STAMP: spatial–temporal analysis of moving polygons. Journal of Geographical Systems 9, 207-227. Robinson HS, Wielgus RB, Cooley HS, Cooley SW (2008) Sink populations in carivore managment: cougar demography and immigration in a hunted population. Ecological Applications 18, 1028-1037. Rogers LL (1977) Social relationships, movements, and population dynamics of black bears in northeastern Minnesota, University of Minnesota. Rogers LL (1987) Effects of food-supply and kinship on social-behavior, movements, and population-growth of black bears in northeastern Minnesota. Wildlife Monographs, 1-72. Rousset (2000) Genetic differentiation between individuals. Journal of Evolutionary Biology 13, 58-62. Roy J, Yannic G, Côté SD, Bernatchez L (2012) Negative density-dependent dispersal in the American black bear (Ursus americanus) revealed by noninvasive sampling and genotyping. Ecology and Evolution 2, 525-537. Royama T (1992) Analytical Population Dynamics. Chapman and Hall, London. Rueness EK, Jorde PE, Hellborg L, et al. (2003a) Cryptic population structure in a large, mobile mammalian predator: the Scandinavian lynx. Molecular Ecology 12, 2623-2633. 144 Rueness EK, Stenseth NC, O'Donoghue M, et al. (2003b) Ecological and genetic spatial structuring in the Canadian lynx. Nature 425, 69-72. Ryman N (1997) Minimizing adverse effects of fish culture: understanding the genetics of populations with overlapping generations. ICES J. Mar. Sci. 54, 1149-1159. Saura S, Pascual-Hortal L (2007) A new habitat availability index to integrate connectivity in landscape conservation planning: Comparison with existing indices and application to a case study. Landscape and Urban Planning 83, 91-103. Saura S, Torne J (2009) Conefor Sensinode 2.2: A software package for quantifying the importance of habitat patches for landscape connectivity. Environmental Modelling & Software 24, 135-139. Schwartz CC, Franzmann AW (1992) Dispersal and survival of subadult black bears from the Kenai Peninsula, Alaska. Journal of Wildlife Management 56, 426-431. Schwartz MK, Luikart G, Waples RS (2007) Genetic monitoring as a promising tool for conservation and management. Trends in Ecology & Evolution 22, 25-33. Segelbacher G, Cushman SA, Epperson BK, et al. (2010) Applications of landscape genetics in conservation biology: concepts and challenges. Conservation Genetics 11, 375-385. Shih CC, Huang CC, Li SH, Hwang MH, Lee LL (2009) Ten novel tetranucleotide microsatellite DNA markers from Asiatic black bear, Ursus thibetanus. Conservation Genetics 10, 1845-1847. Short Bull R, Cushman S, Mace R, et al. (2011) Why replication is important in landscape genetics: American black bear in the Rocky Mountains. Molecular Ecology 20, 10921107. Smouse PE, Long JC, Sokal RR (1986) Multiple regression and correlation extensions of the Mantel test of matrix correspondence. Systematic Zoology, 627-632. Solberg KH, Bellemain E, Drageset OM, Taberlet P, Swenson JE (2006) An evaluation of field and non-invasive genetic methods to estimate brown bear (Ursus arctos) population size. Biological Conservation 128, 158-168. Spear SF, Balkenhol N, Fortin MJ, McRae BH, Scribner K (2010) Use of resistance surfaces for landscape genetic studies: considerations for parameterization and analysis. Molecular Ecology 19, 3576-3591. Spear SF, Storfer A (2010) Anthropogenic and natural disturbance lead to differing patterns of gene flow in the Rocky Mountain tailed frog, Ascaphus montanus. Biological Conservation 143, 778-786. Storfer A, Murphy MA, Evans JS, et al. (2007) Putting the 'landscape' in landscape genetics. Heredity 98, 128-142. 145 Storfer A, Murphy MA, Spear SF, Holderegger R, Waits LP (2010) Landscape genetics: where are we now? Molecular Ecology 19, 3496-3514. Taberlet P, Camarra JJ, Griffin S, et al. (1997) Noninvasive genetic tracking of the endangered Pyrenean brown bear population. Molecular Ecology 6, 869-876. Tammeleht E, Remm J, Korsten M, et al. (2010) Genetic structure in large, continuous mammal populations: the example of brown bears in northwestern Eurasia. Molecular Ecology 19, 5359-5370. Thomlinson JR, Bolstad PV, Cohen WB (1999) Coordinating methodologies for scaling landcover classifications from site-specific to global: Steps toward validating global map products. Remote Sensing of Environment 70, 16-28. Thompson LM, van Manen FT, King TL (2005) Geostatistical analysis of allele presence patterns among American black bears in eastern North Carolina. Ursus 16, 59-69. Urban D, Keitt T (2001) Landscape connectivity: A graph-theoretic perspective. Ecology 82, 1205-1218. Urban DL, Minor ES, Treml EA, Schick RS (2009) Graph models of habitat mosaics. Ecology Letters 12, 260-273. Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P (2004) MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes 4, 535-538. van Strien MJ, Keller D, Holderegger R, et al. (2014) Landscape genetics as a tool for conservation planning: predicting the effects of landscape change on gene flow. Ecological Applications 24, 327-339. Vandergast A, Inman R, Barr K, et al. (2013) Evolutionary Hotspots in the Mojave Desert. Diversity 5, 293-319. Vandergast AG, Perry WM, Lugo RV, Hathaway SA (2011) Genetic landscapes GIS Toolbox: tools to map patterns of genetic divergence and diversity. Molecular Ecology Resources 11, 158-161. Wagner AP, Creel S, Kalinowski ST (2006) Estimating relatedness and relationships using microsatellite loci with null alleles. Heredity 97, 336-345. Wandeler P, Hoeck PE, Keller LF (2007) Back to the future: museum specimens in population genetics. Trends in Ecology & Evolution 22, 634-642. Waples RS (1989) A generalized approach for estimating effective population size from temporal changes in allele frequency. Genetics 121, 379-391. Ward MD, Gleditsch KS (2008) Spatial regression models Sage. Wasserman TN, Cushman SA, Littell JS, Shirk AJ, Landguth EL (2013) Population connectivity and genetic diversity of American marten (Martes americana) in the United States 146 northern Rocky Mountains in a climate change context. Conservation Genetics 14, 529541. Wright S (1931) Evolution in Mendelian populations. Genetics 16, 0097-0159. Wright S (1943) Isolation by distance. Genetics 28, 114-138. Wu H, Zhang SN, Wei FW (2010) Twelve novel polymorphic microsatellite loci developed from the Asiatic black bear (Ursus thibetanus). Conservation Genetics 11, 1215-1217. 147