Prelude to a panzootic : Gene flow and immunogenetic variation in northern little brown myotis vulnerable to bat white-nose syndrome

The fungus that causes bat white-nose syndrome (WNS) recently leaped from eastern North America to the Pacific Coast. The pathogen’s spread is associated with the genetic population structure of a host (Myotis lucifugus). To understand the fine-scale neutral and immunogenetic variation among northern populations ofM. lucifugus, we sampled 1142 individuals across the species’ northern range. We used genotypes at 11 microsatellite loci to reveal the genetic structure of, and directional gene flow among, populations to predict the likely future spread of the pathogen in the northwest and to estimate effective population size (Ne). We also pyrosequenced the DRB1-like exon 2 of the class II major histocompatibility complex (MHC) in 160 individuals to explore immunogenetic selection by WNS. We identified three major neutral genetic clusters: Eastern, Montane Cordillera (and adjacent sampling areas), and Haida Gwaii, with admixture at intermediate areas and significant substructure west of the prairies. Estimates of Ne were unexpectedly low (289–16 000). Haida Gwaii may provide temporary refuge from WNS, but the western mountain ranges are not barriers to its dispersal in M. lucifugus and are unlikely to slow its spread. Our major histocompatibility complex (MHC) data suggest potential selection by WNS on the MHC, but gene duplication limited the immunogenetic analyses.


