Spatiotemporal analyses suggest the role of glacial history and the ice‐free corridor in shaping American badger population genetic variation

Abstract Recurring glacial cycles through the Quaternary period drastically altered the size and distribution of natural populations of North American flora and fauna. The “southerly refugia model” has been the longstanding framework for testing the effects of glaciation on contemporary genetic patterns; however, insights from ancient DNA have contributed to the reconstruction of more complex histories for some species. The American badger, Taxidea taxus, provides an interesting species for exploring the genetic legacy of glacial history, having been hypothesized to have postglacially emerged from a single, southerly refugium to recolonize northern latitudes. However, previous studies have lacked genetic sampling from areas where distinct glacial refugia have been hypothesized, including the Pacific Northwest and American Far North (Yukon, Alaska). In order to further investigate the phylogeographic history of American badgers, we collected mitochondrial DNA sequence data from ancient subfossil material collected within the historical range (Alaska, Yukon) and combined them with new and previously published data from across the species' contemporary distribution (n = 1,207). We reconstructed a mostly unresolved phylogenetic tree and star‐like haplotype network indicative of emergence from a largely panmictic glacial refugium and recent population expansion, the latter further punctuated by significantly negative Tajima's D and Fu's Fs values. Although directionality of migration cannot be unequivocally inferred, the moderate to high levels of genetic variation exhibited by American badgers, alongside the low frequency of haplotypes with indels in the Midwest, suggest a potential recolonization into central North America after the hypothesized ice‐free corridor reopened ~13,000 years ago. Overall, the expanded reconstruction of phylogeographic history of American badgers offers a broader understanding of contemporary range‐wide patterns and identifies unique genetic units that can likely be used to inform conservation of at‐risk populations at the northern periphery.


