Phylogeography of a widespread small carnivore, the western spotted skunk (Spilogale gracilis) reveals temporally variable signatures of isolation across western North America

Abstract We analyzed phylogeographic patterns in the western spotted skunk, Spilogale gracilis Merriam, 1890 (Carnivora: Mephitidae) in relation to historical events associated with Pre‐Pleistocene Divergence (PPD) and Quaternary climate change (QCC) using mitochondrial DNA from 97 individuals distributed across Western North America. Divergence times were generated using BEAST to estimate when isolation in putative refugia occurred. Patterns and timing of demographic expansion was performed using Bayesian skyline plot. Putative climatic refugia resulting from Quaternary climate change were identified using paleoecological niche modeling and divergence dates compared to major vicariant events associated with Pre‐Pleistocene conditions. We recovered three major mitochondrial clades corresponding to western North America (California, Baja, and across the Great Basin), east‐central North America (Texas, central Mexico, New Mexico), and southwestern Arizona/northwestern Mexico. Time to most recent common ancestor for S. gracilis occurred ~1.36 Ma. Divergence times for each major clade occurred between 0.25 and 0.12 Ma, with signature of population expansion occurring 0.15 and 0.10 Ma. Ecological niche models identified three potential climatic refugia during the Last Interglacial, (1) west coast of California and Oregon, (2) northwestern Mexico, and (3) southern Texas/northeastern Mexico as well as two refugia during the Last Glacial Maximum, (1) western USA and (2) southern Texas/northeastern Mexico. This study supports PPD in shaping species‐level diversity compared to QCC‐driven changes at the intraspecific level for Spilogale, similar to the patterns reported for other small mammals (e.g., rodents and bats). Phylogeographic patterns also appear to have been shaped by both habitat and river vicariance, especially across the desert southwest. Further, continuing climate change during the Holocene coupled with anthropogenic modifications during the Anthropocene appears to be removing both of these barriers to current dispersal of western spotted skunks.


| INTRODUCTION
Recent phylogeographic studies of a diverse array of North American plants and animals have resulted in identification of broadscale patterns of genetic variation attributable to both Pre-Pleistocene Divergence (PPD) resulting from vicariance and Quaternary climate change (QCC) associated with expansion and contraction of glaciers during the Pleistocene-Holocene (Hewitt, 2004). Although the relative roles of each of these forces in shaping patterns of speciation have been debated, with some arguing QCC has little if anything to do with species-level evolution (Bennett, 2004) while others support lineage-specific impacts of QCC on vertebrate speciation (Peterson & Ammann, 2013), newer methodologies including ecological niche modeling (ENM) and molecular dating of divergence times has revolutionized our ability to characterize species-specific responses to both PPD and QCC (Alvarado-Serrano & Knowles, 2014;Hickerson et al., 2010). For western North America, a region of intense phylogeographic study (Swenson & Howard, 2005), PPD tends to be characterized by vicariant, tectonic events of the Neogene such as mountain orogensis or sea level inundation, although species currently inhabiting the region display mixed evolutionary responses to such barriers (Graham, Hendrixson, Hamilton, & Bond, 2015). In contrast, the impact of QCC on species evolution is characterized by refugia creation via range fragmentation followed by subsequent expansion in response to changing climates, often resulting in regionally defined responses by animals adapted to particular habitats (e.g., rodents of the warm deserts, Riddle & Hafner, 2006). Additional lines of evidence continue to refine locations of such Pleistocene refugia, enhancing our ability to interpret phylogeographic patterns for wellstudied regions (Holmgren et al., 2014). Thus, species-level studies continue to add value to investigations of shared phylogeographic patterns related to PPD and QCC across broad geographic areas, even for taxa with previously published phylogeographic studies (Bryson, Jaeger, Lemos-Espinal, & Lazcano, 2012;Jaeger, Riddle, & Bradford, 2005).