Introduction
Humans are mediating the global dispersal of species at unprecedented rates, and the resulting biological invasions can have rapid, substantial effects on native species (Ricciardi 2007;Ricciardi et al. 2013).Contact among previously isolated species is a major driver of emerging infectious diseases in humans, plants, and animals, which threaten public health, food security, and biodiversity (Fisher et al. 2012;Voyles et al. 2015).Conservation genomics and conservation biogeography represent key tools in the fight against emerging diseases.
Neutral genetic structure of endangered host species informs the delineation of meaningful management units for threat mitigation.Genetic structure can sometimes predict the spread of host-dispersed pathogens by estimating the magnitude and direction of dispersal among populations, or identifying barriers to dispersal that could limit pathogen spread (Blanchong et al. 2008;Biek and Real 2010;Wilder et al. 2015).Delimiting neutral population structure is also a prerequisite for assessing local adaptation because it allows the effects of genetic drift and selection to be separated (Ekblom et al. 2007;Spurgin and Richardson 2010;Rico et al. 2016).
Functional genetic markers further inform the management of endangered species by identifying populations with strong local adaptation.In populations under pressure from virulent pathogens, immunogenetic local adaptation may be particularly relevant.Pathogens exert powerful selective pressures on naïve host populations (Sironi et al. 2015), and genetic signatures of infection-related selective sweeps may be detectable for generations (de Groot et al. 2002;Di Gaspero et al. 2012;Deschamps et al. 2016).Detection of selective sweeps on particular genes or gene families can confirm or exclude potential mechanisms of susceptibility (in the host) or virulence (in the pathogen), informing the development of treatments (Campbell and Tishkoff 2008).
The recently emerged fungal pathogen Pseudogymnoascus destructans causes white-nose syndrome (WNS) in hibernating bats (Cryan et al. 2013).The fungus is endemic to Eurasia (Leopardi et al. 2015) and was introduced to North America, where it was first observed in New York State in 2006 (Blehert et al. 2009).Since 2006, P. destructans has spread rapidly, causing local extinctions and widespread population declines in some Nearctic species of bat, including the little brown myotis (Myotis lucifugus LeConte, 1831), northern myotis (Myotis septentrionalis), and Indiana myotis (Myotis sodalis) (Langwig et al. 2012;Frick et al. 2015).Bats are exposed when they enter hibernacula containing P. destructans.Fungal lesions may develop on the wings during hibernation and disrupt homeostasis, triggering a cascade of behavioural and physiological responses that result in high mortality rates in susceptible species of bats (reviewed in Willis 2015).
The spread of WNS is associated in part with genetic structure among populations of M. lucifugus (Miller-Butterworth et al. 2014;Davy et al. 2015;Wilder et al. 2015;Vonhof et al. 2016).Genetic, ecological, and environmental data have been used to model the direction and speed with which bats may facilitate dispersal of P. destructans towards the Pacific Coast of North America (Maher et al. 2012;O'Regan et al. 2015).However, models cannot fully predict pathogen dispersal, as illustrated by the recent, unexpected emergence of WNS in western Washington State in March 2016 (Lorch et al. 2016).Previous, large-scale population genetics studies (Vonhof et al. 2015;Wilder et al. 2015) suggest that genetic structure in M. lucifugus may be higher in the western portion of its range.If these patterns can be clarified with more targeted sampling, they may predict the path along which P. destructans is most likely to spread from its new "western front".
Increasing structure in the west (Lausen et al. 2008;Vonhof et al. 2015;Wilder et al. 2015) corresponds approximately with the Rocky Mountains-Great Plains transition.Topographic features such as the Front Range of the Rocky Mountains were proposed as potential drivers of genetic structure in western North America (Vonhof et al. 2015;Wilder et al. 2015).However, the direction of gene flow in the west and the role of potential barriers to dispersal such as the Rocky Mountains or Hecate Strait have not been explicitly tested.
Tolerance of WNS in Eurasian bats is attributed to coevolution with P. destructans (Zukal et al. 2016), and a similar evolutionary process may currently be underway in North American bats (Frick et al. 2017).In four years, following the initial, precipitous declines from WNS, estimated annual adult survivorship of M. lucifugus in some populations increased from 0.68 to 0.75 (for males), and from 0.65 to 0.70 (for females; Maslo et al. 2015).Local immunogenetic adaptation to P. destructans provides one potential explanation for this trend, and the ubiquity of M. lucifugus and its transcontinental range allow extensive sampling across a large geographic area.Thus, M. lucifugus provides an excellent model for investigating local immunogenetic adaptation in bats, and ultimately for testing whether WNS exerts selection on immune genes.
The major histocompatibility complex (MHC) is a family of immune genes involved in initiation of the adaptive immune response in vertebrates (de Groot et al. 2002;Spurgin and Richardson 2010).These genes are affected by pathogen-mediated selection in a wide range of vertebrate species (Wegner et al. 2003;Piertney and Oliver 2006;Savage and Zamudio 2011), and spatial patterns of variation in the MHC can serve as an immunogenetic proxy for the potential of a declining population to adapt to shifting selection by pathogens (Hawley and Fleischer 2012;Kyle et al. 2014).The MHC has been studied in several bats, including a suite of Neotropical phyllostomid species (Schad et al. 2012b;Real-Monroy et al. 2014;Salmier et al. 2016), the sac-winged bat (Saccopteryx bilineata; Mayer and Brunner 2007), and two North American Myotis species (M.velifer, M. velesi; Richman et al. 2010).MHC variants are correlated with ectoparasite load in the lesser bulldog bat (Noctilio albiventris; Schad et al. 2012a).Baseline MHC variation in North American bat populations and the role of MHC in the survival of bats with WNS are unknown.
Following the recent western introduction of P. destructans (Lorch et al. 2016), management biologists need detailed data on directional dispersal in M. lucifugus in northwestern North America to prioritize areas for surveillance and mitigation efforts.The range-wide approach of previous studies (Vonhof et al. 2015;Wilder et al. 2015) provides a valuable perspective on the contemporary genetic structure of this endangered species, but small sample sizes from northwestern locations (Manitoba to British Colombia, and northward to Alaska) limit their resolution.Here, we investigated fine-scale contemporary neutral genetic population structure in 1142 M. lucifugus sampled across the species' northern range.We tested for genetic population structure and used Bayesian estimates of gene flow to explore the most likely direction of spread from the new source in Washington State.We also estimated the effective population size for each identified population.Finally, we used 454 pyrosequencing of the DRB1-like exon 2 of the MHC to explore immunogenetic adaptation of M. lucifugus to P. destructans.We hypothesized that populations of M. lucifugus sampled prior to the arrival of WNS were locally adapted to existing selective pressures on the MHC, predicting discordant patterns of structure between neutral and immunogenetic markers.Finally, we tested the hypothesis that WNS exerts direct selection on the MHC by comparing DRB1 diversity in eastern Canada before and after the arrival of WNS.

Sample collection and DNA extraction
We isolated genomic DNA from non-harmful wing biopsies (Worthington-Wilmer and Barratt 1996) collected from free-ranging M. lucifugus following approved animal care authorizations from the were released unharmed at their capture location.Samples were stored in RNAlater ® (Applied   Biosystems, Carlsbad, California), lysis buffer, or 95% ethanol prior to analysis.To supplement these, we included tissue samples from M. lucifugus carcasses stored in the collection of the New Brunswick Museum, which were collected in Ontario between 1991 and2006, andin Prince Edward Island (PEI) in 1989 (Brown et al. 2007).We also sampled bat carcasses collected under wind turbines in southern Ontario by environmental consultants and submitted to the Canadian Wildlife Health Cooperative (CWHC).Mixing of M. lucifugus samples collected from hibernacula, maternity roosts, and foraging sites does not skew analyses of genetic population structure across the landscape (Wilder et al. 2015), and we, therefore, considered our samples as a single data set.Tissues from bat carcasses stored at −20 or −80 °C were sampled with flame-sterilized surgical scissors.We extracted DNA using the Qiagen DNeasy Blood and Tissue Kit protocol (Qiagen, Valencia, California).
Our samples represent 132 distinct locations, including maternity colonies, hibernacula, and foraging sites, spanning the Canadian range of M. lucifugus, and portions of Montana (MT) and Alaska (AK), USA (Fig. S1).

Microsatellite genotyping
We amplified eleven microsatellite loci for 1142 individuals using fluorescently labelled primers in four multiplex reactions (Burns et al. 2012;Johnson et al. 2014; Table S1).Polymerase chain reaction (PCR) products were analyzed on an ABI 3730 Genetic Analyzer (Applied Biosystems) using ROX-500 size standard (Applied Biosystems).We independently extracted and re-genotyped approximately 4% of our samples to assess potential genotyping error (Pompanon et al. 2005).Re-scoring was blinded (i.e., samples were not identifiable as duplicates at the time of scoring).We scored genotypes manually using GeneMarker v.2.4.2, and checked for long-allele dropout, null alleles, and genotyping stutter using MICROCHECKER (van Oosterhout et al. 2004).We tested for deviations from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) using GENEPOP (Rousset 2008).

