Beyond the snapshot: Landscape genetic analysis of time series data reveal responses of American black bears to landscape change

Abstract Landscape genetic studies typically focus on the evolutionary processes that give rise to spatial patterns that are quantified 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) over time with concurrent changes in landscape pattern. We introduce a novel approach to analyze landscape genetic data through time. We demonstrate this approach using genotyped samples (n = 569) from a large black bear (Ursus americanus) population in Michigan (USA) that were harvested during 3 years (2002, 2006, and 2010). We identified areas that were consistently occupied over this 9‐year period and quantified temporal variation in SGS. Then, we evaluated alternative hypotheses about effects of changes in landscape features (e.g., deforestation or crop conversion) on fine‐scale SGS among years using spatial autoregressive modeling and model selection. Relative measures of landscape change such as magnitude of landscape change (i.e., number of patches changing from suitable to unsuitable states or vice versa), and during later periods, measures of fragmentation (i.e., patch aggregation and cohesion) were associated with change in SGS. Our results stress the importance of conducting time series studies for the conservation and management of wildlife inhabiting rapidly changing landscapes.

and ecological processes that shape SGS is often limited because sampling is typically conducted at a single point in time (Goetze, Andrews, Peijnenburg, Portner, & Norton, 2015). These "snapshots" of data only measure populations in their current state (Anderson et al., 2010;Martensen, Saura, & Fortin, 2017;Schwartz, Luikart, & Waples, 2007). Therefore, inferences can be problematic when studying long-lived and iteroparous species such as black bears (Ursus americanus) that are sensitive to landscape features (i.e., land cover) that are expected to undergo future modification (Cushman, McKelvey, Hayden, & Schwartz, 2006;Cushman, Wasserman, Landguth, & Shirk, 2013). Thus, incorporating a time series approach becomes necessary to make inferences about how a species responds to landscape or environmental change (Sun & Hedgecock, 2017).
Time series data are increasingly being 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 contemporary genetic diversity (Wandeler, Hoeck, & Keller, 2007). Similarly, there has been increasing interest to include genetic monitoring as an important component of management programs (De Barba et al., 2010;Rudnick, Katzner, Bragin, Rhodes, & Dewoody, 2005;Steele et al., 2013). Assessments of changes in genetic diversity provide a means to evaluate trends in connectivity, to infer demographic histories of populations, and to gauge loss of genetic diversity (Schwartz et al., 2007). However, implementing such studies is challenging because sampling the same population(s) at regular intervals is difficult despite the potential value for conservation and management.
There is increased interest in examining how landscape configuration assessed at different times has influenced contemporary SGS Pavlacky, Goldizen, Prentis, Nicholls, & Lowe, 2009;Zellmer & Knowles, 2009). Yet, few if any, landscape genetic studies have applied a time series approach using multiple genetic and landscape data sets collected at the same time points.
Such studies are limited in part by the rarity of multiple samples of the same population from different time points, but are further impeded by the necessity of obtaining complementary time series landscape data. Despite these challenges, time series landscape genetic studies are valuable as natural ecosystems are spatially heterogeneous and landscape composition/configurations can change over time, sometimes drastically so, due to natural and anthropogenic factors (Bolliger et al., 2010;Spear et al., 2010). Furthermore, there may be a discord between the time genetic processes occur (over generations) and the landscape change that reflect contemporary landscape effects (Anderson et al., 2010). Explicitly adding a temporal component in landscape genetics analyses may provide valuable additional resolution on directional trends in gene flow (Martensen et al., 2017;Wagner & Fortin, 2013). By quantifying concurrent changes in SGS and landscape structure over time, we only improve the long-term predictive power of the effects of ongoing or future landscape change on genetic connectivity (Zellmer & Knowles, 2009).
In this study, we propose a novel approach to quantify changes in SGS through time and relate those changes to underlying landscape change for a population of American black bears in the Northern Lower Peninsula (NLP) of Michigan, USA. The NLP black bear population is bounded to the south by expansive areas of unsuitable human-dominated urban and agricultural landscapes. The population therefore experiences little to no emigration and immigration.
The NLP black bears provide a unique opportunity to apply time series analyses in a landscape context. Indeed, we capitalized on harvest samples acquired over a large extent and long-term data spanning a 9-year period. 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 rapid population turnover, and therefore rapid genetic change.
Black bears are affected by forest loss and Michigan's NLP is a heterogeneous landscape that has experienced past and current landscape change (primarily deforestation and agricultural conversion). Hence, we used remotely-sensed land cover data (Coastal Change Analysis Program, CCAP; https://www.coast.noaa.gov/ digitalcoast) that are publicly available every 5 years (2001, 2006, and 2011) for the whole of the NLP during the period over which our genetic samples were collected.
Using a times series approach in a landscape genetic framework, we (i) develop a set of landscape resistance models that incorporate a suite of resistance surfaces representing alternative hypotheses concerning the associations between interindividual genetic distance and landscape resistance distance; (ii) identify and compare the best performing landscape resistance models among years; (iii) evaluate whether local SGS patterns changed over time; and (iv) quantify landscape change and use spatial autoregressive model selection to test for associations between changes in SGS and landscape change.