Phylogeographic patterns and processes across western North
America have been well-characterized for three regions in particular, (1) northwestern North America (as defined by Shafer, Cullingham, Cote, & Coltman, 2010), (2) the warm deserts (i.e., Chihuahuan, Mojave, and Sonoran Deserts as defined by Riddle & Hafner, 2006), and (3) the Great Basin Desert (as defined by Riddle, Jezkova, Hornsby, & Matocq, 2014). For northwestern North America, animals were originally thought to have been isolated between two major refugia during glacial maximum, Berengia and the Pacific Northwest (Pielou, 1991). However, a more complex phylogeographic pattern has been proposed, indicating, for example, that both PPD events such as the orogeny of the Cascade-Sierra mountain chain 5-2 Ma and QCCdriven refugia have played a role in shaping genetic patterns of the regions flora and fauna (Shafer et al., 2010). Similarly for North America's warm deserts, both PPD and QCC may have contributed to patterns of genetic subdivision, with PPD-associated late Miocene-Pliocene vicariant events such as the formation of the Bouse Embayment potentially leading to speciation via allopatry followed by subsequent intraspecies-level divergence attributable to QCC events such as the climatically driven expansion and contraction of the Deming Plains (Hafner & Riddle, 2011;Riddle & Hafner, 2006). Finally, mammals and other biota in the Great Basin Desert have phylogeographic structure that appears to be first shaped by PPD and then followed by dramatic alterations to species' ranges attributable to QCC (Riddle et al., 2014).
These shed light on how wide-ranging, vertebrate taxa with broad ecological niches respond to biogeographic events associated with different historical signatures. Not surprisingly, larger, more mobile species tend to maintain weaker levels of genetic subdivisions (Latch et al., 2014) than smaller species (Mantooth et al., 2013) but see (Nava-García et al., 2016) across these regions. Patterns associated with mid-sized mammals remain obscure for much of this region, including important yet often neglected members of these ecological communities, such as small carnivores (Prugh et al., 2009;Roemer, Gompper, & Van Valkenburgh, 2009).
One notable group of carnivores that remains understudied, especially in western North America, despite their diverse ecologies and broad distributions across important biogeographical areas is the skunks (Carnivora, Mephitidae, Dragoo, 2009). One genus of skunks, Spilogale, is represented by four extant taxa (Dragoo, 2009): the western spotted skunk (S. gracilis), eastern spotted skunk (S. putorius), southern spotted skunk (S. angustifrons), and pygmy spotted skunk (S. pygmaea). Members of Spilogale have existed in North America since the early Pliocene (Van Gelder, 1959), with the oldest fossils of the extinct Rexroad spotted skunks, S. rexroadi from Kansas and Texas dating back to c. 3.0-3.5 Ma (Dalquest, 1971;Hibbard, 1941a,b). The western spotted skunk, S. gracilis, is a small (200-800 g) species currently distributed from central Mexico north to British Columbia, with an east-west distribution reaching from the California coast to the central Great Plains (Verts, Carraway, & Kinlaw, 2001).
The broad geographic distribution of S. gracilis across three wellcharacterized regions of western North America (Figure 1) indicates that the evolutionary history of this species could have been shaped by a number of PPD events, including but not limited to mountain orogeny (e.g., Rocky Mountains, Sierra Nevada, and Sierra Madre Occidental), river drainage alterations (e.g., Rio Grande, Rio Conchos, and Rio Nazas), and sea level changes (e.g., inundation along the Colorado River forming the Bouse Embayment). Characterization of these PPD events for each respective region and their potential impacts on phylogeographic patterns in a wide array of plants and animals provides useful hypotheses for testing within S. gracilis. In addition, habitats across the range of S. gracilis are characterized by changes in extent and distribution associated with fluctuating climates of the Quaternary. Thus, the catholic choice of habitat maintained by S. gracilis also makes it well-suited for examining patterns of genetic subdivision associated with QCC. This broad ecological niche combined with the estimated age and distribution of this lineage across a well-studied but complex biogeographic region makes this species well-suited for analyses of the relative influence of PPD and QCC in shaping phylogeographic patterns of small carnivores across western North America.
Using mitochondrial DNA (mtDNA) and ecological niche modeling, we explored patterns of genetic variation and divergence time estimates for S. gracilis and other Spilogale species in relation to wellcharacterized events associated with PPD as well as fluctuating climates attributable to QCC across western North America. We predict that S. gracilis will show phylogeographic patterns intermediate to those of larger, more mobile species and smaller, less mobile species, with species-level divergence occurring around well-defined phylogeographic breaks associated with PPD followed by population-level divergence associated with QCC.

