Mitochondrial Genomes of the United States Distribution of Gray Fox (Urocyon cinereoargenteus) Reveal a Major Phylogeographic Break at the Great Plains Suture Zone

Reding, Dawn M.; Castañeda-Rico, Susette; Shirazi, Sabrina; Hofman, Courtney A.; Cancellare, Imogene A.; Lance, Stacey L.; Beringer, Jeff; Clark, William R.; Maldonado, Jesus E.

doi:10.3389/fevo.2021.666800

ORIGINAL RESEARCH article

Front. Ecol. Evol., 08 June 2021

Sec. Evolutionary and Population Genetics

Volume 9 - 2021 | https://doi.org/10.3389/fevo.2021.666800

This article is part of the Research Topic Genetics Applied to Mammal Conservation View all 9 articles

Mitochondrial Genomes of the United States Distribution of Gray Fox (Urocyon cinereoargenteus) Reveal a Major Phylogeographic Break at the Great Plains Suture Zone

$\r\nDawn M. Reding*$ Dawn M. Reding^1*

Susette Castañeda-Rico^2,3,4

Sabrina Shirazi^2†

Courtney A. Hofman^2†

Imogene A. Cancellare⁵

Stacey L. Lance⁶

Jeff Beringer⁷

William R. Clark⁸

Jesus E. Maldonado^2,3

¹Department of Biology, Luther College, Decorah, IA, United States
²Center for Conservation Genomics, Smithsonian Conservation Biology Institute, National Zoological Park, Washington, DC, United States
³Department of Biology, George Mason University, Fairfax, VA, United States
⁴Smithsonian-Mason School of Conservation, Front Royal, VA, United States
⁵Department of Entomology and Wildlife Ecology, University of Delaware, Newark, DE, United States
⁶Savannah River Ecology Laboratory, University of Georgia, Aiken, SC, United States
⁷Missouri Department of Conservation, Columbia, MO, United States
⁸Department of Ecology, Evolution, and Organismal Biology, Iowa State University, Ames, IA, United States

We examined phylogeographic structure in gray fox (Urocyon cinereoargenteus) across the United States to identify the location of secondary contact zone(s) between eastern and western lineages and investigate the possibility of additional cryptic intraspecific divergences. We generated and analyzed complete mitochondrial genome sequence data from 75 samples and partial control region mitochondrial DNA sequences from 378 samples to investigate levels of genetic diversity and structure through population- and individual-based analyses including estimates of divergence (F_ST and SAMOVA), median joining networks, and phylogenies. We used complete mitochondrial genomes to infer phylogenetic relationships and date divergence times of major lineages of Urocyon in the United States. Despite broad-scale sampling, we did not recover additional major lineages of Urocyon within the United States, but identified a deep east-west split (∼0.8 million years) with secondary contact at the Great Plains Suture Zone and confirmed the Channel Island fox (Urocyon littoralis) is nested within U. cinereoargenteus. Genetic diversity declined at northern latitudes in the eastern United States, a pattern concordant with post-glacial recolonization and range expansion. Beyond the east-west divergence, morphologically-based subspecies did not form monophyletic groups, though unique haplotypes were often geographically limited. Gray foxes in the United States displayed a deep, cryptic divergence suggesting taxonomic revision is needed. Secondary contact at a common phylogeographic break, the Great Plains Suture Zone, where environmental variables show a sharp cline, suggests ongoing evolutionary processes may reinforce this divergence. Follow-up study with nuclear markers should investigate whether hybridization is occurring along the suture zone and characterize contemporary population structure to help identify conservation units. Comparative work on other wide-ranging carnivores in the region should test whether similar evolutionary patterns and processes are occurring.

Introduction

Past climatic fluctuations have been a major force driving lineage differentiation and contributing to biodiversity patterns (Hewitt, 1996, 2000). Like many regions, North America experienced dramatic changes in climate as repeating cycles of glacial expansion and retreat occurred throughout the Quaternary Period (i.e., the past 2.6 million years). The resulting shifts in vegetation and food resources caused the distributions of temperate vertebrates to contract and expand, often leaving phylogeographic signatures of population decline, fragmentation, and divergence in isolated glacial refugia contrasted with growth, recolonization, and admixture during interglacial periods (Lessa et al., 2003; Hewitt, 2004). For habitat specialists or organisms with low vagility, these expansion-contraction cycles and intermittent biogeographic barriers set populations on divergent evolutionary trajectories that can be recognized today as separate species or subspecies (Stone et al., 2002; Johnson and Cicero, 2004; Kerhoulas and Arbogast, 2010; McDonough et al., 2020). In other cases, flexible ecological requirements and/or high vagility seem to have either prevented genetic structure from forming or largely erased the structure following contact and homogenization (Smith et al., 2011; Koblmüller et al., 2012; Kierepka and Latch, 2016). Yet recent studies have revealed cryptic phylogeographic structure in mobile taxa with broad, seemingly continuous distributions (Barton and Wisely, 2012; Harding and Dragoo, 2012; Reding et al., 2012; Goddard et al., 2015; Puckett et al., 2015), indicating species responses to past environmental change can be complex (Graham et al., 1996; Hofreiter and Stewart, 2009). Such cryptic genetic diversity presents a challenge for assessing and predicting the effects of past and future global climate change on biodiversity (Pauls et al., 2013), evaluating speciation hypotheses, understanding the relative roles of glacial vicariance and ongoing isolating mechanisms at secondary contact zones (Swenson, 2006), and directing conservation efforts to preserve evolutionarily significant units and thus the genetic legacy of species (Coates et al., 2018).

