Range-Wide Genetic Analysis of Little Brown Bat (Myotis lucifugus) Populations: Estimating the Risk of Spread of White-Nose Syndrome

The little brown bat (Myotis lucifugus) is one of the most widespread bat species in North America and is experiencing severe population declines because of an emerging fungal disease, white-nose syndrome (WNS). To manage and conserve this species effectively it is important to understand patterns of gene flow and population connectivity to identify possible barriers to disease transmission. However, little is known about the population genetic structure of little brown bats, and to date, no studies have investigated population structure across their entire range. We examined mitochondrial DNA and nuclear microsatellites in 637 little brown bats (including all currently recognized subspecific lineages) from 29 locations across North America, to assess levels of genetic variation and population differentiation across the range of the species, including areas affected by WNS and those currently unaffected. We identified considerable spatial variation in patterns of female dispersal and significant genetic variation between populations in eastern versus western portions of the range. Overall levels of nuclear genetic differentiation were low, and there is no evidence for any major barriers to gene flow across their range. However, patterns of mtDNA differentiation are highly variable, with high ΦST values between most sample pairs (including between all western samples, between western and eastern samples, and between some eastern samples), while low mitochondrial differentiation was observed within two groups of samples found in central and eastern regions of North America. Furthermore, the Alaskan population was highly differentiated from all others, and western populations were characterized by isolation by distance while eastern populations were not. These data raise the possibility that the current patterns of spread of WNS observed in eastern North America may not apply to the entire range and that there may be broad-scale spatial variation in the risk of WNS transmission and occurrence if the disease continues to spread west.

637 little brown bats (including all currently recognized subspecific lineages) from 29 locations across North America, to assess levels of genetic variation and population differentiation across the range of the species, including areas affected by WNS and those currently unaffected. We identified considerable spatial variation in patterns of female dispersal and significant genetic variation between populations in eastern versus western portions of the range. Overall levels of nuclear genetic differentiation were low, and there is no evidence for any major barriers to gene flow across their range. However, patterns of mtDNA differentiation are highly variable, with high Φ ST values between most sample pairs (including between all western samples, between western and eastern samples, and between some eastern samples), while low mitochondrial differentiation was observed within two groups of samples found in central and eastern regions of North America. Furthermore, the Alaskan population was highly differentiated from all others, and western populations were characterized by isolation by distance while eastern populations were not. These data raise the possibility that the current patterns of spread of WNS observed in eastern North America may not apply to the entire range and that there may be broad-scale spatial variation in the risk of WNS transmission and occurrence if the disease continues to spread west.