Microsatellite analysis
Detailed methods for microsatellite analysis, including settings used in each program, are presented in Table S2.In brief, we conducted analyses both with and without prior information on the geographic origin of the samples.Some analyses of genetic diversity require individuals to be assigned to prior groups (rarified allelic richness and expected heterozygosity), and some can provide greater resolution when given prior information on likely group membership (i.e., STRUCTURE, Pritchard et al. 2000).For these analyses, we clustered samples collected at nearby sampling locations into sampling areas, aiming to (1) maximize the number of distinct sampling areas, and (2) ensure that each defined group of samples contained >25 individuals where possible, to avoid biasing frequency-based analyses with inadequate sample sizes (Hale et al. 2012).In total, we defined 23 sampling areas from 13 ecozones for analysis based on geographic clustering and on physiographic variations across our study area (Table 1; Fig. S1).We quantified genetic variation within and among sampling areas, characterized genetically differentiated clusters, and applied spatially explicit analyses of genetic structure (summarized in Table S2).
Clustering analyses in STRUCTURE (see Table S2 for methods) were run both with and without the LOCPRIOR function.STRUCTURE analyses with the LOCPRIOR function and discriminant analyses of principal components (DAPC; Jombart et al. 2010) were run twice: first using the sampling areas (Table 1) as priors, and then using the ecozones where samples were collected.Running the clustering analyses allowed us to ensure that our results were not biased by the priors used  (Puechmaille 2016).Each cluster recovered by STRUCTURE was re-analyzed independently with STRUCTURE and DAPC to test for further substructure.
To directly test the impact of the Rocky Mountains on population structure in M. lucifugus (Wilder et al. 2015), we also repeated analyses on a subset of samples from the Montane Cordillera, including samples collected directly west and east of the continental divide, in the Okanagan Valley, and on Vancouver Island, which is in the Pacific Maritime ecozone but has montane features (sampling areas 8-10, 12-14; Table 1).
We were particularly interested in quantifying directional gene flow across the landscape near the recent, western record of P. destructans, and we used BAYESASS v.3.0 to estimate proportional migration among genetic clusters identified by STRUCTURE focusing on samples collected west of Alberta (sampling areas 1-10; 13-15; Table S2).BAYESASS estimates gene flow over the last "several generations" (Wilson and Rannala 2003).This is typically interpreted as 1-5 generations (Chiucchi and Gibbs 2010;van der Meer et al. 2013), or 3-50 years for M. lucifugus (COSEWIC 2014).We also used the LD method implemented in NeEstimator v2.01 (Do et al. 2014) to estimate the contemporary effective population size (N e ) for each genetic cluster identified by STRUCTURE.Estimating N e based on LD is appropriate when migration among the tested groups is low (<5%-10%, Waples and England 2011) and is more accurate than other single sample estimators in this scenario (Gilbert and Whitlock 2015), excluding sibship assignment (Jones and Wang 2010).The sibship assignment method cannot be implemented for most studies of wild bats, as this approach requires samples from multiple cohorts, whereas bats can only be aged accurately as juvenile (<1 year) or adult (2-30+ years).In contrast, NeEstimator V2 (Do et al. 2014) does not assume that all samples are from a single cohort, and has been previously used to estimate N e in bats (Plecotus auritus; Razgour et al. 2014).We, therefore, chose this method as the most appropriate for our data set, recognizing that methods for estimating N e each have their own strengths and weaknesses (Wang et al. 2016;Waples 2016).
Rapid, extreme bottlenecks (such as those caused by panzootic events) are associated with the stochastic loss of alleles (Luikart et al. 1998) that do not necessarily indicate selection.Therefore, we also used a discriminant analysis of principal components to compare neutral genetic diversity and population structure in bats from eastern Canada sampled before the arrival of WNS, and bats from the same areas that died soon after the arrival of WNS.This test allowed us to determine whether any observed shift in MHC allelic diversity was likely stochastic (i.e., if a stochastic shift in neutral allele frequencies occurs, the same process could also impact MHC diversity), or whether shifts in MHC allelic diversity might indicate selective pressure from WNS (see below).