Although the magnitude and location of phylogeographic breaks vary depending on an organism’s degree of habitat specialization, dispersal ability, and potential for development of reproductive isolating barriers, comparative work has documented some geographic hot spots in North America where a variety of sister taxa meet at areas of secondary contact called suture zones (Remington, 1968; Swenson and Howard, 2004, 2005; Swenson, 2006; Rissler and Smith, 2010). The Great Plains Suture Zone, located around the 100th meridian in the central United States, is one of the strongly supported secondary contact zones, particularly for forest birds (Swenson and Howard, 2004, 2005; Swenson, 2006). Indeed, many avian genera include sister species pairs represented by eastern and western counterparts (Mengel, 1970; Rising, 1983; Lovette, 2005). Although debate has centered on whether this species-level diversification resulted primarily from Pleistocene vicariance (Johnson and Cicero, 2004; Weir and Schluter, 2004) or much earlier Pliocene events (Klicka and Zink, 1997; Zink et al., 2004; Barnosky, 2005), mounting evidence supports a role for Pleistocene events causing species-level divergence in many cases (Lovette, 2005; Peterson and Ammann, 2013; McDonough et al., 2020). A leading hypothesis for east-west divergences is that populations were separated into two or more disjunct Pleistocene forest refugia south of the glacial extent, an idea supported by climatic modeling of ecoregions (Hargrove and Hoffman, 2005) and ecological niche modeling of historical species distributions (Peterson and Ammann, 2013; Puckett et al., 2015; Loveless et al., 2016; Ferguson et al., 2017). These populations then diverged in allopatry, potentially with the eastern counterpart adapting to wetter and cooler environments and the western to drier and warmer environments (Rising, 1970, 1983; Webb and Bartlein, 1992; Swenson, 2006). Changing climatic conditions at the end of the Pleistocene (∼12,000 ybp) permitted refugia to expand into secondary contact, with exogenous (i.e., environmental) factors such as temperature (Rising, 1969; Swenson, 2006), precipitation (Moore and Price, 1993), and/or vegetation (Moore, 1977) gradients at the Great Plains likely playing pivotal roles in maintaining species boundaries along the hybrid zone, perhaps in combination with endogenous (i.e., genetic-based) factors (Bronson et al., 2003). Such a major biogeographic pattern should be common across diverse vertebrate taxa, but evidence for this phylogeographic break and species-level east-west divergence remains more equivocal for mammals (Barnosky, 2005). Particularly for generalist mammals, phylogeographic divergence is typically absent or shallow, and although east-west breaks have been found within mammalian species at the Great Plains (Reding et al., 2012), they may instead be located at physical barriers such as the Rocky Mountains Continental Divide (Rueness et al., 2003), or Mississippi River (Brant and Ortí, 2002; Cullingham et al., 2008; Barton and Wisely, 2012; Harding and Dragoo, 2012).

Recent genetic research involving the gray fox (Urocyon cinereoargenteus) has indicated populations in the western United States are surprisingly divergent from those in the eastern United States (Goddard et al., 2015; Hofman et al., 2015), though many gaps in our characterization of this pattern remain. Gray fox range across much of North America and into Central America and northern South America (Figure 1A; Hall, 1981). Generally associated with forest, woodland, and brushland, it is an opportunistic carnivore with a variable diet (small mammals, invertebrates, and fruits and other plant matter; Fritzell and Haroldson, 1982) and male-biased dispersal with recorded distances up to 135 km (Banfield, 1974). Across their range, 16 subspecies have been described primarily from body size, pelage color, and geographic differences; 7 subspecies of U. cinereoargenteus and the island fox endemic to the Channel Islands (U. littoralis) occur in the United States (Figure 1A; Hall, 1981). Gray foxes are legally harvested for fur throughout much of the range (Deems and Pursley, 1978). Despite the ecological and economic importance of gray fox, few studies have examined intraspecific genetic variation or evaluated genetic support for subspecies designations, with existing studies focused at local or regional scales (Bozarth et al., 2011; Goddard et al., 2015; Hofman et al., 2015). Bozarth et al. (2011) sequenced a portion of the mitochondrial DNA control region from gray fox from 15 states within the ranges of the 3 subspecies along the United States East Coast and found a lack of support for differentiation across this broad area, though the Northeast samples did show lower levels of genetic diversity indicative of recent range expansion. In mtDNA analyses focused on island foxes and mainland gray foxes in California, both Hofman et al. (2015) and Goddard et al. (2015) included a small number of gray fox samples from the eastern United States (Virginia and Georgia, respectively) for comparison and found (1) island foxes are nested within the gray fox phylogeny (i.e., gray fox is paraphyletic) and (2) eastern and western gray foxes show a deep divergence on par with that seen between some currently recognized species within the Canidae. However, in both of these studies, sampling was geographically limited and disjunct, making it uncertain where the suture zone lies and if additional cryptic structure may exist. In addition, sampling gaps and limited genetic sequence information also presents challenges for accurately estimating divergence time to identify whether the separation occurred during or prior to Pleistocene glaciation events. Using cytochrome b sequence data, Goddard et al. (2015) estimated a split time of approximately 500,000 ybp between the eastern and western lineages. Studies have shown that entire mitogenomes can sometimes call into question inferences made using short mitochondrial sequences (Knaus et al., 2011; Feutry et al., 2014; Hofman et al., 2015), so additional sequence data would help to clarify the date of this split.

FIGURE 1