Introduction
Understanding how host movement patterns influence the transmission of pathogens is critical to the development of effective prevention and control strategies, and to the conservation and management of host populations during and after disease outbreaks. However, for many host species, data on individual movements and contact rates are difficult or impossible to collect because of cryptic behavior, the geographic scale of movements, or methodological considerations that limit our ability to follow individuals through time and space. Evidence from empirical studies employing population and landscape genetic approaches has demonstrated that landscape features, such as mountains and rivers that limit host gene flow, often represent barriers to disease transmission [1][2][3][4][5][6], although alternative mechanisms of pathogen dispersal, including humans and other highly mobile intermediate hosts, may override the influence of primary host population genetic structure [1]. Nevertheless, where they exist, such barriers to host gene flow can have a dramatic impact on initial disease establishment, the rate and direction of disease spread, spatial patterns of host resistance, and dynamics and genetic structure of pathogen populations [1][2][3][4][5][6]. Assuming that rates of contact among individuals leading to gene flow are indicative of contacts that could result in disease transmission, genetic methods provide a useful alternative to traditional demographic approaches as a means of examining host movements and their impact on disease transmission [1].
White-nose syndrome (WNS) is an emerging fungal disease causing high levels of mortality in hibernating North American bats [7][8][9]. The causative agent, Pseudogymnoascus destructans (hereafter Pd), is a cold-loving fungus that affects bats during hibernation and subsequent arousal. Pd causes characteristic cup-like erosions of the epidermis of the wings and muzzle, and may invade sebaceous glands and hair follicles [10]. Mortality of bats likely occurs through the loss of physiological homeostasis [11], possibly associated with dehydration and electrolyte depletion [12,13], leading to more frequent arousal behavior and premature loss of fat reserves [14,15]. Since it was first discovered in New York State during the winter of 2006-07, WNS has since spread to 27 additional states and five Canadian provinces, and is known to affect at least seven species of hibernating bats [16]. Mortality rates vary considerably among species but can be very high (>90% for little brown bats, Myotis lucifugus, and northern long-eared bats, M. septentrionalis [9]), and cumulative mortality of all affected bat species has been estimated at 5.7 to 6.7 million individuals as of January 2012 [17]. The rapid emergence, and the geographic and taxonomic spread of the disease have raised serious concerns about the long-term survival of hibernating bat species in eastern North America, and have highlighted our lack of knowledge of the factors that may influence WNS transmission and spread to currently unaffected regions.
Little brown bats were among the first species to be diagnosed with WNS [7], and population models indicate that if mortality rates stay constant, this species could be extirpated from the northeastern United States within 16 years [18]. Hibernating populations of all sizes have been affected by WNS, but the probability of infection increases with increasing colony size [19,20], although mortality within populations is density-independent and characterized by frequency-dependent transmission [21]. Thus, there is a high probability that little brown bat populations in currently affected regions will be highly reduced and possibly be extirpated in coming decades. Additionally, WNS may pose a threat to the entire species if the disease continues to spread across the species' range. Further, because little brown bats are one of two affected bat species whose geographic ranges span temperate North America, they may drive transmission of Pd to a novel suite of western North American hibernating bat species that might otherwise remain geographically isolated from the disease. Therefore, it is crucial that we understand the probability of Pd transmission across the range of little brown bats, and whether there are barriers to gene flow that could restrict the geographic spread of WNS.
Here we apply genetic approaches to understand levels of gene flow and population connectivity in the little brown bat. This small (6-10 g) insectivorous bat species is among the most widespread (Fig 1) and well-studied in North America [22,23]. During the summer, reproductive females form maternity colonies in buildings, trees, or crevices where parturition and postnatal care take place, while males and non-reproductive females typically roost solitarily [22]. In winter, both sexes congregate in hibernacula, and mating takes place during the pre-hibernation swarming period, or during hibernation itself [24,25]. The size of hibernating populations may vary considerably, on the order of 10's to 100,000's, but most of the larger known hibernacula occur in karst regions of eastern North America, and very little is known about the distribution or size of hibernacula in western North America. Because individuals from many breeding groups come together at swarming or hibernation sites with males that may or may not have originated from the same breeding group [26,27], these sites have been suggested to represent 'hot spots' of gene flow for temperate bats [28][29][30]. Thus, patterns of gene flow will represent the interplay of movements of individuals between summer and/or winter populations, and levels and spatial patterns of connectivity among summer and winter populations that determine the composition of mating aggregations.
There are currently five recognized subspecies of little brown bats (M. l. alascensis, M. l. carissima, M. l. lucifugus, M. l. pernox, and M. l. relictus [22,31]; see Fig 1) based on morphology, but the extent to which these subspecies diverge genetically is unclear. Coalescent analyses of nuclear DNA (nucDNA) and mitochondrial DNA (mtDNA) suggest that some subspecies may represent independent evolutionary lineages, but that M. lucifugus may be paraphyletic with respect to the western long-eared bat, M. evotis [32]. In addition, two distinct mtDNA lineages (corresponding to M. l. carissima and M. l. lucifugus) co-occur in southern Alberta and northcentral Montana, but these two groups of bats are not differentiated based on nuclear microsatellite DNA or morphology, suggesting that the subspecies in question may interbreed and exchange genes [33]. A single mitochondrial lineage corresponding to M. l. lucifugus was observed in the Minnesota populations, and there was a strong signal of population expansion dating to 18 kya [34]. Environmental niche modeling based on conditions during the Last Glacial Maximum (LGM) indicated the presence of a single large refugium extending across the southeastern and south-central United States, and more fragmented refugia in the southern portion of the mountainous western United States [34], suggesting a possible mechanism for lineage differentiation within this species where separation into disjunct glacial refugia was followed by subsequent post-glacial range expansion and secondary contact.
Few studies have examined genetic variation in little brown bats, and there has been no comprehensive range-wide population genetic analysis of this species. Fine-scale genetic studies in Minnesota described high levels of mtDNA structure and weak but significant nucDNA structure among maternity colonies [35], a pattern consistent with many other temperate bat species characterized by female philopatry to summer breeding habitat and extensive gene flow via mating at swarming sites during the autumn [36][37][38][39][40][41]. In southestern Canada, high levels of genetic connectivity were identified among swarming sites, however minimal structuring at both mtDNA and nucDNA markers suggested dispersal may occur in both sexes, although it may be male-biased [42]. In Pennsylvania, no significant nuclear structure was identified among hibernating populations, but these populations were structured matrilineally. This mtDNA structure was correlated with local topography, which may have delayed the spread of WNS to western parts of the state [6].
The rapid spread of WNS through eastern North American populations of little brown bats (and other affected species) suggests that few barriers to transmission exist within the current range of the disease. Here we utilize mtDNA sequence and nucDNA microsatellite variation from a large sample of little brown bats collected across the range of the species to address the following objectives: 1) assess levels of genetic variation in little brown bat populations, including areas affected by WNS and those currently unaffected; 2) quantify genetic differentiation among populations sampled across the range of the species, including populations in eastern North America within the current range of WNS, as well as additional populations situated both east and west of the transition between the Great Plains and Rocky Mountains; and 3) assess the current geographic distribution of and levels of genetic differentiation among currently-recognized subspecific lineages.
There are few physiographic barriers that would limit movement of highly vagile organisms east of the Rocky Mountains. Phylogeographic studies of widespread bats and birds in North was not included in mtDNA analyses. Population abbreviations are detailed in Table 1 [58]), and the presence of three distinct clades within M. lucifugus, based on maximum likelihood analysis of partial COI sequences in PhyML 3.0. A member of the Neotropical Myotis clade (M. austroriparius) was included as the outgroup. Leaves are collapsed to highlight well-supported clades, and the vertical dimension of the triangles is proportional to the number of samples included. SH-like branch support values are provided for all major clades. Clades containing M. lucifugus are designated by the specific abbreviation followed by the subspecies name (e.g., M. l. lucifugus refers to the nominal subspecies). Note that one clade (including M. l. carissima) also contains members of other species (including M. evotis, M. keenii, and M. thysanodes) as previously described [32]. America typically report little differentiation among populations within eastern and central portions of North America, significant differentiation among eastern and western populations, and higher levels of differentiation among populations within the mountainous west [41,[43][44][45][46]. We predict that the Rocky Mountains will represent a barrier to gene flow, and that we will therefore observe genetic differentiation between sample sites east versus west of the Great Plains-Rocky Mountains transition. Further, because of higher topographic variability, we predict higher levels of genetic differentiation among sample sites in the mountainous west than in eastern North America. If subspecific lineages represent reproductively-isolated units that arose during the LGM, then we predict that patterns of differentiation at both nucDNA and mtDNA markers will match the described geographic distribution of subspecies. Our study provides valuable data on population connectivity and hence opportunities for WNS transmission across the range of little brown bats that may be used to inform the management and conservation of affected species.

Sample collection
Tissue samples were obtained during the summer (between May and August) from 637 individuals at 29 locations across the range of little brown bats (Table 1, S1 Table, and Fig 1). Two 3 mm biopsy punches, one from each wing, were taken from each bat and stored in 5 M NaCl with 20% DMSO [47]. The bats were released after sampling. The majority of population samples were collected at maternity colonies (N = 16) or single or several closely-spaced (<10 km) netting sites (N = 12). However, the Idaho sample constituted bats collected in 8 different counties in the southeastern portion of the state. When samples came from more than one capture location, centroids were calculated and used as approximate sample locations.

Ethics statement
One author (MJV) collected samples for this study from one population in Michigan. The samples were collected under permission granted by the State of Michigan Department of Natural Resources (Permit SC 1257), and the methods were approved by the Western Michigan University Institutional Animal Care and Use Committee Protocol Number 08-05-03. All other samples were collected by university and government researchers performing other research who were required to have appropriate permits and other necessary permission to undertake their work.

Mitochondrial DNA sequencing and analysis
Total genomic DNA was isolated using DNeasy Tissue Kits (Qiagen, Valencia CA). We amplified and sequenced a 636 bp fragment of the mitochondrial cytochrome c oxidase subunit I (COI) gene using primers HCO2198 and LCO1490 [48] or primers VF1 and VR1 [49]. Bats from all sample sites except Michigan were sequenced, for a total of 617 individuals. PCRs were conducted in 25 μl volumes containing 0.4 μM of each primer and 20-50 ng of DNA template, using Illustra PuReTaq Ready-To-Go PCR beads (GE Healthcare Life Sciences, Pittsburgh PA). When reconstituted to 25 μl with water, these beads contained 2.5 units PuReTaq DNA polymerase, 200 μM each dNTP in 10 mM Tris-HCl (pH 9), 50 mM KCl, 1.5 mM MgCl 2 , and an unspecified concentration of bovine serum albumin (BSA). No other additives were added to the solution. Cycling conditions consisted of one cycle of 5 min at 94°C, 30 cycles of 30 sec at 94°C, 45 sec at 68°C and 1 min at 72°C, and a final cycle of 2 min at 72°C. PCR products were purified by digestion with exonuclease I and shrimp alkaline phosphatase (EXOSAP), and were sequenced in both directions, using the amplification primers, at the University of Arizona Genetics Core Facility. Sequences were edited using CodonCode Aligner 3.0 (Gene Codes Corp.) and aligned using the default settings in MAFFT [50].

Microsatellite genotyping and analysis
We genotyped individuals at eleven highly variable microsatellite loci using primers previously developed for other vespertilionid bats (IBat CA5, CA11, CA43, CA47, and M23 [51]; MS3D02 and MS3F05 [52]; E24 and G9 [53]; Cora_F11_C04 [54]; Coto_G02F_H10R [55]). We did not genotype population samples with 15 individuals (CA-Ma, CA-Mo, CA-Sh, MB-1, MB-2; Table 1), with the exception of BC-N that had been genotyped for another study. Based on preliminary results of low levels of differentiation in eastern North America and to reduce costs, we did not genotype individuals from three additional sites east of the Rocky Mountains (ON-2, WV, and NJ-Sa). The total number of individuals genotyped was 510. Amplifications were carried out in four multiplex reactions and two single-locus amplifications (see S2 Table), and were subsequently pooled into three different loads for fragment analysis on an ABI 3130 sequencer. The basic cycling conditions consisted of 1 min at 94°C, three cycles of 30 sec at 94°C, 20 sec at T a (54 or 60°C), and 5 sec at 72°C, 33 cycles of 15 sec at 94°C, 20 sec at T a , and 10 sec at 72°C, followed by a final extension at 72°C for 30 min. Some amplifications required additional cycles or the removal of the final extension step (S2 Table). Fragments were analyzed and scored using GeneMarker software (SoftGenetics LLC, State College, PA).

Mitochondrial DNA analysis
To describe overall levels of mtDNA diversity within populations, we calculated haplotype (h) and nucleotide (π) diversities in DNASP v.5.10.1 [56], and determined the number of private haplotypes in each site after collapsing sequences from the entire dataset to unique haplotypes using FABOX [57]. We used several population genetic approaches to establish whether current patterns of variation are indicative of the presence of distinct genetic clusters. We calculated pairwise F ST values between sites and tested for significance with 10,000 permutations in Arlequin v.3.11 [58] to identify pairs of sites that were genetically distinct. We also performed an analysis of molecular variance (AMOVA [59]) to describe the relative amount of genetic variation within and among sites. Based on initial pairwise F ST results, we then performed nested AMOVAs to identify natural groups of sites. Sites were initially grouped together if they had low pairwise F ST values, and the analysis was rerun. Any ambiguous sites (sites that had low F ST values with sites in more than one group) were sequentially moved between groups and the analysis was rerun. All logical combinations were tested to identify the grouping that minimized among-site/within-group variation and maximized between-group variation.
To test the significance of defined subspecific lineages within M. lucifugus using our nationwide dataset, we used a maximum likelihood phylogenetic approach implemented in PHYML  S3 Table for list of specimens). We used the best fit model of sequence evolution (HKY +G) as determined using Mega v.5.0 [62], with the gamma distribution of variability of rates among sites calculated empirically from the data, SPR moves to explore tree space, and SH-Like procedure to assess branch supports [60]. The proportion of each sampled population falling within each subspecific clade was then calculated and plotted on a map produced in ARC-GIS v.10.1 to visualize the geographic distribution of the clades.