| Study area
Our study area covers the northern two-thirds of the lower peninsula of Michigan (47,739 km 2 ) ( Figure 1). The NLP is a fragmented mosaic of different land cover and land-use types including development, agriculture, upland nonforested openings, northern hardwood and mixed hardwood, oak, aspen, pine, forested wetland, and nonforested 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 Annually, MDNR requires all harvested bears to be registered at hunter check stations. During registration, a premolar tooth is extracted for aging and DNA extraction (Coy & Garshelis, 1992).
Hunters report the bear's harvest location and sex. Locations of bear harvest samples were recorded as township, range, and section, which were then georeferenced to the section centroid and converted into UTM coordinates.

| Defining temporal sampling locations
To remove potential sampling error that can arise from differences in dispersion and extent of sampling in opportunistically collected annual harvest samples, we identified areas that were consistently occupied in both space and time. First, for each year (2002, 2006, and 2010), we created Voronoi polygons for all samples within a year. Voronoi polygons are created by partitioning the study site with n sampling points into convex polygons, in which a polygon contains only one sampling point and every location in a given polygon is closer to its generating sampling point than any other point (Longley, Goodchild, Maguire, & Rhind, 2010). We then overlaid the three annual Voronoi polygons and used spatial-temporal analysis of moving polygons in ArcGIS 10.1 to define areas where the three polygon layers overlapped and contained at least one bear from each year ( Figure S1). We defined these as "areas of consistent sampling" ( Figure S2) over the entire (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) sampling period. Voronoi polygons are used as a way to ensure that throughout the years that subareas (here Voronoi polygons) contain enough samples so that comparisons from 1 year to another is possible. Hence, the Voronoi polygons are not used to estimated/mapped density but rather to delineate subareas such that changes are attributed to more or less the local conditions that the bears were exposed to. Samples within the overlap area of consistent sampling (2002, n = 204; 2006, n = 199; 2010, n = 166) were used in all downstream population and landscape genetic analyses.