Figure 1. Analysis of 400 bp mtDNA control region sequence data from samples of gray fox (Urocyon spp.). (A) Map shows the ranges of 7 morphologically-based gray fox (U. cinereoargenteus) subspecies (Hall, 1981) and the island fox (U. littoralis); the locations of sample groups (n = 34) and samples (n = 781) used in the analysis, including newly generated sequences (circles), and previously published sequences (squares; Bozarth et al., 2011; Goddard et al., 2015; Hofman et al., 2015); and 3 possible locations of a hypothesized secondary contact zone between deeply divergent eastern and western gray fox lineages (Goddard et al., 2015). Samples are colored according to whether its haplotype fell into the East (white) or West (black) clade. (B) Median joining haplotype network. Each circle represents a unique haplotype, with size proportional to frequency, as well as a pie chart with color indicating the morphologically-based taxonomic assignment of the corresponding samples. Black circles reflect inferred haplotypes not present in the dataset, and hash marks on the lines connecting haplotypes represent the number of mutations separating them. A minimum of 13 steps separates East and West haplotypes. (C) Principal coordinates plot from population pairwise Fst (K2P-corrected) values, with sample groups colored according to morphologically-based taxa shown in the map. The first axis explained 47.8% of the variation and clearly divided western from eastern sample groups at the Great Plains, the second axis explained 11.3% of the variation and separated sample groups within the East (note only scores for East sample groups are shown since West sample groups had scores near 0 for the second axis), and the third axis explained 7.9% of the variation and separated sample groups within the West (note only scores for West sample groups are shown since East sample groups had scores near 0 for the third axis).

Identifying the location and strength of genetic boundaries in gray fox is relevant for the conservation and management of this species. In particular, a recent petition to list the prairie gray fox (U. c. ocythous; Figure 1A) under the Endangered Species Act has stimulated the United States Fish and Wildlife Service to initiate a status review to determine if listing is warranted (Department of the Interior, 2012). In the petition, Wade and Alton (2012) argued that the prairie gray fox has experienced dramatic declines in states like Iowa, Arkansas, Missouri, and Minnesota due to the loss of early successional habitats from intensified agriculture or forest maturation, competition with expanding coyote populations, and hunting/trapping pressure. It is unclear, however, how much the population has declined relative to historic abundance or whether the prairie gray fox is a genetically distinct segment of the contiguous gray fox range in the eastern United States. Recent studies of mammals have shown that genetic patterns often do not match previously described subspecies (Culver et al., 2000; Cullingham et al., 2008; Reding et al., 2012). Patterns of genetic variation, however, should be concordant with subspecific boundaries if morphologically-based subspecies are to be considered valid.

We aimed to characterize gray fox mtDNA phylogeographic structure across the United States to: (1) identify the location of suture zone(s) between eastern and western lineages and investigate the possibility of additional cryptic divergences; (2) verify the date of divergence between the eastern and western lineages through additional sampling and data from full mitogenomes, and (3) assess population genetic structure within these major maternal lineages to evaluate support for morphologically-described subspecies.

Materials and Methods

Sampling

We collected gray fox DNA samples in 2003–2017 from legally harvested and road-killed animals in 26 United States via assistance from state wildlife agencies and furharvesters (Figure 1A). Tissue sources included muscle, skin or dry pelt, and toe pad stored at room temperature in sterile vials filled with silica desiccant. We collected associated geographic location information for each sample, generally at least to the county level. We used DNeasy Blood and Tissue centrifugation kits (Qiagen, Valencia, CA, United States) to extract DNA from 387 samples, selected to provide broad geographic coverage (Supplementary Appendix 1).

mtDNA Control Region Sequencing and Analysis

To compare results with the eastern samples analyzed in Bozarth et al. (2011) and western samples from Goddard et al. (2015), we sequenced approximately the same ∼426–431 bp segment of the mtDNA control region using universal primers H16498 (CCTGAACTAGGAACCAGATG) and L15910 (GAATTCCCCGGTCTTGTAAACC; Kocher et al., 1989). We used 15 μL PCR reactions containing: 0.5 U IDProof High Fidelity DNA Polymerase (IDLabs, London, ON, Canada), 1X IDProof buffer containing 2 mM MgSO₄, 0.2 mM dNTPs, 0.4 μM each primer, and 20 ng template DNA. Cycle conditions were: 95°C for 2 min, 30 cycles of 95°C for 30 s, 50°C for 30 s, 68°C for 30 s, followed by a final 10 min elongation at 68°C. We cleaned the PCR products using the ExoSAP method (Werle et al., 1994) and submitted them to the Iowa State University DNA Facility for cycle sequencing and analysis on an ABI 3730xl DNA Analyzer. Both directions were sequenced with the same primers used for PCR. We used Sequencher 5.2.3 (Gene Codes Corporation) to trim and align forward and reverse reads. We obtained full sequence coverage from 362 samples (including 7 samples containing 1 or more ambiguities), partial sequence coverage from 18 samples, and failed to obtain sequence data from 7 samples (Supplementary Appendix 1). Sequences are available in GenBank as accession numbers: MW597742 – MW598121.

In addition to the 380 sequences we generated, we downloaded from GenBank an additional 403 control region sequences (Supplementary Appendix 1): 229 samples with full sequence information from Bozarth et al. (2011), 134 samples from mainland (California and Georgia) gray fox from Goddard et al. (2015), as well as 26 mainland samples (25 California and 1 Virginia) and a single representative sample from each of the 14 Channel Island gray fox mitogenome haplotypes from Hofman et al. (2015). We aligned sequences using MAFFT in Geneious 10.0.7¹, trimmed to the same area analyzed in Goddard et al. (2015) resulting in a 408 bp alignment, and deleted an 8 bp region of homopolymer C to end with a 400 bp alignment. For downstream analyses, we included samples with at least 90% coverage (363 bp or more) of control region sequence, which resulted in omission of 2 of our sequenced samples. We ordered the remaining sequences (n = 781) with full sequences first, followed by sequences containing an ambiguity, and lastly partial sequences. We exported the aligned sequences from Geneious in FASTA format and manually replaced ambiguities with Ns. We used DnaSP 5.10.01 (Librado and Rozas, 2009) to export the data as a Roehl Data File with gaps/missing considered. Because it collapses ambiguous and missing data using the most common state found in the closest sequences, we then used Network 10.1 (Bandelt et al., 1999) to collapse sequences into unique haplotypes and to generate a median joining haplotype network. For improved visualization, we used PopART 1.7 (Leigh and Bryant, 2015) to draw a median joining network based on the haplotypes identified by Network.