| Sampling and laboratory methods
Genetic data from recent field expeditions and museum tissue collections were obtained from 97 individual S. gracilis representing 78 unique localities (see Appendix S1, Table S1.1). All genetic sequences were deposited onto GenBank together with associated voucher specimen information, when available (KY679911-KY680207).
Cycle sequencing was performed using BiGDyE version 3.1 (Applied Biosystems, Foster City, CA, USA) following the manufacturer's protocols. Sequences were electrophoresed on an ABI-3100 or 3130 (Applied Biosystems). Models of evolution for each gene were selected using mEGa 6 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013) based on BIC and AICc scores.
F I G U R E 1 Phylogeographic patterns estimated using Bayesian inference for Spilogale gracilis based on Cytb mitochondrial gene. Nodal support for the major clades depicted in the open circle (Bayesian posterior probability above and maximum-likelihood bootstrap below). Shading indicates regions with elevations >1,700 m. Numbers refer to historical barriers associated with phylogeographic subdivisions in warm desert taxa: (1) Boues Embayment/Colorado River, (2) Deming Plains, (3) Rio Grande-Rio Conchos, and (4) Sierra Barabampo-Rio Fuerte adapted from Riddle and Hafner (2006). Map of clades projected in Mollweide with a WGS84 datum. Colored clades correspond to respective mitochondrial DNA clades referred to in text as Arizona (Blue, diamonds), East-Central (Green, circles) and West (Orange, squares). Photograph of a western spotted skunk from San Angelo, Texas. Photograph Credit: Robert C. Dowler

| Phylogenetic estimates
Phylogenetic trees were estimated using Bayesian Inference with mR-BayES 3.2.2 (Ronquist et al., 2012) within the CIPRES Science Gateway (Miller, Pfeiffer, & Schwartz, 2010). Two MCMC chains were run for 10 million generations, sampling trees every 100 iterations and using a burn-in period of 10,000 trees. Base frequencies were assigned an equal Dirichlet distribution (1.0, 1.0, 1.0, 1.0). A maximum-likelihood phylogeny was estimated using phyml version 3.0 (Guindon et al., 2010). The BIONJ option was selected for the starting tree (Gascuel, 1997), and 1,000 bootstrap replicates were performed to estimate nodal support. Proportion of invariable sites and the shape of the gamma distribution were estimated within the maximum-likelihood framework. Phylogenetic trees were estimated for individual mitochondrial genes for Cytb and D-loop only and a concatenated dataset of Cytb and D-loop combined.

