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 of M. 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”.
Myotis lucifugus form a large, panmictic cluster in northeastern North America (Burns et al. 2014; Miller-Butterworth et al. 2014), but begin to show genetic structure along the transition from the eastern boreal forests to the prairies of central Canada (Davy et al. 2015; Wilder et al. 2015). 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.
Materials and methods
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 Ontario Ministry of Natural Resources and Forestry and the University of Manitoba. Sampled bats 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 and 2006, and in 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).
Analysis of neutral genetic structure
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).
Table 1.
Table 1. Neutral genetic diversity in Myotis lucifugus (LeConte, 1831) sampled across 23 sampling areas and genotyped at 11 microsatellite loci.
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; HO, observed heterozygosity; HE, expected heterozygosity; F, fixation index [(HE − HO)/HE = 1 − (HO/HE)].
a
Samples from areas 13 and 14 were collected from bats immediately west and east of the Continental Divide, in Glacier National Park, Montana.
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 (Ne) for each genetic cluster identified by STRUCTURE. Estimating Ne 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 Ne 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 Ne 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).
Analysis of immunogenetic variation
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 105 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.
Table 2.
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 Table 1.
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.
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.0 Microbat (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.
Results
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 DEST ranged from 0.363 (Haida Gwaii and Kirkland Lake, Ontario) to 0.0001 between adjacent areas in Montana and Alberta (Table S3; FST values also shown in Table S4). Mantel tests detected no correlation between genetic distance and Euclidean distance (permutation test with 9999 replicates; R2 = 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 (R2 = 0.0141; p = 0.041). AMOVA indicated low but significant neutral genetic variation (global FST = 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).
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).
Fig. 2.
Fig. 2. Myotis lucifugus in northwestern North America sampled prior to the introduction of white-nose syndrome (WNS) exhibit complex fine-scale population structure based on microsatellite genotypes (11 loci; n = 478, including sampling areas 1–10; 13–15). (A) STRUCTURE barplots showing individual q-values for K = 2 and K = 4. (B) STRUCTURE population q-values for each sample area (corresponding to Table 1, with sampling area 9 divided into samples from Vancouver Island (left) and samples from the Okanagan Valley (right). (C) BAYESASS analysis of samples collected north and west of Alberta (n = 478) detected recent gene flow from Yukon west into Alaska, and from the western Yukon west to the cluster concentrated on the Montane Cordillera and surrounding areas (6–10, 13–15), but did not detect substantial recent gene flow from Haida Gwaii to the mainland. Arrow thickness is proportional to estimated gene flow. (D) Discriminant analysis of principal components, using sampling areas as prior groups and (E) spatial analysis of principal components generally support the clusters identified by STRUCTURE. Map in (E) shows interpolation of the first lagged score only, which by definition can divide the samples into no more than two groups; note that bats from western Yukon are not strongly assigned to either of these most differentiated groups. 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).
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).
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 Table 1. Map created using Natural Earth (naturalearthdata.com).
Estimates of contemporary Ne 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 Ne estimates for western populations were >1600.
Table 3.
Table 3. Effective population size (mean Ne, 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.
Lowest allele frequency used
Haida Gwaii, n = 57
Alaska/Yukon West, n = 198
Montane Cordillera, n = 295
Prairie (admixed), n = 150
Eastern Canada, n = 442
0.050
289.5 (123.9–∞)
498.2 (308.9–1116.6)
1519.2 (582.1–∞)
479.4 (248.9–2672.8)
1682.9 (794.0–261 726.6)
0.020
606.1 (223.1–∞)
845.9 (566.3–1592.2)
822.4 (628.8–1165.9)
573.7 (414.4–907.3)
16 062.4 (3465.0–∞)
0.010
305.8 (182.2–848.9)
998.7 (718.0–1604.0)
457.1 (407.2–518.7)
585.3 (471.0–766.1)
27 778.4 (5519.3–∞)
0+
175.6 (130.1–263.2)
1642.2 (1179.3–2665.1)
733.2 (669.6–808.7)
943.1 (759.0–1238.1)
1034.3 (964.1–1114.5)
Note: Estimates using the lowest allele frequencies are italicized to indicate the higher probability that sampling bias will affect the estimates.
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 (r2 = 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.
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.
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 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 Yukon and Alaska may experience the same delay in P. destructans arrival that has 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 Ne. 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 Ne 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 Ne 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; C.L. Lausen and T.S. Jung, unpublished data). Future research can test potential behavioural or morphological differences 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 mortality 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.
Acknowledgements
This study was supported in part by a Liber Ero Fellowship to CMD. CLL received funding for field sampling of M. lucifugus from Waterton-Glacier Peace Park; BC Habitat Conservation Trust Foundation; Fish and Wildlife Compensation Program Columbia; Flathead Wild; Alaska Department of Fish and Game. MHC analysis of MHC was supported by the Government of Ontario (SARRFO grant to CJK/CMD). We are grateful to the following for providing DNA samples from M. lucifugus: Lesley Hale, Amy Cameron, and Heather Riddell (OMNRF), Lenny Shirose (CWHC), Jesika Reimer, and Joanna Wilson (Government of the Northwest Territories), the Royal Alberta Museum, and the Alaska Museum of the North. Sample processing at the New Brunswick Museum was provided (to DFM/KJV) from the Canadian Wildlife Federation, New Brunswick Wildlife Trust Fund, New Brunswick Environmental Trust Fund, New Brunswick Department of Natural Resources, and Parks Canada. Erica Newton assisted with figure production. Finally, we thank the many field assistants and volunteers who made it possible to capture and sample the bats used in this study.
References
Bataille A, Cashins SD, Grogan L, Skerratt LF, Hunter D, McFadden M, et al. 2015. Susceptibility of amphibians to chytridiomycosis is associated with MHC class II conformation. Proceedings of the Royal Society of London B: Biological Sciences, 282(1805): 20143127.
Blanchong JA, Samuel MD, Scribner KT, Weckworth BV, Langenberg JA, and Filcek KB. 2008. Landscape genetics and the spatial distribution of chronic wasting disease. Biology Letters, 4(1): 130–133.
Blehert DS, Hicks AC, Behr M, Meteyer CU, Berlowski-Zier BM, and Buc EL. 2009. Bat white-nose syndrome: an emerging fungal pathogen? Science, 323(5911): 227.
Blejwas K, Michael K, and Beard L. 2015. Little brown bats in southeast Alaska hibernate in holes: implications for the spread of white-nose syndrome. Bat Research News, 56(4): 70.
Bower SM, and Woo PTK. 1981. Two new species of trypanosomes (subgenus Schizotrypanum) in bats from southern Ontario. Canadian Journal of Zoology, 59(3): 530–545.
Brown JA, McAlpine DF, and Curley R. 2007. Northern Long-eared Bat, Myotis septentrionalis (Chiroptera: Vespertilionidae), on Prince Edward Island: first records of occurrence and over-wintering. Canadian Field-Naturalist, 121: 208–209.
Burns LE, Broders HG, and Frasier TR. 2012. Characterization of 11 tetranucleotide microsatellite loci for the little brown bat (Myotis lucifugus) based on in silico genome sequences. Conservation Genetics Resources, 4(3): 653–655.
Burns LE, Frasier TR, and Broders HG. 2014. Genetic connectivity among swarming sites in the wide ranging and recently declining little brown bat (Myotis lucifugus). Ecology and Evolution, 4(21): 4130–4149.
Byun SA, Koop BF, and Reimchen TE. 1997. North American black bear mtDNA phylogeography: implications for morphology and the Haida Gwaii glacial refugium controversy. Evolution, 51(5): 1647–1653.
Campbell MC, and Tishkoff SA. 2008. African genetic diversity: implications for human demographic history, modern human origins, and complex disease mapping. Annual Review of Genomics and Human Genetics, 9: 403–433.
Chiucchi JE, and Gibbs HL. 2010. Similarity of contemporary and historical gene flow among highly fragmented populations of an endangered rattlesnake. Molecular Ecology, 19: 5345–5358.
Davy CM, Martinez-Nunez F, Willis CKR, and Good SV. 2015. Spatial genetic structure among bat hibernacula along the leading edge of a rapidly spreading pathogen. Conservation Genetics, 16: 1013–1024.
Davy CM, Mastromonaco GF, Riley JL, Baxter-Gilbert JH, Mayberry H, and Willis CKR. 2017. Physiological carry-over effects in bats exposed to an emerging pathogen implicate an overlooked potential driver of the extinction vortex. Conservation Biology, 31: 615–624.
de Groot NG, Otting N, Doxiadis GGM, Balla-Jhagjhoorsingh SS, Heeney JL, van Rood JJ, et al. 2002. Evidence for an ancient selective sweep in the MHC class I gene repertoire of chimpanzees. Proceedings of the National Academy of Sciences of the United States of America, 99(18): 11748–11753.
Deschamps M, Laval G, Fagny M, Itan Y, Abel L, Casanova J-L, et al. 2016. Genomic signatures of selective pressures and introgression from archaic hominins at human innate immunity genes. The American Journal of Human Genetics, 98(1): 5–21.
Di Gaspero G, Copetti D, Coleman C, Castellarin SD, Eibach R, Kozma P, et al. 2012. Selective sweep at the Rpv3 locus during grapevine breeding for downy mildew resistance. Theoretical and Applied Genetics, 124(2): 277–286.
Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, and Ovenden JR. 2014. NeEstimator v2: re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Molecular Ecology Resources, 14: 209–214.
Ekblom R, Saether SA, Jacobsson P, Fiske P, Sahlman T, Grahn M, et al. 2007. Spatial pattern of MHC class II variation in the great snipe (Gallinago media). Molecular Ecology, 16: 1439–1451.
Evanno G, Regnaut S, and Goudet J. 2005. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Molecular Ecology, 14(8): 2611–2620.
Excoffier L, Laval G, and Schneider S. 2005. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evolutionary Bioinformatics Online, 1: 47–50.
Falush D, Stephens M, and Pritchard JK. 2003. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics, 164: 1567–1587.
Fenton MB. 1969. Summer activity of Myotis lucifugus (Chiroptera: Vespertilionidae) at hibernacula in Ontario and Quebec. Canadian Journal of Zoology, 47: 597–602.
Francl K, Ford WM, Sparks DW, and Brack V Jr. 2012. Capture and reproductive trends in summer bat communities in West Virginia: assessing the impact of white-nose syndrome. Journal of Fish and Wildlife Management, 3: 33–42.
Frick WF, Puechmaille SJ, Hoyt JR, Nickel BA, Langwig KE, Foster JT, et al. 2015. Disease alters macroecological patterns of North American bats. Global Ecology and Biogeography, 24(7): 741–749.
Gilbert KJ, and Whitlock MC. 2015. Evaluating methods for estimating local effective population size with and without migration. Evolution, 69(8): 2154–2166.
Hale ML, Burg TM, and Steeves TE. 2012. Sampling for microsatellite-based population genetic studies: 25 to 30 individuals per population is enough to accurately estimate allele frequencies. PLoS ONE, 7(9): e45170.
Hawley DM, and Fleischer RC. 2012. Contrasting epidemic histories reveal pathogen-mediated balancing selection on class II MHC diversity in a wild songbird. PLoS ONE, 7(1): e30222.
Henry M, Thomas DW, Vaudry R, and Carrier M. 2002. Foraging distances and home range of pregnant and lactating little brown bats (Myotis lucifugus). Journal of Mammalogy, 83(3): 767–774.
Jakobsson M, and Rosenberg NA. 2007. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics, 23: 1801–1806.
Johnson JB, Roberts JH, King TL, Edwards JW, Ford WM, and Ray DA. 2014. Genetic structuring of northern myotis (Myotis septentrionalis) at multiple spatial scales. Acta Theriologica, 59(2): 223–231.
Jombart T, Devillard S, and Balloux F. 2010. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genetics, 11: 94.
Jones OR, and Wang J. 2010. colony: a program for parentage and sibship inference from multilocus genotype data. Molecular Ecology Resources, 10(3): 551–555.
Jung TS, Blejwas KM, Lausen CL, Wilson JM, and Olson LE. 2014. Concluding remarks: what do we need to know about bats in northwestern North America? Northwestern Naturalist, 95: 318–330.
Keever CC, Sunday J, Puritz JB, Addison JA, Toonen RJ, Grosberg RK, et al. 2009. Discordant distribution of populations and genetic variation in a sea star with high dispersal potential. Evolution, 63(12): 3214–3227.
Kyle CJ, Rico Y, Castillo S, Srithayakumar V, Cullingham CI, White BN, et al. 2014. Spatial patterns of neutral and functional genetic variations reveal patterns of local adaptation in raccoon (Procyon lotor) populations exposed to raccoon rabies. Molecular Ecology, 23(9): 2287–2298.
Lausen CL, Delisle I, Barclay RMR, and Strobeck C. 2008. Beyond mtDNA: nuclear gene flow suggests taxonomic oversplitting in the little brown bat (Myotis lucifugus). Canadian Journal of Zoology, 86(7): 700–713.
Li H. 2013. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv: 1303.3997 [online]: Available from arxiv.org/abs/1303.3997.
Lorch JM, Palmer JM, Lindner DL, Ballmann AE, George KG, Griffin K, et al. 2016. First detection of bat white-nose syndrome in western North America. mSphere, 1(4): e00148-16.
Maher SP, Kramer AM, Pulliam JT, Zokan MA, Bowden SE, Barton HD, et al. 2012. Spread of white-nose syndrome on a network regulated by geography and climate. Nature Communications, 3(1306): 1–8.
Maslo B, Valent M, Gumbs JF, and Frick WF. 2015. Conservation implications of ameliorating survival of little brown bats with white-nose syndrome. Ecological Applications, 25(7): 1832–1840.
Mayer F, and Brunner A. 2007. Non-neutral evolution of the major histocompatibility complex (MHC) class II gene DRB1 in the sac-winged bat Saccopteryx bilineata. Heredity, 99: 257–264.
Middleton D, Johnson KO, Rosatte RC, Hobbs JL, Moore SR, Rosella L, et al. 2015. Human rabies post-exposure prophylaxis and animal rabies in Ontario, Canada, 2001–2012. Zoonoses and Public Health, 62(5): 356–364.
Miller-Butterworth CM, Vonhof MJ, Rosenstern J, Turner GG, and Russell AL. 2014. Genetic structure of little brown bats (Myotis lucifugus) corresponds with spread of white-nose syndrome among hibernacula. Journal of Heredity, 105(3): 354–364.
Misra V, Dumonceaux T, Dubois J, Willis C, Nadin-Davis S, Severini A, et al. 2009. Detection of polyoma and corona viruses in bats of Canada. Journal of General Virology, 90(Pt 8): 2015–2022.
Norquay KJO, Martinez-Nuñez F, Dubois JE, Monson KM, and Willis CKR. 2013. Long-distance movements of little brown bats (Myotis lucifugus). Journal of Mammalogy, 94(2): 506–515.
O’Regan SM, Magori K, Pulliam JT, Zokan MA, Kaul RB, Barton HD, et al. 2015. Multi-scale model of epidemic fade-out: will local extirpation events inhibit the spread of white-nose syndrome? Ecological Applications, 25(3): 621–633.
Peakall R, and Smouse PE. 2012. GenALEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics, 28(19): 2537–2539.
Pompanon F, Bonin A, Bellemain E, and Taberlet P. 2005. Genotyping errors: causes, consequences and solutions. Nature Reviews Genetics, 6(11): 847–846.
Pond BA, Brown GS, Wilson KS, and Schaefer JA. 2016. Drawing lines: spatial behaviours reveal two ecotypes of woodland caribou. Biological Conservation, 194: 139–148.
Pruett CL, Topp CM, Maley JM, McCracken KG, Rohwer S, Birks S, et al. 2013. Evidence from the genetics of landbirds for a forested pleistocene glacial refugium in the Haida Gwaii area. The Condor, 115(4): 725–737.
Puechmaille SJ. 2016. The program structure does not reliably recover the correct population structure when sampling is uneven: subsampling and new estimators alleviate the problem. Molecular Ecology Resources, 16: 608–627.
Razgour O, Rebelo H, Puechmaille SJ, Juste J, Ibáñez C, Kiefer A, et al. 2014. Scale-dependent effects of landscape variables on gene flow and population structure in bats. Diversity and Distributions, 20: 1173–1185.
Real-Monroy MD, Martínez-Méndez N, and Ortega J. 2014. MHC-DRB Exon 2 diversity of the Jamaican fruit-eating bat (Artibeus jamaicensis) from Mexico. Acta Chiropterologica, 16(2): 301–314.
Reimer JP, Lausen CL, Barclay RMR, Irwin S, and Vassal MK. 2014. Bat activity and use of hibernacula in Wood Buffalo National Park, Alberta. Northwestern Naturalist, 95(3): 277–288.
Richman AD, Herrera MLG, Ortega-García S, Flores-Martínez JJ, Arroyo-Cabrales J, and Morales-Malacara JB. 2010. Class II DRB polymorphism and sequence diversity in two vesper bats in the genus Myotis. International Journal of Immunogenetics, 37(5): 401–405.
Rico Y, Ethier DM, Davy CM, Sayers J, Weir RD, Swanson BJ, et al. 2016. Spatial patterns of immunogenetic and neutral variation underscore the conservation value of small, isolated American badger populations. Evolutionary Applications, 9(10): 1271–1284.
Rollins LA, Moles AT, Lam S, Buitenwerf R, Buswell JM, Brandenburger CR, et al. 2013. High genetic diversity is not essential for successful introduction. Ecology and Evolution, 3(13): 4501–4517.
Salmier A, de Thoisy B, Crouau-Roy B, Lacoste V, and Lavergne A. 2016. Spatial pattern of genetic diversity and selection in the MHC class II DRB of three Neotropical bat species. BMC Evolutionary Biology, 16(1): 229.
Savage AE, and Zamudio KR. 2011. MHC genotypes associate with resistance to a frog-killing fungus. Proceedings of the National Academy of Sciences, 108(40): 16705–16710.
Schad J, Dechmann DK, Voigt CC, and Sommer S. 2012a. Evidence for the ‘good genes’ model: association of MHC class II DRB alleles with ectoparasitism and reproductive state in the neotropical Lesser Bulldog Bat, Noctilio albiventris. PLoS ONE, 7(5): e37101.
Schad J, Voigt CC, Greiner S, Dechmann DK, and Sommer S. 2012b. Independent evolution of functional MHC class II DRB genes in New World bat species. Immunogenetics, 64: 535–547.
Schowalter DB. 1980. Swarming, reproduction, and early hibernation of Myotis lucifugus and M. volans in Alberta, Canada. Journal of Mammalogy, 61(2): 350–354.
Sironi M, Cagliani R, Forni D, and Clerici M. 2015. Evolutionary insights into host-pathogen interactions from mammalian sequence data. Nature Reviews Genetics, 16(4): 224–236.
Sovic MG, Carstens BC, and Gibbs HL. 2016. Genetic diversity in migratory bats: results from RADseq data for three tree bat species at an Ohio windfarm. PeerJ, 4: e1647.
Spurgin LG, and Richardson DS. 2010. How pathogens drive genetic diversity: MHC, mechanisms and misunderstandings. Proceedings of the Royal Society B: Biological Sciences, 277: 979–988.
Talerico JM. 2008. The behaviour, diet and morphology of the Little Brown Bat (Myotis lucifugus) near the northern extent of its range in Yukon, Canada. M.Sc. thesis. University of Calgary, Calgary, Alberta. 104 p.
van der Meer MH, Horne JB, Gardner MG, Hobbs J-PA, Pratchett M, and van Herwerden L. 2013. Limited contemporary gene flow and high self-replenishment drives peripheral isolation in an endemic coral reef fish. Ecology and Evolution, 3: 1653–1666.
van Oosterhout C, Hutchinson WFD, Wills DPM, and Shipley P. 2004. micro-checker: software for identifying and correcting genotyping errors in microsatellite data. Molecular Ecology Notes, 4: 535–538.
Veselka N, McGuire LP, Dzal YA, Hooton LA, and Fenton MB. 2013. Spatial variation in the echolocation calls of the little brown bat (Myotis lucifugus). Canadian Journal of Zoology, 91: 795–801.
Vonhof MJ, Russell AL, and Miller-Butterworth CM. 2015. Range-wide genetic analysis of little brown bat (Myotis lucifugus) populations: estimating the risk of spread of white-nose syndrome. PLoS ONE, 10(7): e0128713.
Vonhof MJ, Amelon SK, Currie RR, and McCracken GF. 2016. Genetic structure of winter populations of the endangered Indiana bat (Myotis sodalis) prior to the white nose syndrome epidemic: implications for the risk of disease spread. Conservation Genetics, 17(5): 1025–1040.
Waples RS, and England PR. 2011. Estimating contemporary effective population size on the basis of linkage disequilibrium in the face of migration in the face of migration. Genetics, 189(2): 633–644.
Webber QMR, McGuire LP, Smith SB, and Willis CKR. 2015. Host behaviour, age and sex correlate with ectoparasite prevalence and intensity in a colonial mammal, the little brown bat. Behaviour, 152(1): 83–105.
Wegner KM, Reusch TBH, and Kalbe M. 2003. Multiple parasites are driving major histocompatibility complex (MHC) polymorphism in the wild. Journal of Evolutionary Biology, 16: 224–232.
Wilder AP, Kunz TH, and Sorenson MD. 2015. Population genetic structure of a common host predicts the spread of white-nose syndrome, an emerging infectious disease in bats. Molecular Ecology, 24: 5495–5506.
Wright SK, and Moran JR. 2011. Ocean-going vessels: a possible conduit for the introduction of white-nose syndrome fungus (Geomyces destructans) into bats in Alaska. Northwestern Naturalist, 92(2): 133–135.
Zhang L, Wu Q, Hu Y, Wu H, and Wei F. 2015. Major histocompatibility complex alleles associated with parasite susceptibility in wild giant pandas. Heredity, 114(1): 85–93.
Zukal J, Bandouchova H, Brichta J, Cmokova A, Jaron KS, Kolarik M, et al. 2016. White-nose syndrome without borders: Pseudogymnoascus destructans infection tolerated in Europe and Palearctic Asia but not in North America. Scientific Reports, 6(19829): 1–13.
All relevant data are within the paper, in the Supplementary Material, and in the Dryad data repository. The genotype data file containing microsatellite genotypes for 11 loci (Myotis lucifugus) is accessible at doi:https://doi.org/10.5061/dryad.h7n25.
Canadian Wildlife Health Cooperative, Atlantic Region, Atlantic Veterinary College, University of Prince Edward Island, 550 University Avenue, Charlottetown, PEI C1A 4P3, Canada
Canadian Wildlife Health Cooperative, Atlantic Region, Atlantic Veterinary College, University of Prince Edward Island, 550 University Avenue, Charlottetown, PEI C1A 4P3, Canada
CMD, MED, YR, CLL, KD, KR, CKRW, DWB, TSJ, SM, AP, DFM, and KJV performed the experiments/collected the data.
CMD, MED, YR, KD, and KR analyzed and interpreted the data.
CMD, CLL, CKRW, SM, DFM, and CJK contributed resources.
All drafted or revised the manuscript.
Competing Interests
The authors have declared that no competing interests exist.
Metrics & Citations
Metrics
Other Metrics
Citations
Cite As
Christina M. Davy, Michael E. Donaldson, Yessica Rico, Cori L. Lausen, Kathleen Dogantzis, Kyle Ritchie, Craig K.R. Willis, Douglas W. Burles, Thomas S. Jung, Scott McBurney, Allysia Park, Donald F. McAlpine, Karen J. Vanderwolf, and Christopher J. Kyle. 2017. Prelude to a panzootic: Gene flow and immunogenetic variation in northern little brown myotis vulnerable to bat white-nose syndrome. FACETS.
2(): 690-714. https://doi.org/10.1139/facets-2017-0022
Export Citations
If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.
Cited by
1. COULD WHITE-NOSE SYNDROME MANIFEST DIFFERENTLY IN MYOTIS LUCIFUGUS IN WESTERN VERSUS EASTERN REGIONS OF NORTH AMERICA? A REVIEW OF FACTORS
2. Effects of capture on stress-axis measures in endangered little brown bats (
Myotis lucifugus
)
3. Foraging habitat drives the distribution of an endangered bat in an urbanizing boreal landscape
4. Major histocompatibility complex variation is similar in little brown bats before and after white‐nose syndrome outbreak
5. Genome-Wide Changes in Genetic Diversity in a Population of
Myotis lucifugus
Affected by White-Nose Syndrome
6. Bats in the changing boreal forest: response to a megafire by endangered little brown bats (
Myotis lucifugus
)
7. Life in a northern town: rural villages in the boreal forest are islands of habitat for an endangered bat
8. White-nose syndrome detected in bats over an extensive area of Russia
