High genomic diversity and candidate genes under selection associated with range expansion in eastern coyote (Canis latrans) populations

Abstract Range expansion is a widespread biological process, with well‐described theoretical expectations associated with the colonization of novel ranges. However, comparatively few empirical studies address the genomic outcomes accompanying the genome‐wide consequences associated with the range expansion process, particularly in recent or ongoing expansions. Here, we assess two recent and distinct eastward expansion fronts of a highly mobile carnivore, the coyote (Canis latrans), to investigate patterns of genomic diversity and identify variants that may have been under selection during range expansion. Using a restriction‐associated DNA sequencing (RADseq), we genotyped 394 coyotes at 22,935 SNPs and found that overall population structure corresponded to their 19th century historical range and two distinct populations that expanded during the 20th century. Counter to theoretical expectations for populations to bottleneck during range expansions, we observed minimal evidence for decreased genomic diversity across coyotes sampled along either expansion front, which is likely due to hybridization with other Canis species. Furthermore, we identified 12 SNPs, located either within genes or putative regulatory regions, that were consistently associated with range expansion. Of these 12 genes, three (CACNA1C, ALK, and EPHA6) have putative functions related to dispersal, including habituation to novel environments and spatial learning, consistent with the expectations for traits under selection during range expansion. Although coyote colonization of eastern North America is well‐publicized, this study provides novel insights by identifying genes associated with dispersal capabilities in coyotes on the two eastern expansion fronts.