To permit population-level analyses, we grouped samples into 34 a priori populations based on geographic proximity, state boundaries, biogeographic patterns, and population information from Bozarth et al. (2011) and Goddard et al. (2015; Figure 1A and Table 1). We used ArcGIS 10 (ESRI, Redlands, CA, United States) to calculate the geographic mean center of individuals assigned to each sample group, which we used as the spatial coordinates for the group. We used these sample groups and the unique haplotypes identified in Network (Supplementary Alignment 1) to generate Arlequin input files. We used Arlequin 3.5 (Excoffier and Lischer, 2010) to calculate population genetic diversity estimates (haplotype number, haplotype diversity, mean number of pairwise differences, and nucleotide diversity) and to examine population structure with pairwise F_ST and AMOVA. For these, we used the Kimura 2 parameter (K2P; Kimura, 1980) model with gamma = 0.948, since it was found to be the best model of sequence evolution with BIC criteria in jModeltest 2.1 (Darriba et al., 2012). Because the haplotype network revealed a clear separation between eastern and western United States samples (see section “Results”), we performed AMOVA with sample groups split into these two sets (East had 552 individuals in 19 sample groups; West had 229 individuals in 15 sample groups; Table 1). We used the REG procedure in SAS OnDemand (SAS Institute Inc., Cary, NC, United States) to perform a regression between diversity estimates and latitude of the mean center for each sample group, using sample size and area of the minimum convex polygon as covariates. We used GenAlEx 6.5 (Peakall and Smouse, 2006, 2012) to calculate geographic distance between sample groups and to perform a Mantel test (999 permutations) to test for a positive relationship between K2P-corrected genetic (F_ST/[1–F_ST]) and ln-transformed geographic distance, both overall and separately within the east and west sample groups. To see if additional substructure could be revealed within each region, we performed SAMOVA (Dupanloup et al., 2002) from K = 2–10 separately on East and West samples using 500 iterations and the same sequence model specified above. We also conducted a principal coordinates analysis (PCoA) on the Arlequin-generated pairwise F_ST values using GenAlEx 6.5. We removed three populations (WA, CA-SE, CA7) that had fewer than 10 individuals from the Mantel and regression analyses, SAMOVA, and PCoA.

TABLE 1

Table 1. Information and indicators of molecular diversity (400 bp mtDNA control region) for gray fox (Urocyon cinereoargenteus) sample groups, except *CI = Urocyon littoralis.

Mitogenome Sequencing and Analysis

We carried out sequence capture to generate complete mitogenomes from 67 of our samples as well as 28 samples original to the Bozarth et al. (2011) study (Supplementary Appendix 1). Captures were completed with samples pooled in groups of ∼6 libraries following MyBaits protocol v 3. Following capture, libraries were amplified 14 cycles. Pools were sequenced on five partial Illumina Miseq runs (2 × 150 bp) at the Smithsonian Institution Center for Conservation Genomics and one partial lane of a HiSeq (2 × 138 bp) at UCLA. We processed the combined reads from all runs using Prinseq-lite v. 0.20.4 for quality filtering (–min_len 20 –min_qual_mean 15) and mapped filtered reads to MVZ_206290 using the BWA mem algorithm. Duplicates were removed with Samtools v. 1.1 and the bam files were imported into Geneious v. 10.0.8. Consensus sequences were called in Geneious and all samples with less than 150 ambiguities including missing data (n = 75) were aligned with MAFFT v.1.36 as implemented in Geneious. The resulting alignment was visually inspected and indels were assessed by returning to the original bam alignment files. In one case (sample VT-101), a NUMT was identified and removed. Additionally, the genes were translated into amino acids and inspected for stop codons in Geneious. We used Network 10.1 to collapse the mitogenomes into unique haplotypes. Sequences are available in GenBank as accession numbers: MW599994 – MW600068. We then used PopART 1.7 to visualize a median joining network consisting of these 75 mitogenome samples with 40 published Urocyon mitogenomes from Hofman et al. (2015) that included all 26 U. cinereoargenteus mitogenomes (25 from California, 1 from Virginia) and a representative of each of the 14 unique U. littoralis mitogenome haplotypes (Figure 2A).

FIGURE 2

Figure 2. Phylogeographic patterns of United States distribution of Urocyon taxa based on mitogenome sequence data. (A) Map shows distributions of the morphologically-based Urocyon taxa and locations of mitogenome samples (n = 115) used in the analysis, including newly generated sequences (circles) and previously published sequences (squares). Samples are colored according to whether its mitogenome haplotype fell into the East (white), West (black), or Island (gray) clade. (B) Median joining network of mitogenomes from 115 samples, including newly sequenced U. cinereoargenteus samples (n = 75) and a subset of samples from Hofman et al. (2015) downloaded from GenBank (U. cinereoargenteus, n = 26; U. littoralis, n = 14). Each circle represents a unique haplotype, with size proportional to frequency, as well as a pie chart with color indicating the morphologically-based taxonomic assignment of the corresponding samples as shown in the map. Black circles reflect inferred haplotypes not present in the dataset, and hash marks on the lines connecting haplotypes represent the number of mutations separating them. A minimum of 286 steps separates East and West haplotypes. (C) Bayesian tree based on the analysis of 104 unique Urocyon mitogenome haplotypes and two outgroups (Vulpes vulpes and Canis latrans, removed for clarity). Nodal support for clades is depicted by the size of black circles, with most clades showing posterior probability values ≥0.99 (see Supplementary Figure 5 for precise values). Sample codes are colored to reflect the morphologically-based taxonomic assignment.