Microsatellite DNA analysis
Deviations from Hardy-Weinberg equilibrium (HWE) were estimated for each locus, and loci were confirmed to be in linkage equilibrium using FSTAT v.2.9.3 [63]. To test for differences in levels of genetic diversity among sites and regions, several indices of nuclear genetic diversity were estimated, including number of alleles per locus, allelic richness, and the inbreeding coefficient (F IS ) using FSTAT, private allelic richness using HP-RARE v.1.0 [64], and observed and expected heterozygosity using GENODIVE [65]. We then tested for differences among sites (or groups of sites) in allelic richness and F IS in FSTAT, and expected heterozygosity in GEN-ODIVE, using 10,000 permutations. Tests were performed among clusters of sites identified using clustering techniques (see below), and among sites falling within states or provinces that were WNS-positive (KY, OH, TN, WV, ON-2, PA, MD, NY, NJ-Mo, NJ-Sa, QB; Table 1) or WNS-negative (all other sites) as of 2012, although it should be noted that tissues were collected prior to the emergence of WNS at these localities in all cases.
We applied three different approaches to determine the most likely number of distinct genetic clusters independent of original sampling locations, as different clustering algorithms can produce different solutions and concordance among multiple techniques is suggestive of the presence of a strong genetic signal [66]. First, we utilized the model-based Bayesian clustering approach in STRUCTURE v.2.3.3 [67,68] with population membership as a prior [69]. To determine the optimal number of clusters (K), we ran 10 runs per K, for K = 1-10, with a 100,000 MCMC iteration burn-in followed by 400,000 iterations using the admixture model with correlated allele frequencies. The most likely number of clusters was determined using the Evanno et al. [70] method implemented in the program STRUCTURE HARVESTER [71]. The Evanno et al. [70] method is not informative for the highest and lowest K; therefore, if the highest log likelihood value was observed for K = 1 or 10 across all replicates, we accepted that as the K with the highest probability. For the best value of K we used CLUMPP [72] to harmonize individual assignments to clusters.
Second, we applied the K means clustering approach in GENODIVE v.2.0 as outlined by Meirmans [65]. This approach is based on an AMOVA framework and uses a simulated annealing algorithm to minimize the among-populations/within-groups sum of squares through maximization of F CT , the variance among clusters relative to the total variance. We determined the most likely number of clusters using the Pseudo-F summary statistic, which performs better than the alternative Bayesian Information Criterion (BIC) when migration rates are high and mating is random [65].
The third approach was that of Duchesne and Turgeon [73] implemented in the software FLOCK. Samples are randomly partitioned into K clusters (2), allele frequencies are estimated for each of the K clusters, and each genotype is then reallocated to the cluster with the highest likelihood score. Repeated reallocation based on likelihood scores (20 iterations per run) resulted in genetically homogeneous clusters within a run. Fifty runs were carried out for each K, and at the end of each run the software calculated the log likelihood difference (LLOD) score for each genotype (the difference between the log likelihood of the most likely cluster for the genotype and that of its second most likely cluster) and the mean LLOD over all genotypes. Strong consistency among runs (resulting in 'plateaus' of identical mean LLOD scores) is used to indicate the most likely number of clusters [73]. Based on this analysis, we re-reran the iterative reallocation procedure for the most likely K and plotted the mean LLOD score against geographically ordered sites as a means of identifying admixture levels between genetic clusters.
For all three methods of genetic cluster identification, we tested whether cluster assignment was valid by quantifying the number of individuals within each sample that were allocated to each cluster, and then building an r × c contingency table where r is the number of genetic clusters and c is the number of sample sites. We then tested for random allocation to the genetic clusters across empirical samples using a likelihood-ratio test with Williams' correction with the null hypothesis that cluster assignments were random across sampled sites. A rejection of the null hypothesis indicated that that cluster composition was unlikely to be random across the samples, and that cluster assignments were therefore valid [74]. In addition, given that most clustering techniques assume that genotypic proportions within each cluster are in HWE and at linkage equilibrium, we tested identified clusters for compliance with these assumptions as suggested by Guillot et al. [66]. To test whether cluster assignment was independent of subspecific mtDNA clade membership we could not simply test for an association, as cluster assignment was confounded by spatial variation in the distribution of mtDNA clades. Therefore, we compiled cluster and clade membership for individuals in each of the four sites that contained members of more than one mtDNA clade (see Results, Fig 1a), and performed a likelihood-ratio test to determine whether cluster assignments were independent of subspecific clade membership within heterogeneous populations.
The level of genetic differentiation among pre-defined sites and an alternative grouping based on subspecific clade membership, where individuals were classified as belonging to the M. l. alascensis, M. l. carissima, or M. l. lucifugus clades based on the mtDNA phylogenetic analysis, was determined by calculating pairwise distance measures, including F ST [75], and two measures independent of the amount of within-population diversity: Jost's D [76], and G" ST [77]. Differences in the magnitude of pairwise distance measures among groups of sites was tested using 10,000 permutations in GENODIVE (G" ST and Jost's D) and FSTAT (F ST ). As with mtDNA, we performed an AMOVA on microsatellite genotypes using ARLEQUIN, and subsequently performed nested AMOVAs by grouping sites with low F ST values to identify the grouping that maximized among-group variation and minimized among-site/within-group variation.