Broadly, range expansion is expected to result in reduced genome-wide diversity relative to the core range as a consequence of small population sizes and serial founder events (Mayr, 1954;Nei, Maruyama, & Chakraborty, 1975). Strong population structure is also expected along the expansion axis, with recently expanded populations often representing differentiated genetic clusters from those in the core range (Ibrahim and Nichols, 1996). Though demographic factors, such as fecundity and population density, influence population structure and genome-wide diversity of recently expanded populations (Hagen et al., 2015), natural selection of traits (i.e., reproduction, dispersal) associated with range expansions may also play an important role. Theoretically, traits facilitating range expansion should experience differential selection pressures along the axis of expansion (Travis & Dytham, 2002;Phillips et al., 2008;Burton, Phillips, & Travis, 2010). For instance, genes associated with exploratory behavior and dispersal abilities are predicted to be beneficial at the front of the expansion axis given such traits directly facilitate movement into and subsequent colonization of a new habitat (Burton et al., 2010;Hughes, Dytham, & Hill, 2007;Phillips, Brown, Travis, & Shine, 2008;Travis & Dytham, 2002). Reproductive traits are also predicted to be under selection, as reduced competition and smaller population sizes at the front of the expansion may favor increased reproductive effort (Burton et al., 2010).
Although predictions for adaptive evolution during range expansion are well described in theory, it is a challenge in practice to identify loci under selection at the range periphery for several reasons.
First, a stochastic phenomenon known as "allele surfing", a consequence of serial founder events and drift at the expansion front, may drive even deleterious alleles to high frequencies along the expansion axis. This process has a strong theoretical basis (Edmonds, Lillie, & Cavalli-Sforza, 2004;Klopfstein, Currat, & Excoffier, 2006) and been suggested in empirical studies for a range of taxa (Hofer, Ray, Wegmann, & Excoffier, 2009;Gralka et al., 2016;Streicher et al., 2016). Therefore, identifying genomic variants with substantial changes in frequency along the expansion axis alone is not a sufficient evidence of recent selection. Additionally, variation in allele frequency may be driven by environmental factors (e.g., novel habitats and food resources) that occur in the expanding range but are independent of traits (e.g., dispersal and reproductive capabilities) that directly facilitate range expansion.
While the effects of allele surfing and environmental factors cannot be completely accounted for when studying range expansion, replicate expansion fronts across distinct environments can help disentangle the relative impact of these forces (Swaegers et al., 2015;White, Perkins, Heckel, & Searle, 2013). As allele surfing is a stochastic process, it is less likely that the same genomic variant would undergo a frequency shift in the same direction relative to the historical range along multiple independent expansion axes. Similarly, when species traverse distinct environments, increases or decreases in frequency at the same loci are less likely to be driven by local adaptation. Therefore, genomic variants that undergo similar frequency shifts across multiple independent axes of expansion are reasonable candidates for range expansion genes.
Coyotes (Canis latrans) provide a tractable system to address questions related to range expansion genomics. Confined to the western and central regions of North America prior to 1900 (Nowak, 1979(Nowak, , 2002Young & Jackson, 1951), hereafter referred to as the coyote historical range, coyotes have substantially expanded their geographic range over the last century to occupy every continental US state and Canadian province (Hody & Kays, 2018). Here, we focus on the eastward expansion across the midwestern US and southeastern Canada, culminating along the eastern seaboard.
This expansion began in the early 20th century and followed two broad expansion routes across distinct environments. In the northeast, coyotes expanded across the Great Lakes region of the United States and Canada into New England, New York, and Pennsylvania.
The southeastern expansion occurred at a slower rate and followed an approximate trajectory through Louisiana, Alabama, and Georgia, with initial reports of coyotes in the Carolinas as recently as the 1980s (DeBow, Webster, & Sumner, 1998). Though fine-scale variation in expansion routes has been documented for the northeast (Kays, Curtis, & Kirchman, 2010;Wheeldon et al., 2010), a recent genetic survey (Heppenheimer et al., 2018) supports two genetically distinct eastern coyote populations across the eastern seaboard that correspond to these broadly described northeastern and southeastern expansion routes, suggesting that fine-scale expansion routes have likely converged.
In addition to geographic isolation, each expansion front represents distinct ecoregions in North America. For example, the northeastern expansion front is primarily northern forests and eastern temperate forests, which is further divided primarily into mixed wood shield, Atlantic Highlands, and mixed wood plains (Omernik & Griffith, 2014). In contrast, the southeastern expansion front is almost entirely eastern temperate forests, but transitions to tropical wet forest and great plains designations along the Gulf of Mexico (Omernik & Griffith, 2014). Furthermore, each expansion front also differs in the presence and abundance of closely related Canis species. Generally, under a range expansion scenario, hybridization between closely related and previously isolated species may occur as a result of low population density of the expanding species along the | 12643 HEPPENHEIMER Et al. range periphery (Seehausen, 2004). As such, coyote hybridization with remnant populations of eastern (C. lycaon) and/or gray wolves (C. lupus) has been documented along the northeastern expansion route (Kays et al., 2010;Rutledge, Garroway, Loveless, & Patterson, 2010;vonHoldt et al., 2011;vonHoldt, Kays, Pollinger, & Wayne, 2016), as well as with red wolves (C. rufus) in the southeastern expansion front (Nowak, 2002). In particular, red wolves are believed to be extirpated outside of the North Carolina recovery area, but hybridization between red wolves and coyotes is well documented within that area (Bohling et al., 2016;Hinton, Gittleman, Manen, & Chamberlain, 2018). Further, several previous studies have also shown that eastern coyote populations have interbred with domestic dogs (Adams, Leonard, & Waits, 2003;Wilson, Rutledge, Wheeldon, Patterson, & White, 2012;Wheeldon, Rutledge, Patterson, White, & Wilson, 2013;Monzõn, Kays, & Dykhuizen, 2014). While there is evidence that these hybridization events have been adaptive (vonHoldt et al., 2016), it is important to note that the full genome-wide consequences and the geographic extent of interspecies hybridization have not been documented throughout the entire eastern range.
Overall, our objectives were to quantify genomic structure and diversity across the historical coyote range and the two recently expanded eastern coyote populations. We then identify outlier loci that may have been under selection in both populations as a result of range expansion. We predict that population structure will correspond to the known demographic history of North American coyotes, that is, a historical range population and two distinct recently expanded groups. In accordance with theoretical assumptions, we expect reduced genomic diversity in the two recently expanded eastern populations relative to the historical range. However, hybridization with other Canis species may result in deviations from this expectation. Finally, we expect genomic variants that underwent frequency shifts in the same direction in both groups to have putative functions related to range expansion, such as dispersal and reproduction.

| Sample collection
We obtained coyote blood and tissue (e.g., liver, kidney, tongue) from state management programs (Princeton IACUC #1961A-13), government organization archives (e.g., Florida Fish and Wildlife, US Department of Agriculture, Ontario Ministry of Natural Resources and Forestry), or museum archives (New York State Museum, Oklahoma Museum of Natural History). In all cases, state or province of origin was documented (Figure 1), and in many cases, sex, approximate age, and fine-scale geographic data were also known. Samples were collected between 1998 and 2017, with the majority collected within the last 10 years (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017).
Samples with unknown collection dates were either known or assumed to fall within the approximate time period (Supporting information Table S1).
For downstream analyses, we considered samples collected from AZ, CA, ID, MN, MO, NE, NM, NV, OK, SK, TX, WA, and WY, to be part of the historical range (i.e., pre-1900; Figure 1) as described by Hody and Kays (2018). Additionally, samples collected from ME, NB, NJ, ON, and PA were considered part of the northeast expansion, and samples collected from AL, FL, GA, KY, LA, NC, SC, TN, and VA were considered part of the southeastern expansion.

| Sampling and DNA extraction
High molecular weight genomic DNA was extracted from all samples with either the Qiagen DNeasy Blood and Tissue Kit or the BioSprint 96 DNA Blood Kit in conjunction with a Thermo Scientific KingFisher Flex Purification platform, in both cases following instructions provided by the manufacturer. We quantified DNA concentration with either PicoGreen or Qubit 2.0 fluorometry and standardized samples to 5 ng/μl. Only high-quality DNA samples, as determined by the presence of a high molecular weight band on a 1% agarose gel, were retained for sequencing.

| RADsequencing and bioinformatics processing
We prepared genomic libraries following a modified restriction-associated DNA sequencing (RADseq) protocol described by Ali et al. (2016). Briefly, samples were digested with sbfI, and a unique 8 bp barcoded biotinylated adapter was ligated to the resulting fragments. Samples were then pooled (96-153 samples/pool) and randomly sheared to 400 bp in a Covaris LE220. Following shearing, we used a Dynabeads M-280 streptavidin bead binding assay to enrich for adapter-ligated fragments. Final sequencing libraries were then prepared using either the NEBNext Ultra DNA Library Prep Kit or the NEBNext UltraII DNA Library Prep Kit. Size selection was made for 300-400 bp insert with Agencourt AMPure XP magnetic beads, which were also used for library purification. Libraries were standardized to 10 nM and sequenced (2X150nt) on two lanes on the Illumina HiSeq 2500.
As this RADseq protocol is unique in that the barcode may be on either the forward or reverse read, data processing was required prior to variant calling. Accordingly, forward and reverse raw sequencing reads were processed such that any read containing the remnant sbfI cut site and one of the possible barcodes were aligned in a single file, while the matching read pairs that lacked the cut site were aligned in a separate file, and all remaining reads were discarded. This was accomplished using a custom Perl script (flip_trim_ sbfI_170601.pl, see Supporting information).
First, reads were demultiplexed using process_radtags, allowing a 2 bp mismatch for barcode rescue and discarding reads with either uncalled bases or a low-quality score (<10) within a sliding window of 0.15. Next, PCR duplicates were removed with the paired-end sequencing filtering option in clone_filter. We excluded all samples with low read count (<500,000) after removal of PCR duplicates from further analysis. Remaining samples were then mapped to the dog genome CanFam3.1 assembly (Lindblad-Toh et al., 2005) in Stampy v 1.0.21 (Lunter & Goodson, 2011). To reduces biases associated with incorrectly mapped loci (e.g., paralogs), we additionally filtered mapped reads for a minimum MAPQ of 96 and converted to bam format in Samtools v 0.1.18 (Li et al., 2009).
We completed SNP calling in STACKS following the recommended pipeline for reference mapped data (i.e., pstacks → cstacks → sstacks → populations; Catchen et al., 2013). In pstacks, we required a minimum depth of coverage of three to report a stack (-m 3). Cstacks and sstacks were run as recommended by Catchen et al. (2013). The populations module was run twice to optimize final sample selection to reduce both missing data and biases resulting from uneven sampling across locations. First, we allowed only the first SNP per locus (--write_single_snp) to be reported but did not apply any missing data thresholds. We then evaluated the per-individual genotyping success, as measured by missingness per individual in Plink v 1.90b3i (Purcell et al., 2007), and removed individuals with a high level of missing data (>80% missing). We also removed samples from locations where n = 1.
In our second run of populations, we included only this reduced set of samples, required that reported loci be genotyped in 90% of individuals (−r = 0.9), and again restricted analysis to only the first SNP per locus. Only this latter dataset was used for subsequent analyses. Following SNP calling, we filtered for statistical linkage disequilibrium in Plink with the argument --indep-pairwise 50 5 0.5.
All SNPs were annotated as genic (intron or exon), within a promoter (i.e., within 2 Kb of transcription start site following von-Holdt, Heppenheimer, Petrenko, Croonquist, and Rutledge (2017)), or intergenic using an in-house python script (chr_site.py; See supporting information). All intergenic SNPs were compiled in a second genotype dataset and filtered for Hardy-Weinberg Equilibrium (HWE) in Plink with the argument-hwe 0.001. These intergenic, HWE-filtered genotypes were presumed neutral in downstream analyses (hereafter, putatively neutral loci). Additionally, we compiled a third dataset of putatively functional loci consisting of all SNPs annotated as genic or within 2 Kb of a transcription start site (hereafter, genic loci).

| Population structure analysis
To visualize clustering in our data and identify strong outliers, we conducted a Principal Component Analysis using our full SNP dataset with flashPCA (Abraham & Inouye, 2014). We identified one strong outlier originating from Ontario, which may be a misidentified eastern or gray wolf (C. lycaon or C. lupus). This individual was removed from further analyses.

NM NM
Following the removal of strong outliers (e.g., putatively misidentified wolves), we determined the most likely number of genomic clusters represented by the data, by conducting an analysis of population structure in ADMIXTURE v1.3. (Alexander, Novembre, & Lange, 2009) with the cross-validation flag. We evaluated K = 1-10, with the K value with the lowest cross-validation (cv) score indicative of the best fit K. ADMIXTURE is similar in principle to the classic Bayesian software STRUCTURE (Pritchard, Stephens, & Donnelly, 2000), but uses a maximum likelihood framework and is more computationally efficient for SNP data.

| Identification of loci associated with range expansion
To identify loci as candidates for selection during range expansion, we restricted analyses to individuals with high cluster assignments as identified in the ADMIXTURE analysis (Q ≥ 0.8). Additionally, to prevent biases resulting from long-distance dispersers, we removed three individuals that had high cluster assignments to different populations from which they were sampled (e.g., sampled from Louisiana but clustered with the northeast group). Though restricting the analyses to individuals with high cluster assignments may inflate the genetic distinction between the historical range and the recently expanded coyote populations, we believe the inferred historical range population to be the best available representation of pre-range expansion allele frequencies, as there is likely ongoing contemporary gene flow between coyotes in the historical range and both recently expanded populations.
As F ST outlier-type approaches to identify loci under selection are prone to high rates of false positives, especially among populations with complex demographic histories, we used two distinct approaches to identify loci putatively under selection. First, a Bayesian framework to detect outlier loci was implemented in BAYENV2 (Coop, Witonsky, Rienzo, & Pritchard, 2010;Günther & Coop, 2013 PCadapt is less prone to type II error than alternative methods (e.g., BayeScan) and PCadapt is expected to perform well under a variety of complex demographic scenarios, including range expansion (Luu et al., 2017).
In the PCadapt analysis, the northeastern and southeastern expansion fronts were analyzed separately. In both cases, we identified outliers with respect to underlying population structure given by PC1. We chose to retain only PC1, rather than selecting the optimal number of PCs based on conventional methods (e.g., scree plot), as PC1 primarily captured the major axis of range expansion and therefore corresponded to the level of population structure relevant for this study (See Figure 2). To evaluate significance, p-values for each SNP were transformed into q-values and SNPs with q-values <0.05 were retained, therefore controlling for a false discovery rate (FDR) of 5%. This was implemented in the R package qvalue v2.6 (Bass, Dabney, & Robinson, 2018). Again, we only considered SNPs that were significant in both historical and northeast and historical and southeast comparisons as candidates for loci under selection during range expansion.
Gene functions and gene ontology biological process annotations of all outlier SNPs identified by both the analyses methods were inferred using Ensembl (release 91; Zerbino et al., 2017) and AmiGO 2 accessed in February 2018 (Carbon et al., 2009). We conducted a gene ontology (GO) biological process overrepresentation enrichment analysis on outlier sites located in functional genomic regions (i.e., intron, exon, or promoter) using WebGestalt (Zhang, Kirov, & Snoddy, 2005;Wang et al., 2013). We used our genic SNP dataset as the reference set for the enrichment analysis, and significance was evaluated using an FDR threshold of 5%. Additionally, we used the Ensembl Variant Effect Predictor (McLaren et al., 2016) to predict the functional effect of all outlier SNPs. (Supporting information Figure S1). Average per-individual missingness was highest in the northeast (mean missing = 10.5%) and

| Genotyping
slightly lower in both the southeast (6.0%) and the historical range (5.4%; Supporting information Figure S2). Average depth of coverage across all individuals was 11.333 with a standard deviation of 6.626 and was similar among the three regions surveyed (historical range: 11.420, stdev = 5.647; northeast: 11.634, stdev = 9.991; southeast: 11.096, stdev = 4.988; Supporting information Figure   S3). Furthermore, overall allele balance for all heterozygotes (i.e., Overall, we primarily captured rare variation, with an average global minor allele frequency of 1.73%. Further, within each of the three regions sampled, minor allele frequencies were typically ≤5% (Supporting information Figure S4). Approximately, half of the sites were intergenic (n intergenic = 12,676), with 14,108 SNPs found within genes (n intron = 12,024; n exon = 1,532; n promoter = 552). These annotations sum to >22,935 as SNPs may have multiple annotations (e.g., promoter and intron). Additionally, our putatively neutral dataset, which consisted of intergenic SNPs in HWE, retained 11,518 SNPs.
Our genic data included 10,259 SNPs within introns, exons, and promoters.

| Population structure corresponds to expansion axis
Our PCA divided sampling locations as predicted, with PC1 (   Figure 4a). Private allele counts (Table 1) (Figure 5a). These SNPs were primarily intronic (n intron = 45), with six exonic SNPs, and two located in putative promoter regions (Table 2; Supporting information Table S5).

| Outlier loci associated with range expansion
With the PCadapt approach, 59 SNPs were significant outliers (FDR 5%) in both the northeast expansion and historical range and southeast expansion and historical range analyses (Figure 5a). Of these, 22 SNPs were within genes (n intron = 20; n promoter = 2) and 37 were intergenic (Table 2; Supporting information Table S5). Though these two analyses methods to identify outlier SNPs are similar, BAYENV2 is based on a priori defined groups (i.e., sampling location), while PCadapt is based on principal component scores without predefined groups. However, as PC1 was highly correlated with the expansion axis (see Figure 2b), there was no discordance between placement along PC1 and the categorization of samples as either historical or recently expanded. Accordingly, we consider this analysis to be directly comparable and focus our results and discussion on SNPs and genes identified by both analyses (Lotterhos & Whitlock, 2015).
A total of twelve SNPs of 22,935 were outliers in both the outlier analyses (Table 2; Figure 5a). In all cases, the change in allele frequency relative to the historical range was in the same direction for both recently expanded eastern populations (Figures 5b   and 6). For nine of these outlier loci, the less common allele overall (i.e., the global minor allele), decreased in the recently expanded populations, and for the remaining three loci, the minor allele increased in frequency (Figures 5b and 6). In one case, WDR17, the minor allele was lost in both of the recently expanded eastern populations ( Figure 6). For three additional outlier loci (PAX5, EPHA6 and CARMIL1), the minor allele was lost in one of the recently expanded populations and substantially reduced in frequency in the other (Figure 5b). Furthermore, there was only one locus (ZDHHC16) where the minor allele was absent from the historical range, but present at appreciable frequencies, in both recently expanded populations (q Northeast = 12.12%; q Southeast = 24.09%; Figure 6).
In our GO enrichment analysis, these outlier SNPs were not significantly enriched for any biological process (Supporting information Table S6) after applying an FDR threshold of 5%. Furthermore, all sites were annotated as "modifiers" in the VEP analysis, which are defined as variants with no predictable functional effects on coding regions (i.e., noncoding variants).

| D ISCUSS I ON
Recently expanded coyote populations in eastern North America provide a unique opportunity to explore how range expansion shapes genomic diversity at neutral and putatively adaptive loci.
Here, we identified three genetic groups of coyotes, which largely correspond to the historical range and two distinct expansion fronts. Instances of discordance between sampling location and cluster assignments were relatively rare, and likely due to recent shared ancestry as well as ongoing gene flow. In particular, coyotes from OK, NE, and MN exhibited intermediate assignments to the southeastern cluster, and coyotes from MO exhibited intermediate assignments to both the southeastern and northeastern clusters. With the exception of MN, this midwestern region has previously been suggested to represent the source population for the southeastern expansion (Nowak, 1979;vonHoldt et al., 2011).
Additionally, as coyotes are highly mobile, there is likely ongoing gene flow among the recently expanded populations and those in the historical range, particularly along the eastern extreme.
To directly address the relative impacts of shared ancestry and ongoing gene flow in the eastern historical range, a more extensive sampling scheme throughout this region is needed and pre-range expansion samples from prior to 1900 should also be included. Sample size (g)

Private allelic richness
Historical range Southeast expansion Northeast expansion relative to their source populations Swaegers et al., 2015). In contrast, we observed approximately equivalent heterozygosity between coyote populations in the historical and recently established northeastern ranges, and more intriguingly, we observed that coyotes in the southeastern US had slightly greater heterozygosity than those in the historical range. As this trend is consistent across various subsections of the genome (i.e., genic and putatively neutral regions), it is not immediately clear from this study what is driving this lack of reduction in genomic diversity in the recently established eastern populations. However, the observed trends for private allele counts and private allelic richness in coyotes are consistent with a bottleneck scenario, as the historical range was observed to have the highest number of private alleles and the highest private allelic richness. As the loss of rare alleles will have a greater impact on allelic richness than heterozygosity, (Greenbaum et al., 2014), allelic richness has been suggested to be a better indicator of a bottleneck, particularly following range expansion.  (Kabitzke et al., 2017). Additionally, ALK knockout mice exhibit enhanced performance in novel object recognition tests, suggesting that this gene plays a role in the ability of animals to explore a novel environment (Bilsland et al., 2008) Notes. Chr, Chromosome; BF, Bayes Factor; prom, promoter.
While these finds are intriguing, the implications for how mutations in these outlier loci impact dispersal capabilities in coyotes should be interpreted with extreme caution. None of these outlier SNPs were found within an exon, and it is unclear from this data whether the genotyped variants directly impact gene expression, or whether these variants are in linkage disequilibrium with one or more nonsynonymous mutations in coding regions. Further, evidence for a dispersal-related function of these genes is entirely based on mouse studies, and it is unknown if these functions are conserved across mammals. We also did not incorporate any behavioral data for coyotes in this study, and it is unclear if eastern coyotes exhibit differences in exploratory behavior relative to coyotes from the historical range, much less if this behavior is correlated with genotype. Future studies should target the coding sequence of these genes in coyotes to further elucidate how these outlier loci may influence phenotypic traits in expanding coyote populations.
It is additionally possible that our outlier detection approach identified loci that are systematically different between both expansion fronts and the historical range as a result of parallel adaptations to similar environments rather than the range expansion process. While the northern and southern expansion fronts are distinct ecoregions (Omernik & Griffith, 2014), environmental similarities between the regions do exist, most notably in the high abundance of deer. However, diet studies reveal that deer consumption varies widely among eastern coyote populations (Kilgo, Ray, Ruth, & Miller, 2010;Mastro, 2011;Robinson, Diefenbach, Fuller, Hurst, & Rosenberry, 2014) and that deer consumption is also reasonably common throughout the historical range (Ballard, Lutz, Keegan, Carpenter, & deVos Jr, 2001;Carrera et al., 2008;Gese & Grothe, 1995). As such, selection associated with the range expansion process is perhaps more likely than adaptation to deer rich environments, though the possibility remains that outlier SNP frequencies are driven by selection associated with unmeasured environmental variables rather than by range expansion.
One additional outlier locus, ZDHHC16, which has putative functions related to eye development, cellular response to DNA damage, heart development, and protein palmitoylation, was monomorphic in the historical population yet polymorphic in both recently expanded populations. There are three general explanations for the origin of this variant in the recently expanded eastern coyote populations: (a) The mutation was present in the historical range, but at an extremely low frequency that was not captured by our sampling, (b) a de novo mutation occurred, either convergently along both expansions fronts, or along one front and then was transferred via intraspecific gene flow, or (c) this allele introgressed from a closely related Canis species following a hybridization event and was subsequently transferred via gene flow. As discussed above, interspecies hybridization occurs among Canis species, and there is evidence that this has been adaptive in the context of coyote range expansion (Thornton & Murray, 2014;vonHoldt et al., 2016). It is therefore conceivable that ZDHHC16 represents an additional case of adaptive introgression. However, the data presented here are not sufficient to address questions related to the origin of genomic variants, though this remains an interesting question for future studies regarding the role of hybridization in facilitating range expansion.
Taken together, we present a comprehensive genome-wide survey of coyote populations across much of the contiguous US as well as southeastern Canada. Despite pronounced geographic structuring among the historical and two recently expanded eastern coyote populations, we did not observe a strong decline in genomic diversity that is characteristic of a range expansion bottleneck, suggesting that coyote range expansion dynamics are more complex than those described in theoretical  and other empirical studies (e.g., White et al., 2013;Swaegers et al., 2015), and is likely attributable to interspecies hybridization.
Further, we identify several genomic variants that are candidates for gene regions under selection during range expansion, which Fish and Wildlife Service.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.
F I G U R E 6 Allele frequencies for all twelve outlier SNPs across sampling locations. Allele one was arbitrarily defined as the most common allele overall (i.e., the major allele)