In phylogenetic analyses, we used an alignment consisting of 106 mitogenomes: 68 unique mitogenomes we recovered from our 75 U. cinereoargenteus assemblies; 36 unique published Urocyon mitogenomes from Hofman et al. (2015) which included 14 U. littoralis and 22 U. cinereoargenteus mitogenomes; and two outgroups from GenBank: Canis latrans (NC_008093.1) and Vulpes vulpes (NC_008434.1). We used Partition Finder 2.1.1 (Lanfear et al., 2016) to test model and partition selection using linked, corrected Akaike Information Criterion and greedy parameters. The Partition Finder analysis detected 25 partitions that we incorporated in phylogenetic reconstruction. We performed a Bayesian Inference analysis (BI) in MrBayes 3.2.6 (Huelsenbeck and Ronquist, 2001; Ronquist and Huelsenbeck, 2003) using 50,000,000 generations sampling every 1,000 generations. We visualized output parameters using Tracer v1.7.1 (Rambaut et al., 2018) to check for convergence between runs and we discarded the first 25% of the trees as burn-in. The final tree and support values were visualized using FigTree v1.4.4.

To estimate molecular dates of divergence, we used partial mitochondrial genomes (no D-loop included). The alignment (Supplementary Alignment 2) included 104 unique haplotypes of Urocyon and three outgroups from GenBank: Canis latrans (NC_008093.1), Vulpes vulpes (NC_008434.1), and Ursus americanus (NC_003426.1). We used BEAST v2.5.2 (Bouckaert et al., 2014) under an uncorrelated lognormal strict molecular clock model and the Yule speciation processes model with 25 partitions. We performed Bayesian Markov chain Monte Carlo searches to obtain the time to the most recent common ancestor for the main lineages. We sampled trees and divergence dates for all nodes every 10,000 iterations for 100,000,000 generations. Divergence time estimates were based on three fossil calibration points. The first was based on the early Pliocene Urocyon fossil dated to 5.33–2.558 million years ago [Ma; mean = 3.5 (in real space), stdev = 0.5; McKenna and Bell, 1997]. The second point was for the root age of Canidae (offset = 10 Ma, mean = 0.7 [12 Ma], stdev = 0.415 [14 Ma]), as used by Slater et al. (2012). The third point was for the split of U. cinereoargenteus and U. littoralis (offset = 0.007, mean [in real space] = 0.704 [0.5 Ma], stdev = 0.36 [1.2 Ma]), as used by Slater et al. (2012). We checked convergence statistics for effective sample sizes using Tracer v1.7.1 (Rambaut et al., 2018) and TreeAnnotator v2.5.2 (available in the BEAST package) to get a consensus tree with node height distribution after elimination of 25% of trees as burn-in. We visualized BEAST results using FigTree v1.4.4. Partition Finder, MrBayes and BEAST analyses were performed on the Smithsonian Institution High Performance Computer Cluster².

Results

mtDNA Control Region

Among 781 individuals successfully sequenced or downloaded from Genbank, we identified 114 unique haplotypes, encompassing 63 substitutions and 2 indels observed at 61 polymorphic sites across the 400 bp control region segment. Overall, haplotype diversity = 0.961, nucleotide diversity = 0.0255, and mean number of pairwise differences = 10.181. At the population-level, all three diversity metrics decreased at northern latitudes in the East (P < 0.05), but there was no significant relationship in the West (P > 0.05; Figure 3, Supplementary Figure 1, and Table 1).

FIGURE 3

Figure 3. Estimates of haplotype diversity for 31 gray fox sample groups (n > 10) based on 400 bp of mtDNA control region. (A) Circles are located at the geographic mean center of each sample group and size reflects the quantile of the haplotype diversity value. (B) Linear regression of haplotype diversity by latitude for East and West sample groups.

The median-joining network (Figure 1B) strongly supported two highly diverged clades (East and West) separated at the Great Plains rather than the Mississippi River or Continental Divide (Figure 1A). Haplotypes from these two clades on the network are separated by at least 13 mutations (considering gaps). All individuals within sample groups had strictly East or West haplotypes with one exception: one sample in Oklahoma (KS-OK sample group) had a West haplotype when all others had East haplotypes (Figure 1A).

Estimates of pairwise F_ST values between the sample groups ranged from effectively 0 to 1, with the most differentiation occurring between East and West sample groups (Supplementary Figure 2). PCoA showed the first axis explained 47.8% of the variation and clearly divided West from East sample groups at the Great Plains (Figure 1C). The second axis explained 11.3% of the variation and separated sample groups in the East, with the lower-diversity northern populations in the Upper Midwest (MN-WI-ND, IA-NE) and New England (RI-CT, NH-VT-MA) clustering away from each other and the rest of the sample groups. The third axis explained 7.9% of the variation and separated sample groups in the West, with the interior sample groups clustering from those in California. Although there was a general pattern of geographic concordance, the sample groups did not clearly cluster by morphologically-defined subspecies (Figure 1C). Similarly, the AMOVA showed East and West sample groups were significantly differentiated (F_CT = 0.741, P < 0.001). Results of the SAMOVA indicated that the grouping with the highest F_CT value in the East differentiated the Upper Midwest (MN-WI-ND + IA-NE), New England (RI-CT + NH-VT-MA), LA, and all other sample groups (K = 4; F_CT = 0.2738, P < 0.001; Supplementary Table 1 and Supplementary Figure 3). The pattern was more complicated in the West, where F_CT values continued to incrementally increase with increasing K (Supplementary Figure 3). However, the largest increase in F_CT occurred when moving to K = 4 (F_CT = 0.3307, P < 0.001), which differentiated two interior groups ([CO_NM + TX] and [AZ_NM + NV + UT]), CA2, and all other California sample groups (Supplementary Table 1). The data demonstrated a significant pattern of isolation by distance overall (r = 0.376, P = 0.001), and within the East sample groups (r = 0.351, P = 0.001) but not the West sample groups (r = 0.175, P = 0.159; Supplementary Figure 4).

