Profiling the immunome of little brown myotis provides a yardstick for measuring the genetic response to white‐nose syndrome

Abstract White‐nose syndrome (WNS) has devastated populations of hibernating bats in eastern North America, leading to emergency conservation listings for several species including the previously ubiquitous little brown myotis (Myotis lucifugus). However, some bat populations near the epicenter of the WNS panzootic appear to be stabilizing after initial precipitous declines, which could reflect a selective immunogenetic sweep. To investigate the hypothesis that WNS exerts significant selection on the immunome of affected bat populations, we developed a novel, high‐throughput sequence capture assay targeting 138 adaptive, intrinsic, and innate immunity genes of putative adaptive significance, as well as their respective regulatory regions (~370 kbp of genomic sequence/individual). We used the assay to explore baseline immunogenetic variation in M. lucifugus and to investigate whether particular immune genes/variants are associated with WNS susceptibility. We also used our assay to detect 1,038 putatively neutral single nucleotide polymorphisms and characterize contemporary population structure, providing context for the identification of local immunogenetic adaptation. Sequence capture provided a cost‐effective, “all‐in‐one” assay to test for neutral genetic and immunogenetic structure and revealed fine‐scale, baseline immunogenetic differentiation between sampling sites <600 km apart. We identified functional immunogenetic variants in M. lucifugus associated with WNS susceptibility. This study lays the foundations for future investigations of rangewide immunogenetic adaptation to WNS in M. lucifugus and provides a blueprint for studies of evolutionary rescue in other host–pathogen systems.