DRB1-like exon 2 amplicon library preparation and 454 pyrosequencing
We amplified a 213 bp fragment of MHC class II DRB1-like exon 2 (hereafter DRB1) in a subset of 160 individuals from our data set (Table 2) using a two-step PCR universal tail design.The final amplicons were gel extracted, quantified, diluted, and pooled to construct amplicon libraries of 10 5 molecules/μL.Emulsion PCR (using 0.25 molecules per bead) and 454 sequencing were conducted according to the manufacturer's recommended protocol using a PicoTitre plate on a GS Junior system (Roche Diagnostics, Laval, Quebec, Canada; Supplementary Material 1).We sequenced approximately 5% of our samples on multiple runs to control for genotyping error (Pompanon et al. 2005).
Our initial intent was to sequence a much larger set of samples to characterize MHC diversity across the study area.However, analyses of the first sequenced samples detected signatures of gene duplication (see below) that confounded analyses of genetic structure and diversity.We, therefore, restricted the MHC analysis to considering allelic diversity in the first 160 samples sequenced.

DRB1 sequence analysis
We discovered extremely high polymorphism at our target DRB1 locus, and we developed the following approach for genotyping these sequences.We used CLC Genomics Workbench v7.5.1 (Qiagen, Toronto, Ontario, Canada) to trim adapter sequences from the sequence data and aligned the sequences to the Myoluc2.0Microbat (M.lucifugus) genome assembly (release version 79) from Ensembl (Cunningham et al. 2015) using the Burrows-Wheeler Aligner software package v0.7.12 (Li 2013).We isolated sequences that uniquely mapped to a putative DRB1-like gene (ENSMLUG00000007755), and specifically to Myluc2.0 scaffold AAPE02060444.We used CLC Genomics Workbench to detect SNPs in the resulting sequences and remove pyrosequencing artifacts (detailed in Supplementary Material 1).The mapped variants and aligned sequences were compiled as "track lists" in CLC Genomics Workbench, which we visually inspected to identify homozygous or heterozygous allele sequences in each sample for DRB1.Each sequence was independently assessed twice.We also converted the nucleic acid sequences to amino acid sequences, while maintaining the correct reading frame for the resulting DRB1-like protein (see Supplementary Material 1).
We estimated the mean number of pairwise differences (k) for each population in ARLEQUIN v3.11 (Excoffier et al. 2005).We also calculated the relative frequency of MHC alleles by counting the number of individuals carrying a particular allele divided by the total number of alleles in each population (Ekblom et al. 2007).

Detection of local adaptation at MHC prior to WNS
Due to the duplication of the MHC locus we targeted, we could not reliably assign alleles to specific loci.Therefore, we did not calculate frequency-based diversity measures or use clustering analyses  that rely on the assumptions of HWE.However, we did quantify allelic richness and frequency of the detected MHC alleles and rare alleles (<1% frequency) to qualitatively compare neutral and MHC diversity.We also compared MHC allelic diversity in bats sampled from a WNS-naïve, panmictic population (n = 42) with diversity in bats that had survived several years before succumbing to WNS in the same population (n = 22), to explore the hypothesis that WNS exerts immunogenetic selection for or against particular allelic variants.