| Divergence time estimates
Divergence times were estimated using the three-gene dataset (Cytb,ND5, in BEaSt version 2.0 (Drummond, Suchard, Xie, & Rambaut, 2012). A likelihood ratio test for a molecular clock (Felsenstein, 1988) performed in mEGa6 was rejected, and therefore, an uncorrelated log-normal relaxed clock (representing heterogeneous evolutionary rates) was employed. Both birth-death (Gernhard, 2008) and Yule (Yule, 1925) speciation priors were independently evaluated using Bayes factor tests following Suchard, Weiss, and Sinsheimer (2001). Substitution models were unlinked and models of evolution selected in the phylogenetic analysis were used as priors. Taxon sampling included three individuals from each major S. gracilis clade recovered from the phylogenetic analysis as well as the following outgroup taxa:
Divergence time estimates were based on three priors, one fossil, and two molecular. The fossil prior was used to inform time to most recent common ancestor (TMRCA) for all Spilogale (offset of 3 Ma), based on the oldest Spilogale fossil (Hibbard, 1941a,b;Dalquest, 1971) with an exponential distribution with a mean of 1.63 such that 97.5% of the prior probability was below 9 Ma, a date which corresponds to the origin of Mephitidae in North America (Wang, Carranza-Castañeda, & Aranda-Gómez, 2014;Wang, Whistler, & Takeuchi, 2005). The second prior was based on a molecular estimate of TMRCA for the three extant genera of New World mephitids (Conepatus, Mephitis, and Spilogale) and was set using a normal distribution (mean = 9.2, SD = 1.7) to place the 2.5 and 97.5 quantiles between 6 and 13 Ma, the highest posterior density (HPD) interval associated with this molecular date (Eizirik et al., 2010). Finally, tree height was given a prior of 20 Ma based on TMRCA between Mydaus and extant New World skunks (Eizirik et al., 2010).
For S. gracilis, a BSP was estimated using the three-gene dataset for 13 individuals. Mutation rates estimated from previous BEAST analysis were used to inform prior distribution in BSP using a log-normal relaxed clock with a normally distributed tree height prior corresponding to MRCA for all S. gracilis (mean = 0.98, SD = 0.2). The analysis consisted of 10 million iterations sampled every 1,000 iterations. For all BEAST analyses, effective sample sizes of posterior parameters were evaluated using tRacER 1.5 (Rambaut, Suchard, & Drummond, 2009).

| Present-day models
Ecological niche models (ENM) were generated following the Biotic, Abiotic, and Movement (BAM) conceptual framework (Soberón & Peterson, 2005). Estimation of abiotically suitable conditions under this framework implies delimitation of the Grinnellian or scenopoetic fundamental niche (Soberón, 2007) and requires explicit assumptions regarding model input and interpretation (Soberón & Peterson, 2005), thereby providing a more transparent and repeatable methodological approach to ENM development (McDonough et al., 2015).
For S. gracilis, we initially defined "M" or the region of the world accessible to this species since its origin (Soberón & Peterson, 2005), as a 500-km buffer around all localities (See Appendix S3, Fig. S3.1). Fossil records of S. gracilis from Alberta, Canada (Kurten & Anderson, 1980), were used to demarcate the northeastern extent of the training area while the eastern edge was limited to the western most boundaries of the congener, S. putorius. The southern edge was truncated along the northern border of the distribution of the southern spotted skunk (S. angustifrons). All range boundaries were interpreted from expert generated distribution maps downloaded from the IUCN Redlist (www. iucnredlist.org/).
Occurrence records were obtained from museum voucher specimens via the Mammal Networked Information System (MaNIS; http://manisnet.org/) or from tag information collected from visited collections or provided by curators and collection managers (see Acknowledgments and Appendix S3, Table S3.1). All textual locality information was georeferenced following best-practices guidelines (Chapman & Wieczorek, 2006) and only specimens with a coordinate uncertainty less than 5 km were used in model development. To reduce potential biases associated with clustered sampling localities, occurrence data were spatially thinned by selecting a random subset of >100 occurrence records with a minimum distance between records of 0.6111 decimal degrees (~50 km) in ArcGIS 10.1 (ESRI, Redland, CA). This reduction in data was performed three times, and each subset was used to generate its own unique ENM. However, due to the disproportionate number of specimens from the 255 unique occurrence records found across the distribution of S. gracilis, spatially thinned subsets resulted in models that tended to over-predict suitability in regions associated with more records (e.g., California, Texas, see Appendix S3, Fig. S3.1). Thus, a manual thinning process whereby we selected a subset of records evenly distributed across the range of S. gracilis was used to generate final ENMs. This process yielded 84 unique occurrence records that were used to generate final models of the abiotic niche of S. gracilis (see Appendix S3, Table S3.1).
Present  environmental data in the form of bioclimatic variables were accessed from WorldClim (www.worldclim.org; Hijmans, Cameron, Parra, Jones, & Jarvis, 2005) at a spatial resolution of 2.5′ (~5 × 5 km). This coarser resolution was chosen over the available 30 arc-second data to better match the 5-km uncertainties associated with georeferenced specimen-locality coordinates. Although utilizing all 19 variables without addressing issues associated with correlated variables has been identified as a concern in ENM ( Guisan & Thuiller, 2005), we chose to include all 19 variables in model development following reasons outlined in McDonough et al. (2015).
Models were generated using the maximum entropy algorithm in maxENt 3.3.3 (Phillips, Anderson, & Schapire, 2006). Maxent models were generated under the "autofeatures"' option, allowing for linear, quadratic, product, threshold, and hinge features to be used to describe relationships between locations and environmental conditions (Merow, Smith, & Silander, 2013) with both clamping and extrapolating options deselected (Owens et al., 2013). Although raw outputs were recently suggested as preferable to logistic output within maxENt, given assumptions associated with τ, logistic output was chosen to allow for comparisons across species with inherently different model calibration regions (Ferguson, 2014;Merow et al., 2013). A total of 10 bootstrap model replicates were used to generate predictions of suitable conditions with final models representing an average of these 10 replicates. All raster processing and assessments were performed using a combination of aRcGIS 10.1 and the R statistical software environment (http://www.r-project.org/).

| Historic models
The resulting predictive models developed for present-day conditions were used to project suitable environmental conditions onto two Pleistocene climatic coverages, representing the last glacial maximum (LGM; ~21,000 years BP) and the last interglacial (LIG; ~120,000-140,000 years BP) periods, respectively. Environmental layers for both LGM and LIG were accessed from the WorldClim database (http://www.worldclim.org/past). LGM coverages are provided at 2.5 arc-minutes resolution for two climate models: the Community Climate System Model (CCSM, http://www.ccsm.ucar.edu/; (Kiehl & Gent, 2004) et al. (2006) are provided at 30 arc-seconds and thus were rescaled to 2.5 arc-minutes in aRcGIS 10.1 to match the present and LGM resolutions of 2.5 arc-minutes (~5 km). Results generated for both CCSM and MIROC circulation models were qualitatively similar (data not shown); thus, only CCSM was used to develop subsequent models.

| Mitochondrial phylogeny and geographic distribution of genetic variation
Phylogenetic analyses of the Cytb dataset (n = 97) recovered three well-supported clades within S. gracilis, subsequently referred to as West, East-Central, and Arizona ( Figure 1). These clades differed from one another with interclade Kimura-2 parameter (K2P) genetic distances ranging from 0.032 to 0.043, compared to 0.074-0.189 for outgroup taxa (Table 1). Individual gene trees for Cytb and D-loop recovered the three clades as did the concatenated 2-gene (Cytb + Dloop) dataset. ND5 was not included in phylogenetic estimates due to low sample size and minimal geographic coverage (see Appendix S1, Table S1.1). Given the depth of sampling for the Cytb tree and the congruency seen among D-loop and the concatenated dataset, we chose to present only the Cytb phylogeny.   (Hafner & Riddle, 2011;Riddle & Hafner, 2006

| Historical demography
Demographic reconstruction using BSP indicated a population increase for all clades starting between c. 0.15 and 0.10 Ma with continued population increase into the present (Figure 2).

| Present-day models
Modern climate Maxent models for S. gracilis indicated high suitability throughout the known distribution of this species with the exception of the "semidesert grassland" and "Chihuahua desert scrub" habitats (Gori & Enquist, 2003) of northern Mexico and New Mexico ( Figure 3).
However, these regions contain few records in our modeling occurrence dataset (Figure 3; Appendix S3,

| Pre-Pleistocene divergence
Divergence estimates for Spilogale indicated that the TMRCA occurred between late Miocene and early Pliocene (c. 6.53 Ma with a 95% HPD of 8.64-4.34 Ma), a date that is concordant with estimates based on both molecular data (Eizirik et al., 2010) and the fossil record (Wang et al., 2014). Similar divergence was observed among antelope squirrels Ammospermophilus across the Miocene-Pliocene boundary, where their three major lineages diverged across the warm deserts of North America (Mantooth et al., 2013 (Bell, Hafner, Leitner, & Matocq, 2010), and diversification of Thomomys across parts of the Great Basin (Belfiore, Liu, & Moritz, 2008). Thus, although PPD appears to have played a role in early stages of evolution among Spilogale, its role in shaping genetic diversity within species appears minimal, at least for S. gracilis.
Despite the lack of a strong signature of PPD on the genetic structure of S. gracilis as a whole, the genetic isolation of the Arizona clade does appear to partially correspond with well-known barriers that predate the Pleistocene (Figure 1, Numbers 1-4). Specifically, it appears that S. gracilis individuals within the Arizona clade were isolated from the West clade by the Bouse Embayment, a Late Miocene to Late

| Quaternary climate change
Although patterns associated with the Arizona clade appear to support the potential impacts of PPD on S. gracilis genetic structure, phylogeographic patterns at the intraspecific level appear to have resulted more from QCC than PPD according to data from divergence estimates, de-  (Rzedowski, 1978

| Signatures of secondary contact and dispersal across the desert southwest
One of the most important barriers separating fauna of the Chihuahuan and Sonoran Deserts is the Cochise filter barrier, a 1,500 m high plain spanning approximately 200 km (east-west) characterized by cool, transitional desert grasslands, known as the Deming Plains (Morafka, 1977). This barrier resulted in separation of desert lowland fauna throughout the Pleistocene, although oscillating climate cycles often led to expansion and contraction of grassland habitat which could explain why some species and not others show phylogeographic signatures associated with this barrier (Morafka, 1977). Phylogeographic breaks along this barrier have been recovered for bats (Weyandt & Van Den Bussche, 2007), birds (Klicka et al., 2016;Zink, Kessen, Line, & Blackwell-Rago, 2001), small mammals, plants, and reptiles (Hafner & Riddle, 2011;Riddle & Hafner, 2006 (Figure 3).

| Biogeography of small mammals in western North America
Recent studies on small mammals with similar distributions to S.

| Biogeography of small Carnivores in western North America
In comparison with similar, broadly distributed and generalist species of small carnivores, including other mephitid species (Barton & Wisely, 2012), S. gracilis displays a pattern more consistent with small mammals like rodents than similar or even smaller-sized carnivores (Aubry et al., 2009;Dawson, Hope, Talbot, & Cook, 2014;Harding & Dragoo, 2012). For example, studies of striped skunks, M.
mephitis (Barton & Wisely, 2012), ermine, Mustela erminea (Dawson et al., 2014), and red foxes, Vulpes vulpes (Aubry et al., 2009;Volkmann et al., 2015) all indicate recent divergences (<400 Ka) of intraspecific lineages following glaciation patterns across North America, but see Harding and Dragoo (2012). However, intraclade divergence within a western North America clade of the long-tailed weasel, Mustela frenata appears more in line with patterns seen in S. gracilis, with major lineage divergence occurring around 1 Ma (Harding & Dragoo, 2012). Although none of these small carnivore studies used ENM to identify putative Pleistocene refugia, inferred refugial locations also differed from those identified for S. gracilis.
However, striped skunks appear to have shared a similar refugium to S. gracilis individuals in the East-Central clade in southern Texas and northeastern Mexico (Barton & Wisely, 2012). Comparisons with these studies support the hypothesis that S. gracilis responded to QCC differently than other, similar sized and related small carnivores, indicating a unique evolutionary history for this species across western North America.

| CONCLUSIONS
Major phylogeographic patterns within S. gracilis appear to have resulted more from QCC than PPD. This pattern differs from other codistributed small carnivores, highlighting that individual species may respond differently to PPD and QCC, even when closely related (e.g., striped and spotted skunks). Our multifaceted approach highlights the complex response of S. gracilis populations to changing climates across the Quaternary, indicating that it was not glacial maxima alone that shaped genetic structure in this species, but instead a unique combination of both physiographic features (e.g., major rivers but not major mountain ranges) and temporally variable changes in habitat suitability associated with changing climates that best explain the observed variation. These patterns would have