| Microsatellite genotyping
We extracted DNA from bear teeth using Qiagen DNEasy Tissue Kits (Qiagen Inc., Valencia, CA, USA) following manufacturer protocols. We quantified DNA using a Nanodrop spectrophotometer (Thermo Scientific, Waltham, MA, USA) and diluted samples to a 20 ng/μl working concentration. We used PCR to amplify 12 vari- We amplified DNA according to conditions outlined in Moore, Draheim, Etter, Winterstein, and Scribner (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, USA). All alleles were scored independently by two experienced laboratory personnel using SAGA genotyping software (LI-COR Inc.). To assess genotyping error, 10% of samples were randomly selected and genotyped twice to yield a genotyping error rate of <2%.
We used program MICRO-CHECKER (Van Oosterhout, Hutchinson, Wills, & Shipley, 2004) to test for the presence of null alleles and allelic dropout. We tested for deviations from Hardy-Weinberg equilibrium based on exact tests χ 2 using program GENEPOP (Version 3.1d; Raymond & Rousset, 1995) and used sequential Bonferroni tests to correct for multiple comparisons (Rice, 1989). We used Bonferroni corrections (Goudet, 1995(Goudet, , 2001 to test for linkage disequilibrium (evaluated by permutation tests, based on 660,000 randomizations) in program FSTAT 2.93.

| Landscape genetic analysis: isolation by distance versus isolation by resistance
We characterized land cover for the three sampling periods using NOAA CCAP Land Cover digital coverage maps, derived from Landsat TM imagery for 2001, 2006, and 2011. To test for isolation by resistance (IBR; to determine 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.
We weighted each land cover class according to positive or negative associations with black bear presence based on habitat suitability indices (HSI) developed independently from telemetry data by Carter, Brown, Etter, and Visser (2010). Landscape resistance 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). Models included resistance surfaces using a land cover classification from CCAP Land Cover datasets (resolution = 150 m, 25 classes). We reclassified datasets into seven land cover classifications according to bear habitat suitability (HSI; cost = 1-100, least to most resistant to bear movement): mixed deciduous forest (MF = 1), forested wetland (FW = 1), evergreen forest (EF = 10), nonforested upland (NFU = 20), agriculture (AG = 50), nonforested 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), because they represent putative physical barriers to dispersal. In addition, we created a 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 evaluated 11 landscape hypotheses of IBR (Table 1).
Using the sum function of raster algebra in Spatial Analyst ArcGIS We used Mantel (Mantel, 1967) and partial Mantel tests (Smouse, Long, & Sokal, 1986) to correlate genetic and geographic/resistance distances. Legendre and Fortin (2010)  and of those neighboring a sampling grid cell, respectively. Also, because deforestation is the predominant type of landscape change, we also described landscape change in terms of proportion of area within a sampling grid cell with forest loss (FL).

| Defining genetic change
Individual-based landscape genetic analyses are often limited to link-based methods. These type of analyses relate pairwise genetic distance between individuals to their landscape distance in which hypotheses are explicitly defined in terms of distances. However, to move beyond link-based methods we must transform individual-based pairwise genetic distances into a single Y vector. Landscape genetic neighborhood-based approaches can use information from pairwise links (genetic distances) to create a node-based data structure to relate GD patterns to local landscape predictors (James, Coltman, Murray, Hamelin, & Sperling, 2011;Wagner & Fortin, 2013). Here we introduce a novel analytical framework to test whether SGS changed over time.
We generated our response variable by producing spatial genetic data layers using the Genetic Landscapes Toolbox (

| Spatial modeling analyses
We compared multiple spatial autoregressive lag models to determine whether landscape change was associated with genetic change. These models are linear regressions that account for spatial nonindependence in the response variable (genetic change) layer with an additional spatial lag term as an explanatory variable and parameter estimation assessed by Maximum Likelihood Estimation (Ward & Gleditsch, 2008). Models without significant regression coefficients (p < .05) were discarded. Because of the inherent relationships between landscape change magnitude metrics NP, PLAND, DEG, FL, and DMFL (Table 2), we performed five univariate analyses with each of these explanatory variables and compared them 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 change magnitude metric of the single regression models was then used in multiple regression analyses.
All tests were performed in the spdep package in R v. 3.0.1 (R Core Team, 2013).

| Population genetic structure
After accounting for spatial differences in sampling among years,  Figure S3).

| Isolation by resistance testing
Landscape genetic analyses revealed discrepancies in landscape resistance model performance among years. Statistically significant landscape resistance models (p < .05) were only identified in 2006 and 2010 after partialling out the effects of geographic distance (

| Genetic change and landscape change
The cumulative genetic change surfaces show areas of increased and decreased GD over the study area ( Figure 4) Figure 4). Therefore, NP was used in multiple regression analyses. Multivariate models with significant (p < .05) coefficients are presented in Table 3. All temporal comparisons indicated that genetic change was significantly influenced by habitat alteration calculated as the number of landscape change patches (NP ; Table 3; Figure 4). However, for all comparisons that included 2010, genetic change was better predicted by inclusion of landscape change heterogeneity (Table 3) Table 3).