Neutral genetic variation and population structure
Microsatellite genotypes showed no evidence of long-allele dropout.Potential signatures of stutter or null alleles were inconsistent among areas and loci (Table S3), as were deviations from HWE and LD.Genotyping error for replicated samples was 0%.We, therefore, used the entire data set for subsequent analyses.
Microsatellite heterozygosity of sampled areas ranged from 0.62 (Haida Gwaii) to 0.83 (Northwestern Ontario; Table 1).Pairwise D EST ranged from 0.363 (Haida Gwaii and Kirkland Lake, Ontario) to 0.0001 between adjacent areas in Montana and Alberta (Table S3; F ST values also shown in Table S4).Mantel tests detected no correlation between genetic distance and Euclidean distance (permutation test with 9999 replicates; R 2 = 0.001; p = 0.5227).However, a significant correlation occurred when testing only the samples from Haida Gwaii, mainland BC, Alaska, the Yukon, and Montana (R 2 = 0.0141; p = 0.041).AMOVA indicated low but significant neutral genetic variation (global F ST = 0.012; p < 0.001), with the greatest genetic variation occurring within populations (79%), followed by variation within individuals (11%) and among populations (9%).STRUCTURE identified three major genetic clusters with substantial admixture at intermediate locations (Figs.1A, 1B).These clusters were recovered by STRUCTURE analyses with and without the LOCPRIOR function (Fig. S2), by sPCA (Fig. 1C), and by DAPC based on sampling areas or ecozones (Figs.1D, 1E).Bats from PEI to Manitoba (n = 442) formed a single cluster.Independent STRUCTURE and DAPC analyses of this cluster revealed no further substructure (Fig. S2).Lack of structure in eastern samples, including between samples collected in 1989 and samples collected in 2015, confirmed that neither drift nor selection had significantly changed the allele frequencies of this panmictic population, indicating that year of collection was not a concern when interpreting other results from clustering analyses.Bats captured west of Manitoba were clustered into a strongly assigned group of bats on Haida Gwaii, and a cluster of genetically similar bats sampled along the Montane Cordillera, which also included samples from Vancouver Island and the Okanagan Valley.Subsampling of the westernmost areas detected further differentiation between Haida Gwaii and continental samples (Figs. 2, S5), and suggested that the western Yukon bats may form a further, distinct cluster, with bats from Alaska showing Yukon and Haida Gwaii coancestry (Fig. 2).This pattern was further supported by DAPC, which resolves the same three clusters as the STRUCTURE analyses (Fig. 2D), and sPCA, in which the first axis divided Haida Gwaii samples from the other western samples (Fig. 2E).
Subsampling of bats from the Montane Cordillera and prairies indicated that genetic partitioning was associated with ecozone boundaries, not the position of the mountain ranges themselves.This result was consistent regardless of whether sampling area or ecozones were used as priors, and when no LOCPRIOR was indicated (Fig. S2).The exception occurred when bats from the Montana prairies (sampling area 15) were assigned strongly to the Montane Cordillera cluster in the analyses of western sites (Figs.2A, 2B), likely because of the smaller sample size for this area (n = 16).When analyzed alongside the other prairie samples (sampling area 16), the prairie areas either grouped together or  showed the same level of admixture between the two adjacent clusters (Figs.1A, S5).Further analyses of the prairie samples to investigate potential structure did not provide further resolution (e.g., Figs.S2, S5).
In our analyses of samples collected in and directly adjacent to the Montane Cordillera (sampling areas 8-10; 12-14; n = 228), sPCA and DAPC both indicated that the mountains do not act as a barrier to gene flow, and supported the STRUCTURE analyses suggesting that bats in and directly adjacent to the Montane Cordillera represent a biologically meaningful genetic cluster (Fig. 3).
Estimates of contemporary N e based on allele frequencies >0.02 ranged from 606.1 (Haida Gwaii) to 16 062.4(eastern Canada).Despite large confidence intervals and substantial variation among estimates based on different cut-offs for the lowest allele frequency considered (Table 3), none of the mean N e estimates for western populations were >1600.
We used a final STRUCTURE analysis of bats collected north and west of Alberta (Haida Gwaii, Alaska, British Colombia, the Yukon, and Montana; sampling areas 1-10; 13-15; Figs.2A, 2B) to group sampling areas with population q values >0.7 to a particular genetic cluster.Samples from Alaska were strongly admixed, and we, therefore, considered Alaska as a separate group for the BAYESASS estimates of directional gene flow (Fig. 2C).We detected concordant genetic structure with DACP and sPCA analyses (Figs.2D, 2E).BAYEASS detected substantial migration from the western Yukon cluster to Alaska, with migrants accounting for approximately 27% of Alaskan ancestry per generation.Migrants from the western Yukon cluster accounted for approximately 13% per generation of the cluster containing samples from the Montane Cordillera and adjacent areas (Fig. 2C).Other estimated migration rates (range: 4.3%-1.4%)had confidence intervals that overlapped zero, indicating that these migration rates were not significantly different from zero (Table S3).

Immunogenetic variation and population structure
Our DRB1 analysis mapped an average of 1118 reads/sample to the bat genome.Coverage averaged 277× for reads mapped to ENSMLUG00000007755 (putative DRB1-like exon 2).All technical repeats of MHC genotypes were identical, providing no evidence for genotyping error.During analysis, it became clear that the target locus is highly polymorphic in M. lucifugus.Allelic richness was strongly correlated with sample size (r 2 = 0.827; Table 2; Fig. S3).We, therefore, report MHC allelic diversity, but we do not consider this value a statistically appropriate estimate of heterozygosity because this value can be confounded by duplicated loci, even with our intensive approach of mapping reads directly to the genome.