| INTRODUC TI ON
Understanding how historical processes shape the genetic variation of natural populations has long been of interest to ecologists and evolutionary biologists. Glaciation, in particular, has been extensively studied in the context of population genetic structure, with numerous studies finding concordance between patterns of genetic variation and the historical locations of ice sheets (Hewitt, 1999(Hewitt, , 2000(Hewitt, , 2004. Early studies, based on fossil evidence, suggested that natural populations survived the harsh environments imposed by glaciation by residing in refugia that constituted geographic regions hospitable for flora and fauna, primarily south of glacial extents (Petit et al., 2003). Indeed, genetic data have corroborated this hypothesis for numerous species, with high levels of genetic diversity at lower latitudes decreasing in a clinal fashion toward northern latitudes, coinciding with the recolonization events that followed glacial retreat (Hewitt, 2004).
This "southerly refugia model" (Bennett, Tzedakis, & Willis, 1991) has been the longstanding framework for testing the effects of glaciation on contemporary genetic patterns.
Recently, studies of greater breadth and depth, as well as advances in ancient DNA technologies, have depicted more detailed and oftentimes complex histories for North American populations (Shafer, Cullingham, Côté, & Coltman, 2010;Soltis, Morris, McLachlan, Manos, & Soltis, 2006). Cryptic refugia, which were semihospitable environments within ice sheets (Provan & Bennett, 2008), and a northern refugium in ancient Beringia (Tremblay & Schoen, 1999) have both been suggested to account for the reconstructed patterns of North American populations (Shafer et al., 2010;Soltis et al., 2006). Indeed, genetic studies have identified "hot spots" of genetic diversity and discrete genetic variation for populations at the northern extent of species' ranges, suggesting a refugial history (Rowe, Heske, Brown, & Paige, 2004;Shafer, Côté, & Coltman, 2011). However, patterns of genetic diversity indicating potential glacial refugia can be obscured by complex scenarios, such as admixture between lineages from separate refugia (Petit et al., 2003), genetic structure within refugia (i.e., refugia within refugia; Gómez & Lunt, 2007), and deeper historical associations (Lovette & Bermingham, 1999). Identifying the locations of glacial refugia and accounting for these alternative scenarios can also directly inform conservation management by identifying lineages of distinct ecological and evolutionary significance (Bhagwat & Willis, 2008;Hampe & Petit, 2005).
One interesting species to explore the genetic legacy of glacial history is the American badger (Taxidea taxus). As a semifossorial mustelid, American badgers are typically restricted to dry grasslandshrubland habitats, where they can burrow in silty, sandy-loam soils and catch fossorial prey (Committee on the Status of Endangered Wildlife in Canada, 2012). Despite its habitat specialization, the American badger's high dispersal capabilities have enabled colonization to diverse areas throughout much of central and western North America (Committee on the Status of Endangered Wildlife in Canada, 2012;Kierepka & Latch, 2016). Thus, while habitat specialization suggests that badgers should have resided in one or multiple refugia, their dispersal capacity and extensive species' range strongly suggest the potential for the latter (Shafer et al., 2010). It has been hypothesized that glaciation displaced badgers into a single, southerly refugium from which the species recolonized northern latitudes, with isolation by major geographic barriers leading to the recognition of four subspecies (T. t. jacksoni, T. t. taxus, T. t. berlandieri, and T. t. jeffersonii;Kierepka & Latch, 2016;Long, 1972). Phylogeographic studies have partially supported this hypothesis, providing evidence for fragmented units with low genetic diversity at the northern periphery (Ethier, Lafleche, Swanson, Nocera, & Kyle, 2012), and connected populations with high levels of gene flow at the center of the species' range (Kierepka & Latch, 2016). However, a lack of genetic sampling from areas where distinct glacial refugia have been hypothesized (e.g., the Pacific Northwest and Alaska, the latter a location where badgers resided historically) has hindered a complete picture of how postglacial colonization has affected contemporary genetic patterns. Understanding the relative impact of postglacial colonization is essential, especially for the management of endan- Here, we aimed to gain a broader perspective of the historical influences that shaped the genetic variation of American badger populations across the entire species' range. To do so, we collected mitochondrial DNA sequence data from ancient subfossil material collected within the historical range of American badgers in Alaska (USA) and the Yukon (Canada), and paired these with new data from the contemporary range in the American West and Pacific Northwest where coverage has been previously thin or lacking. We combined these data with those previously published from across the species distribution (Ethier et al., 2012;Ford, Weir, Lewis, Larsen, & Russello, 2019;Kierepka & Latch, 2016), producing the most comprehensive dataset to date, covering the American badger's historical and contemporary ranges. We used these data to reconstruct the phylogeographic history of the American badger to specifically test alternative hypotheses of (1) a single refugium, which would be generally evidenced by a star-like topology showing little differentiation among a few common haplotypes surrounded by rare haplotypes, or (2) multiple refugia, which would be partially supported by well-defined geographic groups with highly divergent haplotypes.

K E Y W O R D S
ancient DNA, carnivore, peripheral populations, phylogeography, Pleistocene glaciation, Taxidea taxus 2 | ME THODS

| Sample collection
We sampled ancient subfossils (n = 13) with estimated ages dating to the late Wisconsin (>12,000 ybp) with provenances in Alaska, USA, and Yukon, Canada, from the American Museum of Natural History and the Canadian Museum of Nature (Figure 1; Appendix S1; Table S1).
Methods to prepare subfossil bones for DNA extraction followed those described in Lippold et al. (2011). Briefly, a Dremel ® tool fitted with a cutoff wheel was used to remove surface contaminants from bones. A 0.25 g sample was then removed from each specimen, taking caution not to disturb bone processes used for species identification.
Samples were covered with aluminum foil and pulverized with a mortar and pestle before DNA extraction. All surfaces and equipment were thoroughly cleaned with bleach solutions between each sample. Columbia (n = 31), Washington (n = 20), Idaho (n = 9), and Montana (n = 14) (Appendix S1). Samples were collected following methods described in Casas-Marce, Revilla, and Godoy (2010). Briefly, a Dremel ® bit tool (1.5-2mm diameter) instrument was used to drill powder from the base of a single foreclaw. Surfaces and tools were cleaned with bleach between collection of each sample. Powder samples were stored dry in 2 ml centrifuge tubes until DNA extraction.

| Data collection
Genomic DNA was extracted from contemporary samples using the Nucleospin Tissue Kit (Macherey-Nagel). DNA extractions from historical samples were conducted using a MinElute Kit (Qiagen) and the protocol in Jensen et al. (2018) within the dedicated Historical DNA Analysis Facility at the University of British Columbia Okanagan, Canada. Genomic DNA was extracted from subfossil specimens using a phenol-chloroform method (Losey et al., 2013) in the ancient DNA laboratory at the Estación Biológica de Doñana, Spain. All extractions from historical and ancient specimens were conducted in batches of 16 tubes including three or four negatives that were carried through all steps, including PCR amplification, when an additional two negatives were added per batch. Different species were included in ancient DNA extraction batches such that two samples of the same species were never in adjacent tubes.
Haplotypic data were collected from a fragment of the D-loop from the mitochondrial genome so as to make use of the expansive range-wide data from several previously published studies (Table 1; Ethier et al., 2012;Ford et al., 2019;Kierepka & Latch, 2016). We used the primer set designed in

| Phylogenetic analysis
All newly acquired and previously published sequences from Ford et al. (2019) and Ethier et al. (2012) were trimmed at the 5′ end to align with sequences from Kierepka and Latch (2016), and 3′ terminal gaps were coded as missing data. All sequences were aligned using default parameters in MUSCLE TA B L E 1 Sample distribution and diversity indices based on a fragment of the mtDNA D-loop for American badgers across the species' range

Number of Haplotypes
To reconstruct relationships between mtDNA haplotypes, Bayesian phylogenetic analysis was conducted using MrBayes 3.2 (Ronquist & Huelsenbeck, 2003). A previously identified 26-bp indel was coded as a fifth character state following Kierepka and Latch (2016). The best fit model of nucleotide substitution (GTR + G + I) was determined using jmodeltest2 (Darriba, Taboada, Doallo, & Posada, 2012) and the Bayesian information criterion (BIC). Overlapping D-loop sequence from the previously published European badger (Meles meles) mitochondrial genome was used as an outgroup. The analysis ran four simultaneous chains for 2.0 × 10 6 total generations, each using a random tree as a starting point, the default heating scheme, and saving a tree every 200 generations for a total 10,000 trees. The first 1,000 trees were discarded as burn-in samples and the remaining 9,000 trees were used to construct a consensus tree and derive posterior probability values.
The dissimilarity matrix from the TCS output was used to visualize the haplotype network using Hapstar 0.7 (Teacher & Griffiths, 2011).

| Haplotypic variation and geographic structure
To were used for the geographic coordinates of sampling units following Kierepka and Latch (2016). Estimates of haplotype diversity (H d ), nucleotide diversity (π), and pairwise difference were calculated for each sampling unit in Arlequin 3.5 (Excoffier, Laval, & Schneider, 2005 To measure differentiation between sampling units across the species' distribution, we conducted a spatial analysis of molecular variance (SAMOVA 1.0; Dupanloup, Schneider, & Excoffier, 2002).
SAMOVA defines groups of populations by identifying the configuration that maximizes the amount of genetic variance explained among differentiated groupings (Dupanloup et al., 2002). The optimal configuration of groups was tested from K = 2 to K = 15, using 100 simulated annealing steps, as suggested in the manual (Dupanloup et al., 2002). We chose the optimal configuration based on the grouping that maximized the total amount of genetic variance between groups (ϕ CT ).

| Demographic history
To test for evidence of population expansion, we used two neutrality tests, Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997), for each of the SAMOVA groupings and for the entire contemporary dataset.
Tajima's D and Fu's Fs are expected to be significantly negative for populations that have undergone a recent expansion. We tested for significance of the two summary statistics using 1,000 random permutations implemented in Arlequin 3.5 (Excoffier et al., 2005).
To infer changes in effective population size over time, a Bayesian skyline plot analysis (Drummond, Rambaut, Shapiro, & Pybus, 2005) was implemented in BEAST 2.5.2 (Bouckaert et al., 2014) using the best fit model of nucleotide substitution (GTR + G + I) as described above. We ran three chains for 1 × 10 7 Markov chain Monte Carlo   Table 1;   Table S1).

| Haplotypic variation and geographic structure
Variation in haplotype diversity was large, ranging from 0.000 ± 0.000 in BC-Cariboo to 1.00 ± 0.096 in Colorado (Table 1).
Haplotype diversity was significantly associated with distance to center, decreasing as sampling units were closer to the periphery (R 2 = .44; p < .001; Figure 4). Haplotype diversity generally decreased with an increase in latitude, although not significantly (R 2 = .10; p = .08, data not shown).
Haplotypes recovered from ancient badgers from Alaska and the Yukon were exact matches to haplotypes in the Canadian "Prairie  (Table 2), with 39.26% of variation explained among the two groups. The percent of variation among groups and ϕ CT gradually decreased as the number of groups considered (K) increased. The next configuration that explained the most variation between groups further separated BC-Cariboo and BC-Nicola units into their own groups at K = 3 (Table 2).
Consistent with results from Kierepka and Latch (2016), we found that the Lower Peninsula, Ohio, and Ontario samples grouped together at K = 5 relative to all other sampling units in North America, besides those identified in British Columbia (Table 2).

| Demographic history
We found significant evidence for population expansion using Tajima's D and Fu's F neutrality statistics consistent with a previous hypothesis of a large, panmictic refugium in the center of the species' range (Kierepka & Latch, 2016;Long, 1972); however, no significant association was detected for the differentiated units (SAMOVA Group 1:

BC-Cariboo and BC-Nicola) at the northwestern periphery in British
Columbia (Table 3). The timing of the demographic expansion for the majority of the species' distribution appears to have occurred following the last glacial maximum, approximately 20,000 years ago, based on our Bayesian skyline plot analysis ( Figure 6).

| D ISCUSS I ON
The American badger has had a long history in North America, with the earliest fossil specimens dating back to the early Pliocene, ~3.5 million years ago (Long, 1972). American badgers are thought to have survived the Pleistocene ice ages by residing in a single refugium, south of glacial extents, from which they recolonized northern latitudes following glacial retreat (Kierepka & Latch, 2016). By expanding on previous studies with newly acquired genetic data, we shed light on more complex patterns across the species' range. Ontario, and Ohio. Moreover, haplotype diversity generally decreased with increasing peripherality. In addition, the largely unresolved tree and star-like haplotype network both displayed little discrete geographic structuring, indicative of a largely panmictic glacial refugium and recent population expansion (Figure 4) (Slatkin & Hudson, 1991). Furthermore, coalescent-based analyses suggested a population expansion following the last glacial maximum. These patterns are consistent with another highly mobile generalist carnivore, the coyote (Canis latrans), which also lacks distinct phylogeographic structure and faced population size changes coinciding with historic changes in climate (Koblmüller et al., 2016).
On a more regional level, we expected to see shared or connected haplotypes from ancient populations in eastern Beringia with those from British Columbia, Canada, based on previous studies of North American postglacial colonization in other species (Shafer et al., 2010). Contrary to these expectations, we discovered that haplotypes recovered from Alaska and the Yukon were identical to haplotypes found in the Prairie provinces, the American Midwest, and the American West. Other species' populations have displayed this connectivity (Heintzman et al., 2016), owing to a hypothesized ice-free corridor between the Cordilleran and Laurentide ice sheets through which animals migrated during the last glacial maximum, just east of present-day British Columbia (Stalker, 1977). Recent evidence, however, suggests that the Cordilleran and Laurentide glaciers coalesced between 23,000 and 13,000 years ago, completely blocking any migration during this time period (Heintzman et al., 2016). Therefore, the shared haplotypes between ancient and contemporary populations suggest one of three possibilities: (1) badgers migrated to eastern Beringia from the grasslands of central North America when the ice-free corridor reopened 13,000 years ago; (2) badgers from central North America migrated north more than 23,000 years ago prior to the closing of the ice-free corridor; or (3) badgers from eastern Beringia colonized central North America after the corridor reopened 13,000 years ago, with badgers being extirpated from eastern Beringia following each scenario.
Unfortunately, we did not have radiocarbon dates for the ancient specimens from which full sequences were recovered, but they are suggested to be late Wisconsin in age, with radiocarbon dates for other species ranging between 12,600 and 59,000 ybp at the Alaska site (Anderson, 1977), and between 20,900 and 39,000 ybp at the Yukon site (Fox-Dobbs, Leonard, & Koch, 2008). Furthermore, radiocarbon dating for unsampled Yukon badger specimens date to 37,930 and 15,190 years ago (Harington & Clulow, 1973). These lines of evidence suggest that scenarios 2 or 3 are more plausible, as badgers were present in eastern Beringia before the corridor reopened.

ACK N OWLED G M ENTS
We would like to thank the Royal British Columbia Museum,  , NI | TH, OK, WA, OR, EK, AL, ID, MT, SK,  MB, WY, UT, CO, NM, ND, SD, NE, OL, KS,  MN, IA, MO, WI, IL, IN, ON,  Note: Acronyms as in Table 1.
TA B L E 2 SAMOVA groupings for K = 2 to K = 10 TA B L E 3 Neutrality statistics for a fragment of the mtDNA D-loop for the entire dataset and the two SAMOVA groupings at K = 2