Development of a genotype‐by‐sequencing immunogenetic assay as exemplified by screening for variation in red fox with and without endemic rabies exposure

Abstract Pathogens are recognized as major drivers of local adaptation in wildlife systems. By determining which gene variants are favored in local interactions among populations with and without disease, spatially explicit adaptive responses to pathogens can be elucidated. Much of our current understanding of host responses to disease comes from a small number of genes associated with an immune response. High‐throughput sequencing (HTS) technologies, such as genotype‐by‐sequencing (GBS), facilitate expanded explorations of genomic variation among populations. Hybridization‐based GBS techniques can be leveraged in systems not well characterized for specific variants associated with disease outcome to “capture” specific genes and regulatory regions known to influence expression and disease outcome. We developed a multiplexed, sequence capture assay for red foxes to simultaneously assess ~300‐kbp of genomic sequence from 116 adaptive, intrinsic, and innate immunity genes of predicted adaptive significance and their putative upstream regulatory regions along with 23 neutral microsatellite regions to control for demographic effects. The assay was applied to 45 fox DNA samples from Alaska, where three arctic rabies strains are geographically restricted and endemic to coastal tundra regions, yet absent from the boreal interior. The assay provided 61.5% on‐target enrichment with relatively even sequence coverage across all targeted loci and samples (mean = 50×), which allowed us to elucidate genetic variation across introns, exons, and potential regulatory regions (4,819 SNPs). Challenges remained in accurately describing microsatellite variation using this technique; however, longer‐read HTS technologies should overcome these issues. We used these data to conduct preliminary analyses and detected genetic structure in a subset of red fox immune‐related genes between regions with and without endemic arctic rabies. This assay provides a template to assess immunogenetic variation in wildlife disease systems.