| INTRODUCTION
Host-pathogen dynamics are changing at an unprecedented rate as climate change and human-mediated transport expand the range of pathogens into previously inhospitable/inaccessible environments (Fisher et al., 2012). As pathogen ranges shift, disease-related population declines in naïve wildlife populations often threaten population persistence, as evidenced by several emerging wildlife diseases Puechmaille, Wibbelt, et al., 2011). Controlled experiments with captive bats show that identical strains of P. destructans cause mortality in North American little brown myotis (Myotis lucifugus) but not in European greater mouse-eared bats (Myotis myotis; , suggesting a genetic basis for immunotolerance or immunoprotection. There are several promising leads for the development of treatments for WNS (e.g., Cheng et al., 2016;Cornelison, Gabriel, Barlament, & Crow, 2014;Wilcox & Willis, 2016), but no effective mitigation or treatment protocols are currently available. However, some populations near the epicenter of WNS may be stabilizing following their initial, precipitous declines (Dobony et al., 2011;Langwig et al., 2012Langwig et al., , 2017. Persistence of these populations does not seem to be associated with immigration (Maslo, Valent, Gumbs, & Frick, 2015), but may indicate evolution of resistance or tolerance to the disease (Langwig et al., 2017). Thus, adaptation and evolutionary rescue may be the best hope for recovery of bat populations affected by WNS (Maslo & Fefferman, 2015). Understanding patterns of immunogenetic adaptation to WNS is therefore critical to determining disease management strategies and recovery programs for the affected populations.
Immune genes mediate the initial response of individuals to pathogens and in many cases, the acquisition of immunity. At the population level, genetic diversity of immune genes influences resistance or tolerance to disease via pathogen-mediated balancing selection (Eizaguirre et al., 2012;Rico et al., 2015). Studies of wildlife populations generally focus on adaptive immunity, which has often been assessed by using genetic diversity in the major histocompatibility complex (MHC).
Diversity at the MHC provides a proxy for potential to adapt to shifting pathogen pressures, due to the role of MHC in pathogen recognition and pathogen susceptibility (Acevedo-Whitehouse & Cunningham, 2006;Eizaguirre et al., 2012;Kyle et al., 2014). Some studies of immunogenetic diversity also include receptor genes associated with innate immunity (e.g., Toll-like receptors and interleukins) and these markers have revealed spatial patterns of resistance to emerging infectious diseases such as chytridiomycosis and mycoplasmosis (Bonneaud, Balenger, Zhang, Edwards, & Hill, 2012;Savage & Zamudio, 2011).
The MHC Drb1 locus in M. lucifugus is among the most polymorphic recorded in mammals to date . Pyrosequencing of 160 individuals sampled across Canada suggests that balancing selection has maintained similar MHC diversity among genetically differentiated subpopulations, which may be disrupted by WNSmediated immunogenetic selection . However, the extreme observed polymorphism of the Drb1 locus in M. lucifugus is due partly to multiple gene duplications, which limits the use of these data. Furthermore, susceptible bats infected with P. destructans upregulate multiple, complementary immune responses (Field et al., 2015;Lilley et al., 2017;Moore et al., 2013;Rapin et al., 2014), so immunogenetic selection by WNS cannot be fully captured by experimental designs that target single, candidate genes. No other population-level immunogenetic analyses exist for M. lucifugus, or for any other North American species of bats threatened by WNS. Fortunately, new molecular tools allow more comprehensive investigation of immunogenetic adaptation (Harrisson, Pavlova, Telonis-Scott, & Sunnucks, 2014).
Genotype-by-sequencing (GBS) assays have emerged as a costeffective method for obtaining population-level assessments of neutral and functional genetic variation, and identifying local adaptation (Tiffin & Ross-Ibarra, 2014). GBS assays involve enriching for genomic subsets of DNA via restriction enzyme-, amplicon-, or hybridizationbased methods (Jones & Good, 2016), conducting high-throughput sequencing and identifying single nucleotide polymorphisms (SNPs).
Targeted approaches, including amplicon-and hybridization-based GBS, have been used in wildlife studies to identify SNPs in specific coding and regulatory regions of immune genes, collectively called the "immunome." Targeted GBS can identify population-level immunogenetic shifts in response to pathogens, and has been applied to a range of species, including the Tasmanian devil (Sarcophilus harrisii; Morris, Wright, Grueber, Hogg, & Belov, 2015), turkey (Meleagris gallopavo; Reed, Mendoza, & Settlage, 2016), gray wolf (Canis lupus ;Schweizer et al., 2016), thinhorn sheep (Ovis dalli; Roffler et al., 2016), and red fox (Vulpes vulpes; Donaldson et al., unpublished). GBS is an attractive option for understanding the impacts of WNS on immunogenetic diversity in bat populations, because it allows accurate characterization of diversity at duplicated loci, and cost-effective targeting of multiple, relevant genes. Regardless of the genomic coverage of highthroughput sequencing methods, population genetic analyses still rely on adequate sample sizes to detect the genetic signature of selection by pathogens or other selective pressures, reinforcing the importance of a cost-effective approach.
We developed a novel hybridization-based GBS assay to characterize the genetic diversity of the M. lucifugus immunome. Our assay includes 170 loci, including 120 immune genes and their regulatory areas, 18 Drb1-like exon 2 regions, and 32 neutral loci to allow characterization of neutral population structure, against which hypotheses of local adaptation can be tested. We applied this assay to test the hypothesis that WNS exerts significant selection on the immunome of affected bat populations. Controlling for neutral genetic population structure, we predicted immunogenetic divergence would be detectable between WNS-naïve populations and populations affected by WNS. This study provides a foundation for future investigations of rangewide immunogenetic adaptation to WNS in M. lucifugus and other affected species of bats.

| Microsatellite markers, immune genes and probe development
We developed our assay for primary application to M. lucifugus because this species' genome is publicly available (Myoluc2.0 genome assembly, Ensembl release version 81; Cunningham et al., 2015), and recent research has identified putative "WNS-response" genes for this species (Rapin et al., 2014), which informed our selection of target genes for sequence capture.
To assess functional immunogenetic variation, we assembled a list of 120 candidate genes related to immune system processes based on (i) the Human Innate & Adaptive Immune Responses RT 2 Profiler PCR Array (Qiagen); (ii) a review of innate and adaptive immunity, development, and signaling (Knight, 2013); (iii) a study of gene expression in M. lucifugus following infection with P. destructans (Rapin et al., 2014); and (iv) a gene ontology (GO) term search in the M. lucifugus Ensembl database for GO records related to fungi (including cellular response to molecule of fungal origin, defense response to fungus, and neutrophil-mediated killing of fungus). We used this candidate gene list to query the M. lucifugus Myoluc2.0 genome assembly and created a BED-formatted file containing coordinates for all exons. Additionally, we targeted potential regulatory regions by including coordinates for the 1,500-bp region upstream from the first exon for each gene.
Finally, we added exon 2 coordinates for 18 Drb1-like genes identified in Ensembl that putatively encode functional full-length proteins.
Using these primer sets, we added coordinates for these markers to the BED-formatted file. In total, the final BED-formatted file contained coordinates for 170 loci. Descriptions for protein-coding and microsatellite regions are provided in the Supporting Information (Tables S1-S2).
Custom NimbleGen SeqCap EZ probes (Roche) were produced for "primary targets" using the BED-formatted file and the M. lucifugus Myoluc2.0 genome assembly as a reference. We added 100-bp "padding" to each target to increase the efficiency of the sequence capture, and we used a "relaxed" probe design that allowed up to 20 close matches to the M. lucifugus reference genome. We compared our probes to the M. lucifugus reference genome to ensure that our assay had a low likelihood for "off-target" sequence capture: 91% of the probes matched only their target sequence, and 99% had five or fewer matches to the M. lucifugus reference genome.

| Sample collection, DNA extraction, and quantification
All work was conducted under approved animal care protocols from the University of Winnipeg and the Ontario Ministry of Natural Resources and Forestry. To test the relative impacts of geographic location and exposure to P. destructans on neutral and immunogenetic population structure in M. lucifugus, we assigned bat samples collected in eastern Canada to three post hoc groups (Table S3).
The first group included bats that were nonharmfully sampled at a hibernaculum in Manitoba, Canada (MB, n = 28), that did not contain P. destructans at the time of sampling. The second group contained bats from two hibernacula near Thunder Bay, Ontario, which were also sampled before the arrival of P. destructans (ON, n = 36).
Wing biopsies from these bats were immediately stored in RNAlater (Qiagen) following sampling. These two groups represent our "pre-WNS treatment." The third group came from populations of bats in Atlantic Canada that had been exposed to WNS for at least one winter, but were found moribund or dead in the winter of 2014 in the Atlantic provinces of Nova Scotia and Prince Edward Island (ATL, "post-WNS treatment," n = 28). These bats were submitted to the Canadian Wildlife Health Cooperative (CWHC), Atlantic Region for necropsy. Of these post-WNS bats, 15 were diagnosed as positive for WNS, 12 were diagnosed as suspect D for WNS, and 1 was negative for WNS using the approved diagnostic categories for WNS found in the Canadian Bat WNS Necropsy Protocol (CWHC, 2014), and they were assumed to not be immunotolerant nor immunocompetent to WNS. Wing tissue was collected from the left dactylopatagium major during these necropsies and stored in lysis buffer (4 M urea, 0.2 M NaCl, 0.5% n-lauroyl sarcosine, 10 mM 1,2-cyclohexanediaminetetr aacetic acid, 0.1 mM Tris-HCl pH 8.0) until analysis. We dissolved all tissue samples in lysis buffer containing 600 U/ml proteinase K at 56°C for 2 hr. We extracted DNA using either the automated 96-well MagneSil Blood Genomic Max Yield System (Promega) or the DNeasy Blood and Tissue Kit (Qiagen). We then quantified all DNA extractions using the Quant-iTPicoGreen dsDNA Assay Kit (ThermoFisher Scientific).
To investigate the possibility that our assay could also be used to investigate immunogenetic variation and adaptation in other species affected by WNS, we also isolated DNA from "post-WNS" Eptesicus fuscus (n = 2) from New Brunswick, both suspect B for WNS (CWHC, 2014) and M. septentrionalis (n = 2) from Nova Scotia and Prince Edward Island, positive and suspect D for WNS, respectively (CWHC, 2014), and included these samples in the assay.

| DNA library preparation, sequence capture, and high-throughput sequencing
We prepared DNA libraries using the KAPA HTP Lib Prep Kit (Roche) and performed the sequence capture using the NimbleGen SeqCap EZ Developer Library kit v5.1 (Roche) with the following modifications to the manufacturer's protocol. Each DNA library preparation used 150 ng total DNA. TruSeq HT Dual-Index Adapters (Integrated DNA Technologies) resuspended in Nuclease Free Duplex Buffer (Integrated DNA technologies) were used at a final concentration of 600 nM instead of the SeqCap Adapter Kits A and B (Roche) during adapter ligation. We performed 11 cycles during the LM-PCR, and initial DNA library quality was assessed by ethidium bromide-stained gel electrophoresis using a 2% E-Gel (ThermoFisher Scientific). We

| Analysis of targeted microsatellites
We used two different approaches to assign microsatellite genotypes. The first method (gatk) relied on sequence alignment to the M. lucifugus genome. We identified a single INDEL to represent each microsatellite by selecting the short tandem repeat that yielded the highest: (i) percentage of heterozygotes; (ii) gatk "quality" score; or (iii) number of alleles. For each of these three scenarios, we used the gatk to calculate the number of heterozygotes for each marker using a subset of our data that included only the 36 ON samples and the 28 MB samples. These 64 samples were previously genotyped based on traditional PCR amplification and sequencing of 11 microsatellite markers . We calculated the number of heterozygotes for each of these markers, to assess whether our sequence capture assay could be used to build on previous microsatellite-based studies. Our second genotyping method used the Galaxy platform (Afgan et al., 2016) to run STR-FM (Galaxy Version 1.0.0 ;Fungtammasan et al., 2015) and identify di-and tetra-nucleotide STRs from the raw Illumina.fastq data, without genome alignment.

| Analysis of functional loci and identification of novel, putatively neutral SNPs
We used gatk to assemble a master variant call format file (.vcf) that included SNPs with a maximum missing genotype frequency of 5% and a minimum minor allele frequency of 2%. We then used gatk to generate subdatasets of SNPs from specific categories (exon, intron, regulatory region, and Drb1-like exon 2). For the "off-target" SNPs, we used the Ensembl variant effect predictor tool to determine the bp distance between a SNP and the closest gene and generated a list of putatively neutral SNPs that were at least 100,000 bp from a gene (e.g., Kawakami et al., 2014), which we considered to be in linkage disequilibrium. We "binned" these SNPs based on the minor allele frequency, and tested for genetic structure (see below) using the SNPs with minor allele frequency values of 2% and 25%. All.vcf files were reformatted using pgdspider v2.0.9.2 (Lischer & Excoffier, 2012) for downstream analyses.
To explore variation in functional regions, we ran two lositan analyses (Antao, Lopes, Lopes, Beja-Pereira, & Luikart, 2008;Beaumont & Nichols, 1996) to identify F ST outliers that are putatively under selection. lositan parameters included 1,000,000 iterations, a 99.5% confidence interval, a false discovery rate (FDR) threshold of 0.05, and a stepwise mutation model. We enabled the "Neutral mean F ST " and the "Force mean F ST " options. The first analysis used population priors based on geography (MB, ON, and ATL) and the second considered exposure to WNS (pre-WNS, post-WNS). We extracted the subset of directional F ST outliers identified in each analysis with vcftools v0.1.14 (Danecek et al., 2011) and used them to explore immunogenetic population structure (see below).
SNPs that alter amino acids or affect splicing regions can have major effects on the function of encoded proteins. We considered F ST outliers that had these particular consequences as the most likely signals of either local adaptation to pre-occurring pathogens (in the geographic comparison) or alleles disproportionately selected against by WNS (in the pre-and post-WNS comparison). Mutations in regulatory regions can also influence gene expression and ultimately affect disease outcome (Fraser, 2013) so we also identified SNPs within regulatory regions, although the functional results of these mutations cannot be inferred from sequence capture data alone.

| Characterization of neutral and immunogenetic population structure
We used two a priori groupings to test for neutral genetic and immunogenetic population structure: (i) geographic grouping (MN, ON, and ATL), or (ii) grouping by exposure to WNS (pre-and post-WNS). We explored genetic structure based on the different SNP datasets using structure v2.3.4 (Pritchard, Stephens, & Donnel l y, 2000) and paral l el ized the runs using the strauto v0.3.1 script (Chhatre & Emerson, 2017). We ran structure with a burn-in length of 50,000 followed by 200,000 iterations for K = 1 through 4, and each run was performed 20 times. We used structure harvester web v0.6.94 (Earl & VonHoldt, 2012) to calculate the ΔK statistic (Evanno, Regnaut, & Goudet, 2005). Multiple structure runs were combined with clumpp v1.1.2 (Jakobsson & Rosenberg, 2007) using the Greedy option (10,000 repeats), and we visualized the results using distruct v1.1 (Rosenberg, 2004). We also performed principal component analysis (PCA) 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).  Table 2). To visualize the variation in depth of coverage across samples and the primary targets, we plotted the mean depth of coverage for primary targets across all samples (Figure 1) and the depth of coverage obtained from each sample for the primary targets ( Figure 2).

