Landscape heterogeneity in landform and land use provides functional resistance to gene flow in continuous Asian black bear populations

Abstract Context Genetic diversity is one of the most important facets of biological diversity, and changes in the spatial pattern of habitats, often modified by human activity, are believed to have affected the genetic diversity of resident natural populations. Objectives We undertook a landscape genetic analysis in order to determine which landscape features influence gene flow within Asian black bear populations and to identify the underlying processes. Methods In our evaluation of gene flow, we estimated four parameters of resistance with regard to landscape elevation: the mean, the difference between the highest and lowest, the standard deviation, and the coefficient of variation of elevation among individuals. We then examined the resistance effect of different land use types. Results With the exception of mean elevation, we found that all parameters showed a significant relationship with genetic distance, indicating that unevenness in elevation provides functional resistance to gene flow. Although we found no evidence of landscape barriers (isolation‐by‐barrier), there was an indication of landscape resistance (isolation‐by‐resistance). Urban area and farmland are suggested to be the strong factors contributing to the resistance to gene flow, even though isolation‐by‐distance was also detected. When we examined gene flow for pairs of males and pairs of females, both isolation‐by‐distance and isolation‐by‐resistance were stronger in order of female pairs, male pairs, all individual pairs. Conclusions We conclude that landscape resistance was detectable with a high contrast in landscape heterogeneity and they are more influential on females than males. OPEN PRACTICES This article has been awarded Open Data badge. All materials and data are publicly accessible via the Open Science Framework at https://doi.org/10.5061/dryad.gn0qf16. Learn more about the Open Practices badges from the Center for Open Science: https://osf.io/tvyxz/wiki.


| INTRODUC TI ON
Genetic diversity is one of the most important elements of biological diversity because it determines the fitness and survival of individuals, the viability of populations, and the ability of species to adapt to environmental change (Balkenhol, Cushman, Storfer, & Waits, 2016). The spatial pattern of habitats influences organism perception and behavior, which drive the higher-level processes of population dynamics, gene flow, and adaptive evolution (Cushman, McRae, & McGarigal, 2016). However, human activity has, to varying degrees, modified the spatial pattern of habitats, which is believed to have affected the genetic diversity of resident natural populations.
Although many population genetic studies over the past few decades have revealed patterns in the genetic variation in mammal populations, few have statistically examined the underlying processes (Cushman, Shirk, & Landguth, 2013).
Population boundaries of the American black bear (Ursus americanus), which is known as a wide-ranging omnivore that is dependent on forest habitat, are genetically obscure when the forest habitat is contiguous and the species is widely distributed. In the Great Lakes region in Canada, clinal structuring of distribution is induced by isolation-by-distance (IBD) at 550 km and anthropogenic effects are rarely observed (Pelletier et al., 2012). In contrast, the findings of two studies have indicated that landscape features have a greater influence on the genetic structure of black bear populations in mountains on the border between the states of Montana and Idaho than does Euclidean distance (Cushman, McKelvey, Hayden, & Schwartz, 2006;Short Bull et al. 2011). On the basis of an examination of the effect of land use, no landscape barriers (IBB: isolation-by-barrier) were detected, whereas there was an indication of landscape resistance (IBR: isolation-by-resistance). Moreover, resistance due to elevation has been suggested to be lowest in the midelevation range (Cushman et al., 2006). The authors of a further study concluded that the effect of resistance was more detectable in landscapes with high heterogeneity in elevation and land use (Short Bull et al. 2011).
Simulation studies have clarified that a higher contrast in landscape resistance between habitat and nonhabitat increases the ability to detect the influence of landscape resistance on gene flow (Cushman et al., 2013;Jaquiery, Broquet, Hirzel, Yearsley, & Perrin, 2011). It has also been suggested that genetic differentiation is independent of Euclidean distance and significantly related to landscape structure when habitats are highly fragmented (Cushman et al., 2013). It has, however, been argued that it is difficult to identify the landscape features that influence genetic variation without reference to Euclidean distance, because in most real landscapes, the landscape features are generally aggregated into a few patches.
In the present study, we undertook a landscape genetic analysis in order to determine which landscape features influence gene flow within Asian black bear populations in the northern part of Japan and to identify the underlying processes. We considered that northern Japan would be a suitable area to examine these issues because it is characterized by a range of different topographical features and has a mix of both artificial and natural landscapes.
The present study was conducted in three stages. Initially, we determined the genetic structure of continuous populations of the black bear using Bayesian clustering analysis. We then proceeded to assess gene flow resistance in terms of four parameters of elevation: the mean, the difference, the standard deviation, and the coefficient of variation of elevation among individuals. Finally, we examined the resistance effect of land use by comparing IBD, IBB, and IBR. We conducted these analyses for all sample pairs, pairs of males, and pairs of females. Due to male-biased dispersal, we predicted that the relationships between genetic distance and each parameter would be detected in the following order, female pairs >all sample pairs >male pairs.

