Climate connectivity of the bobcat in the Great Lakes region

Abstract The Great Lakes and the St. Lawrence River are imposing barriers for wildlife, and the additive effect of urban and agricultural development that dominates the lower Great Lakes region likely further reduces functional connectivity for many terrestrial species. As the climate warms, species will need to track climate across these barriers. It is important therefore to investigate land cover and bioclimatic hypotheses that may explain the northward expansion of species through the Great Lakes. We investigated the functional connectivity of a vagile generalist, the bobcat, as a representative generalist forest species common to the region. We genotyped tissue samples collected across the region at 14 microsatellite loci and compared different landscape hypotheses that might explain the observed gene flow or functional connectivity. We found that the Great Lakes and the additive influence of forest stands with either low or high canopy cover and deep lake‐effect snow have disrupted gene flow, whereas intermediate forest cover has facilitated gene flow. Functional connectivity in southern Ontario is relatively low and was limited in part by the low amount of forest cover. Pathways across the Great Lakes were through the Niagara region and through the Lower Peninsula of Michigan over the Straits of Mackinac and the St. Marys River. These pathways are important routes for bobcat range expansion north of the Great Lakes and are also likely pathways that many other mobile habitat generalists must navigate to track the changing climate. The extent to which species can navigate these routes will be important for determining the future biodiversity of areas north of the Great Lakes.