| Microsatellite genotyping via INDEL detection
When processing the 32 microsatellite loci included in our assay, gatk analysis identified 400 INDELs for the 32 loci, demonstrating that unique microsatellite regions contained multiple INDEL calls.
However, the relatively short read length obtained with our sequencing method failed to reliably capture entire short tandem repeat (STR) regions. Thus, microsatellite genotypes could not be recovered for all samples. As a result, heterozygous genotypes scored from our highthroughput sequencing differed from the previous results obtained using traditional PCR methods  by −25% to −42%.
The STR-FM analysis, which does not rely on aligning reads to the genome, was also unable to generate genotypes for more than two microsatellite markers using a subset of our samples (data not shown). Therefore, we did not conduct further analyses with the microsatellite data.
T A B L E 1 High-throughput sequencing and read alignment summary statistics

| Analysis of neutral genetic structure
The gatk analysis identified 16,115 "off-target" SNPs. The Ensembl variant effect predictor tool found that 1,038 of these SNPs were located >100,000 bp from a neighboring gene (Table 3; Table   S5). The putatively neutral SNPs map to 111 different clusters (>100,000 bp from the next cluster) on 88 different scaffolds of the M. lucifugus genome sequence assembly. We found no evidence for genetic structure based on these putatively neutral SNPs in the adegenet-or structure-derived plots based on geography (Figure 3a), or based on the presence of WNS in those areas (Figure 3b), regardless of the minor allele frequency cutoff used in the analysis (data not shown).