Mitogenomes

We identified 104 unique mitogenome haplotypes among 115 Urocyon samples. The median-joining network supported a clear east-west division, with a minimum of 286 steps (considering gaps) separating East (n = 53) and West (n = 51) haplotypes (Figure 2B). The island fox (U. littoralis) mitogenome haplotypes formed a unique cluster in the West and were not shared with any U. cinereoargenteus samples (in contrast to the control region haplotypes, Figure 1B), but overall, there was a lack of phylogeographic cohesiveness of morphologically-identified taxa (Figure 2B). Bayesian inference also strongly supported monophyletic East and West clades (BPP = 1), and a monophyletic clade (BPP = 1) of island foxes was nested within the West gray fox clade (Figure 2C). Although phylogeographic structure was somewhat more apparent within the West, with some subclades restricted to California and others more broadly distributed across the western interior (Supplementary Figure 5A), neither the West nor the East (Supplementary Figure 5B) showed support for geographically-restricted subclades consistent with recognized subspecies. The BEAST analysis estimated a mean divergence time of 0.796 Ma (95% HPD = 0.498–1.153 Ma) for the East and West lineages (Figure 4). Time to most recent common ancestor (TMRCA) for the East lineage was estimated at 0.132 Ma (95% HPD = 0.080–0.194), TMRCA for the West lineage was 0.180 Ma (95% HPD = 0.113–0.263 Ma), and TMRCA for the island fox was 0.087 Ma (95% HPD 0.053–0.125 Ma).

FIGURE 4

Figure 4. Estimated divergence time of eastern and western lineages of Urocyon based on BEAST analysis of mitochondrial genomes (D-loop omitted) for 104 Urocyon samples and 3 outgroups. (A) Full tree with East and West Urocyon clades collapsed and (B) Urocyon subtree. Blue bars on nodes and corresponding values in parentheses indicate the 95% highest probability density (HPD) interval for the mean time estimates.

Discussion

The control region and mitogenome data strongly support two major lineages of gray foxes in the United States: East and West clades that meet at the Great Plains. The pattern was consistently recovered in haplotype networks, phylogenetic reconstructions, and population-based analyses. These findings are concordant with the hypothesis of divergence in disjunct eastern and western North American forest refugia south of the ice sheets during glacial episodes. Phylogeographic support for eastern and western Pleistocene refugia has been reported in other North American carnivores, including American black bears (Ursus americanus, Wooding and Ward, 1997; Puckett et al., 2015), American and Pacific martens (Martes americanus and M. caurina, Stone et al., 2002), bobcats (Lynx rufus, Reding et al., 2012), red fox (Vulpes vulpes, Aubry et al., 2009), striped skunk (Mephitis mephitis, Barton and Wisely, 2012), spotted skunk (Spilogale spp., McDonough et al., 2020), and long-tailed weasel (Mustela frenata, Harding and Dragoo, 2012).

The Great Plains Suture Zone is recognized as a hotspot of secondary contact for diverse taxa (Swenson and Howard, 2004, 2005; Lovette, 2005; Swenson, 2006), but it is not as well-documented for wide-ranging carnivores. Many phylogeographic studies have focused on montane or boreal species, where the major divergence occurs between Pacific coastal and interior continental populations (e.g., V. vulpes, Aubry et al., 2009; M. caurina/americana, Stone et al., 2002; U. americanus, Wooding and Ward, 1997; Puckett et al., 2015), but with the latter sometimes showing shallow subdivision into eastern and western lineages at or west of the Great Lakes (Aubry et al., 2009; Pelletier et al., 2011; Puckett et al., 2015). Studies have also examined open habitat or grassland-associated species, which typically lack phylogeographic structure consistent with a single mid-continent Pleistocene refugium (e.g., black-footed ferret [Mustela nigripes], Wisely et al., 2008; American badger [Taxidea taxus], Kierepka and Latch, 2016; and coyote [Canis latrans], Koblmüller et al., 2012). Large carnivores such as the cougar (Puma concolor, Culver et al., 2000) and gray wolf (Canis lupus, Koblmüller et al., 2016) also lack phylogeographic structure due to recent (∼10 ka) recolonization of North America. Because they are often not of conservation concern, relatively few studies have examined continental patterns in common, temperate generalist and forest-associated carnivores, where a mid-continent phylogeographic break may be most likely (Barton and Wisely, 2012; Reding et al., 2012; McDonough et al., 2020).