Preliminary investigation of the impacts of WNS on genetic variation
The DAPC of the microsatellite data confirmed that bats from the eastern cluster sampled before or shortly after the arrival of WNS were drawn from a genetically homogeneous population, with similar levels of heterozygosity (Fig. S4).However, the distribution of DRB1 alleles shifted in samples collected from bats that had survived several winters with the pathogen (Fig. 4).The most common allele prior to WNS (H01) remained the most common allele detected in multi-year survivors, but we observed a higher frequency of rare alleles (i.e., n = 21 alleles with frequency <1% in the entire data set) in these samples, as well as a general qualitative increase in MHC allelic diversity in multi-year survivors.Further work would be required to test the statistical significance of this trend.We did not observe the emergence of any new "common" alleles.bats on Haida Gwaii appear to be effectively isolated.Our results predict that Haida Gwaii may act as a temporary refugium from the disease, whereas M. lucifugus may carry P. destructans more quickly through southern BC, across the Rockies, and into Montana and Alberta.Southward gene flow from the Yukon opposes the direction of the oncoming pathogen, providing some hope that bats in the  Note: Estimates using the lowest allele frequencies are italicized to indicate the higher probability that sampling bias will affect the estimates.Yukon and Alaska may experience the same delay in P. destructans arrival that occurred in central Canada (Davy et al. 2015;Wilder et al. 2015).We detected unexpected fine-scale population structure in the Yukon, with differentiation between the eastern and western Yukon sampling areas (Figs.2A, 2B).Preliminary immunogenetic analysis suggests that baseline patterns of MHC may be altered by the impacts of WNS.This could indicate that the MHC is under balancing selection that is disrupted by the arrival of WNS.However, our immunogenetic analyses were confounded by high gene duplication in this species (see also Palmer et al. (2016), who detected 24 DRB1-like alleles in only 15 individuals), and new sequencing methods will be required to disentangle the immunogenetic impacts of WNS on bats.We also relied on samples from bats that died in the first waves of mortality following WNS.Although our results are suggestive of selection by WNS on the MHC, we were not able to confirm or refute this hypothesis.
Estimates of effective population size are surprisingly low given M. lucifugus' ubiquitous distribution and high genetic diversity, which are often associated with extremely high N e .Our estimates were 10-1000× lower than estimates for three migratory bat species (Lasionycteris noctivagans, Lasiurus cinereus, and Lasiurus borealis), based on samples from a single wind farm (Sovic et al. 2016).Estimates of N e based on single samples can be affected by many factors, and the values presented here should not be misinterpreted as an estimate of the absolute number of breeding individuals.Rather, the estimates indicate relative effective population sizes among genetic clusters.The low values estimated here suggest that despite its historical ubiquity, M. lucifugus may be unexpectedly vulnerable to population declines from additive mortality, with more reproductive variation among individuals than previously thought.The particularly low N e of western populations may leave them more vulnerable to declines once the pathogen arrives, although the genetic structure of M. lucifugus in the west may slow the spread of WNS.
Our results expand on previous studies of M. lucifugus (Burns et al. 2014;Davy et al. 2015;Vonhof et al. 2015;Wilder et al. 2015) by clarifying the locations of genetic transitions among western bat populations, where previous studies already detected increased genetic variation, and revealing further population structure in the west, particularly on Haida Gwaii.Haida Gwaii bats apparently exhibit low dispersal across the Hecate Strait, and dispersal towards the island is negligible.Haida Gwaii acted as a glacial refugium during the Quaternary glaciation, as did parts of Alaska and the Yukon, and the genetic endemism of bats on Haida Gwaii mirrors that found in other taxa, including a sea star (Patiria miniata), the black bear (Ursus americanus), and several land birds (Byun et al. 1997;Keever et al. 2009;Pruett et al. 2013).Although M. lucifugus undertake migrations >500 km (Fenton 1969;Norquay et al. 2013) and disperse to and from Vancouver Island, Newfoundland, and PEI (Burns et al. 2014;this study), our results suggest that the longer flight distance to Haida Gwaii from adjacent coastal islands (≥48 km) and the prevailing westerly winds may deter bats from crossing.Nevertheless, accidental translocation of bats on boats could move individuals across the Hecate Strait (e.g., Constantine 2003;Wright and Moran 2010).We strongly recommend that stringent preventative measures and public outreach are used to avoid human dispersal of P. destructans to the archipelago and delay the pathogen's arrival for as long as possible.
Previous studies suggested that mountain ranges may limit population connectivity in M. lucifugus (Vonhof et al. 2015;Wilder et al. 2015), but our sampling design revealed that the mountain ranges are not barriers to dispersal in this species.In our study area, the closely related M. californicus and M. yumanensis occur only on the western side of the Continental Divide (C.L. Lausen, unpublished data), indicating that the mountains act as a barrier.In contrast, genetic clusters of M. lucifugus align with ecozones, suggesting possible "montane", "prairie", and "coastal" ecotypes.Morphological differences in M. lucifugus including pinnae length and forearm length support this hypothesis, also corresponding to ecozone boundaries (Lausen et al. 2008 data).Future research can test potential behavioural or morphological among putative ecotypes that could explain their continued differentiation such as differences in movement, communication, hibernation, or foraging behaviour (Veselka et al. 2013;Jung et al. 2014;Pond et al. 2016).
In the western Yukon, M. lucifugus take advantage of insect-rich summer feeding grounds, but apparently migrate over the Coast Mountains to hibernate in talus slopes and tree-root wads of coastal Alaska, where winters are milder (Jung et al. 2014;Blejwas et al. 2015).Movement of Yukon bats into Alaska is also supported by our data, and by anecdotal reports of bats frozen into glaciers of the Coast Mountains (Slough and Jung 2008).In the eastern Yukon, karst formations provide potential hibernacula and remove the need for energetically costly winter migrations.Thus, different availability of overwintering sites may drive the observed differentiation between bats in western and eastern Yukon.As above, morphological data (forearm length) also differentiate bats from the western and eastern Yukon (Lausen et al. 2008;Talerico 2008;C.L. Lausen et al., unpublished data).
The Canadian prairies appear to represent a zone of admixture between M. lucifugus from the Eastern and Montane Cordillera clusters; a pattern that may be maintained by annual migrations among maternity, swarming, and overwintering sites, or that may represent the presence of a further, distinct genetic cluster that was not detected by our analyses.Maternity colonies of M. lucifugus occur in anthropogenic structures on the Great Plains during the summer.Our results suggest that many of the M. lucifugus using prairie maternity colonies must migrate to the edges to swarm (mate) and hibernate, restricting genetic exchange to the periphery of the prairies and preventing differentiation of the prairie population.No M. lucifugus hibernacula are known in Saskatchewan or from the Alberta prairies, although hibernacula are known to the north in the Alberta boreal (Reimer et al. 2014;C.L. Lausen et al., unpublished data) and to the west in the Alberta Rocky Mountains (Schowalter 1980).Collecting genetic samples from these sites in future may provide further resolution on the ancestry of prairie bats and their relationship with surrounding genetic clusters.
Despite genetic population structure and high MHC polymorphism, M. lucifugus exhibit low spatial immunogenetic structure across the northern extent of their range.The lack of MHC differentiation observed despite low levels of gene flow may reflect balancing selection prior to the introduction of WNS, from immunogenetic selective pressures on M. lucifugus that may have been relatively constant among locations.However, gene duplication at the DRB1 locus limits the use of this data set.Specific MHC alleles have been linked to disease susceptibility in other host-pathogen systems, including amphibian resistance to the pathogenic fungus Batrachochytridium dendrobatidis and the susceptibility of giant pandas (Ailuropoda melanoleuca) to intestinal parasitic infections (Savage and Zamudio 2011;Bataille et al. 2015;Zhang et al. 2015).Sequencing approaches that can distinguish between alleles from duplicated genes are required to further test the hypothesis that WNS exerts immunogenetic selection.
The idea that bats in western North America may harbour greater potential for the evolution of disease resistance based on standing genetic variation (Wilder et al. 2015) is extremely attractive, especially given our current inability to prevent or cure WNS.The link between standing genetic variation within a population and adaptive potential with regard to shifting selective pressures is a central tenet of conservation genomics, but unfortunately, it does not always apply in practice (Rollins et al. 2013;Mittell et al. 2015).Bat WNS provides an excellent illustration of high standing genetic variation providing little protection in the face of emerging infectious diseases.Myotis lucifugus was one of the most common, widespread and abundant mammal species in North America prior to the arrival of WNS.Heterozygosity of neutral markers and allelic diversity of at least one functional marker are high across the species' range.Nevertheless, an emerging pathogen has brought this species to the brink of regional extinction in less than a decade.The new selective pressure to which northeastern populations of M. lucifugus were exposed (i.e., WNS) is simply not one to which those populations could quickly adapt; at least not without high individual rates.Although the conservation of genetic diversity is critical to long-term species conservation, it is not a guarantee of adaptive potential.Our preliminary results also indicate that although M. lucifugus in western North America exhibit neutral genetic differences compared with the eastern population, this is not reflected in a functional gene expected to experience direct selection from WNS.
Selective sweeps from emerging diseases such as WNS could erode immunogenetic adaptations to other preexisting or future pathogens.The lack of regional MHC differentiation in M. lucifugus across its northern range suggests that disease-related selective pressures on this species may be similar among genetically differentiated clusters, but testing this hypothesis requires more robust genetic tools that can differentiate among duplicated loci, and that can simultaneously target multiple immune genes.Bats in our sampling area coexist with a variety of pathogens including several strains of bat rabies (Middleton et al. 2015), corona and polyoma viruses (Misra et al. 2009), endemic fungal pathogens such as Trichophyton redelli (Lorch et al. 2015), internal parasites such as Trypanosoma myoti (Bower and Woo 1981), and ectoparasites such as fleas and mites (Webber et al. 2015).Most of these are relatively widespread, affecting bats across North America.If WNS causes a significant selective sweep at MHC in M. lucifugus across the continent, bats that survive WNS may experience a trade-off in susceptibility to other pathogens; a carryover effect of the disease that extends its influence past the stage of acute infection and recovery.Examples of other potential carryover effects from WNS include delayed parturition (Francl et al. 2012) and increased chronic stress following recovery (Davy et al. 2017).
Emerging infectious diseases provide excellent opportunities to document selective sweeps in real time, and to test hypotheses about the impacts of disease on genetic and phenotypic diversity.This requires a backdrop of robust data on neutral genetic diversity against which immunogenetic diversity can be contrasted (Ekblom et al. 2007;Rico et al. 2016).However, studies of immune response to infectious diseases rarely consider more than one pathogen at a time and often focus on particular immune genes (as we have done here).In reality, the immune system responds to a myriad of pathogens and other stressors simultaneously, and no free-ranging population is likely to coevolve in isolation with a single pathogen.The impacts of emerging infectious diseases on immunogenetic variation should, therefore, be considered in the context in which they emerge-as novel selective pressures on host populations that already have long evolutionary histories with a complex set of other pathogens, parasites, and commensal microflora.