| Immunome SNP detection and analyses
We identified 17,495 SNPs within the primary target loci, located  Tables 5 and S6. Focusing on SNPs most likely to cause major functional changes, we found that 23 outlier SNPs in 19 genes in the geographic comparison resulted in an amino acid change, as did 28 SNPs in 21 genes in the WNS comparison (Table 6).
In the WNS comparison, an outlier SNP in the intron region of HLA-DPB1 resulted in likely modification of the splice donor sequence. We also detected 194 SNPS in the regulatory regions of 78 genes (Table S6), of which 11 were identified in both comparisons (Table 4).
Analyses of genetic structure in adegenet and structure did not identify geography-or WNS-associated genetic structure using the entire primary target locus, exon, intron, or regulatory region SNP datasets (data not shown). However, using the lositan-predicted outlier SNPs, we observed subtle immunogenetic structure based on geography or the presence of a brief period of co-occurrence with P. destructans (Figure 4). The analysis based on geography iden- ber of loci has many advantages in elucidating the genetic variation from known targets; however, in this instance, we chose a sequence capture approach to also pull down large segments of the immunome that included upstream regulatory regions of the genes of interest.
As such, the candidate gene GBS approach employed in this study provided several advantages over other means in obtaining immunogenetic information that is likely to be influenced by the selective pressures from disease such as that caused by P. destructans (Table 6; Figure 4).
Here, we found target enrichment led to even sequencing uniformity/coverage, which has been noted by other research groups (Powell, Amish, Haynes, Luikart, & Latch, 2016;Samorodnitsky et al., 2015;Schweizer et al., 2016). The assay we developed provided a high on target means to obtain moderate to high coverage of each target (26-463X; Table 1) that was relatively even across samples and loci (Figures 1 and 2). The assay also worked across other species (M. septentrionalis and E. fuscus) that are known to also be impacted by WNS to varying degrees (Frank et al., 2014;Frick et al., 2015;Langwig et al., 2012). As such, this assay sets the stage for cross-species analyses to further our understanding of the variable immune responses to this disease. The assay, however, yielded far too high a percentage of duplicates (54%-92%; Table 1) that compromised the level of coverage. The duplicates were likely a matter of too many PCR cycles at the adapter ligation stage during DNA library preparations. In the future, we would decrease from 11 to 6-8 cycles in the LM-PCR step. One aspect of the assay that did not meet expectations was the amplification of microsatellite loci, largely as a matter of the sequencing technology used (HiSeq 2500 rapid run using 2 × 100-bp reads on a single flow cell). The