Our study adds to a small but growing list of phylogeographic research on wide-ranging, common carnivores revealing east-west divergence. In gray fox, we found the location of secondary contact to center on the southern Great Plains (e.g., between western Texas and eastern Oklahoma). Bobcats show a nearly identical pattern, although the break is more diffuse and extends north into the northern Great Plains where bobcats are more common than gray fox (Reding et al., 2012). Divergence in the spotted skunk complex (Spilogale spp.) is also well characterized in the region, with the western lineage S. leucoparia inhabiting western Texas and S. interrupta in eastern Texas (McDonough et al., 2020). The Great Plains is also a region of admixture for southern and western lineages of striped skunk, though eastern lineages appear not to have expanded west of the Mississippi River drainage, a putative historical and contemporary geographic barrier (Barton and Wisely, 2012). Similarly, long-tailed weasels show diverged eastern and western lineages that appear to split at the Mississippi River (Harding and Dragoo, 2012). However, sampling gaps in the central United States lead to uncertainty in the location of secondary contact for these two species. Although all striped skunks sampled in Louisiana showed western haplotypes, consistent with a Mississippi River break, no samples were collected from areas just west of the river such as Arkansas, Missouri, Iowa, and Minnesota (Barton and Wisely, 2012). These states were also unsampled in long-tailed weasels, and eastern haplotypes were found in Louisiana and Texas, suggesting the break could exist farther west (Harding and Dragoo, 2012). Although we lacked samples from eastern Texas to pinpoint the precise location of secondary contact in gray foxes, all individuals sampled in states just west of the Mississippi River showed eastern haplotypes, providing no support for the Mississippi River as a major phylogeographic barrier in this species, unlike other carnivores (Cullingham et al., 2008; Barton and Wisely, 2012; Shaffer et al., 2018). Additional sampling and comparative work in the Great Plains region could help reveal concordant patterns in these and other carnivore species. For example, American mink (Neovison vison, García et al., 2017) and northern raccoon (Procyon lotor, Cullingham et al., 2008; Louppe et al., 2020) show promising trends of an east-west division.

Using mitogenomes, we estimated the divergence between eastern and western gray fox lineages to date to ∼0.8 Ma, a surprisingly deep division slightly older than the ∼0.5 Ma estimate reported by Goddard et al. (2015) based on mtDNA cytochrome b gene but still within the Pleistocene. For mobile carnivores, the east-west divergence is often shallow, typically dating to the Last Glacial Maximum at ∼23 ka (Aubry et al., 2009; Reding et al., 2012) or near the Penultimate Glacial Maximum at ∼140 ka (Puckett et al., 2015). Populations of species with high dispersal and flexible ecological requirements may have repeatedly expanded and interbred during interglacial periods when habitat was more connected, largely erasing signatures of divergence during earlier glacial extents (Barnosky, 2005). Compared to many carnivores, gray foxes have more intermediate dispersal distances (Whitmee and Orme, 2013), and with their unique tree-climbing behavior for foraging and predator avoidance, they are more strongly associated with forests (Fritzell and Haroldson, 1982). Fossils of U. cinereoargenteus are recorded at over 40 different Pleistocene sites in North America, predominantly in southwestern (e.g., California, Arizona, New Mexico, and Texas) and southeastern (Florida, Georgia, South Carolina, and Tennessee) United States (Kurtén and Anderson, 1980; Williams et al., 2018). Although records extend to Pennsylvania, Indiana, and Missouri, fossils appear absent from the Great Plains region prior to the Holocene (Williams et al., 2018), consistent with the area forming an ecological barrier for forest species (Mengel, 1970). Given the degree of divergence, eastern and western gray fox lineages have been isolated for multiple glacial-interglacial cycles. The timing is just shy of east-west divergence dates estimated for smaller carnivores such as spotted skunk and long-tailed weasel at around 1.5 and 2.2 Ma, respectively, (Harding and Dragoo, 2012; McDonough et al., 2020), which may warrant species-level recognition (McDonough et al., 2020). However, sampling of Urocyon taxa from across its range and the addition of nuclear markers will be needed to verify the divergence date estimated here.

In addition to the east-west divergence, Quaternary climate change influenced patterns of genetic diversity within each of these lineages. For example, low genetic diversity in northern areas in the East (e.g., Great Lakes and New England), where glacial ice persisted until ∼10–15 ka, is consistent with founder effects associated with rapid, recent recolonization (Hewitt, 1996). A similar pattern was seen with bobcats (Reding et al., 2012), indicating contemporary gene flow has not erased this pattern of recent expansion (Smith et al., 2011). Neither a pattern of declining diversity with increasing latitude, nor isolation by distance, was seen in the West. Such results may be indicative of a more complex scenario of multiple glacial refugia in the western United States, as indicated in studies of other carnivore taxa (Harding and Dragoo, 2012; Ferguson et al., 2017). To provide a more comprehensive view of the biogeography of Urocyon, additional sampling to cover the full range (e.g., from Mexico and Central and South America) of gray foxes will be needed, and ecological niche modeling could help reconstruct the locations of glacial refugia and recolonization routes (Puckett et al., 2015; Ferguson et al., 2017).

Our results provide some insight into the validity of current taxonomic classifications for Urocyon. Despite broad-scale sampling, we did not recover any additional, cryptic lineages of gray foxes in the United States. The east-west split does coincide with some noted phenotypic differences, with pelts from eastern foxes tending to be darker, duller, and less silvery (Obbard, 1987), and western forms generally having a more slender body with a longer tail and longer and more pointed ears (Mearns, 1891). As argued by Goddard et al. (2015), the divergence between eastern and western gray fox lineages is on par with some species-level distinctions in Canidae. Within the gray fox east-west lineages, we did not recover any clear phylogeographic structure besides a mitogenome monophyletic group of island foxes. Although the island fox is currently recognized as a separate species (U. littoralis), it is nested within the western gray fox lineage. A similar pattern occurs with the co-occurring island spotted skunk (Spilogale gracilis amphialus), which is considered a subspecies of the western spotted skunk (McDonough et al., 2020). Mitogenomes for the prairie gray fox (U. c. ocythous) showed no evidence of monophyly, but some haplotypes were present only in the Midwest and population-level analyses (PCoA and SAMOVA) tended to differentiate this subregion. Overall, haplotypes from major geographic regions and putative subspecies do not form discrete lineages, but structure within the eastern and western lineages is likely to be more recent and better characterized with nuclear markers. Our findings suggest taxonomic revision will be needed, and we recommend a comprehensive range-wide analysis considering phenotypic and molecular (mtDNA and nuclear markers) data to better inform pending conservation decisions (Department of the Interior, 2012).