Fig. 1 .
Fig. 1.Myotis lucifugus (LeConte, 1831) fall into three genetic clusters across the northern limits of their range with substantial admixture at intermediate sites.(A) STRUCTURE analysis of 1142 M. lucifugus genotyped at 11 polymorphic microsatellite loci shows clustering of samples into three genetically differentiated clusters.(B) Ancestry estimates for each sampled area.Interpolation of the first lagged score of a spatial analysis of principal components supports the STRUCTURE assignment of bats into a Haida Gwaii and Eastern cluster, with admixture around each, and shows that bats from the Montane Cordillera and adjacent sites are not strongly assigned to either of these clusters.(D) Discriminant analysis of principal components using sampling areas as prior groups, or (E) ecozones as prior groups generally support the clusters identified by STRUCTURE.Ecozones: PM, Pacific Maritime; BC, Boreal Cordillera; SAP, Semi-arid Plateau; MC, Montane Cordillera; PR, Prairies; TP, Taiga Plains; BP, Boreal Plains; BS, Boreal Shield; GL, Great Lakes/Mixedwood Plains; AM, Atlantic Maritime.Numbers for sample areas correspond to Table 1; areas marked with "*" in (A) are already impacted by white-nose syndrome (WNS).Red stars on the maps indicate the most westerly current record of WNS in North America (March 2016).Map created using Natural Earth (naturalearthdata.com).
Davy et al.FACETS | 2017 | 2: 690-714 | DOI: 10.1139/facets-2017-0022 700 facetsjournal.comFACETS Downloaded from www.facetsjournal.comby Canadian Science Publishing on 05/01/18 Discussion Neutral genetic structure in northern populations of M. lucifugus generally increases from east to west, and the spatial distribution of neutral genetic diversity across the northwestern range of M. lucifugus is intriguing: neither the Coastal nor Rocky Mountains act as barriers to gene flow, yet