& Joseph, 2009), but other species that face geographic barriers and cannot keep pace with the velocity of climate change might risk extinction (Thomas et al., 2004).
Many highly mobile species might indeed be able to track contemporary climate change across natural landforms, but the addition of cities, highways, roads, and agricultural crops can potentially hinder mobility in an additive fashion (Epps et al., 2005;Riley et al., 2006;Robillard et al., 2015). Furthermore, suitable habitat might be found several hundreds of kilometers north of a species' current range, but the environmental characteristics of the interstitial landscape could be well outside of its niche (Early & Sax, 2011). In such a case, the environment might impede or block dispersal and the colonization of newly available habitat (McRae, 2006;Wang & Bradburd, 2014). Ultimately, future biodiversity across the globe will depend on the ability of species to rapidly disperse throughout continental-scale habitat networks to keep up with changing environmental conditions. Many animals will need to migrate across human-dominated and highly modified landscapes to colonize new habitats. It is therefore necessary to understand the effect that natural and anthropogenic barriers have on a species dispersal ability and to understand how these barriers influence connectivity across entire regions.
The Great Lakes are visibly the largest natural barrier to terrestrial species migration in eastern North America. Currently, we do not know for many species to what extent the Great Lakes have influenced movement and consequently gene flow, although the impacts are likely profound; much smaller barriers such as canals, highways, mountains, rivers, roads, sea lochs, and urban development have been shown to restrict the gene flow of many terrestrial vagile species (Blanchong et al., 2008;Coulon et al., 2006;Cushman & Lewis, 2010;Epps et al., 2005;Koen et al., 2015;Kuehn et al., 2007;Pérez-Espona et al., 2008;Proctor, McLellan, Strobeck, & Barclay, 2005;Riley et al., 2006;Robinson, Samuel, Lopez, & Shelton, 2012;Robinson, Samuel, Rolley, & Shelton, 2013;Vander Wal, Paquet, & Andraés, 2012). Unfortunately, the additive influence of anthropogenic disturbance between and within the vicinity of the Great Lakes will likely further restrict gene flow through these large natural barriers for many species.
The bobcat (Lynx rufus) is the most widely distributed feline species in North America, and there is evidence that it was more abundant across the continent before European colonization and during the Pleistocene (Deems & Pursley, 1983;Graham & Lundelius, 2010;Lariviere & Walton, 1997). It is generally thought that intensive trapping and land clearing led to the extirpation of the species in the midwestern United States and many parts of the Great Lakes region and this may also have caused the apparent absence of the species in many areas of the Midwest United States (rode Deems & Pursley, 1978;Deems & Pursley, 1983;Vos, 1964;Woolf, Nielsen, & Gibbs-Kieninger, 2000).
In recent decades, bobcat sightings, road kills, and individuals incidentally harvested by trappers have become more common in the region (Marrotte, Bowman, & Morin In Review;Roberts & Crimmins, 2010;Woolf & Hubert, 1998). There is evidence that bobcat populations are spreading into areas where they were thought to be absent (Linde, Roberts, Gosselink, & Clark, 2012;Woolf & Hubert, 1998 Bobcats are not spreading to the same extent into southern Ontario however, even though they once inhabited this landscape (de Vos, 1964). Landscape configuration could be playing an important role in structuring the recolonization of the bobcat in the Great Lakes region. For instance, the Great Lakes and the St. Lawrence River are imposing barriers to movement for wildlife. In addition, urban development and agricultural development dominate southern Ontario and may be impeding bobcats from colonizing this range frontier, over and above the barrier effect of the Great Lakes. There is evidence, however, that the bobcat can cope with an anthropogenic environment (Lee et al., 2012;Riley et al., 2003;Tigas, Vuren, & Sauvajot, 2002;Woolf et al., 2000). For example, in Illinois, the bobcat occupies landscapes with intensive agriculture (Woolf et al., 2000). It also seems capable of occupying areas surrounded by transportation infrastructure and urban development (Lee et al., 2012;Riley et al., 2003;Tigas et al., 2002). However, urban land cover and major highways have caused reduced gene flow in bobcat populations in California (Kozakiewicz et al., 2019;Lee et al., 2012). Snow is also considered by some to be a limiting factor to bobcat expansion north of its range, as many researchers have suggested that the species has high foot loading and cannot efficiently travel and hunt in deep snow (Hoving, Joseph, & Krohn, 2003;Marston, 1942;McCord, 1974;Parker, Maxwell, Morton, & Smith, 1983). For example, McCord (1974) found that the bobcat had a difficult time traveling through areas that had a sinking depth exceeding 15 cm. Also, Parker et al. (1983) suggested that the reason the bobcat did not invade the highlands of Cape Breton was because of the deeper snow. In addition, snow clearing and compaction near human settlements may mediate the influence of snow on colonization and may promote bobcats from occupying areas north of their range (Marrotte et al., In Review).
We investigated several land cover and bioclimatic hypotheses that may explain the northward expansion of the bobcat throughout the Great Lakes region in North America, because any restrictions imposed on a highly mobile species would likely be even more perilous for less vagile species. We hypothesized that barriers to northward expansion would hinder gene flow. We considered that there are 3 scenarios that may describe range expansion in the Great Lakes region: • H 0 : Panmixia: Natural and anthropogenic barriers have no effect on gene flow; thus, individuals are panmictic.
• H 1 : Isolation by distance (IBD): Natural and anthropogenic barriers have no effect on gene flow, but gene flow decays over geographic distance.
• H 2 : Isolation by resistance (IBR): Natural and anthropogenic barriers constrict gene flow; thus, flow percolates through land bridges between the lakes.
The bobcat is an ideal study species to test our hypotheses of range expansion in the context of anthropogenic change, because it is a vagile habitat generalist that is currently expanding its range and demonstrates some limitations to human disturbance and climate.
We predicted that gene flow of the bobcat is obstructed naturally by the Great Lakes and deep snow but also hindered by low forest cover and by the transportation infrastructure. This model most closely follows the isolation-by-resistance hypothesis (H 2 ) previously described.
Therefore, we predicted that gene flow is constricted through certain pathways that connect individuals throughout the region ( Figure 1). We predicted that gene flow in southern Ontario originated mostly from the east from the province of Quebec and New York State, since flow is limited through the Lower Peninsula of Michigan (LPM) and between Lake Ontario and Erie, because of the high road density and low forest cover of these regions. On the other hand, northern Ontario is connected to the south by a more natural landscape with high forest cover and less human disturbance. Consequently, gene flow to northern Ontario is facilitated by the Upper Peninsula of Michigan (UPM) and the largely forested area to the west of Lake Superior. Our rationale is that range expansion in this region is restricted by the additive effect of natural and anthropogenic barriers. Gene flow should be constricted and forced to pass through land in between and around the Great Lakes, while deep snow should reduce the capability of flow northwards and cause gene flow to deviate around areas that receive high annual snowfall caused by the lake effect (Norton & Bolsenga, 1993).
The upper Great Lakes are periodically hit by frequent lake-effect snowfall or snow squalls with over 15 cm of snow accumulation in a single day (Baijnath-Rodino & Duguay, 2018). In addition, gene flow should be hindered by agricultural areas with low cover such as the corn belt areas of the Midwest and areas with high density of roads such as urban areas. Resources and Forestry, researchers, and trappers ( Figure 2). We sampled bobcat pelts found within an area around the Great Lakes defined by the maximum dispersal distance of the bobcat of 300 km (Johnson, Walker, & Hudson, 2010). We sampled on both sides of the international border between Canada and the United States.
We then explored the spatial structure of these data using a spatial principal component analysis (sPCA ;Jombart, Devillard, Dufour, & Pontier, 2008) and tested for patterns of spatial autocorrelation.
We used a distance-based nearest neighbor approach, where individuals within 300 km were assumed to be neighbors, since bobcat have been observed to disperse up to 288 km from their natal range (Johnson et al., 2010;Knick & Bailey, 2006). We also tested different neighborhood approaches and found similar spatial patterns.
We used a Monte Carlo test to determine whether there were any global spatial structures (broadscale clusters or clines) or local spatial structures (fine-scale disparity among neighbors) worth investigating (Jombart, 2008). We permuted the alleles scores 9,999 times to test the significance of the spatial structure. If there was no spatial structure, then the panmixia hypothesis (H 0 ) would be concluded, because our isolation-by-distance (H 1 ) and isolation-by-resistance (H 2 ) hypotheses were inherently spatial. After exploring the spatial structure, we then used the proportion of shared alleles as a metric of genetic similarity between individuals and tested our spatial hypotheses.

| Gene flow covariates
We built 4 different landscape maps, which we thought could explain bobcat gene flow in the Great Lakes region and which would allow us to test our isolation-by-resistance hypothesis ( Figure 3). To establish the spatial extent of our analysis, we used the minimum convex hull that contained the Great Lakes with a 400 km buffer to leave 100 km between the edge of the map and any bobcat samples (Koen, Garroway, Wilson, & Bowman, 2010). The Great Lakes spatial layer we used to create this boundary and the St. Lawrence river layer that we later used were gathered from the Lakes and Rivers, 2009, spatial layers freely available by the Commission for Environmental Cooperation (CEC; cec.org).
To create our Great Lakes landscape layer, we assigned values of 1 to areas where there was either the Great Lakes or the St.
Lawrence and the value 2 to the land ( Figure 3a).
We used the forest cover layer provided by the University of Maryland, which has a continuous forest cover field where each pixel has an assigned value that represented the percentage of tree cover as the data source for our second map (glcf.umd.edu/data/ treecover; DeFries, Hansen, Townshend, Janetos, & Loveland, 2000). This forest cover layer extended across North America at a resolution of 1 km 2 with values that ranged from 10% to 80% forest cover. However, there were values of 254 and 255, which represented areas that were nonvegetated and areas where tree cover was less than 10%. We assigned a value of 1 to areas that were nonvegetated and a value of 5% to areas that were less than 10% vegetated (Figure 3b). We were not able to find tree cover at such a fine resolution that matched the temporal resolution of our bobcat data, but there has not been much recent forest loss in our study area; from 2000 to 2017, the largest forest loss was in Minnesota, which had a decrease in forest cover by 6.9% (globalforestwatch.org).
We used the freely available road layer provided by Natural Earth (naturalearthdata.com) to create our next landscape. We calculated the road density within a radius of 1 km from the center of a pixel with a resolution of 1 km ( Figure 3c).
For our snow layer, we gathered annual data from the Global Historical Climatology Network (GHCN; ncdc.noaa.gov/snowand-ice), which provides data tables for climate stations found across the world. We chose stations found within Canada and the United States that had more than 15 years of data between 1980 and 2011 and calculated the annual snowfall mean of each climate station. We then interpolated these data using ordinary spherical kriging ( Figure 3d).
All surfaces were aggregated using the mean for continuous surfaces and the mode for discreet surfaces to a resolution of 2 km to reduce computation time. We masked out the ocean using a landform layer also provided by Natural Earth. All our layers were projected to North America Lambert Conformal Conic (https ://epsg. io/102009).

| Statistical framework
Our hypotheses were that bobcat gene flow across the Great Lakes region is panmictic (H 0 ), a function of distance (H 1 ), or a function of resistance (H 2 ). Therefore, we first investigated 2 null models that simply tested for panmixia and isolation by distance, and we then tested 8 isolation-by-resistance models that combined landscape features ( Table 1). The Great Lakes barrier was present in all 8 landscape models. We did not find it logical to test isolation-by-resistance models in the absence of the Great Lakes, because forest cover, road density, and snow cover are additive effects on gene flow and not solitary effects. We reasoned that bobcats cannot inhabit a lake but could occupy an area where there are no forest cover, high road density, and deep snow. Also, these isolation-by-resistance models included the influence of isolation by distance, because of the nature of resistance distance (McRae, 2006).
For each hypothesis, we fit the proportion of shared alleles using generalized linear mixed models with a normal error structure and a "log" link, since we reasoned that the influence of landscape features on gene flow decays exponentially across space. We also found higher variance explained using an exponential model than a linear model. We used a maximum-likelihood population-effects covariance structure to account for the nonindependence of the pairwise nature of the data (Clarke, Rothery, & Raybould, 2002). For our isolation-by -resistance models, we were interested in estimating the resistance of the landscape; consequently, we used circuit theory and landscape resistance optimization (Marrotte, Gonzalez, & Millien, 2014;McRae, 2006;Peterman, 2018). We optimized the resistance of each landscape map with the functions provided in the "ResistanceGA" package (Peterman, 2018)  We also accounted for uneven sampling intensity, because it could lead to sampling artifacts (Kierepka & Latch, 2016) and eventually spurious conclusions if not overlooked (Balkenhol & Fortin, 2015). Consequently, we resampled the individuals with replacement into 999 sets of individuals that were at least 100 km apart across the study area. We chose a minimum distance of 100 km, because it gave us the ability to homogenize sampling intensity across the study area and gave us a sufficient number of individuals to investigate. In contrast, a larger distance would have left us with less than 30 individuals and a smaller distance would have left us with quite variable sampling intensity across the study area ( Figure 2). In addition, resampling gave us the ability to replicate our models on 999 different combinations of data and therefore gave us the ability to measure the consistency of our results. Given the computation time required to optimize the resistance surfaces and the large number of replicates (7,992), we fit our models using several computer clusters (Cedar, Graham and Orca; computecanada.ca).
We finally ranked each model within its set of replicates with AICc. We also calculated the overall average rank, AICc, ΔAICc, and AICcW t to compare each model. We then used a branch-and-bound algorithm in R to find the consensus median ranking of all 10 models using the "ConsRank" package (D'Ambrosio, Amodio, & Mazzeo,

| Functional connectivity
By Ohm's law, circuit resistance is reciprocal to current and fundamental work by McRae (2006) demonstrated that current density is proportional to gene flow (McRae, 2006). To help answer our research question concerning how bobcat populations are connected through the Great Lakes region, we created an omnidirectional current density map. We gathered all 999 optimized resistance surfaces of the top landscape model, and for each surface, we standardized the optimized resistance to the mean. We then calculated the average resistance of each pixel to produce a map of standard average resistance (resistance from this point forward). We then used circuit theory in Circuitscape version 4.05 to produce a current density map (Shah & McRae, 2008). We generally followed the methods of Koen, Bowman, Sadowski, and Walpole (2014) to produce an omnidirectional current density map. However, resistance values cannot be negative; therefore, we first scaled the values between 1 and 100.
We then regularly placed 100 nodes at the periphery of the map and simulated current passing between all pairs and summed the total current passing through the Great Lakes region. In conjunction with the optimized average standard resistance surface, this current map gave us an idea where current or gene flow was being impeded.

F I G U R E 3
Landscape maps used to test isolation-by-resistance hypotheses of bobcat gene flow across the Great Lakes region. (a) Great Lake land barrier, (b) forest cover, (c) road density, and (d) annual snowfall. In total, 8 isolation-by-resistance models were tested and included only combinations of the Great Lakes with all other three maps. We also compared these models to a null model of panmixia (H 0 ) and isolation by distance (H 1 )

| RE SULTS
We used 240 samples after removing those with missing alleles that were not accurately georeferenced or that were outside our deline- After performing an sPCA on the allele scores, we found 1 significant global pattern (Observation: 0.020, p-value: 0.001, Figure 2).
There was a NE to NW pattern in bobcat alleles scores across the Great Lakes region, where individuals in both northern corners of the Great Lakes region were found at opposite ends of this gradient.
The proportion of shared alleles ranged from 0.07 to 0.63. All 999 sets of samples had on average 37 individuals and ranged from 31 to 43 individuals. Due to the way we bootstrapped our samples, some samples were selected more often than others. On average, individuals were sampled 154 times and this ranged from 2 to 999 times. There were only 3 samples (from Indiana, Illinois, and Ontario) that were sampled in every set. The 3 individuals were sampled each time, because they were isolated in an area more than 100 km from the nearest other individuals. We checked whether the 3 samples might have driven the optimization of the resistance surface, and we found that the resistance values within 100 km of each of these sites were near the average range compared with other areas on the map (0.17, 0.38 and 0.23). Only 1 individual was sampled as little as twice; this sample was in the highly sampled area in western Minnesota. The composite landscape model that generally ranked first with AICc within its set using the consensus branch-and-bound algorithm included the Great Lakes, forest cover, and annual snowfall (Table 1).
In addition, this composite landscape model had the lowest average rank, AICc and ΔAICc. The composite model ranked first 19.2%, second 22.1%, and third 16.5% of the time, and had the largest proportion of its replicates in the top 3 ranks, but the Great Lakes and Forest cover model did have more replicates that ranked first (22.2%; Table 1 Summary statistics from the resistance optimization algorithm allowed us to determine how much each covariate contributed to the optimized resistance surface (Peterman, 2018). We found that the contribution of the forest cover was the highest (µ = 58.5%), followed by snow (µ = 37.8%) and then the Great Lakes barrier (µ = 3.7%). The mean slope of this model was −0.079, with only 6 of 999 iterations having a positive slope. Thus, in most cases, the effective resistance was positively correlated with genetic distance.
After scaling all replicates of the optimized surfaces of the best combined landscape model and taking the average through each pixel, we found that the Great Lakes had a higher resistance than other features on the landscape (Figures 4,5). In fact, the re-  Note: Consensus rank was determined using the AICc between the 10 models within each 999 set of replicates. Values in bold font are the best value of each metric. All models except the panmixia (H 0 ) and isolation by distance (H 1 ) were isolation-by-resistance models (H 2 ).

*Great Lakes.
There were generally spatial patterns of resistance and current density over our study area. For instance, resistance was high in the lower Great Lake states, but there was a zone of low resistance that overlapped Pennsylvania and Ohio ( Figure 5). Conversely, resistance was lower in the upper Great Lake States and farther south. Though resistance seemed high in New York State, there was an "L"-shaped corridor of low resistance and high current that followed the border between Vermont and New York State from the border of Quebec and turned west from the tristate boundary and continued all the way to the Canada and US border between Lake Ontario and Lake Erie. This area with low resistance connected QC, VT, and NY State with a square corridor around the Adirondack (Figure 6). This corridor also connected to southern Ontario and jumped the St. Lawrence

River from the Thousand Islands Archipelago between Canada and
the United States into Ontario.
On the Lower Peninsula of Michigan (LPM), resistance was low and current was high compared with the Upper Peninsula. The high resistance and low current area of the UPM also overlapped to some extent northern Wisconsin and northeastern Minnesota. In southern Ontario, Canada, resistance was high and decreased into central Ontario but increased once again toward northeastern Ontario.
The current mimicked this pattern but was also amplified on the Niagara Escarpment. On the shores of Lake, Huron resistance was low and current high, and this was also true for Manitoulin Island.
Comparatively, northwestern Ontario had low resistance and consequently high current. Other important corridors were the high current areas that followed Mississippi River into Manitoba and northern Ontario and the pinch point that connected Michigan to northern Ontario from the Straits of Mackinac.

| D ISCUSS I ON
We originally hypothesized that gene flow percolated between the Great Lakes and deep snow areas but was also hindered by low forest cover and by the transportation infrastructure. We found significant spatial structure where gene flow was constricted by the Great Lakes and areas with low and high forest cover with deep lake-effect snow impeded gene flow, while intermediate forest cover facilitated gene flow in the Great Lakes region. Although we did not anticipate the bimodal effect of forest cover, our findings were consistent with our isolation-by-resistance hypothesis (H 2 ).
The warming climate will only aid in the expansion of vagile habitat generalists, since areas of deep lake-effect snow may eventually disappear. However, the Great Lakes and areas with low forest cover limited gene flow more than snow; therefore, we can only predict that less mobile and generalist species will have a more difficult time spreading northward across this landscape as they track climate. In

F I G U R E 5
The average standard resistance from 999 replicates of the top model that was fit using resistance surface optimization of a landscape model that included the additive effect of the Great Lake, forest cover, and annual snowfall. These models were fit to the genetic similarity of bobcat samples across the study area. Labels are as follows: F I G U R E 6 Current density for gene flow through the Great Lakes region. Current density was estimated from the pairwise current of 100 nodes placed on the extremity of all sides of the study area. We used the average standard resistance surface of the top model and rescaled the values from 1 to 100 and used Circuitscape to calculate the cumulative current density of the pairwise iteration of the 100 nodes. We than natural log-transformed, standardized, and scaled the current density to the mean. A value of 0 indicates areas that have average log-transformed current density, and value below and above indicate below and above average log-transformed current density We also predicted that certain pathways connected populations throughout the region (Figure 1). We wrongly expected the Upper Peninsula of Michigan to be a main pathway to northern Ontario but found instead that the Lower Peninsula was more probable.
This would mean that species would be forced to cross the Straits of Mackinac and the St. Marys River. The route from the Lower Peninsula is likely difficult for many species, since it requires the crossing of a 5-km stretch of water, or ice in the winter months. The latter will become less common as the climate warms. It is also likely that the UPM will be favored as climate warms, since deep lake-effect snow will become less likely and this route also does not require crossing a 5-km stretch of water. However, both the LPM and UPM require the crossing of the St. Marys River. We also incorrectly expected that the Niagara region between Lake Ontario and Erie would be an unimportant route. However, the Niagara region might only be traversable by mobile species such as the bobcat that are resilient to anthropogenic disturbances, since the area has a high density of roads compared with other pathways through the Great Lakes ( Figure 3c). For example, gene flow of the highly adaptable raccoon (Procyon lotor) was restricted through this route between Canada and the United States and this also matched the pattern of raccoon rabies incidences at the time (Cullingham, Kyle, Pond, Rees, & White, 2009).

| Natural barriers
We found that the Great Lakes, on average, had high resistance compared with any other feature on the landscape (Figures 4 and   5). Although at first glance, we did find that the Great Lake barrier itself did not contribute much to the optimized resistance values of the top models, but this was due to the snow layer that created a spatial trend within the Great Lakes ( Figure 5). The pattern within the lakes was caused by the large quantity of annual snow received due to the lake effect ( Figure 3d). The variability in snow found within the Great Lakes seemed to be important, since models that included snow ranked better than the Great Lakes and forest cover models (Table 1). For example, in Keweenaw Bay in Lake Superior (KB; Figure 5), resistance was higher than the rest of the lake and this was due to the high amount of snow that the bay received annually due to the lake effect (Figure 3d).
In all, even if the Great Lakes and annual snow were confounded, the outcome was the same, the Great Lakes were without a doubt a barrier to gene flow whether it was caused by the lakes themselves as a barrier or the deep snow that accumulated on them in winter when they freeze or both. On land, snow alone did not seem to be quite important overall, but there were a few areas where the lake-effect snow was quite important, this was the UPM, the Bruce Peninsula, and the area to the east of Lake Ontario that intersected some parts of the Adirondacks (Figure 3d). However, only the UPM and the east side of Lake Ontario had higher resistance ( Figure 5).
Our prediction that snow restricted gene flow northward over our study area did not hold, this could have been due to the low number of samples in the more northern areas of the bobcat's distribution where snowfall is much higher (Figure 3d).

| Gene flow into the northern range limit
Gene flow in southern Ontario was more likely through the Niagara region from New York State, since we found that the land that connected both areas was more resistant to gene flow and this was due to low forest cover (Figures 3b and 4). From this point, gene flow was possible into central Ontario, since resistance decreased and current increased northward. Bobcat have been reported in central Ontario in the past, and occurrences are more common than in southern Ontario. In fact, bobcats were once common on the Bruce Peninsula (de Vos, 1964).
Gene flow to northern Ontario was not facilitated by the UPM and the forests of western Ontario as we previously hypothesized.
We found that gene flow into northern Ontario through the UPM is less likely and had more likely occurred through the LPM ( Figures   5,6). This was a surprising result, since the UPM previously seemed more likely because of the high amount of forest cover and because the LPM was an area where the bobcat was slowly recolonizing after it was extirpated ( Figure 2). However, considering our results the LPM is more appropriate because of the intermediate amount of forest cover which seemed to amplify gene flow (Figures 4 and 5).
In addition, snowfall on the UPM was much higher than the LPM.
In fact, the average annual snowfall on some areas of the UPM exceed well over 6 meters (Figure 3d). In contrast, the average annual snow on the LPM was shallower with depths not exceeding 4 meters, but compared to the UPM, these snowy areas occupied less of the land. Furthermore, the major thruway from the LPM to northern Ontario was across the 5.6-km-long Straits of Mackinac, which was only feasible in the winter when the lake was frozen. From our own experience tracking bobcats, it is not uncommon to observe bobcat that cross large bodies of water in winter. A bobcat with a GPS collar from our study crossed the North Channel to the Grant Islands, a distance > 5 km across Lake Huron ice (unpublished data).

| The importance of intermediate forest cover
Even if the Great Lakes had the highest resistance, forest cover had the highest average contribution to restricting gene flow. Though we previously thought that high tree cover would amplify gene flow, we found that intermediate amount of forest amplified gene flow.
Also, forests with 80% tree cover hindered gene flow more than areas with no tree cover (Figure 4). One common assumption is that the bobcat is a habitat generalist (Anderson & Lovallo, 2003); consequently, it does not specialize on any specific habitat type across its range, and therefore, it may perform better in environments with an average amount of forest cover compared with area that are 0 or 80% forest cover. Areas with low forest cover are either urban centers or areas that are predominantly used for agricultural purposes.
These areas are mostly found in the Midwest corn belt of the United States where tree cover was low over our study area (Figure 3b).
The corn belt was first cleared for agriculture in the 1850s and since then has been an area with low biodiversity and intensive agriculture use (Jenkins, Houtan, Pimm, & Sexton, 2015;Nassauer, Santelmann, & Scavia, 2007). Compared with more forested areas, the corn belt may have lower abundance and diversity of prey species that may not be able to sustain the bobcat.
On the other extreme, areas with high forest cover were areas where bobcat gene flow was obstructed, these forests were generally found in the upper Great Lakes region and were also found in the Appalachian corridor and the Adirondacks. Some of these forests with high amount of canopy cover also had high annual snowfall compared with other areas of the bobcat range in our study area (Figure 3b and d). Bobcat may not be able to effectively hunt in dense forests with high annual snowfall compared with forest with similar snow and intermediate forest cover.
The bobcat needs some forest cover to stalk prey (McCord, 1974), but perhaps the forest cover cannot be so dense as to reduce visibility and muffle sound, since the bobcat relies heavily on sight and sound to hunt (McCord & Cardoza, 1982); therefore, forests with an intermediate amount of cover might be more preferred by bobcat.
The ability to see and catch prey is a function of forest cover, but also, there is an interplay with snow depth, since deep snow reduces their ability to hunt (McCord, 1974). To some extent, forest cover was also associated with road density, in cases where road density is high, and forest cover was reduced, and this happened near urban areas (Figure 3b-c). In more rural communities, where forest cover is intermediate with road density, bobcat gene flow was amplified.

| Bobcat range expansion
In general, our results suggest that the northward expansion of the bobcat in the Great Lakes region has been facilitated by intermediate forest cover. Therefore, the expansion of the bobcat is in part a response to the decrease in forest cover due to land clearing and forestry in the Great Lakes region. In fact, in the northern Great Lakes US states the area occupied by open land has increased from 12.3% to 41.3% since 1836 (Schulte, Mladenoff, Crow, Merrick, & Cleland, 2007), which is within the range of optimal forest cover for bobcat gene flow that we found in our analysis (Figure 4). This land clearing may have opened previously unavailable habitat in northern Ontario to bobcat, but deep winter snow may have still been a limiting factor until snow depth subsided due to climate change in later years (Dyer & Mote, 2006).
After this point, further north, bobcat began to be harvested by trappers in northwestern Ontario in the early 1900s and in northeastern Ontario in the mid-1900s (de Vos, 1964. The 50-year lag period between these two areas could have been due to the disparity in the amount of annual snow received in both areas (Figure 3d).
Northwestern shores of Lake Superior received less snow than the northern shores of Lake Huron. Currently, on the north shores of Lake Huron intermediate forest cover is still important to bobcat, since bobcat in this area are almost exclusively found in rural communities within 50 km of urban centers (Marrotte et al., In Review).
The increasing density of human disturbances such as roads, rail lines, urban areas, rangeland, and agricultural land would have also further amplified colonization, because road plowing and snow compaction would have become more frequent, which allowed bobcat to move around and hunt more effectively. Therefore, at their northern limit, areas with intermediate forest cover that have an intermediate density of roads may have mediated the colonization of bobcats into areas the bobcat generally would not have occupied, because of deep annual snow.
Overall, land use and cover change and the decreasing snow pack due to climate change will only facilitate the expansion of bobcat and we can only expect to find bobcats farther north each year. However, it is important to note that the landscape of southern Ontario has impeded gene flow and consequently movement of bobcat over the past decades. As the climate continues to warm and species are tracking their bioclimatic niche through the Great Lakes region, we can only expect that less mobile species are less likely to cross southern Ontario. Other routes are already blocked by natural features such as the St. Lawrence River, the 5.6-km-long Straits of Mackinac, and the St Marys River, and the additive effect of human modification will undoubtedly further restrict these routes and reduce future potential biodiversity.

ACK N OWLED G M ENTS
This work was supported through funding from NSERC (a scholarship to RRM and a Discovery Grant to JB), Trent University, and OMNRF.

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

AUTH O R CO NTR I B UTI O N S
RRM and JB conceived the study. RRM, JB, and PJW collected the samples. RRM wrote code, ran simulations, analyzed the data, and wrote the manuscript. JB and PJW critically reviewed the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data and R scripts that support the findings of this study are openly available on Dryad at https ://doi.org/10.5061/dryad.zgmsb cc6d.