The Great Plains may be a common but understudied phylogeographic break and zone of secondary contact in North American carnivores. If so, the finding has important implications for evaluating the importance of Pleistocene climate change as a driver of mammalian evolutionary divergence and speciation (Barnosky, 2005), informing conservation and management of biodiversity (Coates et al., 2018), and providing insight into the potential for recombination and spread of lineage-specific pathogens including zoonotic diseases (Dragoo et al., 2006; Barton and Wisely, 2012). Given the pattern we revealed with mtDNA, follow-up with nuclear markers would help with: (1) characterizing hybrid-zone structure and dynamics, including the extent and direction of admixture; (2) identifying genome regions potentially under divergent selection across the zone; and (3) evaluating population genetic structure within these two lineages. Furthermore, sampling from the gray fox range in Mexico and Central and South America will be necessary to fully understand biogeographic patterns and taxonomic classification of species and subspecies within the Urocyon genus. Finally, additional wide-ranging carnivore taxa such as long-tailed weasel, American mink, North American river otter (Lontra canadensis), and northern raccoon should be genetically sampled in the Great Plains region to examine broader support for a mammalian phylogeographic break at the Great Plains Suture Zone.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI GenBank; MW597742 – MW598121 and MW599994 – MW600068.

Ethics Statement

Ethical review and approval was not required for the animal study because data were derived from samples taken from legally-harvested or road-killed animals, or from previously published studies.

Author Contributions

DR, JB, WC, CH, and JM conceived and designed the study. DR, JB, SL, IC, and JM provided samples. DR, SS, and CH performed the laboratory work and sequence assemblies. DR performed the data analysis and SC-R performed the mitogenome phylogenetic and dating analyses. DR evaluated the results and drafted the manuscript together with CH, SC-R, and JM. DR and JM provided laboratory space for the project. All authors contributed to manuscript writing and approved it for publication.

Funding

This work was supported by Luther College, the Missouri Department of Conservation, Smithsonian Institution, and through grant ISF13-04 from the Iowa Science Foundation administered by the Iowa Academy of Science. This work was partially supported by the Department of Energy Office of Environmental Management under Award Number DE-EM0004391 to the University of Georgia Research Foundation.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The project would not have been possible without the cooperation and support of many fur harvesters, fur buyers, state and federal agencies and biologists in collecting and donating samples and information (see Supplementary Appendix 1). We thank M. Oberfoell, T. May, J. Ash, B. Wollenzien, M. Crotty, J. Hadish, and L. Brondyke for assisting with DNA extractions and N. Rotzel McInerney for lab support at the Center for Conservation Genomics. We also thank F. Quintela and two reviewers for comments on an earlier draft.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.666800/full#supplementary-material

Supplementary Figure 1 | Geographic pattern of nucleotide diversity for 31 Urocyon sample groups (n > 10) based on 400 bp of mtDNA control region. (A) Circles are located at the geographic mean center of each sample group and size reflects the quantile of the nucleotide diversity value. (B) Linear regression of nucleotide diversity by latitude for East and West sample groups.

Supplementary Figure 2 | Matrix of pairwise F_ST values (K2P-corrected) between Urocyon sample groups (n = 34) based on 400 bp of mtDNA control region.

Supplementary Figure 3 | Results of SAMOVA for (A) 19 East and (B) 12 West Urocyon sample groups (n > 10) based on 400 bp of mtDNA control region.

Supplementary Figure 4 | Relationship between pairwise geographic distance and genetic distance for East (n = 19) and West (n = 12) Urocyon sample groups based on 400 bp of mtDNA control region.

Supplementary Figure 5 | MrBayes subtrees for (A) West haplotypes and (B) East haplotypes based on the analysis of 104 unique Urocyon mitogenome haplotypes and two outgroups (Vulpes vulpes and Canis latrans, removed for clarity). Support for clades (BPP) is shown next to each node. Urocyon sample names reflect species (Uc = U. cinereoargenteus; Ul = U. littoralis), morphologically-based subspecies, location (state or island), and sample ID or GenBank accession number, and colors correspond to identified major clades. Maps show the ranges of morphologically-based subspecies, and the sampling locations of the mitogenomes are depicted as colored dots corresponding to major clade (as shown in the tree).

Supplementary Table 1 | Results of SAMOVA for 19 East and 12 West Urocyon sample groups (n > 10) based on 400 bp of mtDNA control region.

Supplementary Appendix 1 | An Excel file containing information about the samples used in this study.

Supplementary Alignment 1 | Fasta alignment file for mtDNA partial control region haplotypes.

Supplementary Alignment 2 | Fasta alignment file for mtDNA complete genome sequences (D-loop omitted) for BEAST analysis.

Footnotes

References

Aubry, K. B., Statham, M. J., Sacks, B. N., Perrine, J. D., and Wisely, S. M. (2009). Phylogeography of the North American red fox: vicariance in Pleistocene forest refugia. Mol. Ecol. 18, 2668–2686. doi: 10.1111/j.1365-294x.2009.04222.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Bandelt, H.-J., Forster, P., and Röhl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. doi: 10.1093/oxfordjournals.molbev.a026036

PubMed Abstract | CrossRef Full Text | Google Scholar

Banfield, A. W. F. (1974). The Mammals of Canada. Toronto: University of Toronto Press.

Google Scholar

Barnosky, A. D. (2005). Effects of Quaternary climate change on speciation in mammals. J. Mamm. Evol. 12, 247–264. doi: 10.1007/s10914-005-4858-8