| INTRODUCTION
Populations are exposed to spatially heterogeneous selective pressures that often lead to localized patterns of adaptation, where local individuals are expected to have higher fitness than those from other populations from different environments (Thompson, 2005). In this context, characterizing patterns of local adaptation illustrates the spatial and temporal interactions of species with their environment, describes processes that generate and maintain biodiversity, and identifies evolutionary forces shaping populations (Carlson, Cunningham, & Westley, 2016;Harrisson, Pavlova, Telonis-Scott, & Sunnucks, 2014).
Host responses to pathogens are largely influenced by immune genes, where genetic diversity influences population-level resistance to disease via pathogen-mediated directional or balancing selection (Eizaguirre et al., 2012;Rico et al., 2015). Studies of wild populations generally focus on adaptive immunity, assessed via a small number of major histocompatibility complex (MHC) loci, due to their role in pathogen recognition and pathogen susceptibility (Acevedo-Whitehouse & Cunningham, 2006;Eizaguirre et al., 2012;Kyle et al., 2014). A growing number of studies also examine key receptor genes associated with innate immunity (e.g., toll-like receptors and interleukins) to clarify patterns of resistance to emerging infectious diseases (e.g., Bonneaud et al., 2012). The aforementioned loci, however, provide insight into only a small fraction of the genes associated with mounting an immune response, disease resistance, and adaptation. Overall, a more holistic assessment of immunogenetic variation is required to understand the capacity for adaptation to disease (Harrisson et al., 2014;Morris, Wright, Grueber, Hogg, & Belov, 2015;Ng et al., 2009).
High-throughput sequencing (HTS) technologies allow for expanded explorations of genomic variation associated with local adaptation; however, for population-level assessments of genetic variation in wildlife with larger genomes, full-genome processing is generally not yet feasible (Ekblom & Wolf, 2014). For most non-model organisms, large chip arrays of single-nucleotide polymorphisms (SNP) are not available, and those are maybe prone to ascertainment bias when applied to new populations (Albrechtsen, Nielsen, & Nielsen, 2010;Lachance & Tishkoff, 2013). As an alternative, genotype-bysequencing (GBS) assays can be employed to obtain genomic subsets of population variation. GBS largely falls into three categories: restriction enzyme-, amplicon-, and hybridization-based models, each coming with important considerations for implementation (Jones & Good, 2016). For example, restriction site associated DNA (RAD) markers are well-suited for detecting neutral sequence variation across the genome to develop models of genetic population structure (Catchen et al., 2017), but there are some controversies as to their applicability in detecting patterns of local adaptation (Lowry et al., 2017). In more targeted GBS approaches, amplicon-based and target hybridizationbased methods have been shown to be highly effective in enriching for subsets of genomic sequence. Amplicon-based methods rely on the ability to PCR amplify large numbers of targeted, overlapping short-sequences, using multiple sets of PCR primers in a single PCR (Samorodnitsky et al., 2015). Target capture utilizes overlapping biotinylated DNA or RNA probes, which bind to complementary targeted regions of DNA that are then selectively pulled down by magnetic streptavidin beads to enrich for target DNA (Chilamakuri et al., 2014;Ng et al., 2009). Both techniques, when coupled with HTS, enable the identification of genetic variation across specific loci in the genome.
Targeted sequence capture approaches have only recently been applied to wildlife systems. For example, Morris et al. (2015) used genome-level information from ten Tasmanian devils to identify SNPs within the "immunome," consisting of the coding and regulatory regions of 167 immune genes. Using this information, they developed an amplicon-based assay for nine immune genes with nonsynonymous SNPs, and then, genotypes were generated for 220 Tasmanian devils in the context of a transmissible facial cancer decimating the species.
Similarly, Schweizer, Robinson et al. (2016) employed an RNA-bait version of target capture to enrich for 1,040 candidate genes and their regulatory regions, as well as 5,000 1-kbp nongenic neutral regions, to determine genetic variation in 107 gray wolves from diverse ecotypes. In light of recent technological advances, multiple GBS-based strategies have emerged with the potential to study nonmodel wildlife species under various selective pressures.
Rabies viruses cause fatal encephalopathies in mammals with strains adapted to infect different primary hosts (Jackson & Wunner, 2007). In Alaska, arctic rabies (AR) is typically restricted to arctic coastal areas following arctic fox (Vulpes lagopus) distributions, where AR variants have discrete spatial distributions (Kuzmin, Hughes, Botvinkin, Gribencha, & Rupprecht, 2008;Nadin-Davis, Sheen, & Wandeler, 2012). Distinct phylogeographic patterns have been observed between AR variants (Kuzmin et al., 2008;Nadin-Davis et al., 2012), suggesting variable disease resistance may influence the spread of rabies. Three AR variants occur in Alaska that are isolated to regions with tundra red foxes (Vulpes vulpes) and arctic foxes, with nonoverlapping, temporally stable distributions, although there are no data regarding the relative virulence of these strains. The AR variants occur on the North Slope (NS; AR variant 3), Seward Peninsula (SP; AR variant 2), and southwestern (SW; AR variant 4) regions of Alaska ( Figure 1).
Conversely, in interior and southcentral (SC) Alaska, only boreal red foxes are found, and AR is not endemic (Goldsmith et al., 2016). Red fox, a highly susceptible AR host, has been associated with AR spreading from Arctic regions to southern Canada (Kuzmin et al., 2008). In the Canadian Arctic, one AR variant predominates, including in regions where arctic foxes are absent and red fox are present. Therefore, it is unclear why rabies has not spread via red fox into central Alaska given the patterns observed across Canada (Mørk & Prestrud, 2004;Nadin-Davis et al., 2012). In red fox, recent experiments detected weak genetic structure in Alaska between tundra and boreal regions consistent with spatial distributions of AR variant presence and gene flow between those populations (Goldsmith et al., 2016). Questions remain, as to how DNA sequence diversity at functionally important immune-related loci for red fox is spatially structured in these regions, and if geographic patterns of diversity are correlated with AR strain distribution.
Herein, we describe the development of an immunogenetic assay designed to test the hypothesis that selective pressure from disease alters the spatial distribution of genetic variants among populations.
Our goals were to (1) develop an assay capable of elucidating genetic variation at red fox immunogenic loci and (2) use the resulting SNPs to apply preliminary tests for genetic structure and F ST outliers. We predicted that even with low sample sizes, we would be able to detect specific genetic variants based on AR presence/absence because AR exerts a strong selective pressure on fox populations. To develop the assay, we annotated immune-related red fox gene sequences and then applied high-throughput targeted sequence capture to enrich for 300-kbp of preselected genomic regions (per individual) that could be used to identify SNPs in a wide range of key immune genes and their regulatory regions, which can influence the expression of immune system genes and disease outcome (Fraser, 2013). To test how robust our assay was, we assessed whether 1 ng DNA/sample could be used for target capture-based GBS to evaluate the potential of this technique for noninvasive samples with less DNA (such as from hair or feces), and whether the red fox hybridization probes could also be used to enrich for similar targets from an arctic fox DNA sample for future AR research. To address our hypothesis and elucidate patterns of local adaptation, we required an understanding of the underlying patterns of gene flow and genetic drift that may influence the presence/absence of genes beyond selective pressures. As such, we also included 23 neutral microsatellite loci and amelogenin (for sex determination) in this assay to provide insight into the demographic processes of the populations under investigation, thus providing an "all-in-one" test to generate data for our research aims. This study provides a template for future experiments aimed at identifying immunogenetic variants and genetic patterns of local adaptation to disease.