| Immunogenetic diversity and structure
We developed a GBS sequence capture assay to cost-effectively and rapidly reveal genetic diversity in the immunome of endangered M. lucifugus threatened by mass die-offs from WNS. The assay characterized neutral population structure to control for stochastic immunogenetic differentiation among sampled areas (Table S5;  F I G U R E 3 Visualizing lack of genetic structure using 1,038 putatively neutral SNPs (>100 kpb from nearest gene) with max-missing genotype of 5% and minor allele frequency of 2% for Myotis lucifugus (n = 92). Samples were grouped based on (a) geographic location or ( (Figures 3b and 4b). If similar selection is occurring in M. lucifugus that are surviving in WNS-impacted areas, there may be potential for rapid local adaptation to WNS, raising the possibility of evolutionary rescue (Carlson et al., 2016;Maslo & Fefferman, 2015). Conversely, immunogenetic selection by WNS may disrupt previously adaptive patterns of immunogenetic variation as P. destructans continues to spread, further complicating the recovery of M. lucifugus. Our interpretation of this data is effected by the possibility that the 28 M. lucifugus we sampled from Atlantic Canada were not exposed to P. destructans during the previous year, and while 27 of these individuals died of WNS in 2014, this might have been their first exposure to an infection with P. destructans. To partially address this concern, we note 14 of 28 M. lucifugus with sample IDs "X2393-14-N" (where N varies; Table S3) came from a hibernaculum in Prince Edward Island where WNS mortality was identified in the previous winter, 1 year prior to these individuals dying of WNS and being collected for this study.
Immunogenetic diversity in M. lucifugus is extremely high.
Previous attempts to quantify variation were complicated by duplication of loci in the MHC of M. lucifugus, which exhibits up to 24 Drb1-like loci Palmer et al., 2016). Our targeted sequence capture assay controls for this gene duplication and allows genotypes to be unambiguously assigned to each individual. We detected functionally significant differentiation in several Drb1-like loci associated with both geography and previous exposure to WNS (Table 6). Exposure to WNS is also associated with a shift in genetic variation at interleukins and Toll-like receptors Field et al., 2015;Lilley et al., 2017;Rapin et al., 2014), consistent with the hypothesis that WNS exerts immunogenetic selective pressure on M. lucifugus. Our research on the interactions between M. lucifugus genetics and P. destructans continues to reinforce the need to take both interindividual and inter-regional variation of both the host and pathogen into account when interpreting genetic data. In this study, bats collected from sites <600 km apart in Manitoba and Ontario belong to a panmictic population based on neutral molecular markers, but exhibit local variation in the immunome that may result in different expression of immune genes among sites (Table 6). For example, it is possible that local immunogenetic differentiation between these sites result in different survival rates following the introduction of WNS. Variation in the regulatory regions (Table S6) could also alter the expression of integral immune genes among sites.
High immunogenetic variation in M. lucifugus has implications for the interpretation of gene expression studies as well. Bats from different sampling sites may respond differently to immune challenges due to variation in exon and regulatory regions of the immunome. Therefore, experimental gene expression studies related to P. destructans or other pathogens should explicitly control for potential geographic variation. Otherwise, observed differences in gene expression cannot be unambiguously attributed to the effects of the pathogen (or other treatments of interest).
The GBS approach used here provides a basis for real-time investigations of evolutionary rescue in populations of bats that persist following initial declines from WNS (Maslo & Fefferman, 2015). The results of our study were based on analysis of a small number of populations, and in the post-WNS population, the M. lucifugus we examined did not survive infection. Ideally, future studies will compare immunogenetic variation not only among exposed and unexposed sites, but also among time-series samples taken from bats that have survived multiple selective sweeps from one, two, or more winters in hibernacula containing P. destructans. T A B L E 4 Summary of lositan F ST outliers (FDR < 0.05) in targeted immunome features of Myotis lucifugus (n = 92). Directional outliers were used for structure and adegenet analyses (Figure 4) T A B L E 6 lositan-detected F ST outliers from SNP analyses based on a priori grouping by geographic location (Manitoba, Ontario, and Atlantic) or by WNS exposure history (pre-WNS, post-WNS). Only the SNPs that are most likely to have a functional impact by altering amino acids or affecting splice sequences are listed (see Methods and

| CONCLUSION
We developed a cost-effective high-throughput sequence capture assay to test for immunogenetic shifts in M. lucifugus populations fol- Dr. Barry Saville (Trent University) for assistance in analyzing the highthroughput sequence data; and two anonymous reviewers for critique leading to an improved version of our manuscript.

CONFLICTS OF INTEREST
The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS
MED, CMD, CKRW, and CJK conceived and designed the experiments. SM, JS, and CKRW collected bat tissue. SM performed postmortem diagnosis of WNS. MED performed the experiments and analyzed the data. CJK contributed reagents/materials/analysis tools.
MED, CMD, CKRW, CJK, and SM wrote and revised the manuscript. F I G U R E 4 Analysis of immunogenetic population structure based on lositan-detected F ST outliers with max-missing genotype of 5% and minor allele frequency of 2% for Myotis lucifugus (n = 92). Samples were grouped based on (a) geographic location or (b) previous exposure to WNS. Principal component analysis plots were produced using adegenet, and the percentage of variation for each axis and a scatter plot of eigenvalues are included for each analysis; barplot shows results of structure analysis (K = 3). MB = Manitoba (black); ON = Ontario (gray); ATL/post = Atlantic Canada/post-WNS (red); pre = pre-WNS (blue)

SUPPORTING INFORMATION
Additional Supporting Information may be found online in the supporting information tab for this article.