In the absence of obvious barriers to dispersal microsatellite studies of vagile mammalian carnivores frequently find panmictic-like genetic structure over wide scales, whereas high levels of differentiation at much finer scales are detected with mitochondrial DNA (mtDNA). Given the maternal inheritance of mtDNA, these differences are often attributed to male-biased dispersal, remnants of postglacial range expansion, or both. Based on such contrasting results, it is not always clear how to delineate contemporary populations. We investigated the genetic structure of American black bears (Ursus americanus) over a wide geographic area (>1,700 km) that has no obvious physiogeographic barriers to gene flow. We analyzed a 315-base pair fragment of the mtDNA control region from 660 individual bears from 23 regions of Ontario, Canada. Relative to black bear studies based on nuclear data, mitochondrial analyses revealed much stronger patterns of genetic structure among regions (0.09 < FST < 0.44), even at small-scale intervals (<150 km), which likely reflects strong female philopatry combined with male-biased dispersal. The patterns of genetic differentiation among regions were consistent with previously described historical patterns in black bears, specifically the division of the species into 2 phylogeographic clades (coastal and continental). We confirmed that further subdivision of the continental clade occurs in a region where obvious physiogeographic barriers do not exist. We postulate that this small-scale differentiation can be explained by residual patterns from postglacial recolonization routes on either side of the Great Lakes. We suggest that it was maintained through extreme female philopatry due to habitat saturation following the postglacial geographic expansion. Based on our results, we propose that a combination of several molecular markers can be more useful in defining population units for conservation and management decisions than biparentally inherited microsatellites.
Understanding the genetic structuring patterns of species increases our knowledge of their ecology and evolution and helps inform conservation and management strategies directed toward maintaining stable populations. A problem that arises, however, is that research using combinations of neutral molecular marker types that have different rates of evolution and modes of inheritance can reveal contrasting patterns of genetic differentiation (Brito 2007; Flanders et al. 2009; Hellborg et al. 2002; Johnson et al. 2003). For this reason studies that incorporate both biparentally inherited nuclear microsatellites and maternally inherited mitochondrial DNA (mtDNA) are becoming increasingly important in describing population delineations of highly vagile species that have different male and female life histories (Chappell et al. 2004; Tomasik and Cook 2005).
In the absence of physiogeographic features that impede the movement of animals, contrasting levels of differentiation detected between microsatellite and mtDNA analyses have been explained by factors such as low dispersal distances, reduced effective population size of mtDNA, long-term isolation of historical lineages, cryptic boundaries, and sex-biased dispersal (Irwin 2002; Tomasik and Cook 2005). For large and mesocarnivores, for which topographic barriers to dispersal were perceived to be absent, microsatellite analyses have shown panmictic structuring patterns, even across large distances (marten [Martes americana—Kyle and Strobeck 2003], brown bear [Ursus arctos—Paetkau et al. 1998], and lynx [Lynx canadensis—Schwartz et al. 2002]). At large scales genetic differentiation can be explained by factors such as isolation by distance (marten [Broquet et al. 2006], wolf [Canis lupus—Geffen et al. 2004], and puma [Puma concolor—McRae et al. 2005]), or anthropogenic and natural influences acting as barriers to dispersal (wolverine [Gulo gulo—Kyle and Strobeck 2002], puma [McRae et al. 2005], and bobcat [Lynx rufus—Millions and Swanson 2007]). However, even over small distances across which no barriers exist, genetic structure can be observed in wide-ranging species (black bear [Ursus americanus—Peacock et al. 2007] and lynx [Rueness et al. 2003]). This suggests that factors such as the maintenance of historical lineages due to an intermediate level of dispersal (Peacock et al. 2007), or cryptic differentiation (Rueness et al. 2003), also play a role in contemporary structuring patterns that are not necessarily defined by the current landscape.
Mitochondrial DNA studies of North American taxa focusing on the identification of such historical lineages show that many species share similar patterns of genetic structure at the continental scale, reflecting common physiogeographic patterns (Arbogast 1999; Aubry et al. 2009; Byun et al. 1997; Conroy and Cook 2000; Demboski et al. 1999; Demboski and Sullivan 2003; Wooding and Ward 1997). Two main mtDNA clades are identified most often, a geographically restricted coastal clade found along the North Pacific Coast and a widespread continental clade. Because many of the species sharing this pattern of differentiation are associated with forest, the genetic division between these 2 main clades has been attributed to the existence of isolated forest refugia located on opposite sides of the continent during the last glaciation.
Although molecular studies focusing on North American species are numerous, few of them deal with highly vagile species found across extensive sampling areas free of physiogeographic barriers. For this reason we see a lack of comprehensive studies of genetic differentiation focusing on taxa that are both continuously and widely distributed, although such research would provide a base of comparison for studies that identify genetic discontinuities in isolated and fragmented populations. The American black bear is no exception, because studies of this species focus mostly on isolated populations that have arisen from habitat loss and human-caused mortality in the southern part of its range (Vaughan and Pelton 1995), such as Florida (Dixon et al. 2006), Louisiana (Csiki et al. 2003; Larkin et al. 2004), Arkansas (Csiki et al. 2003; Van Den Bussche et al. 2009), or Mexico (Onorato et al. 2004), whereas few data exist on the core population that remains in the northern part of the distribution. Although they are capable of extensive dispersal movements of more than 200 km from their natal site (Lee and Vaughan 2003; Rogers 1987), American black bears show relatively high levels of genetic structuring across their range (Csiki et al. 2003; Onorato et al. 2004; Paetkau and Strobeck 1994). Most of these were small-scale population divisions identified through genetic studies and have been attributed to isolation by physiogeographic barriers (islands [Paetkau and Strobeck 1996] or water bodies [Peacock et al. 2007; Robinson et al. 2007]), geographic distance (Mills 2005), and landscape gradients (Cushman et al. 2006; Dixon et al. 2006).
At the continental scale molecular studies of black bears have focused mostly on populations in the western and southern portions of the current distribution (Byun et al. 1997; Paetkau and Strobeck 1996; Peacock et al. 2007; Stone and Cook 2000; Van Den Bussche et al. 2009; Wooding and Ward 1997), resulting in a lack of data from the eastern part of the range (Fig. 1). Ontario, Canada, comprises a significant proportion of the contemporary range of black bears in eastern North America, with a population size of about 100,000 individuals (M. E. Obbard, pers. obs.). In addition, with the exception of the far southern regions of the province, the landscape is largely continuous and homogeneous, with no obvious physiographic features such as large rivers, mountains, or drastic habitat change that would impede dispersal.
Fig. 1
Map of sampling locations of black bears, including our Ontario sites and sites from other studies across North America. The Coastal phylogeographic clade is represented by dark gray triangles and the Continental clade by circles. The circles representing the 2 continental subclades are black for Continental Western subclade and light gray for Continental Eastern subclade. The circles that are both dark gray and black, located in Montana, represent sites where bears from both the Coastal and Continental clades were found. The circle that is both light gray and black, located in Oklahoma, represents a site where bears from both the Continental Western and Continental Eastern subclades were found. Based on the information provided by the median joining network (Fig. 3), the Ontario samples were attributed to each of the Continental subclades, the Northwest cluster belongs to the Continental Western subclade, and the other Ontario clusters belong to the Continental Eastern subclade.
We analyzed mtDNA sequences of the control region of black bears obtained from hair samples of 660 individuals from 23 locations across Ontario (Fig. 2). Our 1st goal was to assess the mitochondrial genetic structure of black bears across a 1,700-km continuum in a landscape that is largely homogeneous. We hypothesized that within this continuous landscape, due to black bear male-biased dispersal combined with female philopatry (Lee and Vaughan 2003; Rogers 1987) and the maternal mode of inheritance of mtDNA, our results would show strong differentiation among regions, a pattern that would not be observed with biparentally inherited neutral markers (Chappell et al. 2004; Johnson et al. 2003; Tomasik and Cook 2005). Our 2nd goal was to place the mtDNA results into a wider continental context by clarifying how Ontario black bears relate to other North American populations.
Fig. 2
Frequencies of haplotypes represented >5 times in the 23 sampling locations distributed across Ontario (small pies) and grouped into 4 geographic clusters (large pies). Triangles represent the sampling sites in the Bruce cluster, diamonds represent the sites in the Southeast cluster, circles represent the sites in the Central cluster, and squares represent the sites in the Northwest cluster.
These results, in light of other microsatellite studies of black bears, can help determine if the observed differentiation is caused by contemporary factors that promote genetic structuring, such as territoriality and natal philopatry, or if it could be explained by historical events that illustrate the consequences of continental processes, such as long-term isolation and colonization. In addition to providing further insights into the ecology of black bears, our findings have implications regarding how data from several molecular markers with different underlying evolutionary histories can be assimilated and interpreted in the implementation of conservation and management plans for large carnivores.
Materials and Methods
Sample collection and DNA analysis
Between 1997 and 2007 black bear hair samples were collected along trap lines grouped in 23 sampling sites located across Ontario (Fig. 2). Samples were obtained both opportunistically (livetrapping, hunting, or road kills) and from baited barbed wire hair traps (Woods et al. 1999). These procedures were consistent with the animal care guidelines approved by the American Society of Mammalogists (Gannon et al. 2007).
For hair samples collected from 1997 to 2004 DNA extraction was performed using a modified version of the DNeasy tissue extraction protocol (Qiagen, Mississauga, Ontario, Canada). For each individual sample 10–15 hairs with visible roots were suspended in a solution containing 500 µl of 1X lysis buffer (Applied Biosystems, Inc., Streetsville, Ontario, Canada) and incubated in 10 µl proteinase K (Qiagen) at 37°C for 12 h. After incubation, standard Qiagen tissue extraction procedures were followed. Samples collected from 2004 to 2007 were extracted following a MagneSil paramagnetic bead automated DNA extraction procedure (Promega, Nepean, Ontario, Canada) using a P3 Evolution (Perkin Elmer, Woodbridge, Manitoba, Canada) liquid handler, eluting in a final volume of 75 µl.
A 315–base pair (bp) fragment of the mtDNA black bear control region was amplified by polymerase chain reaction using the primers H16498 (Ward et al. 1991) and L15997 (Wooding and Ward 1997). The sequences were obtained from black bears that had been identified individually based on 15 microsatellite loci (C. J. Kyle and M. E. Obbard, pers. obs.; Mills 2005) and sex analyses (primers S47 and S48—Ennis and Gallagher 1994). DNA amplification reactions contained 1X polymerase chain reaction buffer, 0.2 mM of deoxynucleoside triphosphates, 2.0 mM of MgCl2, 0.3 mg/ml of bovine serum albumin, 0.2 µM of each primer, 1 U of Taq DNA polymerase/µl (Invitrogen Corp., Burlington, Ontario, Canada), and 10 ng of DNA extract as a template. Amplification reactions were run on a Dyad Disciple Peltier Thermal Cycler (Bio-Rad Laboratories, Inc., Mississauga, Ontario, Canada) programmed for an intial 5-min denaturation step at 94°C, followed by 35 cycles of the following steps: denaturation at 94°C for 1 min, annealing at 60°C for 1 min, and extension at 72°C for 1.5 min. The extension was completed after a final extension step at 60°C for 45 min. Amplified products were separated and quantified via electrophoresis in 1.5% agarose gels. Polymerase chain reaction products were purified using QIAquick Purification Kit (Qiagen) to remove excess primers and deoxynucleoside triphosphates. Forward and reverse sequences were obtained by using BigDye Terminator Cycle Sequencing Kit version 3.1 (Applied Biosystems, Inc.). Sequencing was performed on an automated DNA sequencer (ABI 3730; Applied Biosystems, Inc.).
To facilitate analyses the 23 sampled sites (Appendix I) also were pooled into 4 geographic clusters (Bruce, Southeast, Central, and Northwest) based on both geographic proximity and microsatellite data that suggested weak genetic divisions between these broad geographic regions (Mills 2005).
Sequence analysis
We profiled 660 individuals for which mtDNA fragments were edited and aligned manually with MEGA 4.1 (Tamura et al. 2007) relative to previously identified haplotypes downloaded from GenBank (Appendix II). Sequences that did not align to previously identified haplotypes in the literature were considered new haplotypes only after resequencing with the reverse primer to confirm the sequence. All the sequences obtained in this study were submitted to GenBank (accession numbers GU724158–GU724193).
Haplotype frequencies were calculated with FaBox (Villesen 2007), and levels of genetic diversity were estimated using Arlequin 3.1 (Excoffier et al. 2005) by calculating haplotype diversity (h, the probability that 2 haplotypes drawn randomly from a population are different), nucleotide diversity (π, the mean number of pairwise differences per site between 2 sequences—Nei and Kumar 2000), and genetic divergence (FST—Weir and Hill 2002). Due to discrepancies in sample sizes between clusters (Bruce: n = 38 individuals; Southeast: n = 321; Central: n = 126; Northwest: n = 175), we conducted a rarefaction analysis with ADZE-1.0 (Szpiech et al. 2008) to standardize the levels of haplotypic diversity.
We tested for departure from the neutral model of evolution and population growth by computing Tajima's D (Tajima 1989), and Fu's FS (Fu 1996) tests using Arlequin 3.1 (1,000 permutations). The D-test compares the number of segregating sites in the sample to the mean number of pairwise differences between haplotypes, whereas the FS-test determines the probability of obtaining the observed number of haplotypes given the observed average number of pairwise differences. To achieve an alpha of P = 0.05 for the rejection of the null hypothesis of neutrality FS must be negative (indication of population expansion) and its P-value < 0.02.
Genetic structure
In our study design, where possible, we selected 30 individuals per sampling site to conduct the mtDNA analyses. We estimated the degree of differentiation among sampled sites, and among the broad geographic clusters, by calculating pairwise FST values (Weir and Cockerham 1984) in Arlequin 3.1 (1,000 permutations, P < 0.05—Excoffier et al. 2005). To evaluate the optimal grouping pattern of the sampled sites without a priori assumptions we conducted a spatial analysis of molecular variance (SAMOVA, 1,000 initial conditions—Dupanloup et al. 2002). For this analysis we grouped Algonquin (n = 50) and Bracebridge (n = 5) together, because the only 2 haplotypes found in Bracebridge (HAP1 and HAP2) were common in Algonquin. Of all our SAMOVA results (our sites divided into K = 2–15 groups), those that had the highest variance among clusters (FCT) were reported (K = 2, 4, and 11). Two of these (K = 2 and 4) were compared to the results of an analysis of molecular variance (AMOVA—Excoffier et al. 1992) based on the geographic clusters determined a priori. Both SAMOVA and AMOVA comparisons examined the partitioning of genetic variation among clusters, among sampled sites within clusters, and within sampled sites.
Pairwise FST values also were used to perform Mantel tests (Mantel 1967) to establish whether the level of genetic differentiation was correlated with geographic distance between sampled sites (Wright 1943). In addition, a partial Mantel test was conducted to model a barrier to gene flow between the a priori defined geographic clusters. Through the Isolation by Distance Web Service version 3.14 (Jensen et al. 2005) we regressed the pairwise FST values between all the sampling sites against pairwise geographic distances (km) using 1,000 randomization steps. Geographic distances between each sampling location were obtained by plotting the samples in ArcGIS version 9.0 (Environmental Systems Research Institute, Inc., Toronto, Ontario, Canada) and by calculating the distance between their centroids.
Phylogenetic analyses
The relationships between the haplotypes found in Ontario were estimated by creating a median joining network (Bandelt et al. 1999) with the haplotypes that were observed in more than 5 individuals (or >0.9%) across the entire data set, using the software Network version 4.5 (Network 2008). The cluster differentiation found in Ontario black bears then was assessed with respect to the phylogeographic structure identified at the continental scale (Wooding and Ward 1997) by integrating sequences from Ontario with all available black bear haplotypes (Appendix II). As intraspecific haplotype differences can be low (e.g., only 1 nucleotide substitution), phylogenetic relationships can be represented accurately by a haplotype network (Posada and Crandall 2001). Therefore, we constructed a 2nd median joining network that included both the Ontario and all other North American sequences to clarify the relationships among haplotypes at the continental scale.
Results
Genetic diversity in Ontario
The analysis of the 315-bp fragment of the mtDNA control region obtained from the 660 black bear samples identified 36 haplotypes and 26 variable sites (GenBank accession numbers GU724158–GU724193; Appendix II). Of these haplotypes, 11 were observed previously and 25 were newly identified, of which 14 were identified only in a single individual (Appendixes I and III). Eight haplotypes were observed in >1 but <6 individuals (relative frequency < 0.9%), and 5 (HAP1, HAP2, HAP5, HAP6, and HAP15) had a relative frequency > 10%.
The neutral model of evolution could not be rejected for any of the sampled sites using either D or FS. The values for D ranged from 0.000 for the Bruce Peninsula National Park (BPNP; P = 1.000) to 20.974 (Bracebridge; P = 0.999), with an overall value of D = 9.594 (P = 0.926). Results from the FS-test ranged from 9.038 (Kenora; P = 0.895) to 0.256 (Borland + Ivanohe; P = 0.178), with an overall value of FS = 3.437 (P = 0.591). Population-specific FST values ranged from 0.281 (Red Lake) to 0.323 (BPNP).
Haplotypic diversity within sampled sites ranged from 0.419 (Fort Frances A) to 0.893 (Midhurst), with an overall haplotypic diversity of 0.691. Nucleotide diversity (π) within sampled sites ranged from 0.002 (BPNP) to 0.026 (Red Lake), with an overall nucleotide diversity of π = 0.015 (Appendix III).
At the cluster level a high genetic diversity was detected in the Southeast (0.832) compared to the other clusters (Bruce: 0.501; Central: 0.680; Northwest: 0.753). The low haplotypic diversity found in the Bruce cluster was not a consequence of smaller sample size, because when standardized, its value remained low compared to the other clusters (Table 1).
Table 1
Comparison of observed and standardized haplotypic diversity (obtained with the rarefaction analysis conducted with ADZE-1.0—Szpiech et al. 2008) in each cluster.
Distribution of haplotypes in Ontario
Haplotypic distribution varied both among and within black bear clusters. Among clusters, strong differences in the frequency of the most common haplotypes were observed. The 2 predominant haplotypes in the Southeast cluster (HAP1 and HAP2, cluster frequencies = 26%) differed from the haplotype most frequently observed in the Central cluster (HAP5, cluster frequency = 54%) and from the one that was predominant in the Northwest cluster (HAP15, cluster frequency = 46%; Fig. 2). Only 2 haplotypes were found in the isolated Bruce cluster. One of them was predominant in Central Ontario (HAP5) but had lower frequencies in the other main clusters, whereas the other was common in Southeast Ontario (HAP2) but found in very low frequencies in the Northwest and Central clusters (Appendix I; Figs. 2 and 3). Finally, among the 8 rare haplotypes other than singletons found in Ontario (relative frequency < 0.9%), 100% (n = 8) were restricted geographically to their respective cluster, and 62.5% (n = 5) were restricted to 1 sampling site within a cluster.
Fig. 3
a) Median joining network including sequences from sites located across North America, obtained from GenBank (see Appendix II for citations), and the Ontario sequences that appeared >5 times in the data set. Each color represents the sites where the respective sequences were found (states or provinces), and the slash marks represent the genetic distance (number of base pair differences) between each haplotype. Circle size is not proportional to the frequency of the haplotypes. The haplotypes within the black circle belong to the Continental Western subclade identified in black bears (Wooding and Ward 1997), whereas the rest of the sequences belong to its Continental Eastern subclade, with the exception of H8, H16, H17, H18, and H19, which belong to the Coastal clade of black bears. The red dot is a median vector, which is a potential common ancestor between ≥2 haplotypes. b) Ontario median joining network showing the relationships between the most frequent haplotypes from our study area, based on 315-bp fragments of the mitochondrial control region of black bears. Circle size is proportional to the total number of individuals sharing each haplotype, and slices are proportional to the number of individuals per cluster carrying a particular haplotype. The Bruce cluster is shown in white, the Southeast cluster in light gray, the Central cluster in black, and the Northwest cluster in dark gray. The most frequent haplotypes (HAP15 and HAP24) in the Northwest cluster belong to the Western Continental clade subdivision, whereas all the other Ontario sequences are found in the Eastern Continental clade subdivision.
Genetic structure
Various levels of differentiation were observed among sampling sites, with pairwise FST values (Weir and Cockerham 1984) ranging from 0.849 (BPNP/Fort Frances A; P = 0.000) to −0.078 (Midhurst/Sioux Lookout; P = 0.504). The results indicated that BPNP bears showed a higher degree of genetic differentiation when compared to bears from the other sampling localities.
All of the geographic clusters were highly differentiated from each other, with the lowest level of divergence found between the Central and Southeast clusters (FST = 0.120; P = 0.000) and the remainder of the values ranging from 0.419 (Bruce/Northwest; P = 0.000) to 0.210 (Bruce/Central; P = 0.000; Table 2). In addition, the Northwest cluster was more differentiated from the Central cluster (FST = 0.301; P = 0.000) than the Bruce cluster.
Table 2
Pairwise comparison of mitochondrial DNA (mtDNA) genetic differentiation (FST values) for American black bears (Ursus americanus) between the geographical regions of Ontario. Pairwise FST values are located below the diagonal (in italic type) for microsatellites (Mills 2005) and above the diagonal for mtDNA. Significant (P < 0.05) values are indicated with an asterisk (*). All P-values for the mtDNA data were P = 0.000. Sampling sites are mapped in Fig. 2.
The SAMOVA result that had the highest variance among groups was K = 2 (FCT = 0.343, P = 0.000), which separated Dryden, Fort Frances A, and Thunder Bay (all included in the Northwest geographic cluster) from all the other sampled sites in Ontario. K = 4 also had a high variance among groups (FCT = 0.311, P = 0.000), and grouped Dryden/Fort Frances A/Thunder Bay, Kenora/Fort Frances B (also included in the Northwest geographic cluster) together, Timmins (Southeast cluster) on its own, and the rest of the Ontario localities (Table 3). At K = 11 (FCT = 0.306, P = 0.000), all the sampled sites that were grouped into the Central cluster stayed together (Borland + Ivanhoe/Chapleau Crown Game Preserve/Hearst/Nipigon), as did most of the sites grouped into the Southeast (Sudbury/Bancroft/North Bay/Sault Ste. Marie, Midhurst/Parry Sound, and Algonquin–Bracebridge/Pembroke) and Northwest (Dryden/Fort Frances A/Thunder Bay, and Kenora/Fort Frances B) clusters (Table 3). The corresponding AMOVA based on geographic clusters (K = 2 corresponding to Northwest versus the other clusters, K = 4 to all the geographic clusters) also indicated genetic differentiation, but at lower levels than SAMOVA (FCT = 0.278 and 0.239, respectively). Both SAMOVA and AMOVA demonstrated that a substantial portion of the mtDNA genetic variability was found among groups, whereas the differences among sites within groups accounted for less variation. Variation within sites, on the other hand, accounted for the major part of the observed variation.
Table 3
Results of AMOVA and SAMOVA (in italic type) are indicated for each of the number of groups (K) into which the black bear sampling sites were pooled. Results include degrees of freedom (d.f.), percentage of variance, and fixation indexes. Significant (P < 0.05) values are indicated with an asterisk (*). All P-values were P = 0.000.
The results of the Mantel test showed that the genetic differentiation between sampling sites across Ontario could be explained partly by isolation by distance; the correlation between geographic and genetic distances was significant (r = 0.315, P = 0.002). Isolation by distance was supported more strongly when the BPNP samples were removed from the analysis (r = 0.347, P = 0.002); however, significance decreased when both BPNP and the Northwest samples were removed (r = 0.287, P = 0.005) and was absent when only the Northwest samples were removed (r = 0.09, P = 0.256).
Phylogenetic analyses
The median joining network of the haplotypes found in Ontario showed that the most frequent haplotypes identified were all located in the trunk of the network (Fig. 3). The median joining network of the sequences found across the continent, including the samples we obtained in Ontario, showed 2 genetically distinct groups, 1 largely restricted to the Pacific Northwest region and highly divergent from the other North American haplotypes (10 mutational steps), and a 2nd encompassing the rest of the continent, which corresponded to the 2 clades identified by Wooding and Ward (1997). Within the widespread continental clade we found a geographical distinction between a subclade running along the Eastern Seaboard of North America and another found in western Canadian provinces and American states (Fig. 1). This intraclade divergence was detected because 2 haplotypes that were almost exclusively restricted to the northwestern cluster of Ontario (HAP15 and HAP24) also were found in other western Canadian provinces and American states but not anywhere else along the eastern side of the continent (Fig. 3).
Discussion
Although many studies investigate genetic structuring patterns of fragmented populations to better inform conservation and management initiatives, only a few focus on genetic variation across homogenous landscapes. Such studies are useful because they add context about the state of fragmented populations by showing how genetic variation is distributed in the absence of ecological or anthropogenic disturbance, and they help make inferences about how fast continuously distributed populations can be subjected to extirpation in case of isolation. Our study used an extensive data set to describe black bear mtDNA genetic structure across a presumed continuous landscape. Relative to black bear microsatellite data also obtained in Ontario, which detected FST values illustrating a weak structure for this marker (pairwise FST < 0.02 between the nonisolated geographic clusters—Mills 2005), the values detected by our mtDNA analyses revealed a structure that was defined more strongly for this type of marker. This discrepancy in the levels of structuring detected with mtDNA and microsatellite suggests a male-biased dispersal pattern combined with female philopatry, as seen in other species (Chappell et al. 2004; Johnson et al. 2003; Tomasik and Cook 2005). In addition, integrating the Ontario haplotypes into a network that included other North American sequences showed that historical remnants of phylogeographic isolation were still observed in black bears at restricted spatial scales and were maintained most likely by sex-biased dispersal.
Our mtDNA analyses detected differences in haplotypic composition among Ontario regions and a geographic restriction of haplotypes. The most frequent haplotype found in the Northwest (HAP15) was not shared with any other clusters, and the dominant haplotypes in the Central (HAP5 and HAP6) and Southeast (HAP1 and HAP2) clusters were seldom found in northwestern Ontario (Fig. 2). HAP24, which was close to HAP15 on the network, but less frequent, also was restricted to the northwestern cluster in Ontario. In addition, of the 36 haplotypes found in Ontario, 8 were found in only 1 region, even when we excluded the singletons, illustrating a high proportion of private haplotypes in the province.
The SAMOVA results did not reveal genetic structuring that corresponded completely to the a priori defined geographic clusters. The SAMOVA provided plausible results at K = 2, K = 4, and K = 11, for which the variance among groups was maximized, but the variance among populations within groups was minimized. At K = 2, the SAMOVA grouped only 3 of the 8 Northwest sampled sites together. This was not expected, but could be explained by the fact that these sites are the ones that have the lowest HAP5 frequency compared to the rest of the Northwest populations. At K = 4, 2 additional Northwest sites were grouped together (Kenora/Fort Frances B), both of which had a high HAP15 frequency and a higher HAP5 frequency than the 1st group of localities that was pulled from this cluster. In addition, a site from the Southeast cluster, Timmins, formed a group on its own, which was explained by the fact that the highest HAP18 frequency (36%) was detected at this location. At K = 11, all the sampled sites that were grouped a priori into the defined Central cluster stayed together, and most of the sites grouped into the Southeast and Northwest clusters also grouped together. In addition to Timmins and BPNP (Southeast), the Northwest sites that had the lowest HAP15 frequencies remained separate (Fort Frances C, Sioux Lookout, and Red Lake; Table 3). These differences between the grouping patterns likely illustrate genetic structuring occurring at different geographic scales. At the largest scale (K = 2), sampled sites from the Northwest cluster of Ontario are separating from the others, and at the smallest scale (K = 11), further divisions appear within clusters. Despite this substructuring pattern, the clusters boundaries that were defined a priori are still present overall.
In addition to the SAMOVA results, the separation of the Northwest sites from the rest of Ontario was supported by the Mantel test, whose significance decreased when these populations were excluded. This suggests that the high differentiation between the Northwest and the other clusters might have skewed the results toward supporting isolation by distance across our sampling study. Because the Bruce bears also were highly differentiated from the rest of the clusters and had a level of genetic diversity that was much lower, with only 2 haplotypes detected among 38 individuals (HAP2 and HAP5), we conducted a Mantel test that excluded them in addition to the Northwest sites. With only the Central and Southeast populations, the correlation between geographic and genetic distance was detected only at the threshold of significance, which is concordant with the lower FST values observed between these clusters.
Previous findings that expressed a conservation concern for the Bruce black bears due to their geographic isolation from the rest of the Ontario individuals (Howe et al. 2008) were supported by the genetic results of the present study. Bruce bears were highly differentiated from the others and had a lower haplotypic diversity, a pattern that could be explained by isolation by fire events during the last 150 years (M. E. Obbard, pers. obs.). Because the 2 haplotypes found within the Bruce cluster are common in the rest of the province, and because no unique genetic haplotypes were found in the Bruce black bears, we conclude that they do not form an evolutionary unit. However, the combination of low genetic diversity and strong differences in haplotypic frequencies compared to the other Ontario black bear clusters suggests that the Bruce bears could be defined as a Distinct Management Unit (Committee on the Status of Endangered Wildlife in Canada [COSEWIC] 2005).
In addition to these different levels of structuring, contrasting levels of differentiation were detected between microsatellite (Mills 2005) and mitochondrial analyses. Excluding the Bruce, FST values based on microsatellites illustrated subtle levels of genetic structure (0.008 < FST < 0.140), whereas values based on mtDNA illustrated stronger levels of differentiation (0.120 < FST < 0.419), even at geographic distances as small as 150 km. For example, although the Northwest cluster of Ontario was not isolated geographically from the other clusters (Fig. 2), pairwise FST values showed that it was strongly differentiated from them. Such results, combined with the absence of topographic barrier to dispersal across the sampling area, have been explained by low effective population sizes, low dispersal distances, long-term isolation of lineages, cryptic barriers (Irwin 2002), or sex-biased dispersal (Tomasik and Cook 2005). However, male black bears are known for their long-distance dispersal capabilities (Rogers 1987), and total abundance of black bears is reasonably high in Ontario (approximately 100,000 individuals—M. E. Obbard, pers. obs.), suggesting that black bears may be at equilibrium, and hence we would expect genetic drift to have little impact on them compared to natural selection. Because our results do not support a panmictic structure in Ontario black bears and show discrepancies in differentiation levels between genetic markers, the most likely explanation is a combination of male-mediated gene flow and female natal philopatry, which supports our prediction and previous studies that detected those patterns in black bears (Costello et al. 2008; Onorato et al. 2007; Rogers 1987).
In black bears and other taxa found in North America, such as northern flying squirrels (Glaucomys sabrinus), red foxes (Vulpes vulpes), long-tailed voles (Microtus longicaudus), American pine martens (M. americana), and yellow-pine chipmunks (Tamias amoenus), 2 main historical lineages were identified, a continental one and a coastal one (Arbogast 1999; Aubry et al. 2009; Byun et al. 1997; Conroy and Cook 2000; Demboski et al. 1999; Demboski and Sullivan 2003; Wooding and Ward 1997). Their origin has been suggested to derive from several isolated refugia during the last glacial maximum along the coasts of the North Pacific and East Atlantic; however, the exact locations of these refugia remain unclear. In black bears the continental lineage extends from Alaska southward to New Mexico and eastward to Newfoundland and Florida, and the coastal one extends from Alaska to California and also occurs in British Columbia, Alberta, and Montana (Byun et al. 1997; Peacock et al. 2007; Stone and Cook 2000; Wooding and Ward 1997). In addition to this continental–coastal divergence, Wooding and Ward (1997) found a low east–west genetic differentiation within the continental clade and suggested that the disjunct distribution of these 2 potential subclades was due to a lack of samples from the central part of North America (e.g., no samples from Ontario, Manitoba, or Michigan). This intraclade subdivision is subtle (Wooding and Ward 1997), and it cannot be explained by a prominent physiographic factor such as isolated glacial refugia.
Our samples, collected on a 1,700-km continuum across Ontario, allowed us to fill this sampling gap that existed in the mideastern portion of the black bears' range and subsequently put our results into a broader continental context. Our 2nd network including these sequences from Ontario and sequences from the rest of the North American continent showed that HAP15 and HAP24 (both mostly restricted to the Northwest cluster of Ontario) were restricted to the mideastern to western part of the black bear's range (Alberta, Alaska, British Columbia, Manitoba, Montana, New Mexico, Ontario, and Utah). In contrast, HAP1 and HAP2 (both mostly restricted to the southeastern portion of Ontario) were restricted to the eastern part of their range (Florida, Louisiana, Mexico–Texas, New Brunswick, and Ontario; Fig. 3). In addition, the Ontario Northwest cluster was strongly differentiated from the Central cluster. Because of this geographic restriction of haplotypes, and this high level of genetic differentiation detected at a very small scale in Ontario, we infer that the black bears located on the western (Northwest cluster of Ontario) and eastern (Central and Southeast clusters of Ontario) sides of the province belong to the western and eastern North American continental subclades, respectively. This clade subdivison also was supported by the significant partial Mantel test that modeled a barrier to gene flow between the Northwest and the rest of the Ontario populations (partial r = 0.255, P = 0.007). This pattern of geographic distribution of genetic types shows that the disjunct distribution previously identified by Wooding and Ward (1997) was not due only to a lack of sampling in the central part of the North America, because our results still identify the eastern–western subdivision of the continental clade (Fig. 3) at a very small geographic scale.
For black bears located in the southwestern region of North America, barriers to gene flow were suggested to be driving this type of differentiation (Onorato et al. 2004). The presence of a physiogeographic barrier represented by the Chihuahuan Desert, which restricts gene flow between the sites of Mexico–Texas and New Mexico, could have helped maintain a high level of differentiation between black bear populations (Onorato et al. 2007). However, the absence of topographic barriers to long-distance dispersal on the eastern side of the continent, and at a smaller scale, between the differentiated Central and Northwest clusters of Ontario, seems to rule out a structure linked to long-lasting landscape features. This finding supports the results from Peacock et al. (2007), who suggested that clusters are not necessarily defined by physical barriers. Given the evolutionary rate of mtDNA (∼10−6 substitutions per site per year—Brown et al. 1979), historical factors likely are driving such a differentiation pattern. This continental clade subdivision likely has occurred over a much more restricted length of time than the coastal/continental clade division, because it is not strong enough to suggest isolated glacial refugia on the eastern side of the Rockies. Rather, the shape of our network, with a few ancestral haplotypes (HAP1 and HAP15) having many recent derivatives, suggests range expansion (Avise 2000). Because this east–west subdivision of the continental clade seems to follow the pattern of the retreat of the last ice sheet (Adams and Faure 1997), we suggest that it exists because after departing from an ancestral population located in the main continental refugium during the late Pleistocene, black bears followed 2 opposite recolonization routes on either side of the receding ice sheet. Due to rapid geographic expansion following the melting ice, the 2 subclades met in northern Ontario. The habitat at the contact zone between the 2 subclades likely became saturated, inhibiting future female migration. Thus, the historical genetic structure that arose during the postglacial recolonization of North America, which was 1st due to isolation by distance after the postglacial range expansion, could have been maintained subsequently by female philopatry and male-biased dispersal, resulting in the observed contemporary clusters. That the Mexico–Texas haplotypes are closely related to haplotypes from the eastern North American subclade, whereas those from New Mexico are more closely related to sequences from the western subclade, further confirms our proposition of 2 recolonization routes. It also supports the 2nd long-distance colonization hypothesis proposed by Onorato et al. (2004) suggesting that dispersal of black bears from the eastern United States leads to their current distribution in the Mexico–Texas region.
Our sampling across Ontario allowed us to detect the continental clade subdivision at a small geographic scale (150 km between the 2 subclades). In studies for which samples were collected at longer distance intervals this differentiation was observed at a more intermediate scale, even with markers that have a higher rate of evolution than mtDNA, such as microsatellites (lynx, FST = 0.0622, P = 0.01 [Rueness et al. 2003]; and piping plovers [Charadrius melodus], FST = 0.473, P < 0.000 [Miller et al. 2009]). In these studies it was suggested that this differentiation could be caused by contemporary rather than historical factors, and the structuring patterns were influenced by climate variations through habitat and breeding-site choice. In lynx the cryptic division was suggested to be due to opposite effects of the North Atlantic Oscillation on the snow conditions of different climatic regions, which would affect hunting abilities and habitat choice of lynx, because individuals would stay on 1 specific side of this North Atlantic Oscillation line because of habitat familiarity, which would lead to genetic structuring despite the absence of barriers to gene flow (Stenseth et al. 1999, 2004). For the piping plover (Miller et al. 2009) the genetic structure was explained by differential levels of breeding-site fidelity due to opposite flooding conditions in the neighboring regions. The location of this North Atlantic Oscillation line, which marks the division between the Continental and Atlantic climatic regions (Stenseth et al. 1999), corresponds to where we identified the cryptic genetic subdivision of the wide continental clade in black bears (Fig. 1). We cannot envision how differential climatic conditions could maintain such small-scale differentiation for black bears, but these findings in other species warrant further investigation, at least to verify if males disperse more likely within clusters, as opposed to between them.
In addition to future research aspects, we suggest that future conservation and management decisions for large carnivores, and especially ones that are known to have differential male and female dispersal patterns, are made based on genetic information that uses both microsatellites and mtDNA. As shown here, microsatellites are not fully informative when historical lineages are maintained contemporarily by dispersal patterns, and management decisions solely based on microsatellites can lead to changes in the genetic composition of populations. In Arkansas, for example, the genetic composition of populations that belonged historically to the Continental Eastern subclade changed into a Continental Western subclade type after they received translocated individuals from Manitoba and Minnesota (Van Den Bussche et al. 2009). The mitochondrial genome also has highly functional fragments (Ballard and Whitlock 2004; Rutledge et al. 2010) in addition to the neutral control region. If the variation of neutral fragments reflects that of functional fragments, not accounting for variation in mtDNA could negatively impact management actions that focus on recovery of populations.
To complement this study we suggest gathering more data from potential secondary contact zones between the 2 continental subclades and examining functional markers to look for possible local adaptive responses. The field of ecological genomics, for example, would allow us to identify the genes that are involved in the various responses to differential environmental conditions. Future studies of North American forest species, whose distribution is similar to that of black bears, should focus on explaining the small-scale intralineage diversification on the eastern side of the continent, because it could lead to new findings on the influence of both contemporary and historical forces on the dynamics of species diversification. At the local scale we showed that the genetic structure of Ontario black bears reflects their historical differentiation levels in the absence of barriers to gene flow. Such information can be used as a baseline to quantify the amount of disturbance in the current isolated North American populations of black bears.
Acknowledgments
We thank K. Mills for providing the microsatellite results, E. Howe and many employees from district offices of the Ontario Ministry of Natural Resources across the province for assistance with sampling, and S. Coulson and K. Wozney for laboratory assistance. Funding was provided by the Ontario Ministry of Natural Resources. Finally, we thank C. Wilson, J. Zigouris, L. Finnegan, and 2 anonymous reviewers and the associate editor, whose comments and suggestions on previous drafts greatly improved the quality of this manuscript.
Literature Cited
Appendices
Appendix I
Distribution of the 36 mitochondrial DNA haplotypes at 23 sampling sites of American black bears (Ursus americanus) across Ontario, and measures of their absolute and relative frequencies. Sampling sites are mapped in Fig. 2. No. HAP = number of haplotypes; BPNP = Bruce Peninsula National Park; Bor + Iv = Borland + Ivanhoe; CCGP = Chapleau Crown Game Preserve.
Appendix II
Locations, references, and GenBank accession numbers of all the haplotypes used in this study. AB = Alberta; AK = Alaska; AR = Arkansas; BC = British Columbia; CA = California; FL = Florida; LA = Louisiana; MB = Manitoba; MN = Minnesota; MT = Montana; NB = New Brunswick; NL = Newfoundland; NM = New Mexico; OK = Oklahoma; ON = Ontario; QB = Quebec; TX = Texas, UT = Utah.
Appendix III
Measures of neutrality (Tajima's and Fu's tests), nucleotide (π) and haplotype (h) diversity and their standard deviations (SD π and SD h), and sampling site–specific FST for 23 sampling sites of American black bears (Ursus americanus) across Ontario. Sampling sites are mapped in Figs. 1 and 2. BPNP = Bruce Peninsula National Park; Bor + Iv = Borland + Ivanhoe; CCGP = Chapleau Crown Game Preserve.