| Sample collection, DNA extraction, and quantification
We obtained red fox tissue (n = 33) or extracted DNA samples from rabies-positive red fox (n = 8) as described by Goldsmith et al. (2016) from a variety of organizations and individual trappers. The rabiesnegative red fox tissue samples have voucher specimens in the University of Alaska Museum of the North, and accession numbers are provided, where applicable (Table S1). Rabies-positive red fox tissue was confirmed at the Alaska government health laboratories by serology (ELISA) and then at the Centres for Disease Control by sequencing. The arctic fox DNA sample (n = 1) was obtained from Terry Spraker (Colorado State University, USA). We dissolved all F I G U R E 1 Schematic of Alaska arctic fox and red fox samples analyzed using immunogenetic profiling. Approximate arctic rabies variant (ARV) distribution was modified from Goldsmith et al. (2016). Arctic rabies and arctic rabies-free zones are indicated by blue and white background colors, respectively. Arctic fox and red fox sample locations are denoted by an open-green circle or open-red circles, respectively. Red fox sample size (n) for each location used in the genetic structure analyses is indicated. SW, Southwest; SP, Seward Peninsula; NS, North Slope tissue samples in 1× lysis buffer (4 mol/L urea, 0.2 mol/L NaCl, 0.5% n-lauroyl sarcosine, 10 mmol/L 1,2-cyclohexanediaminetetraacetic acid, 0.1 mmol/L Tris-HCl pH 8.0) containing 600 U/ml proteinase K (time and temperature) at 56°C for two hours and extracted DNA using either the automated 96-well MagneSil Blood Genomic Max Yield System (Promega) or the DNeasy Blood and Tissue Kit (Qiagen).

| Identification of red fox immune genes, neutral markers, and probe development
We compiled a list of candidate immune genes (Table S2) using those listed in The Dog Innate & Adaptive Immune Responses RT 2 Profiler PCR Array (Qiagen) and in the overview of key mediators of innate and adaptive immunity, development, and signaling (Knight, 2013).
Using this candidate gene list, we downloaded genomic sequences (covering an entire locus), individual exon sequences, and protein sequences from the CanFam3.1 dog (Canis lupus familiaris) genome assembly (Hoeppner et al., 2014) using Ensembl release version 79 (Cunningham et al., 2015).
We used the dog genomic sequences to perform blastn (Altschul, Gish, Miller, Myers, & Lipman, 1990) against the V. vulpes draft genome (Kukekova et al., in review) to identify segments of the red fox genome assembly (NCBI BioProject PRJNA378561) with sequence similarity. The red fox genome assembly consists of 2,495,544,672 assembled bp and 676,878 scaffolds (N50 = 11.5 Mb). We extracted red fox immune-gene sequences (2,229,152 bp) from the red fox genome assembly that included upstream and downstream sequences (±1,500 bp). We annotated the red fox immune-gene structures (introns and exons) using dog exon sequences by performing blastn (Altschul et al., 1990) against the red fox DNA sequences with a wordsize parameter set to 7 to include short sequence hits. We input the red fox immune-gene sequences and dog peptide sequences into fgenesh+ (Solovyev, 2007) on the softberry online portal (http://www.softberry. com/berry.phtml) to resolve cases where blastn did not properly annotate the red fox intron/exon boundaries, which led to incomplete open reading frames (start codon to stop codon). We performed blastp (Altschul et al., 1990) against the NCBI refseq_protin database to verify all red fox protein-coding sequences. The red fox gene annotations and predicted red fox immune-related protein sequences are available in Appendices S1 and S2, respectively.
We also generated a list of canine and red fox microsatellite markers from previous studies Fredholm & Winterø, 1995;Mellersh et al., 1997;Ostrander, Sprague, & Rine, 1993). We downloaded the red fox microsatellite region from NCBI when available. Otherwise, we extracted the genomic sequences from the CanFam3.1 dog assembly and used blastn to identify and extract regions of putative sequence similarity in the red fox genome. This yielded another set of red fox DNA sequence targets (Table S3) that were included in the custom NimbleGen SeqCap EZ probe design (Roche).
We added 100-bp "padding" to each of the red fox DNA sequence targets (primary targets) to increase the efficiency of sequence capture. The custom NimbleGen SeqCap EZ probes were produced using the red fox DNA sequence targets. We only had limited access to the red fox genome assembly (Kukekova et al., in review); therefore, the probe set was compared to the dog genome reference sequence to test probe specificity. We finalized a "relaxed" probe design that allowed up to 20 close matches to the reference dog genome; however, 91.5% of the probes had only one match to the reference dog genome sequence, and 95.8% had five or fewer matches (not including zero) to the reference dog genome sequence, indicating a low likelihood for "off-target" sequence capture.