Fig. 3 .
Fig. 3. Genetic differentiation in Myotis lucifugus sampled within and adjacent to the Montane Cordillera (n = 228) does not support the hypothesis that the Rocky or Coastal Mountains are a barrier to dispersal in this species.(A) Interpolated first lagged score of a spatial principal components analysis.(B) Discriminant analysis of principal components.The symbols W and E indicate samples collected west and east of the Continental Divide; numbers in the legend correspond with sampling areas defined in Table1.Map created using Natural Earth (naturalearthdata.com).

Fig. 4 .
Fig. 4. Frequency of rare MHC alleles (<1%) in bats from a panmictic population in eastern Canada (Ontario to Atlantic Canada) shifts following exposure to white-nose syndrome (WNS).Samples from Prince Edward Island were taken before the arrival of Pseudogymnoascus destructans, and a year following the initial observation of WNS (2013), during which time the number of unique alleles detected per individual rose from 0.08 to 0.41.(A) Shows the frequency of occurrence for each allele among the three groups of samples and (B) shows the frequency of each allele within the three sets of samples.
Note: Sampling sites are the number of distinct locations within each area at which bats were collected.N (f, m), total sample size (number of females, number of males); Na, mean number of alleles/locus; Ar, rarified allelic richness; H O , observed heterozygosity; H E , expected heterozygosity;F, fixation index [(H E − H O )/H E = 1 − (H O /H E )].aSamples from areas 13 and 14 were collected from bats immediately west and east of the Continental Divide, in Glacier National Park, Montana.

Table 2 .
Genetic diversity indices for MHC class II DRB1-like exon 2 in little brown myotis (Myotis lucifugus) sampled across Canada before the arrival of Pseudogymnoascus destructans, sampled in Prince Edward Island one winter after the arrival of P. destructans, and sampled following substantial local population declines due to white-nose syndrome (WNS).
Note: Numbers in parentheses indicate the sampling area from which samples were drawn, following Table1.a Samples collected in PEI in 2014 from bats that had survived at least one winter season with WNS but succumbed in winter 2014.These samples were collected from sampling area 23 but are different individuals than those used in the microsatellite analyses.

Table 3 .
Effective population size (mean N e , with 95% confidence interval range in parentheses) estimated for genetic clusters identified by STRUCTURE (Figs.1 and 2) based on genotypes at 11 microsatellite loci.