| DNA samples
DNA samples were obtained from a total of 148 individuals from the northern margins of a large continuous population of black bears in Japan. Among these, 44 bears in Aomori Prefecture were identified F I G U R E 1 An Asian black bear cub with its mother in Iwate Prefecture, Japan. Photo by Y. Sato from bear hairs collected by Yamamoto et al. using noninvasive hairtrapping (Yamamoto et al., 2012). A further 100 tissue samples were obtained in four other prefectures, and four samples were collected from bears in the east of Aomori Prefecture that had been hunted or captured for pest control between 1991 and 2011.

| Bayesian clustering and genetic statistics
The genotype data were analyzed using STRUCTURE Bayesian clustering software (Pritchard, Stephens & Donnelly, 2000) to determine the genetic structure of the black bear population. Initially, we calculated P ID 10 times for each K at 1-10 with 5,000,000 Markov chain Monte Carlo (MCMC) iterations, with a burn-in of 100,000, and independent allele frequencies (lambda 1.0), using an admixture model (alpha inferred, with initial alpha set to 1.0). The deltaK had a single peak when K = 2, thus we set K to 2, and assumed population inferred proportions using 10,000,000 MCMC iterations, with a burn-in of 1,000,000.
After the Bayesian clustering, we estimated genetic statistics of each assumed population. Departure of the observed heterozygosity from the Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) were examined using the Markov chain method with 10,000 permutations (Guo & Thompson, 1992) by using Arlequin 3.5.2.2 (Excoffier, Laval, & Schneider, 2005). The sequential Bonferroni method was used to adjust significance values for all multiple comparisons (Rice, 1989). The presence of null alleles and genotyping errors for each locus and each assumed population were examined using MicroChecker 2.2.3 (Van Oosterhout, 2004

| Analytical unit
We used ~1-km grid cells for spatial analysis. This is one of the standard units used in Japan for gathering statistics relating to populations, land use, and wildlife. For analysis, we demarcated the sites of sample collection in grids.

| The influence of landform on gene flow
We used Mantel tests to investigate the relationship between genetic distance and elevation parameters for all pairs of individuals. We constructed distance networks by initially connecting two points with a straight line for all individual pairs, that is, male-female pairs, all male pairs, and all female pairs. If the two points were located within the same 1-km grid, Euclidean distance was defined as 0. If the two points were located in neighboring grids, Euclidean distance was approximated as 1 km according to the analytical unit. Finally, Mantel tests (Manly, 1997) for genetic distance and candidate elevation parameters were performed with 999 permutations using the vegan 2.5-1 package in R.

| The influence of land use on gene flow
We examined which landscape features have the potential to function as barriers or resistance to gene flow among individuals. For this analysis, we calculated the amounts of different land use types within the network lines. To this end, we used seven land use types (forest, grass field, farmland, urban area, open water, special matrix, and wetland) using the 1-km cell grids of a land use classification map (Ogawa et al., 2013). The special matrix includes natural bare ground, limestone vegetation, volcanic desert, solfatara desert, and cliffs. We collated the data on land use for all network lines and calculated the number of classes in each network. In total, we assessed 138.4 billion grids for all networks (Table 1). Among these, forest represented the largest (64.6%) area, followed by grass field (20.3%).
Special matrix, which comprises areas such as cliffs and bare ground, occurred in 72 million grids (0.052%).
We investigated the resistance effect of land use on gene flow by using simple and partial Mantel tests according to the concept of Ruiz-Gonzalez et al. (2014). Initially, we conducted a simple Mantel test (G × R L ) for genetic distance and the candidate land TA B L E 1 The number of grids of land use types in the study area use class with 999 permutations. In this analysis, we calculated the Mantel r value between genetic distance and the resistance value of each candidate land use class. Here, resistance values were obtained for the product of the number of grids of the candidate land use class within the network and presumed resistance (2, 5, 25, 50, and 100). Assigning resistance values to landscape features is a well-known challenge in landscape genetics (Tucker, Allendorf, Truex, & Schwartz, 2017), and thus, we assigned each landscape feature with a wide range of resistance values. Next, we performed partial Mantel tests between genetic distance and the resistance value (2, 5, 25, 50, and 100) of the candidate land use with partialling out of Euclidean distance of the network (G × R L |Dis) and between genetic distance and Euclidean distance with partialling out the resistance value of candidate land use (G × Dis|R L ). When G × R L |Dis was statistically significant (p < 0.05) and G × Dis|R L was negative or not statistically significant (p > 0.05) in these partial Mantel tests, we regarded that the candidate land use influenced gene flow via resistance according to Ruiz-Gonzalez et al. (2014).
Subsequent to these tests, we additionally estimated the correlation between genetic distance and the combined resistance of land uses by using simple Mantel tests (G × ΣR L ).
We also performed simple and partial Mantel tests to determine whether different land uses function as barriers to gene flow with the same resistance effect. We examined the correlation between genetic distance and the presence/absence of a candidate land use (0 or 1) using a simple Mantel test (G × B L ). Partial Mantel tests were used to assess the relationship between genetic distance and the presence of candidate land use with partialling out of Euclidean distance of the network (G × B L |Dis) and between genetic distance and Euclidean distance with partialling out of the presence of a candidate land use (G × Dis|B L ). When G × B L |Dis was statistically significant (p < 0.05) and G × Dis|B L was negative or not statistically significant (p > 0.05) in these partial Mantel tests, we regarded the candidate land use as representing a functional barrier to gene flow according to Ruiz-Gonzalez et al. (2014).

| Genetic structure
Bayesian clustering analysis inferred five genetic clusters from the 148 genotyped samples. Among these, 126 samples had a high assignment ratio (>0.8) to one cluster when K = 2, and a spatial structuring of clusters was observed ( Figure 2). After the Bayesian clustering analysis, we checked the validity of markers for each assumed population with K = 2. Null alleles were indicated at two loci (G10B and MSUT-7) in larger assumed population, but neither LD nor HWE deficiency was observed in both assumed population.

| The influence of landform on gene flow
The results of Mantel tests, which were used to estimate the relationship between the genetic distance and elevation parameters of all individual networks, are shown in Table 2. We found that the mean elevation within each network was not significantly related to the genetic distance for all individuals, males, and females. In contrast, the difference between the highest and lowest elevation in each network, and the SD and CV within each network all showed higher Mantel r values than the null model (the Euclidean distance) for all three pair categories and male pairs. CV had highest Mantel r values for all three network categories, and only CV showed higher Mantel r values than the Euclidean distance in female pairs. For all parameters, females had larger Mantel r values than males.

| The resistance function of land use on gene flow
We performed simple Mantel tests to determine the relationship between genetic distance and each resistance value (2, 5, 25, 50, and 100) of land use types (G × R L ,), and found that for all sample pairs, farmland, and urban area showed higher Mantel r values for all F I G U R E 2 Maps of the study area and relevant adjacent regions: light shadings show the distribution of the Asian black bear (Ursus thibetanus) in Japan (a); colored maps represent elevation (b) and land use (c) of the study area; sites where bears were captured (b and c); and the different symbols indicate the location of genetic clusters inferred using the Bayesian clustering software STRUCTURE resistance values than the null model (the Euclidean distance: Mantel Table 3). The resistance indicating the highest Mantel r was 25 for both land use types (Mantel r = 0.0779 for farmland and 0.0721 for urban area). Next, we performed partial Mantel tests with partialling out of Euclidean distance of the network (G × R L |Dis) and between genetic distance and Euclidean distance with partialling out the resistance value of the candidate land use (G × Dis|R L ).

| The barrier function of land use on gene flow
Analysis of the correlation between genetic distance and the presence of candidate land use type using a simple Mantel test (G × B L ) revealed that Mantel r values for farmland (0.0787) for all sample pairs and forest (0.2558) for female pairs were higher than those between the Euclidean distance and genetic distance (0.0687 and 0.2373, respectively), whereas no land use type showed higher a Mantel r value than the Euclidean distance for male pairs (Table 4).
When we analyzed the barrier function of farmland area using partial Mantel tests to assess the relationship between genetic distance and the presence of farmland with partialling out of the Euclidean dis-

| D ISCUSS I ON
Gradient concepts such as IBD are sometimes consistent with the genetic structure of continuous populations (Ohnishi, Kobayashi, Nagata, & Yamada, 2017;Pelletier et al., 2012). Although we also detected an IBD effect in the Asian black bear population investigated in the present study, Bayesian clustering analysis suggested a certain degree of genetic clustering. This prominent structuring could be a consequence of gene flow disruption caused by landscape heterogeneity. We determined which landscape factors influence the gene flow producing genetic clusters and concluded that they provided functional resistance to gene flow. Cushman et al. (2013) have explained the IBD patterns determined by landscape heterogeneity based on behavioral ecology.
Animal movement behavior is selected to maximize fitness with respect to resources, such as food, dens, and reproductive partners, and to minimize the risk of predation. Continuous uniform landscapes, by definition, have low heterogeneity, thereby inducing animals to select differential paths. In such landscapes, movement generally conforms to a random walk model, which gives rise to an IBD pattern and does not support IBR. In contrast, in landscapes with high heterogeneity and fragmented habitats, movement paths will be selected primarily to avoid unsuitable conditions. As a consequence, an IBR pattern in path selection will be observed.

Model of elevation
In the present study, we confirmed these assumptions based on our analyses of different landscape features, namely, elevation and land use. The maximum difference and unevenness in elevation between individuals represent functional heterogeneities for the present population. Although we were unable to detect an effect of elevation per se in the black bear population we studied, elevation has been implicated in studies on American black bears (Cushman et al., 2006;Short Bull et al. 2011). The difference in the findings of the American black bears studies and those of the present study may be explained by differences in the elevational ranges in the respective study areas. The elevational ranges of the American black bears study areas are from <2,000 m to over 7,000 m, whereas, the highest mountain in the present study area reaches only 2,236 m. Both Asian and American black bears are omnivores, and the foraging diet of individuals is believed to be learnt from their mothers (Kitamura & Ohnishi, 2011;Mazur & Seher, 2008). Given that vegetation changes according to elevation, American black bears would tend to disperse to areas with an elevation similar to that of their birthplace to locate familiar foods. In contrast, there is little restriction in accessing vegetation for the bears in the present study area due to the relatively low elevation. Consequently, a change in elevation on the route of bear movements would not represent a resistance to gene flow. It has been suggested that elevation does not influence gene flow in American black bear populations in landscapes where the topography is relatively flat and elevation is not highly variable (Short Bull et al. 2011). However, the results of the present study indicate that a topographic influence can be detected in nonflat landscapes, even Euclidean distance were higher for females than for males, and the resistance value for farmland was four times larger for females than for males. This tendency can probably be explained by differences in the dispersal patterns and social systems of males and females.

G*RL G*RL|Dis G*Dis|RL
Natal philopatry in female-and male-biased dispersal have been suggested (Ohnishi & Osawa, 2014), and matrilineal site fidelity has also been detected in black bears (Kozakai et al., 2017). Additionally, a recent behavioral study also suggested that female tended to avoid farmland more than male (Takahata, Takii, & Izumiyama, 2017). An interesting point in this regard is that special matrix and wetland appear to show resistance to females, even though a partial Mantel test did not support the resistance of land use. We assume, however, that farmland has a resistance function because it represents an artificial land use, whereas special matrix and wetland are natural landscapes.
These results tend to indicate that females are strongly dependent on forests and grass fields for their habitat. In contrast, these natural land uses did not show resistance to males, but we do not believe that these males are independent of forests and grass fields. The capacity of long-distance dispersal and larger home ranges would enable males to circumvent these particular landscape types. Moreover, the difference in habitat preference would also cause these results be-  (Tucker et al., 2017). We also predicted that Mantel r values for all sample pairs are lower than those for female pairs but higher than those for male pairs; however, we observed that the Mantel r values for all pairs were actually lower than those for the pairs of each sex. The same tendency in the Mantel r values of IBD has been confirmed in a previously published study (Ohnishi, Yuasa, Morimitsu, & Oi, 2011), in which authors used hair samples collected in different areas using hair traps (−0.113 for females, −0.112 for males, and −0.099 for all sample pairs), even though they obtained negative Mantel r values because of relatedness as genetic distance. Therefore, we assume that this tendency could be a common characteristic of black bear populations in Japan.
Accordingly, it may also be necessary to consider relationships among individuals not only in terms of sex-biased dispersal patterns but also with respect to behavior pattern. Strong matrilineal relationships after natal dispersal is suggested in females (Kozakai et al., 2017) and there are assumed to be certain male-male relationships, although to date these have not been investigated.
Furthermore, difference in habitat selection between sexes in each season would tend to cause lower IBD and IBR. For example, females select subalpine forests during spring and deciduous forests during summer, but males do not prefer to them during both seasons (Takahata et al., 2017). Accordingly, we propose that these male-male and female-male relationships should be studied combined with behavioral investigations.
On the basis of our current findings, we could reject the likelihood of IBB operating in the bear population we studied, which has also been indicated in the case of American black bears (Cushman et al., 2006). Whereas resistance can show degrees of variability, barriers tend to be absolute (i.e., 0 or 1). Such an immoderate habitat would be basically island surrounded by sea or lake for terrestrial animal. We can thus assume that barriers provide extremely strong resistance to gene flow. For example, Saitoh et al. (2001) and Ohnishi et al. (2007) revealed large genetic differentiation among the black bear populations fragmented by rivers in western Japan and concluded that the rivers constitute functional barriers to gene flow among fragmented populations. However, low frequency movements across rivers are also suggested (Ishibashi & Saitoh, 2004;Ohnishi et al., 2009).
Nevertheless, given that residential areas and farmlands tend to concentrate along either side of river courses, the landscape in the vicinity of rivers would have strong functional resistance.

ACK N OWLED G M ENTS
We thank H.B. Tamate for encouragement throughout the study, and also T. Oka and T. Tsubota for providing samples collected in Aomori and Miyagi prefectures, respectively. We also thank the local governments and the branches of the Japan Hunting Association for collecting samples in the following prefectures: Iwate, Akita, and Yamagata. This work was supported by JSPS KAKENHI Grant numbers 20380098 and JP18K06438.

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

DATA ACCE SS I B I LIT Y
We submitted genotype data to the Dryad database on 12 July