| Library preparation, sequence capture, and high-throughput sequencing
We prepared DNA libraries using the Kapa HTP Lib Prep Kit (Roche).
Forty-two DNA libraries were prepared using 1 μg DNA, and an additional three DNA libraries were prepared for the dilution series using 1, 10, and 100 ng DNA from sample KH1354. Special considerations were required during library preparation based on the quantity of DNA put into the reaction (Table S4)

| Sequence alignment and variant annotation
We used the bwa-mem command in the burrows-wheeler aligner v0.7.12 (bwa; Li, 2013) to align the paired-end reads to the red fox immune-gene sequences. We compiled sequence alignment metrics using samtools v1.2 (Li et al., 2009

| SNP/INDEL filtering and analyses
For SNP/INDEL analyses, we removed the arctic fox sample library (KH1527), the three "dilution" sample libraries (KH-1354-1, KH-1354-10, KH-1354-100), and one sample library with <10× coverage (KH989). Further, we removed the eight rabies-positive samples because we did not have enough samples from any one area (3 SP, 1 NS, 4 SW) to provide the sample size required for the preliminary downstream analyses presented in our study. We include the sequencing results from these rabies-positive samples as a future resource for researchers that may be able to use these data. Therefore, SNP/INDEL data from the 32 remaining sample libraries ( Figure 1; Table S1) were used in the preliminary downstream analyses to help inform future experiments. We generated separate SNP and INDEL subdatasets containing genetic variants for those 32 red fox using gatk.
We used the tandem repeat finder trf v4.09b (Benson, 1999) to annotate the INDEL subdataset corresponding to microsatellites and amelogenin and then vcftools to extract the genotypes of each individual for the microsatellite loci and amelogenin. We then used the python script "getgenosfromvcf" (De Wit et al., 2012) to extract microsatellite genotypes with a Phred quality score cutoff of 20, which yields genotypes with a 99% probability of being true. We created a filtered INDEL subdataset containing genetic variants for 15 microsatellite regions by removing ambiguous microsatellite calls with multiple sequence alignments. We also tested for significant departures from  (Narum & Hess, 2011). Therefore, we used three different outlier detection programs to mitigate the number of false positives in the F ST outlier subdataset we used for preliminary analyses. We used bayescan v2.1 (Foll & Gaggiotti, 2008), lositan (Antao, Lopes, Lopes, Beja-Pereira, & Luikart, 2008;Beaumont & Nichols, 1996), and outFLANK (Whitlock & Lotterhos, 2015) to identify F ST outliers putatively under selection from the filtered SNP subdataset. The three outlier detection programs do not rely on a set of "presumed" neutral loci to generate an empirical null distribution of F ST . Rather, they simulate a null distribution of F ST for the sample sizes observed in the dataset and identify departures from neutrality using different analytical approaches (Antao et al., 2008;Beaumont & Nichols, 1996;Foll & Gaggiotti, 2008;Whitlock & Lotterhos, 2015). We ran bayescan using 1:10 prior odds for the neutral model and included 20 pilot runs consisting of 5,000 iterations each, followed by 550,000 iterations with a burn-in of 50,000 iterations. We ran lositan for 1,000,000 iterations within a 99.5% confidence interval with the "Neutral mean F ST " and "Force mean We performed principal component analysis (PCA) on the filtered SNP subdataset, INDEL subdataset, and F ST outlier subdataset using adegenet v2.0.0 (Jombart & Ahmed, 2011). We obtained the required "genlight" objects for the adegenet analysis using a combination of vcftools and plink v1.07 (Purcell et al., 2007) to reformat the .vcf files to plink formatted files (.raw). We also tested for genetic structure using the F ST outlier subdataset in structure. We used the strauto v1.0 script (Chhatre & Emerson, 2017) to run structure over multiple processors at the same time. We ran structure with a burn-in length of 50,000 followed by 200,000 iterations for K = 1 through 8, and each run was performed 20 times. We calculated the ΔK statistic (Evanno, Regnaut, & Goudet, 2005) to help determine the number of inferred genetic clusters using structure harvester web v0.6.94 (Earl & von-Holdt, 2011). We used the LargeKGreedy (10,000 repeats) algorithm in clumpp v1.1.2 (Jakobsson & Rosenberg, 2007) to combine the SNP results from the multiple structure runs, and we visualized those results using distruct v1.1 (Rosenberg, 2004).

| Identification of red fox immune genes, neutral markers, and probe development
We extracted regions of genomic sequence corresponding to 116 immune genes and their upstream regulatory regions (Table S2), 23 microsatellite regions, and a portion of intron 1 in the amelogenin sex-determining marker (Table S3) from the red fox genome sequence assembly (Kukekova et al., unpublished data). Using the red fox immune-gene sequences, we annotated the immunity gene intron/ exon boundaries (Appendix S1). We validated the protein-coding annotations by translating all 116 nucleotide coding sequences to proteins (Appendix S2) and performing BLASTp (Table S2). We generated full-length coding sequence information for 109 immunity-related genes (excluding: CD86, DLA-12, DLA-88, IL18, MAPK1, MAPK8, and MX1). We successfully extracted 1,500-bp of upstream genomic sequence representing putative regulatory regions for 93 genes and partial (<1,500-bp) regulatory region sequences for an additional 15 genes. Our red fox immune-gene sequence database lacked the upstream sequence necessary to target the regulatory regions for CD28, CD86, DLA-12, HLA-DPB1, IKBKG, IL18, MAPK1, and MX1 (Table   S2). Those regulatory regions were either missing from, or spanned, multiple scaffolds in the red fox genome assembly and could not be included in our immune-gene sequence database. Using these annotations, the custom NimbleGen SeqCap EZ Developer Library probe set was synthesized.

| Library preparation, high-throughput sequencing, and sequence alignment
We obtained over 12.4 million paired-end reads from the HTS MiSeq run (

| INDEL detection and microsatellite analysis
Of the 32 libraries processed for downstream analyses, we detected a single bi-allelic INDEL for the amelogenin sex-determining region and detected 25 males and seven females (Table S1). There was a skewed sex ratio and a peak in coverage for the amelogenin locus.
Using independent field information (Table S1)

| SNP detection and analysis
We used the mapped reads from 32 libraries to generate a hard- rabies is present and the boreal southcentral regions where rabies is absent (Fig. S3b). We identified F ST outliers under directional selection using lositan and outFLANK (FDR < 0.05; Figure 2a; Table S8), but zero F ST outliers were identified using bayescan (FDR < 0.05; data not shown). Furthermore, 1,195 SNPs were identified as candidates under balancing selection by both mutation models in lositan (Table   S8). There were 15 F ST outliers predicted by both mutation models in lositan and outFLANK (Table 2). We used the genotypes from these 15 F ST outliers to visualize genetic structure by PCA and structure (Figure 2b,c). Preliminary plots revealed genetic structure of red fox between coastal tundra regions where arctic rabies is present compared to red fox from the boreal southcentral regions where rabies is absent.

| DISCUSSION
In this study, we annotated red fox immune genes and developed a GBS target capture protocol that can be used to elucidate spatial patterns of genetic variation at neutral and functional loci of red fox.
For the functional loci, we targeted genes associated with an innate, intrinsic, or adaptive immune response, as well as their regulatory segments that have been implicated in variable responses to disease exposure. We developed a system that accurately and evenly captured the targeted loci with sufficient coverage (mean = 50×; Fig. S2, Tables S5 and S6) for variant detection. Without a complete red fox genome assembly, we acknowledge the possibility of excess heterozygotes for a small proportion of our SNPs could be due to paralogous loci; however, our probe design tests suggested only one hit for most probes (91.5%). We used these data to perform a preliminary assessment of F ST outliers. The vast majority of outliers were found to be under balancing selection, with a smaller subset under directional selection (Tables 2 and S8). Despite moderate sample sizes, our preliminary analyses found notable differences in the frequencies of the F ST outliers in regions with and without AR. This assay provides a means to elucidate genetic variation from a large portion of the immunome with even coverage across samples and loci while alleviating some of the ascertainment bias of other genotype-by-sequencing approaches. AR. However, their study was limited to neutral marker analyses, and genetic diversity in arctic fox and red fox at functional loci that may be influenced by selective pressure from AR is unknown and the motivation to develop this GBS immunogenetic assay.

| Development of the red fox target capture genotyping-by-sequencing assay
As the cost of HTS continues to decline, strategies have emerged to gain insight into the genetic diversity of regions of the genome that encode proteins, which are presumably functionally relevant for local adaptation of wildlife species. For example, to assess the genetic basis for adaptation of arctic foxes to a cold climate, a comparison of transcriptome sequences from two captive arctic fox and one red fox identified two fat metabolism genes under positive selection in the arctic fox transcriptome (Kumar, Kutschera, Nilsson, & Janke, 2015).
However, isolating quality RNA from wildlife samples is problematic given the many environmental and individual variables associated with RNA expression, and the difficulty in obtaining fresh samples from a broad geographic range in regions that are difficult to access. Recently, GBS assays have emerged as a cost-effective alternative to wholegenome sequencing that aims to detect a genomic basis for population  (Funk et al., 2016). While useful in studies involving nonmodel species, RAD sequencing is not restricted to detecting F ST outliers in protein-coding regions of the genome. Therefore, targeted GBS approaches using either amplicon-or hybridization-based workflows have been employed to identify phenotype-altering mutations in both regulatory and protein-coding regions.
While amplicon-based methods may require more straightforward sample preparation and have the ability to utilize smaller DNA inputs, hybridization-based approaches have performed slightly better in sequencing complexity and uniformity with respect to target enrichment (Samorodnitsky et al., 2015). The decision to use either approach may be a matter of user preference, and both have been successful in wildlife applications. Using the amplicon-based approach, Morris et al. (2015) found low levels of polymorphism in the Tasmanian devil.
However, it is noteworthy that their amplicon design relied on the whole-genome sequence comparison of ten Tasmanian devils, which may be cost-prohibitive as a model for other wildlife investigations.
Moreover, sample bias may have affected their choice in amplicon gene targets. Target capture has been used to identify functional variation in gray wolf ecotypes in North America (Schweizer, Robinson et al., 2016); however, the gene targets were first identified by supporting information from an established Affymetrix v2 Canine SNP array that included ~42K SNPs (Schweizer, Vonholdt et al., 2016). In this study, we designed a custom target capture assay for red fox to screen for genetic variation in functional immune genes (Table S2), including their regulatory regions, and in neutral microsatellite regions (Table S3). The immune-gene annotations and their corresponding protein sequences (Appendices S1 and S2) helped validate a portion of the red fox genome assembly draft and are a resource for future immunome studies in red fox and related fox species. Our goal was to use these probes to screen for genomic signatures of local adaptation between tundra and boreal red fox from different AR regions in Alaska.

| Target capture is successful with low copy number wildlife DNA and across species
For GBS, Samorodnitsky et al. (2015) suggested higher amounts of input DNA required for hybridization-based Nimblegen SeqCap (1 μg), relative to amplicon-based Illumina AmpliSeq (50 ng), is a limiting factor. The difference in the amount of template required for these respective methods is of concern for wildlife studies, which can include T A B L E 2 Directional outliers detected by lositan and outFLANK analyses (FDR < 0.05) potential to generate GBS data using a single MiSeq run; however, optimization of such experiments may require more intense investigations.
We estimated a target enrichment of 61.5% and 59.9% (Table S5) and a mean target coverage of 50.1× and 37.0× (Table S6) for red fox and arctic fox, respectively. Similarly, cross-species exon target capture has been successful for a broad range of cichlids (Ilves & López-Fernández, 2014). Red fox and arctic fox recently diverged from one another ~3.2 mya during the late Pliocene (Kumar et al., 2015), and these data suggest that our existing probe set can be used in future studies aimed at addressing similar questions regarding genetic variation in arctic fox or related fox species.

| Analysis of microsatellites
Genotyping microsatellite regions using PCR-based techniques can be difficult because the polymerase can "slip" at STRs, leading to amplicons that differ in length. Additionally, STR-containing reads from HTS data can be difficult to accurately map to a reference genome due to mismatch/INDEL penalties associated with STR expansion and contraction (Fungtammasan et al., 2015). Therefore, we manually inspected the bi-allelic INDELs of our 23 microsatellite targets for poor alignments, which generated a filtered microsatellite dataset containing 15 loci. Microsatellite genotypes from dinucleotide STRs have been accurately called with a 90% success rate given a 17× depth of coverage (Fungtammasan et al., 2015), and our filtered subset had an average depth of coverage of 48×. We did not detect genetic structure in this microsatellite dataset using structure or adegenet PCA plots (Fig. S3a).
These results did not fully align with genotypes from the same samples generated using traditional microsatellite amplification and profiling using a genetic analyzer by Goldsmith et al. (2016). Our observed heterozygosity is also high compared to reported estimates of expected heterozygosity in Polish red fox (72%, Mullins et al., 2014) and that of other carnivores, including the American badger (81%; Rico et al., 2016) and Canadian black bear (55%-81%; Pelletier et al., 2017).

| Target capture analysis of immunogenetic diversity-considerations
We used lositan, bayescan, and outFLANK to detect candidate outlier SNPs putatively under selection. While lositan and outFLANK detected F ST outliers (Figure 2a; Tables 2 and S8), bayescan did not. Narum and Hess (2011) reported lower type I (false positive) and type II (false negative) error rates for lositan and bayescan compared to arlequin, but false positives for candidate SNPs under directional and balancing selection are abundant with all three approaches. Further, these methods have limited sensitivity in detecting SNPs under weak selection (Narum & Hess, 2011). We also detected 1,195 SNPs under balancing selection dispersed through exons, introns, and regulatory regions of the red fox immunome (Table S8). Signatures of balancing selection in the context of response to pathogens have been found in human immune-related genes, including those associated with the major histocompatibility complex (MHC; Andrés et al., 2009), and in the MHC of wolverine (Gulo gulo; Rico et al., 2015). MHC genes are highly polymorphic and play an important role in the adaptive immune response to pathogens. Therefore, our findings are consistent with previous studies in finding candidate immune-related loci under balancing selection. While these results are promising, they should be viewed with caution as we concede that there is the potential for bias in the outlier analyses because the initial dataset was not a random distribution of genome-wide loci that would include both neutral and adaptive alleles. We also acknowledge that our tests for directional and balancing selection likely contain biases associated with linked loci because all of the SNPs cannot be considered as independent markers.
Outlier SNPs in exons under directional selection that were detected by both lositan and outFLANK included C3, CRP, and TLR7 (Table 2). C3 encodes an activator in the complement system, which is involved in innate and adaptive immune responses and functions to lyse microorganisms, promote phagocytosis, trigger inflammation and aids in immune clearance (www.genecards.org). CRP recognizes foreign pathogens and promotes their elimination (www.genecards. org). TLR7 is a toll-like receptor that recognizes single-stranded RNA and activates the innate immune system (www.genecards.org). The candidate SNP in TLR7 may be an intriguing target for future studies because AR is a single-stranded RNA virus. The preliminary interpretation of our data suggests selective pressure at the molecular level on three immune genes in red fox. Overall, the SNP-based analysis found no evidence of population structure among red fox using all SNPs (Fig. S3b), but indicated population structure based on F ST outliers only, where red fox from coastal tundra regions were distinct from those in the boreal Southcentral regions (Figure 2b,c). To support these preliminary findings, future studies are required using larger sample sizes from Alaska red fox that compare rabies-negative and rabies-positive individuals from the same population and they could help reveal a genetic link to rabies susceptibility/resistance.

| CONCLUSION
Understanding the capacity for local adaptation to disease in wild populations requires expanded genomic assessments of the genetic responses to these selective pressures. Herein, we developed an immunogenetic assay that bridges between genetic and full-genome research and has the potential to generate empirical data that set the basis for predictive models to enhance our ability to anticipate epizootic disease spread and impacts under different climatic scenarios.

ACKNOWLEDGMENTS
This research was funded by discovery grants from the Natural Shafer for critical comments of our manuscript prior to submission; and two anonymous reviewers for critique leading to an improved version of our manuscript.

CONFLICT OF INTEREST
None declared.