| D ISCUSS I ON
The field of landscape genetics can advance beyond characterizing static landscape pattern-genetic process relationships . Incorporating landscape change and quantifying the effects of change on SGS over time is vital to improve understanding of anthropogenic impacts on functional connectivity (Bolliger et al., 2010;Segelbacher et al., 2010;Storfer et al., 2010). Our study illustrates the advantages of joint use of time series genetic and landscape data when using established TA B L E 2 Most supported univariate spatial regression models to predict black bear genetic change using landscape variables that characterize the relative magnitude of landscape change in the Northern Lower Peninsula, Michigan, USA. Abbreviations are as described in Table S1 landscape genetics methods. Using a long-term genetic data set and land cover imagery covering an analogous time period, we found disparities in landscape resistance models during different "snapshots" of time (Table 1). Our results show that if landscape does not appear to be an important component shaping functional connectivity at one time point, it does not mean that attributes will not become important over relatively short time scales in changing landscapes. Our findings have broad implications for conservation and management, and speak to the importance of collecting long-term data.
Black bears in Michigan's NLP exhibit an IBD pattern of gene flow, commonly documented in black bears (Costello, Creel, Kalinowski, Vu, & Quigley, 2008;Moore et al., 2014;Rogers, 1977Rogers, , 1987Schwartz & Franzmann, 1992). Indeed, strong IBD is consistent with studies on bears and other wide-ranging carnivores (Brown, Hull, Updike, Fain, & Ernest, 2009;Paetkau, Waits, Clarkson, Craighead, & Strobeck, 1997;Rueness, Jorde, et al., 2003;Rueness, Stenseth, et al., 2003). However, we found isolation by landscape resistance was more strongly supported than IBD in 2006 and 2010. Land cover was the landscape factor most strongly correlated with interindividual genetic distance in black bears in the NLP after partialling out geographic distance (Table 1). In contrast, roads and rivers alone did not appear to influence genetic distance. However, the model that only included rivers was significantly correlated with GD in 2010. Our results are consistent with a previous study that performed similar landscape genetic analyses on NLP black bears sampled during 2006 but used a much larger sample size (n = 365) and spatial distribution (Draheim, 2015). 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, Weir, Hamilton, & Deal, 2006;Noyce & Garshelis, 2011;Rogers, 1977). Carter et al. (2010) found a negative association between bear location and small and medium roads in the NLP black bear population based on radiotelemetry locations from 1991 to 2000 for 20 males and 35 females. However, our results suggest roads do not unduly influence functional connectivity. Indeed, while previous landscape genetic studies have found roads as a limiting factor to gene flow in black bears, the relative influence of roads varies among populations (Cushman et al., 2006;Short Bull et al., 2011), and in some cases, roads may serve to facilitate rather than impede bear movement ).
The cumulative effects of landscape alteration over increasing lengths of time explain the differences in relative effects of landscape features on SGS in NLP black bears. We found significant effects of landscape change on changes in levels of GD (Tables 2   and 3). However, the amount of variation explained in models that included variables associated with habitat heterogeneity were the best predictors of genetic change for temporal comparisons that included 2010 (Table 3). This discrepancy among years reflects an increasing amount of habitat fragmentation due to deforestation over the 9-year sampling period. Currently, the NLP is a patchy landscape with areas of deciduous, mixed, and coniferous forest frag-  Figure 2). On the whole, F I G U R E 4 Distribution of (a) the cumulative difference in genetic differentiation (from spatial interpolated genetic surfaces) within a grid cell for all temporal comparisons (2002-2006, 2006-2010, 2002-2010) and ( TA B L E 3 Most supported multiple spatial autoregression models to predict black bear genetic change using landscape metrics that characterize the relative magnitude or configuration of landscape change in the Northern Lower Peninsula, Michigan, USA. Abbreviations are as described in Table S1 Model ( 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 in dispersed, complex habitat patches.
Lack of an association between land cover and SGS in bears during 2002 (Table 1) (Anderson et al., 2010;Epps & Keyghobadi, 2015). However, when dispersal rates and distances are large, as exhibited in the NLP black bear population (Draheim, 2015;Draheim et al., 2016;Moore et al., 2014), shorter or no time lags are expected (Epps & Keyghobadi, 2015). Also, legacy effects of historical landscape processes may be reduced using genetic markers with higher mutations rates (i.e., microsatellites) that reach mutation-drift equilibrium quickly and genetic measures that respond rapidly to changes in connectivity (e.g., Dps) (Anderson et al., 2010). A number of simulation studies have found landscape effects on SGS could be detected with relatively short time spans Landguth et al., 2010;Murphy, Evans, Cushman, & Storfer, 2008).
Field sampling of wildlife populations can never be conducted exhaustively nor can spatial coverage be completely replicated during every time period. Thus, stochastic sampling processes, especially in populations where density varies across the landscape, can never really be removed from empirical data sets. Harvest data are opportunistically collected, meaning locations will inevitably vary among temporal periods and may be predisposed to erroneous inferences across years (Schwartz & McKelvey, 2009). In our study system, bear harvest samples do show consistent regional density and distribution patterns among years (Draheim, Lopez, Winterstein, Etter, & Scribner, 2015). We additionally controlled for spatial sampling heterogeneity among years by only including areas where bears were consistently harvested over time. We acknowledge that reducing our sample size in this way may reduce statistical power and add sampling noise, but we weighed the benefit of increased sample size against increased spatial sampling bias. Regardless, our samples sizes after accounting for spatial heterogeneity are respectable (n = 166-204). Also, Graves et al. (2012) showed that samples located at different locations within a home range did not affect the genetic-landscape relationship, we used a comparable approach and assumption here in this study.
There can be a disconnect between what is statistical vs. biologically significant. Thus, we tried to account for this in our study.
First, we based landscape resistance weights on field data from radiotelemetry. Second, we performed a sensitivity analyses on our resistance weights to ensure landscape genetic results were not unduly influenced by the scale of resistance values. Lastly, due to the inherent risk of false positives using Mantel and partial Mantel tests, we implemented an analytical framework proposed by Cushman et al. (2013) that is based on the relative support of each candidate model to separate the true resistance model from a range of erroneous alternative resistance models.
Our finding of genetic change between periods within localized areas over a 9-year sampling period ( Figure 4)  We are unsure as to how genetic drift, which is likely occurring, contributed to our results. If genetic drift was the dominant force contributing to temporal variation in SGS of NLP black bears, we would predict little evidence for IBD/resistance or associations between genetic and landscape change, which is inconsistent with our results. Furthermore, in an graph theory framework, Draheim et al. (2016) found black bears in the NLP exhibit asymmetric gene flow based on areas of high and low net flux indicating source-sink dynamics (Pulliam, 1988). Therefore, if gene flow is mediated by use of corridors among source and sink areas and corridors are disrupted due to landscape change, local SGS is expected to change. Further investigations, for example, using local harvest abundance or other proxies for bear density (e.g., Moore et al., 2014), are needed to address this question.
While our time series approach may be broadly applicable across a range of taxa, the approach needs to be parameterized based on specific life history and landscape characteristics of the species or population of 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 genetic/landscape change 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 SGS were not due to immigration or emigration. Also, it is important that the rate of landscape change be sufficient to influence SGS over relatively short time intervals.
To 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 enable managers to target regions or habitat types that are important for maintaining connectivity across anthropogenically altered habitats and assess the impacts of future landscape change.

ACK N OWLED G EM ENTS
We are grateful to J. Fierke, K. Filcek, B. Dreher, Veronica Lopez, and S. Libants for assisting in data generation. We also thank J.
Messina, B. Epperson, K. Holekamp, and J. Kanefsky, 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.
Financial and logistical support for this project was provided by the

DATA ACCE SS I B I LIT Y
Data used in this study, including bear microsatellite genotypes and harvest locations, are available for download on Dryad (https://doi. org/10.5061/dryad.c61q0).

AUTH O R CO NTR I B UTI O N S
All authors contributed to the theoretical foundation and writing of the manuscript. HD and JM generated data and HD performed analyses.

CO N FLI C T O F I NTE R E S T
None declared.