Isolation by distance
There is considerable evidence to suggest that, regardless of the algorithm employed, clustering methods are confounded by the presence of isolation by distance (IBD), such that consistent clinal genetic variation may be misinterpreted by clustering algorithms as the presence of distinct clusters even though there is no true barrier to gene flow [66,[78][79][80]. We therefore tested for IBD in mitochondrial and nuclear DNA both globally (including all sampled locations) and within identified clusters (for microsatellite data only). We conducted a Mantel test comparing standardized genetic distance [F ST /(1-F ST )] and the natural log of geographic distance [81] using the IBD Web Service [82]. To calculate between-site geographic distances, polylines were constructed from X,Y coordinates in ArcGIS 10.1. The geodesic distance of these polylines was calculated using the "Shape.length@meters" command. For the microsatellites, we followed the recommendations of Guillot et al. [66]: we plotted genetic distance (F ST ) against geographic distance while differentiating between data points for site pairs that belonged to the same genetic cluster and data points for site pairs belonging to different clusters. If clusters are real, then for any given geographic distance, genetic distance between site pairs in different clusters should consistently be greater than distance between site pairs falling within the same cluster (e.g. [83,84]). In addition, as suggested by Meirmans [80], we performed a partial Mantel test to investigate the association between the matrix of genetic distances (F ST and G" ST ) and a matrix of cluster membership for the microsatellite data, with the matrix of geographical distances as a covariate. If there is not a relationship between genetic distance and cluster membership after controlling for geographic distance, then cluster membership is likely confounded by IBD. Given two main clusters (see Results), we coded cluster membership in two ways: a) both members of a sampling site pair belonging to the same cluster (1) or not (2); or b) both members of a sampling site pair belonging to a western cluster (1), both belonging to an eastern cluster (2), or sample sites belong in different clusters (3). Partial Mantel tests were carried out in PASSaGE software and significance was assessed with 10,000 permutations [85].

Genetic diversity
We observed 148 unique haplotypes characterized by 104 segregating sites among the 617 individuals sequenced. The number of haplotypes per site ranged from 2-13 (mean: 7.8), and the number of haplotypes unique to a site ranged from 1-11 (mean: 4.5; Table 1). There was no significant difference in nucleotide diversity (π) between sites east versus west of the Great Plains-Rocky Mountains transition (Mann-Whitney U-Test; Mean East: 0.00, West: 0.01, P = 0.829), or between sites in states that were positive or negative for WNS as of the 2012-2013 winter season (Mann-Whitney U-Test; WNS-Neg: 0.01, WNS-Pos: 0.00, P = 0.132). However, sites west of the Great Plains-Rocky Mountains transition had significantly lower haplotype diversity than those east of the boundary (Mann-Whitney U-Test; East: 0.80, West: 0.56, P = 0.002), and WNS-free sites as of 2012 also had significantly lower haplotype diversity than WNSaffected sites (Mann-Whitney U-Test; WNS-Neg: 0.65, WNS-Pos: 0.81, P = 0.015), although this latter result is likely confounded by the high proportion of western sites in the WNS-Neg group.
Although we originally typed 11 microsatellite loci, two loci (E24 and COTO_G02_H10) had high null allele frequencies and were dropped from further analyses. The remaining nine loci all met HWE expectations and were unlinked. Mean observed and expected heterozygosities were high (0.876 and 0.897, respectively), as was the mean number of alleles per locus (13.9) and allelic richness (11.9), although private allelic richness was low (0.23; Table 1; see S4  Table for diversity statistics for each locus). Comparing sample-level measures of genetic diversity among identified clusters and among samples in WNS-positive and WNS-negative states or provinces revealed no significant differences in observed or expected heterozygosity, allelic richness, or F IS (P > 0.05 in all cases based on permutation tests).

Differentiation of M. lucifugus subspecies
Phylogenetic analysis of the mtDNA COI sequences confirmed the existence of three clades within little brown bats roughly corresponding to previously defined subspecies (Fig 2). Following Carstens and Dewey [32], we refer to these clades as M. l. lucifugus, M. l. carissima, and M. l. alascensis (hereafter lucifugus, carissima, and alascensis, respectively). The carissima clade had three other species (M. evotis, M. keenii, M. thysanodes) nested within it, as previously described [32]. We have focused on M. lucifugus sensu stricto here. Addressing the taxonomic relationships among M. lucifugus, M. evotis, M. keenii, and M. thysanodes is beyond the scope of this paper; therefore we ignored the presence of these additional taxa in the rest of our analyses. For the most part, clades were geographically restricted and followed previously defined geographic ranges (Fig 1a). Almost all of the sampled sites east of the Great Plains-Rocky Mountains transition (with the exception of the Alberta-South site) were composed entirely of individuals classified as lucifugus. The carissima clade was restricted to mountainous but not coastal regions of the west and a portion of the western plains (in southern Alberta). The alascensis clade was found along the Pacific Coast and coastal mountain ranges (Fig 1a). However, individuals in the lucifugus clade were also sampled at locations west of the Great Plains-Rocky Mountains transition, including the British Columbia-South (BC-S), Wyoming (WY), and Washington (WA) samples. The Alberta-South (AB-S) sample east of the Great Plains-Rocky Mountains transition also contained both lucifugus and carissima haplotypes (see also [33]). No alascensis haplotypes were found at the BC-S site, although range maps place this site in that subspecies' range; furthermore, the Alaska (AK) sample clustered with the alascensis clade, even though it fell within the described range of the lucifugus clade.

Spatial patterns of population genetic structure
Mitochondrial DNA. AMOVA analysis considering all samples as a single group revealed high levels of differentiation (F ST = 0.721). We iteratively grouped sites with low pairwise F ST values to determine the best arrangements of sites that maximized 'among-group' and minimized 'among-site/within-group' variation in the AMOVA framework. Most samples were highly divergent from all others (76% of pairwise comparisons had F ST > 0.2, and 62% were > 0.5; S5 Table and Fig 1b), but we identified two groups of sites (one in the central United States and Canada east of the Rocky Mountains, and one in eastern North America) within which divergence was low (F ST = -0.019-0.130; Fig 1b). After grouping these sites together in an AMOVA analysis, among-group variation (F CT ) accounted for 73.5% of variation in haplotype frequencies, and among-site/within-group variation accounted for 1.4%.
Microsatellite DNA. All clustering methods employed [Bayesian clustering (STRUC-TURE), repeated reallocation (FLOCK), and K means clustering (GENODIVE)] identified K = 2 as the most likely number of genetic clusters, roughly corresponding to clusters east versus west of the Great Plains-Rocky Mountains transition, and not corresponding to subspecies affiliation (Figs 3 and 4). However, four sites in the transition, namely British Columbia-North (BC-N), Alberta-North (AB-N), AB-S, and WY had intermediate Q (STRUCTURE) or LLOD (FLOCK) values indicative of admixture, and there was a clear gradation in this region between genetic clusters (Figs 3 and 4). On average, the BC-N and AB-N sites had higher proportional membership with the eastern cluster, and AB-S and WY sites with the western cluster, and therefore they were grouped accordingly in all subsequent analyses. To test the validity of the most likely number of clusters identified using each of the three algorithms, we performed contingency table analyses with the null hypothesis that cluster assignments were random with respect to sampling sites. In all cases we rejected the null hypothesis (STRUCTURE: G = 521.5; FLOCK: G = 422.6; GENODIVE: G = 77.1; df = 20 and P < 0.001 in all cases) and concluded that clusters were valid. However, all clusters identified by the three algorithms failed to meet HWE expectations (P < 0.05 in all cases). Across the four sites (WA, BC-S, AB-S, and WY) containing both carissima and lucifugus mtDNA haplotypes (setting aside the 3 alacensis-type individuals in WA), cluster membership was independent of subspecies affiliation for the clusters defined by FLOCK (Χ 2 = 1.65, P = 0.1990) and GENODIVE (Χ 2 = 0.22, P = 0.6390), but was not for STRUCTURE (Χ 2 = 5.65, P = 0.0175).
AMOVA analysis of microsatellite genotypes indicated weak but significant population structure (global F ST = 0.0161, P < 0.001; proportion of variation within sites = 0.984). The grouping of sites that maximized among-group variation and minimized among-site/withingroup variation included a group containing AK only, a western group of samples (California-Siskiyou (CA-Si), WA, BC-S, Idaho (ID), and WY), and an eastern group containing all other samples (variation among groups = 2.72%, P < 0.001; variation among sites within groups = 0.41%, P < 0.001). Generally, F ST values between Alaska and all other sites were high and significant (0.049-0.089; S6 Table). F ST values between sites in the western and eastern groups ranged from 0.004-0.045, values between sites within the western group ranged from 0.002-0.018, and values between sites within the eastern group ranged from -0.005-0.013 (S6 Table). An AMOVA grouping individuals based on their subspecific membership similarly indicated weak but significant population structure (global F ST = 0.022, P < 0.001, proportion of variation within subspecies = 0.978). The maximum pairwise F ST was between alascensis and lucifugus (0.0270, P < 0.001), the minimum was between alascensis and carissima (0.0141, P < 0.001), and differentiation between carissima and lucifugus was intermediate (0.0200, P < 0.001).
Genetic distance measures were higher among sites within the western group than among sites within the eastern group, and permutation tests approached significance at the P < 0. Isolation by distance. To test for isolation by distance we performed Mantel tests on the logarithm of geographic distance and standardized genetic distance [F ST /(1-F ST )]. There were clear signals of IBD for both mitochondrial (r = 0.346, P < 0.0001; Fig 5a) and microsatellite DNA (r = 0.537, P < 0.0001 ; Fig 5b) across the range of little brown bats. Within identified clusters, there was no signal of IBD in the east (r = -0.307, P = 0.9925 ; Fig 5c), but there was within the west (r = 0.913, P = 0.0411; Fig 5d). To test if IBD in the west was disproportionately driven by the Alaska population, we re-ran the analysis with that site removed, and the signal of IDB remained (r = 0.880, P = 0.0069; S1 Fig).
To assess the validity of clusters given the pattern of isolation by distance, we plotted geographic and genetic distance (F ST ) based on microsatellites according to cluster membership (points identified separately for comparisons within the same cluster vs. in different clusters; Fig 6). There was no clear separation between site pairs in the same vs. different clusters, and low genetic distances for a given geographic distance were observed regularly for site pairs in different clusters. However, partial Mantel tests to examine the association between the matrix Range-Wide Population Genetic Analysis of Little Brown Bats of genetic distances and a matrix of cluster membership for the microsatellite data, with the matrix of geographical distances as a covariate, were significant (P < 0.05 in all cases) regardless of coding scheme (see Methods) or genetic distance measure (F ST or G" ST ) used.

Discussion
The emergence and spread of WNS has decimated bat populations in affected areas and raised the specter of extinction for Myotis lucifugus, M. septentrionalis, and other highly-affected bat species. Therefore, understanding population connectivity and possible barriers to disease transmission is vital to the ongoing management and conservation of affected species and the development of mitigation strategies to limit disease spread and associated mortality. Here, using a combined dataset of mtDNA sequences and nuclear microsatellite genotypes for little brown bats from across their range, we demonstrate considerable spatial variation in patterns of female dispersal and significant genetic variation between sites in eastern versus western portions of the range of little brown bats. Whether the observed variation is representative of discrete genetic clusters rather than isolation by distance is debatable (see below), but overall, it is clear that levels of nuclear genetic differentiation are low, and there is no evidence for any major barriers to nuclear gene flow across the range of little brown bats. However, some key spatial patterns emerge from our analyses, namely (1) patterns of mtDNA differentiation are highly variable, with high F ST values between most sample pairs (including between all western samples, between western and eastern samples, and between some eastern samples), while low mitochondrial differentiation was observed within two groups of samples found in central and eastern regions of North America (shown in AMOVA and pairwise F ST analyses; Fig 1b); (2) the site from Alaska is highly differentiated from all others in our study (shown in AMOVA and pairwise F ST analyses); and (3) western sites are characterized by significant isolation by distance based on microsatellites, while those in the east are not. These data raise the possibility that the current patterns of spread of Pd observed in eastern North America may not apply to the entire range of the little brown bat, and that there may be broad-scale spatial variation in the risk of WNS transmission and occurrence if the disease continues to spread west.
The presence of isolation by distance is a major confounding factor when examining levels of genetic structure among populations of widespread species such as little brown bats, because clinal variation may be interpreted as the presence of discrete clusters even in the absence of barriers to gene flow [66,78,80]. All clustering methods we employed on our microsatellite data identified the presence of two genetic clusters, roughly dividing sites east vs. west of the Great Plains-Rocky Mountains transition. However, we also observed a strong pattern of isolation by distance, indicating that these observed clusters may be an artifact of dispersal limitation and clinal variation across the very large range of this species. Additional analyses to test the validity of clusters provided mixed results. On the one hand, genetic distances between eastern and western sites were relatively low (the best grouping in the AMOVA explained only 2.7% of the variation in microsatellite allele frequencies), a Mantel test showed no clear separation between site pairs in the same versus different clusters (Fig 6), and identified clusters failed to meet HW expectations, suggesting that clusters did not represent genetically panmictic populations (cf. [66]). On the other hand, partial Mantel tests controlling for geographic distance revealed a significant correlation between cluster membership and genetic distance, suggesting that western and eastern clusters were differentiated despite the signal of isolation by distance (cf. [80]).
What is clear from these data is that there is significant genetic variation among samples from east to west across the range of little brown bats (although whether they represent discrete genetic clusters rather than isolation by distance is debatable), but levels of nucDNA genetic differentiation were relatively low, and there was no unambiguous evidence for any major barriers to gene flow that might severely restrict the spread of WNS. Gene flow in temperate bats is mediated through the permanent dispersal of individuals and the exchange of genes at mating congregations during swarming and hibernation. The lack of isolation by distance and low levels of nucDNA differentiation among sites in eastern North America is concordant with the continuous spread of WNS from its origin in New York, and indicates that gene flow via mating has occurred over wide geographic areas. Furthermore, the disease has passed, or is currently passing, through regions in which there are low levels of mtDNA differentiation among  Fig 1b). Most temperate bats are characterized by relatively high levels of female philopatry and male-biased dispersal, resulting in significant matrilineal genetic structuring of populations (e.g. [35,[86][87][88]). However, our data suggest that the exchange of females among populations across large portions of the range of little brown bats is a non-trivial source of gene flow that may be contributing to the spread of WNS, and are consistent with similar inferences of female dispersal among populations over smaller spatial scales in little brown bats [6,33,35,42] and other bat species [38,89] based on mtDNA. In addition, these data are consistent with banding data showing extensive movements by individual little brown bats of both sexes over hundreds of kilometers between summer roosts, swarming sites, and hibernacula within and between years in central Canada [27,42]. Banding studies of little brown bats in other parts of their range would help to resolve whether the observed differences are better explained by historical demography or current gene flow. The high levels of observed mtDNA differentiation in other portions of the range, particularly in the west, suggest important spatial variation in female dispersal patterns, and highlight the need to consider permanent movements of both males and females and incorporate regional variation in dispersal rates and distances in models of WNS transmission dynamics.
Given the lack of major physiographic barriers east of the Rocky Mountains and the high levels of gene flow we inferred, it is likely that WNS will continue to expand its range across eastern North America. Current models of disease spread indicate that WNS exhibits characteristics of an expanding epizootic wherein relatively distant sites have lower infection risk, but over time infection rates increase and the effect of distance diminishes as the disease 'fills in' behind the expansion front [19,20]. However, these models are based on data from eastern portions of the range where gene flow is high. Even within this region the rate of spread may be restricted by consistent spatial variation in above-ground or cave microclimates that may limit the survival and growth of Pd (cf. [20,21]), or by spatial variation in topography or land use that limits movements and dispersal by bats [6]. In Pennsylvania, for example, topographical features such as the Appalachian high plateau and the Allegheny Front escarpment may have influenced seasonal migration patterns of female bats, thereby limiting matrilineal gene flow and disease transmission rates among populations [6]. Indeed, two genetically distinct populations of wintering colonies were observed on either side of the Allegheny Front [6]; hibernating colonies of little brown bats located on the western Appalachian high plateau were infected with WNS 1-2 years later than colonies in the central mountainous and eastern lowland regions of the state [16]. Thus, topographic or climatic variation may slow the spread of WNS through some areas by limiting population connectivity of the host or the survival and growth of Pd, and may explain some of the observed spatial variation to date in the rate and direction of WNS spread through eastern North America.
Our genetic data indicating lower levels of population connectivity in the west suggest that if WNS reaches western populations, the rate of disease spread may decline. The high mtDNA F ST values among western populations (S5 Table) indicate that female movements are highly restricted relative to eastern populations. Overall levels of nucDNA gene flow among western sites were reduced relative to the east, and western sites were characterized by isolation by distance based on microsatellites while eastern sites were not. These results may, in part, be related to the greater topographical and ecological heterogeneity in the west, which includes multiple mountain ranges, plateaus, basins, and coastal lowlands, and which has been implicated in recurrent phylogeographic patterns in a wide variety of other taxa [90]. Hibernation behavior is poorly characterized for western North American little brown bat populations. All known large hibernating populations (>10,000 individuals) are described from eastern North America, and identified hibernacula in the mountainous west typically have lower census sizes than many hibernacula in the east. The high physiographic variation in the mountainous west may limit population connectivity and the scale of bat movements, and the high density of mines and caves in many regions in the west may result in smaller and more diffuse hibernating colonies relative to eastern North America. Comparative data on connectivity between summer and winter sites (as in [27]) are urgently required to quantify spatial and temporal patterns of movement of little brown bats in the western portion of their range and to predict potential rates of WNS transmission. Further, the most distant population we sampled (in Alaska) was by far the most divergent from all other populations, and we require much more dense sampling in the western portion of the range of little brown bats to determine if any other populations are equally or more isolated and hence may have reduced contact rates with other regional populations.
The spatial variation in population connectivity we observed was largely independent of subspecific affiliation. Phylogenetic analysis revealed the presence of three divergent lineages based on mtDNA (corresponding to previously defined subspecies M. l. alascensis, M. l. carissima and M. l. lucifugus; as in [32]), with the notable finding of multiple lineages at the same sampling locations in southern British Columbia, southern Alberta (as in [33]), and Wyoming. Although Carstens and Dewey [32] provided some support from mtDNA and nuclear introns for discrete evolutionary lineages within M. lucifugus, we found that cluster membership based on microsatellites was independent of subspecific affiliation, and we estimated low levels of nucDNA differentiation among subspecies (F ST = 0.022 in AMOVA analysis). Our observed discrepancy between mtDNA and nucDNA signals may be due in part to homoplasy, particularly at rapidly-evolving microsatellite loci. Alternatively, incomplete lineage sorting might be producing discordant patterns among loci, particularly for a species with a relatively large effective population size (N e 400,000 based on Carstens and Dewey's [32] estimate of θ for M. lucifugus) and relatively recent divergence (divergence from western Myotis sp. approximately 1-1.5 Mya [32]). A third alternative is that our finding of lower differentiation at nucDNA compared to mtDNA is that patterns at mtDNA may reflect genetic differentiation that evolved and were reinforced as populations used distinct glacial refuges (as in [34]), while patterns at nucDNA reflect secondary contact particularly mediated via male gene flow.
Our data show extensive spatial variation in levels of connectivity among little brown bat populations and provide valuable information for understanding past and future patterns of WNS spread. However, the use of genetic methods to infer patterns of transmission assumes that patterns of gene flow are indicative of the movement of infectious individuals, and we must recognize that the risk of disease transmission may be higher than genetic data may indicate because there may be more contacts among infected and susceptible individuals, including among members of multiple species, than just those that lead to gene flow. Urgent research is required to determine how and when individual bats may be exposed to Pd spores, and how contacts of varying durations and seasonal timings influence the risk of WNS transmission. Ultimately we need to learn whether brief contacts during mating can result in transfer of spores leading to infection or whether permanent dispersals are driving transmission. The usefulness of our genetic data on little brown bats also rests on the assumption that intraspecific transmission dynamics outweigh the impact of cross-species transmission, given that multiple sympatric bat species are affected by WNS. This assumption may be justified, as post-WNS population declines of affected bat species are not influenced by population sizes of other affected, cohabiting bat species [21], but research is required to assess how often cross-species transmission may take place and how the rate of introduction of infective propagules to environmental reservoirs is influenced by multiple species cohabiting the same hibernaculum.
In conclusion, this study identified high levels of genetic variation among populations of little brown bats across their range, and mitochondrial DNA sequences revealed considerable spatial variation in patterns of female dispersal. Overall levels of nuclear genetic differentiation among M. lucifugus populations are low, and we did not identify any major barriers to gene flow across their range. However, levels of genetic differentiation at both mtDNA and microsatellites are significantly higher among populations to the west of the Great Plains-Rocky Mountains transition, suggesting that the current pattern of spread of WNS and risk of transmission of Pd observed in eastern North America may not apply to the entire range of the little brown bat.