Genetic Structure of Little Brown Bats (Myotis lucifugus) Corresponds with Spread of White-Nose Syndrome among Hibernacula

Until recently, the little brown bat (Myotis lucifugus) was one of the most common bat species in North America. However, this species currently faces a significant threat from the emerging fungal disease white-nose syndrome (WNS). The aims of this study were to examine the population genetic structure of M. lucifugus hibernating colonies in Pennsylvania (PA) and West Virginia (WV), and to determine whether that population structure may have influenced the pattern of spread of WNS. Samples were obtained from 198 individuals from both uninfected and recently infected colonies located at the crest of the disease front. Both mitochondrial (636bp of cytochrome oxidase I) and nuclear (8 microsatellites) loci were examined. Although no substructure was evident from nuclear DNA, female-mediated gene flow was restricted between hibernacula in western PA and the remaining colonies in eastern and central PA and WV. This mitochondrial genetic structure mirrors topographic variation across the region: 3 hibernating colonies located on the western Appalachian plateau were significantly differentiated from colonies located in the central mountainous and eastern lowland regions, suggesting reduced gene flow between these clusters of colonies. Consistent with the hypothesis that WNS is transmitted primarily through bat-to-bat contact, these same 3 hibernating colonies in westernmost PA remained WNS-free for 1–2 years after the disease had swept through the rest of the state, suggesting that female migration patterns may influence the spread of WNS across the landscape.

cytochrome oxidase I, dispersal, microsatellites, sex-biased gene flow, white-nose syndrome

Host dispersal and migration are fundamental parameters in disease ecology. Large-scale movements of host individuals may contribute to the spread of wildlife diseases by introducing a pathogen to naïve populations of the same species or by exposing different species to the pathogen, thereby facilitating host-jumping (reviewed by Altizer et al. 2011). Population and landscape genetic approaches have proven useful for understanding the broad-scale movements of both host and pathogen populations (Archie et al. 2009; Atterby et al. 2010; Biek and Real 2010; Kelly et al. 2010; Cullingham et al. 2011; Vollmer et al. 2011; Gray and Salemi 2012). Assuming that patterns of gene flow in the host influence that of the pathogen (although see Humphrey et al. [2010]), the identification of spatial, geographic, and behavioral factors that influence the dynamics of host dispersal might allow the prediction of likely patterns of disease movement into currently uninfected populations, allow the estimation of relative risk of infection for different populations and regions, and allow managers to strategize management decisions, both spatially and temporally (Biek and Real 2010).

With a range that extends across much of the United States and Canada (Fenton and Barclay 1980; Reid 2006), the little brown bat (Myotis lucifugus) was, until recently, one of the most common bat species in North America. Despite its historical abundance, M. lucifugus now faces significant threats to its continued existence in the eastern United States and Canada, the most severe being white-nose syndrome (WNS; Frick et al. 2010). WNS is an emerging disease of major conservation concern for hibernating bat populations of North America. To date, WNS has been documented in 7 species of bats (Reeder and Turner 2008; US Fish and Wildlife Service 2012a). The etiological agent is the psychrophilic fungus, Pseudogymnoascus destructans (Lorch et al. 2011; Minnis and Lindner 2013), which grows optimally between 5 and 10 °C (Blehert et al. 2009), overlapping with the optimal ambient hibernation temperature for little brown bats (T_a = 7.2±2.6 °C; Brack 2007). The fungus is transmitted primarily through direct contact between individuals (Lorch et al. 2011). WNS has spread rapidly since it was first detected in New York in 2006 (Hicks et al. 2007), particularly along the Appalachian Mountains. It has reached as far south as Alabama, as far west as Oklahoma, and has spread north into Ontario, Québec, and the maritime provinces of Canada (Blehert et al. 2009; US Fish and Wildlife Service 2012b). More than 5.5 million hibernating bats have died from WNS (US Fish and Wildlife Service 2012c), and if current mortality rates continue, there is a 99% probability that M. lucifugus will be extirpated from the northeast within the next 15 years (Frick et al. 2010). Thus, the current status of M. lucifugus as a nonpriority species will need to be reassessed.

In Pennsylvania (PA), winter aggregations of up to 90000 little brown bats have previously been reported (Turner and Butchkoski 2006). However, as in other parts of the northeastern United States, M. lucifugus colonies in PA have suffered high mortality from P. destructans infection, with populations in affected areas experiencing a decline greater than 90% in just a few years (Turner et al. 2011). WNS was first documented in PA in January 2009 (Table 1) and spread predominantly in a southwesterly direction along the Appalachian Mountains and associated ridge and valley regions into Maryland, Virginia, and West Virginia (WV). Interestingly, several colonies in western PA remained WNS-free for 1–2 years after the first infections were confirmed in the eastern and southern parts of the state (Table 1, US Fish and Wildlife Service 2012b). Because WNS is transmitted mostly through direct contact between individual bats (Lorch et al. 2011), movement of bats across the landscape is a determining factor in the spread of the pathogen. By influencing bat dispersal and seasonal migration, it is thus possible that topographic variation across the state may temporarily protect some populations from exposure to the disease.

Table 1.

Location and WNS status of hibernating colonies at the time of sample collection

Hibernaculum	ID	NS/NG	County	WNS status (at collection)	Date WNS infection confirmed
PA (n = 179)
Durham Mine	DM	19/19	Bucks	Positive	Dec 2009
Dunmore Slope	DUN	21/21	Lackawanna	Positive	Feb 2009
Glen Lyon Mine	GL	20/0	Luzerne	Positive	Feb 2009
Lake Mine	LM	6/0	Tioga	Positive	Feb 2010
Shindle Iron Mine	SHIM	20/20	Mifflin	Positive	Jan 2009
US Steel Mine	US	20/20	Armstrong	Negative	Feb 2012
Barton Cave	BC	14/0	Fayette	Negative	Jan 2011
Layton Fire Clay Mine	LAY	20/20	Fayette	Negative	Apr 2010
CSM Mine	CSM	19/20	Lawrence	Negative	Apr 2010
Snowtop Mine	ST	20/20	Lawrence	Negative	Dec 2010
WV (n = 19)
Snedegar’s Cave	SNED	19/19	Pocahantas	Positive	Feb 2010

N_S and N_G represent the number of individuals sequenced and genotyped, respectively.

PA has more than a dozen physiographic regions (Figure 1), ranging from a lowland (sea level) belt in the southeast along the Delaware River, through the hills and lowlands of the Piedmont Plateau and Great Valley, to the ridges and valleys of the Appalachian Mountains in the eastern and central regions, and the Appalachian plateau (AP) in the west and north (Sevon 2000). The Allegheny front forms the escarpment along the eastern edge of the plateau (Figure 1). The initial spread of WNS from its original infection site in New York has been primarily to the north, south, and east, rather than to the west (Swezey and Garrity 2011; US Fish and Wildlife Service 2012b), such that all the PA sites infected with WNS during or prior to 2009 are located to the east of the Allegheny front. We, therefore, hypothesized that the Allegheny front and associated topographic features may act as barriers to bat movement and thus to gene flow across the Allegheny Mountains. This in turn may influence the rate and direction of spread of WNS across PA and into neighboring states.

Figure 1.

Physical map of the northeastern United States showing locations of hibernacula included in this study and the major topographic features of the region (Sevon 2000; Woods et al. 2006). Open circles and filled squares indicate colonies that were WNS-neg or WNS-pos, respectively, at the time of collection. Maps, services, and data are available from US Geological Service, National Geospatial Program.

Little brown bat movement involves seasonal migrations between summer maternity sites, autumn swarming sites, and winter hibernacula. Maternally inherited genetic markers are strongly structured by colony in western portions of its range during the summer, indicating strong female philopatry to maternity colonies (Lausen et al. 2008), which are populated by reproductive females and their young. However, this behavior may be somewhat variable across the species’ range, as weaker levels of female philopatry have recently been reported in M. lucifugus populations in Minnesota (Dixon 2011). In late summer and autumn, bats of all ages and both sexes migrate to swarming sites located at cave or mine entrances. Here, and later at hibernacula, bats from a larger geographic area congregate for mating (McCracken and Wilkinson 2000). This highly promiscuous and indiscriminate mating behavior likely leads to extensive genetic mixing among populations that are spatially and behaviorally isolated at other times of the year (Lausen et al. 2008). As with maternity colonies, little brown bats appear to display high roost fidelity to both swarming sites and hibernacula (Davis and Hitchcock 1965; Humphrey and Cope 1976).

To date, few population genetic studies have been conducted on M. lucifugus (Lausen et al. 2008; Dixon 2011), and none has examined its population structure in the northeastern United States where WNS is most prevalent. Information on the existing population structure of little brown bats in this region is urgently needed in anticipation of WNS recovery because if genetically distinct subpopulations exist, they may need to be managed independently during future WNS recovery programs. Furthermore, if there is differential survivorship at those locations, these should become the primary targets for conservation and management in WNS recovery planning. However, without a detailed understanding of patterns of gene flow and how colonies are linked via seasonal migrations of the bats, it will be difficult to identify and prioritize those subpopulations.

This study utilized highly variable molecular markers to examine the fine-scale population genetic structure of M. lucifugus hibernating colonies in PA and WV, and to correlate these genetic inferences of bat movement with the documented spread of WNS through this region. To compare patterns of migration and/or dispersal between the sexes, we examined both biparentally inherited nuclear markers (microsatellites), as well maternally inherited mitochondrial DNA (mtDNA). The latter marker tracks female movement patterns, whereas the nuclear microsatellites provide information about both sexes; any differences between the markers can thus be interpreted as the result of sex-biased differences in gene flow.

The aims of this study were 1) to determine over what spatial scales and in which directions individuals move among hibernacula; 2) to determine whether local topography influences bat movement, and 3) if population structure does exist, to determine whether it corresponds with the spread of WNS through this region. Specifically, we hypothesized that the Allegheny front and associated topographic features influence dispersal and migration of little brown bats, thereby restricting gene flow among colonies on either side of the Appalachian Mountains. This would not only result in genetically distinct subpopulations but also may have influenced the pattern of spread of WNS through the region.

Materials and methods

Sample Collection

Samples were collected between 2008 and 2010 from 198 M. lucifugus. This included 179 individuals from 10 hibernacula located throughout PA (Table 1, Figure 1), as well as 19 bats from a single colony (Snedegar’s Cave [SNED], Table 1, Figure 1) in neighboring Pocahontas County, WV. The latter was included to extend our sampling area further south, tracking the path of WNS infection along the Appalachian Mountains. To test our hypothesis that local topography influences bat movement, and in turn the spread of WNS, we sampled colonies on either side of the Allegheny front. In anticipation of future infection with WNS, we specifically included sites that were WNS-negative (WNS-neg) at the time of collection, as well as localities already affected by WNS (WNS-pos, Figure 1, Table 1). Tissue samples included wing biopsies (Worthington Wilmer and Barratt 1996) from living bats at WNS-neg sites, and wing or muscle tissue from dead bats at WNS-pos sites. All samples were stored in 500 µL of 90% ethanol or lysis buffer (0.5% sodium dodecyl sulfate, 400mM NaCl, 100mM ethylenediaminetetraacetic acid [EDTA], 50mM Tris–HCl pH 7.5) until DNA extraction. DNA was extracted using a DNeasy kit (Qiagen, Valencia, CA), according to the manufacturer’s instructions, and then stored in 1× Tris–EDTA.

Data Collection

mtDNA Sequencing

Approximately 680bp of the cytochrome oxidase I (COI) gene was amplified from each individual (Table 1) by means of the polymerase chain reaction (PCR), using primers HCO2198 and LCO1490 (Hebert et al. 2003) or primers VF1 and VR1 (Ivanova et al. 2006), and PCR conditions as described in the respective publications. 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.

Microsatellite Genotyping

We genotyped 159 individuals from 8 of the 11 colony locations (Table 1). Preliminary mtDNA indicated no significant substructure within eastern PA. Therefore, to minimize microsatellite genotyping costs, we excluded colonies that would be unlikely to add significant information to the analyses, namely those with small sample size (e.g., Lake Mine [LM], n = 6) and those that were located very close to other colonies in the eastern region (e.g., Glen Lyon Mine [GL] is only 68 km from Dunmore Slope [DUN]). Individuals were genotyped at 8 highly variable microsatellite loci, using primers previously developed for other vespertilionid bats: IBat CA5, CA43, CA47, and M23 from Oyler-McCance and Fike (2011); MS3D02 and MS3F05 from Trujillo and Amelon (2009); E24 from Castella and Ruedi (2000); and Cora_F11_C04 from Piaggio et al. (2009). Amplifications were carried out in 3 multiplexes (two 2-locus and one 3-locus) and 1 single-locus amplification, which were subsequently pooled into 2 different loads for fragment analysis. PCR and cycling conditions generally followed those outlined in Vonhof et al. (2002), with the modification of a 10 s (rather than 1 s) extension step at 72 °C and variable annealing temperatures and number of cycles depending on the multiplex. Information on multiplexes and associated primer concentrations, annealing temperatures, and cycling conditions is available from the authors. Amplified products were sent to the Vanderbilt University DNA Sequencing Facility for genotyping.

Statistical Analyses

COI sequences were manually edited and aligned using CodonCode Aligner (CodonCode Corp., Dedham, MA) and were trimmed to 636bp. Estimates of genetic diversity were obtained using DnaSP version 5.10.01 (Librado and Rozas 2009), and pairwise genetic differences (Φ_ST) between each pair of colonies were calculated in ARLEQUIN version 3.5.1.2 (Excoffier and Lischer 2010) using a Tamura-Nei distance model and a transition:transversion ratio of 26.9:1, as estimated by jModelTest version 0.1.1 (Posada 2008). To examine hypotheses of population substructure, analyses of molecular variance (Amova; Excoffier et al. 1992) were performed. Based on the pairwise Φ_ST results, 2 different Amovas were conducted, using the same model parameters described previously: 1) WNS status at the time of collection (i.e., WNS-neg colonies Snowtop Mine [ST], CSM Mine [CSM], US Steel Mine [US], Layton Fire Clay Mine [LAY], Barton Cave [BC], compared with the rest) and 2) the effect of topography (i.e., colonies on the AP [ST, CSM, US] compared with the rest). The distance (in kilometers) between each pair of hibernacula was estimated from their geographic coordinates using Google Earth version 6.1.0 (Google Inc., Mountain View, CA). A Mantel test (Rousset 1997) was conducted using the Isolation By Distance Web Service version 3.23 (Jensen et al. 2005), testing for a correlation between a matrix of log-transformed geographic distances and a matrix of genetic distances (Slatkin’s linearized F_ST) between the colonies. Significance was determined using 1000 permutations. In fulfillment of data archiving guidelines (Baker 2013), primary data underlying these analyses have been deposited with Dryad.

We used the isolation-with-migration model implemented in the IMa2 version 8.27.12 software package (Hey and Nielsen 2004, 2007) to estimate demographic parameters from the mtDNA sequence data. As implemented here, this analysis estimated 6 parameters: the genetic variability of 2 daughter populations (θ₁ and θ₂); the genetic variability of the ancestor of those daughters (θ_A); directional gene flow between the 2 daughter populations (M₁ and M₂); and the time at which the daughter populations diverged (τ). These coalescent parameters were converted to natural parameters (N_e = effective population size; Nm = number of migrants per generation; T = divergence time) as follows: θ_X = 4N_xμ; Nm_x = (θ_XM_X)/4; and τ = Tμ. Six independent runs were conducted, with parameters drawn from uniform prior distributions set at θ ∈ U[0, 100], m ∈ U[0, 10], and τ ∈ U[0, 50]. Analyses were run for at least 10 million steps, each with 40 Metropolis-Hastings coupled chains and a burn-in of at least 6 million steps. Posterior estimates of scaled parameters (θ, M, and τ) from the coalescent were converted to natural parameters (N_e, m, and T) by assuming a mutation rate of 1.95% per million years (Hickerson et al. 2006) and a generation time of 2 years (Humphrey and Cope 1976). After verifying that all runs converged on similar posterior distributions, we used IMa2 to conduct parameter comparisons and to calculate the joint posterior densities for model parameters based on the 132589 resulting coalescent genealogies from all 6 runs.

Microsatellite fragments were analyzed and scored using GeneMarker software (SoftGenetics LLC, State College, PA). Single base-pair alleles caused by indels in flanking sequences were common, and we converted them to whole-repeat alleles by consistently rounding up or down within a locus. Deviations from Hardy–Weinberg equilibrium (HWE) for each locus were estimated, and loci were confirmed to be in linkage equilibrium using ARLEQUIN version 3.11 software (Excoffier et al. 2005). Null allele frequencies were estimated in CERVUS version 3.0 (Kalinowski et al. 2007). Several indices of nuclear genetic diversity were estimated, including number of alleles per locus, observed heterozygosity using CERVUS, and allelic richness and the inbreeding coefficient, F_IS, using FSTAT version 2.9.3 (Goudet 1995). F_IS is the proportion of the genetic variance in the subpopulation that is contained within an individual (e.g., high F_IS values imply high levels of inbreeding). Mean observed heterozygosity and allelic richness (across loci) were compared between WNS-pos and WNS-neg colonies using Mann–Whitney U-tests.

To determine the most likely number of distinct genetic clusters, we utilized 2 clustering approaches. First, we used the model-based Bayesian clustering approach in STRUCTURE version 2.3.3 software (Pritchard et al. 2000; Falush et al. 2003; Pritchard and Wen 2004). To determine the optimal number of clusters (K), we ran 1000000 MCMC iterations (200000 burn-in) with 10 runs per K, for K = 1–8 using the admixture model with correlated allele frequencies and determined the optimal number of clusters with ΔK, calculated using the criteria outlined in Evanno et al. (2005) in the program STRUCTURE HARVESTER (Earl and Vonholdt 2012). The Evanno et al. (2005) method is not informative for the highest and lowest K, therefore, if the highest log likelihood value was observed for K = 1 or 8 across all replicates, we accepted that value of K as having the highest probability. Second, we used the approach of Duchesne and Turgeon (2012) implemented in the software FLOCK. In this approach, samples are initially randomly partitioned into K clusters (K ≥ 2), allele frequencies are estimated for each of the K clusters, and each genotype is then reallocated to a cluster to maximize the 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) was used to indicate the most likely number of clusters (Duchesne and Turgeon 2012).

The level of genetic differentiation among sampled colonies was determined by calculating pairwise F_ST values (Weir and Cockerham 1984) and by testing for significance using Amovas with 1000 permutations in ARLEQUIN version 3.11 software (Excoffier et al. 2005). We tested all colonies as independent units in the Amovas, as well as the 2 groupings outlined previously for mtDNA. Bayesian clustering methods such as STRUCTURE may overestimate the number of independent genetic units if a pattern of isolation-by-distance (IBD) is present (Frantz et al. 2009). Therefore, we used the Mantel tests implemented in the program IBD Web Service (Jensen et al. 2005) to test whether log-transformed geographic distances were correlated with standardized genetic distance (F_ST/1 − F_ST) as suggested by Rousset (1997).

Results

Genetic Diversity

In total, 37 unique haplotypes were identified among the 198 mitochondrial COI sequences. There were 38 segregating sites, including 23 that were parsimony informative (34 transitions, 4 transversions, no indels). All mutations were synonymous. Hibernating colonies in PA had 36 different haplotypes (n = 179) and Snedegar’s in WV had 6 (n = 19). Mitochondrial genetic diversity in the colonies was moderately high (haplotype diversity: 0.655–0.9421; Table 2). There were no significant differences in haplotype-based diversity statistics between WNS-pos and WNS-neg locations, but we did detect significantly lower nucleotide-based diversity measures at sites that were WNS-pos at the time of sampling (average pairwise difference: Mann–Whitney U = 29, P = 0.009; nucleotide diversity: U = 29, P = 0.009).

Table 2.

Summary of mitochondrial genetic diversity statistics from each hibernating colony

	Number of haplotypes	Haplotype diversity	Mean # pairwise differences (k)	Nucleotide diversity (π)
All colonies (n = 198)	37	0.852	2.4510	0.0042
WNS-pos sites (n = 105)	29	0.8360	1.8020	0.0028
WNS-neg sites (n = 93)	11	0.9024	3.0591	0.0048
PA (n = 179):
DM	6	0.6550±0.1115	1.2217±0.8122	0.0019±0.0014
DUN	13	0.8905±0.0461	2.2657±1.2963	0.0036±0.0023
GL	11	0.8842±0.0535	1.7800±1.0743	0.0028±0.0019
LM	5	0.9333±0.1217	1.6093±1.1004	0.0025±0.0020
SHIM	10	0.8842±0.0479	2.1035±1.2242	0.0033±0.0022
US	12	0.9421±0.0295	3.7552±1.9763	0.0059±0.0035
BC	9	0.9011±0.0624	2.2136±1.2989	0.0035±0.0023
LAY	10	0.8368±0.0757	2.8301±1.5570	0.0044±0.0027
CSM	10	0.9006±0.0390	3.0689±1.6694	0.0048±0.0029
ST	14	0.9316±0.0347	3.4278±1.8283	0.0054±0.0032
WV (n = 19)
SNED	6	0.6023±0.1242	1.3519±0.874	0.0021±0.0015

There was no evidence of linkage among microsatellite loci after applying Bonferroni corrections. Three of the 8 loci exhibited some level of departure from HWE; however, these deviations were limited to a small number of colonies for each locus (1 colony for CORA F11_C04, 1 colony for E24, and 3 colonies for IBAT M23). Nuclear genetic diversity was high, with a mean of 25.6 alleles identified per colony, mean observed heterozygosity of 0.906, and mean allelic richness per colony of 13.77 (Table 3). Similarly, high genetic diversity values were observed on a per locus basis (Supplementary Data). Observed heterozygosity, allelic richness, and F_IS did not differ significantly among WNS-pos and WNS-neg locations (P > 0.05 in all comparisons).

Table 3.

Measures of nuclear (microsatellite) genetic diversity for each hibernating colony for 159 adult Myotis lucifugus captured in 8 colonies in PA and WV

Colony	NA	HO	AR	FIS
DM	13.6	0.906	13.63	–0.009
DUN	15.1	0.920	14.69	0.044
SHIM	15.4	0.920	15.07	0.002
US	13.6	0.900	13.39	0.020
LAY	13.9	0.914	13.63	0.049
CSM	14.5	0.898	14.15	0.075
ST	13.0	0.899	12.80	0.013
SNED	13.4	0.902	13.38	0.052
Overall	25.6	0.906	13.77	0.028

N_A, observed number of alleles; H_O, observed heterozygosity; AR, mean allelic richness per population; F_IS, estimated null allele frequency. Overall values are 8-locus means for number of alleles, heterozygosity, allelic richness, and null allele frequency per population.

Population Structure

Pairwise mtDNA comparisons among hibernating colonies (Table 4) indicated that colonies located north and east of the Allegheny front (Shindle Iron Mine [SHIM], Durham Mine [DM], DUN, LM, and GL), that is, in the central and eastern ridge and valley region of PA, were not significantly differentiated from one another (Φ_ST: −0.0089 to 0.0378; P > 0.05; Table 4). However, WNS-neg colonies located further west on the AP (CSM, US, and ST) were significantly differentiated from those in central and eastern PA (Φ_ST: 0.1480–0.639; P <0.0001–0.0400). The 3 western colonies (CSM, US, and ST) also differed significantly from SNED in WV (Φ_ST: 0.1864–0.2482; P < 0.0001). Two of the 5 PA colonies that were WNS-neg at the time of collection (BC and LAY) yielded mixed results. BC was moderately differentiated from the western WNS-neg colonies (CSM, US, and ST; Φ_ST: 0.0956–0.1427; P < 0.0313), but not from any of the WNS-pos colonies. LAY was weakly differentiated from ST, DUN, and SHIM (Φ_ST: 0.0477–0.0801; P < 0.0469), but not from any other colonies (Table 4). SNED in WV was weakly differentiated from SHIM in central PA (Φ_ST = 0.0488; P = 0.0469), but did not differ significantly from any other central and eastern PA colony (P > 0.0537).

Table 4.

Pairwise comparisons between hibernating colonies, based on mtDNA COI sequences

	DM	DUN	GL	LM	SHIM	US	BC	LAY	CSM	ST	SNED
DM	—	0.0635	0.1699	0.2315	0.1904	<0.0001	0.5068	0.0801	<0.0001	<0.0001	0.7568
DUN	0.0378	—	0.5850	0.9600	0.0908	<0.0001	0.3691	0.0469	<0.0001	<0.0001	0.0537
GL	0.0183	−0.0089	—	0.8936	0.2930	<0.0001	0.3613	0.1006	0.0020	<0.0001	0.1231
LM	0.0266	−0.0700	−0.0576	—	0.5889	0.0244	0.5195	0.2686	0.0400	0.0147	0.1748
SHIM	0.0188	0.0281	0.0058	−0.0211	—	<0.0001	0.1816	0.0440	<0.0001	<0.0001	0.0469
US	0.2452	0.2035	0.1966	0.1691	0.2081	—	0.0244	0.0547	0.3145	0.7998	<0.0001
BC	−0.0036	0.0037	0.0039	−0.0121	0.0230	0.1199	—	0.6660	0.0313	0.0137	0.4395
LAY	0.0473	0.0477	0.0333	0.0223	0.0572	0.0553	−0.0173	—	0.2959	0.0381	0.1045
CSM	0.1781	0.1978	0.1514	0.148	0.1595	0.0058	0.0956	0.0094	—	0.2432	<0.0001
ST	0.2639	0.2298	0.2129	0.1952	0.2183	−0.0275	0.1427	0.0801	0.0115	—	<0.0001
SNED	−0.0251	0.0459	0.0297	0.0466	0.0488	0.2309	−0.0021	0.0410	0.1864	0.2482	—

Φ_ST values are given below the diagonal, associated probability (P) values are above the diagonal. Colonies that are significantly differentiated from one another (P < 0.05) are shaded gray.

Based on the above pairwise Φ_ST comparisons, we conducted 2 different Amovas for the mtDNA data. The first (Amova-geography) compared a subpopulation containing the 3 plateau colonies (CSM, US, and ST) to a subpopulation containing all other colonies (Table 5). This analysis indicated that 19.2% of the genetic variation in the data was attributable to differences between those subpopulations (Φ_CT = 0.1923, P = 0.0054), whereas differences among sites within the subpopulations explained less than 1% of the variation (Φ_SC = 0.0092, P = 0.1387). In the second Amova (Amova-WNS), we grouped sites by their WNS status at the time of sampling (i.e., [DM, DUN, GL, LM, SHIM, SNED] vs. [CSM, ST, US, LAY, BC]). This grouping also loosely corresponded to their locations relative to the Allegheny front (CSM, ST, US, LAY, and BC are all located west of the front), whereas DM, DUN, GL, and SHIM are to the east, LM is northeast, and SNED is slightly southwest of the front. Although there was some weak genetic differentiation between a group including all western colonies (CSM, US, ST, BC, and LAY) and the remaining colonies, the second Amova-WNS suggested that BC and LAY were genetically more closely related to the eastern-central colonies than to the 3 western AP hibernacula. In this latter Amova, 13.24% of the genetic variance was attributable to differences between the 2 subpopulations (Φ_CT = 0.1324, P = 0.0042), compared with 19.23% in the Amova-geography. The proportion of genetic variance (2.82%) distributed among colonies within subpopulations (Φ_SC = 0.0325, P = 0.016) was approximately 3.5 times higher than in the Amova-geography (Table 5), suggesting that the LAY and BC colonies are not genetically homogenous with CSM, ST, and US.

Table 5.

Results of Amova of mtDNA sequences

	Amova-geography			Amova-WNS status
	Fixation index	P-value	% Genetic variation	Fixation index	P-value	% Genetic variation
Between subpopulations (ΦCT)	0.1923	0.0066	19.23	0.1324	0.0042	13.24
Among colonies within subpopulations (ΦSC)	0.0092	0.1339	0.74	0.0325	0.0160	2.82
Within colonies (ΦST)	0.1997	<0.0001	80.03	0.1606	<0.0001	83.94

The Amova-geography analysis was based on geographical features, comparing the 3 colonies on the western AP (CSM, ST, and US) to those in central and eastern PA and WV (SHIM, DM, DUN, LM, GL, BC, LAY, and SNED). The Amova-WNS status analysis categorized colonies by their WNS status at the time of sampling.

A Mantel test including all hibernacula indicated that mtDNA genetic distances between the colonies were correlated with geographic distance (correlation coefficient, r = 0.3223, P = 0.009). This suggests that the mtDNA differentiation identified can partially be explained by IBD.

The program STRUCTURE was utilized to determine the most likely number of independent genetic clusters based on the microsatellite markers. We tested K = 1 – 8 subpopulation clusters, but, in contrast to the mtDNA results, a single cluster (K = 1) scored the highest probability across all replicates (Supplementary Data). Because the lnProb(K) values for K = 1 and K = 2 were very similar, we examined cluster membership after harmonizing individual assignments to clusters using the software CLUMPP (Jakobsson and Rosenberg 2007) for K = 2. All but 2 individuals (n = 157 of 159 individuals genotyped) were assigned with the highest probability to the same cluster, indicating that K = 1 is likely the true number of clusters. Similarly, using the repeated reallocation method in FLOCK, we observed a plateau of mean LLOD values for K = 1, and plateaus were not achieved for any other K.

Amova analyses of the microsatellite data, treating all sampling locations as independent populations, indicated that 99.9% of nuclear genetic variation occurred within populations, and only 0.1% was observed among populations. Grouping locations by WNS status and by geography did not alter this result. Correspondingly, pairwise F_ST values among sampled colonies were consistently low and not significantly different from 0, ranging from −0.0054 to 0.0057 (Supplementary Data). There was no pattern of IBD (P > 0.8).

To determine whether these patterns of population substructure extend beyond PA and WV, north along the Appalachian Mountains, we obtained an additional 30 samples from Aeolus Bat Cave, a hibernating colony in Bennington Co, Vermont, which was first confirmed infected with WNS in 2008. We repeated the above-described analyses for both mtDNA and microsatellites. As mentioned previously, there was no substructure at the nuclear markers. However, the mtDNA marker indicated that, despite the geographic distance from PA and WV, Aeolus was not significantly differentiated from any colonies within the central and eastern subpopulation but was significantly differentiated from the 3 colonies (CSM, ST, and US) on the AP (Supplementary Data).

Historical Demography

We estimated population parameters under the isolation-with-migration model (Hey and Nielsen 2004) for mtDNA data, with daughter populations defined as per the population structuring analyses to include the 3 AP colonies (CSM, ST, and US) in 1 subpopulation and the more eastern colonies (DM, DUN, GL, LM, SHIM, BC, LAY, SNED) in the other subpopulation (Table 6). Six independent runs of at least 10 million steps converged on similar estimates for all parameters (Table 6). Because these analyses were focused on the historical demography of females and therefore based on a single locus, we excluded historical model parameters (θ_MRCA and τ) from the main results in an effort to avoid over-interpreting the analyses (results for these parameters are available from the authors).

Table 6.

Results of isolation-with-migration (IMa2) analyses of mtDNA sequences comparing AP colonies (CSM, ST, and US) to the remainder of colonies (SHIM, DM, DUN, LM, GL, BC, LAY, and SNED)

	AP		East		NmAP	NmEast
	θ	Ne	θ	Ne
Run 1	14.650 (7.55, 47.65)	171145 (88201, 556659)	24.350 (13.65, 89.25)	284463 (159463, 1042640)	13.72 (0.52, 114.42)	7.58 (1.07, 216.32)
Run 2	13.850 (7.05, 46.75)	161799 (82360, 546145)	25.150 (13.85, 90.35)	293808 (161799, 1055491)	10.82 (0.54, 112.73)	7.39 (1.19, 219.44)
Run 3	14.350 (7.65, 46.95)	167640 (89369, 548481)	24.250 (13.55, 89.05)	283294 (158294, 1040304)	3.75 (0.51, 113.68)	10.64 (1.07, 214.50)
Run 4	13.950 (7.45, 46.05)	162967 (87033, 537967)	26.050 (14.35, 90.55)	304322 (167640, 1057827)	4.93 (0.46, 111.04)	7.91 (0.95, 219.47)
Run 5	14.350 (7.65, 48.75)	167640 (89369, 569509)	23.750 (13.65, 89.95)	277453 (159463, 1050818)	5.61 (0.64, 117.67)	6.86 (0.97, 217.79)
Run 6	14.350 (7.65, 47.55)	167640 (89369, 555491)	25.450 (13.75, 89.75)	297313 (160631, 1048481)	3.25 (0.53, 114.89)	7.98 (1.15, 218.20)
Joint posterior	14.350 (7.45, 47.55)	167640 (87033, 555491)	24.850 (13.75, 89.85)	290304 (160631, 1049650)	9.38 (0.53, 114.65)	7.49 (1.08, 218.00)

Migration rate estimates are presented as generational rates of immigration into the population indicated. Error associated with each parameter estimate (shown in parentheses) was calculated as 95% confidence intervals. N_e, effective population size; Nm_x, number of individuals currently migrating into population X.

We found that estimates of θ for the eastern subpopulation were consistently higher than those for the AP subpopulation (θ_East = 24.850, θ_AP = 14.350). Assuming a mutation rate of 1.95% per million years and a generation time of 2 years per generation yielded an effective population size for the eastern subpopulation nearly 2 times greater than that for the AP population (N_e-East = 290304; N_e-AP = 167640). Comparisons of these parameters clearly supported a model in which θ_East is significantly larger than θ_AP [Pr(θ_East > θ_AP = 0.821);Pr(θ_AP > θ_East = 0.179)]. Furthermore, likelihood ratio tests rejected models specifying θ_East = θ_AP (P < 0.01).

Although these analyses did not provide clear point estimates of current migration, 2 significant patterns emerged regarding these parameters. Likelihood ratio tests rejected all models specifying m_x = 0 (P < 0.01), indicating that migration is occurring in both directions. Although point estimates of migration suggest a slightly higher rate of migration from east to west (Nm_AP = 9.38) than in the opposite direction (Nm_East = 7.49), parameter comparisons conducted in IMa2 could not reliably distinguish the 2 posterior distributions [Pr(M_AP > M_East = 0.501); Pr(M_East > M_AP = 0.499)].

Discussion

The results of this study suggest that female philopatry to wintering sites has led to matrilineal genetic structuring of hibernating colonies in northeastern United States. Furthermore, our data suggest that topography plays an important role in limiting female movements and matrilineal gene flow, and thus may represent an important and predictive limiting factor to the spread of WNS as the disease continues to spread toward the mountainous western United States. As we hypothesized, the Allegheny front may influence bat movement to some degree. Colonies CSM, ST, and US, located to the west of the front, were significantly differentiated at the mtDNA marker from colonies to the east and northeast of the front (SHIM, DUN, DM, and LM). However, BC and LAY, which are also located west of the front, showed weaker genetic affiliation with the other western colonies (CSM, ST, and US) than with the remaining colonies located to the east, north, and south of the front (SHIM, DUN, DM, LM, and SNED), suggesting that the Allegheny front is a permeable barrier to gene flow. In contrast, other topographic features such as the AP do appear to influence bat movement significantly. The mitochondrial substructure identified here appears to be driven predominantly by differentiation of the 3 colonies on the AP (CSM, ST, and US), with female-mediated gene flow between the AP colonies and those in the central and eastern Appalachian Mountains and lowland regions (Figure 1) being partially restricted. Approximately, 20% of the population genetic variance in the mitochondrial data can be explained by differences between these 2 clusters of sites, whereas less than 1% occurs among hibernating colonies within each of the clusters.

This pattern of genetic differentiation is correlated with geographic distances among the hibernacula (r = 0.3223, P = 0.009). However, this pattern of population substructure can only partially be explained by IBD. Colonies located within the mountainous Appalachian ridge and valley areas in central and eastern PA are genetically more similar (Φ_ST: −0.0089 to 0.0378, P > 0.05; Table 4) than those found closer to each other, but in different physiographic regions. For example, US (on the AP) and SHIM (ridge and valley region) are significantly differentiated from one another (Φ_ST = 0.2081, P < 0.0001; Table 4), but are separated by only 192 km, which is well within the migration abilities of this species. In contrast to US and SHIM, SNED in WV and DUN in northeastern PA are not genetically distinct from one another (Φ_ST = 0.0459, P = 0.0537; Table 4), despite being 566 km apart, and neither are SNED and Aeolus Bat Cave in Vermont, which are 840 km apart (Supplementary Data). Previous studies (Davis and Hitchcock 1965; Griffin 1970; Humphrey and Cope 1976) have identified M. lucifugus winter hibernacula located up to 300 km from summer maternity colonies, and some may be as far as 1000 km away (Wilson and Ruff 1999). Our results, therefore, suggest that topographical features may play a more significant role in determining patterns of colony connectivity than absolute distances.

Telemetry studies of female M. lucifugus at spring emergence further support our findings of limited female movements between regions (Butchkoski 2009). Individuals captured and radio-tagged at a hibernaculum in the mountains of central PA remained in this region and were tracked to summer colonies in nearby (up to 48 km) central PA counties (Butchkoski 2009). In addition, individual Indiana bats (M. sodalis) banded at maternity sites, fall swarming sites, or at hibernaculum-spring emergence in the mountains of central and southwestern PA were later recovered at hibernacula to the south in WV (Butchkoski 2009). Thus, indirect measures of gene flow from genetic data are largely in agreement with direct measures of movement from tagging studies.

In contrast to the above patterns of female movement, male-mediated gene flow, possibly via mating at swarming and/or hibernation sites, appears to occur throughout the region, with the result that there is no significant differentiation among sites at the nuclear loci (Supplementary Data and Supplementary Data Online). This pattern of female philopatry and male-mediated gene flow is typical for many temperate bat species, with swarming sites acting as hotspots for gene flow (e.g., Kerth et al. 2003; Rivers et al. 2005, 2006), and has been hypothesized for M. lucifugus in other parts of its range (Lausen et al. 2008; Dixon 2011). However, even at the mitochondrial locus, hibernacula on the AP are not completely differentiated from those in the central and eastern regions; there is some genetic similarity. This may be due to low (but nonzero) levels of female gene flow as suggested by the IMa2 analysis (Table 6) and/or may be the result of incomplete lineage sorting.

The pattern of structure observed in our mitochondrial data correlates with the observed spread of WNS infection of the colonies, with the 3 hibernating colonies on the AP (CSM, US, and ST) being among the last to be infected in PA. They were confirmed infected in mid to late 2010, whereas SNED in WV and the majority of other PA hibernacula were infected in the winter of 2009–2010 (Table 1). BC and LAY were also confirmed infected in winter 2010–2011, but they are located in the Monongahela transition zone (Woods et al. 2006) between the AP and the ridge and valley region (Figure 1). This transitional location may partially explain their genetic affiliations. Overall, these 2 colonies are genetically more closely related to hibernating colonies in the central mountainous and lowland regions to the east of the Allegheny front, but they do share more female-mediated gene flow with the AP colonies than do the more eastern sites. We, therefore, propose that this transition zone represents a region of admixture between female lineages from the AP hibernacula and those from hibernacula in the ridge and valley region. To test this hypothesis, we ran 1 additional Amova comparing the AP colonies (ST, CSM, and US) to the rest, but excluding BC and LAY. As expected, the resulting Φ_CT of 0.2264 (P = 0.012) was higher than in our first geography-based Amova, that is, the proportion of genetic variance that can be explained by differences between the subpopulations increased from 19.23% to 22.64%. This suggests that these 2 colonies may indeed represent an area of admixture between the 2 matrilineal subpopulations, and when they are excluded from the analysis, the genetic distinction between the subpopulations increases proportionately.

The pattern and timing of spread of WNS through this region may be partially explained by lowered gene flow between these genetic subpopulations, as female bat movement may have facilitated the initial spread of WNS along the Appalachians through PA, and south into Virginia, Maryland, and WV. In contrast, female philopatry to hibernacula may have delayed the spread of WNS to western PA. It is possible that the regional cluster of AP colonies forms part of a larger genetic subpopulation that extends further west, including hibernating colonies in neighboring states such as Ohio, Indiana, and Michigan. Colonies in Ohio and Indiana were also infected with WNS much later than those along the mountains and in the lowland regions to the east and south, and those in Michigan are as yet uninfected (US Fish and Wildlife Service 2012b). Unfortunately, without additional samples from hibernating colonies in these neighboring states, it is not possible to test this hypothesis.

In conclusion, this study identified significant levels of female-mediated genetic differentiation in PA hibernating colonies of little brown bats, which is correlated with topographical regions in the state. Coalescent-based analyses also revealed differences in the effective population size of these regional subpopulations, with the eastern ridge and valley hibernacula making up a significantly larger population than the western AP colonies. These analyses also indicate a low, but nonzero, rate of migration between the populations, but lack the power to reject a model of equal migration rates. Despite the limitations of available sampling, this study demonstrates the relevance of analyses of population structure in host species to understanding and predicting the spread of wildlife diseases.

Funding

Pennsylvania Game Commission State Wildlife Grant Program (Cooperative Agreement No. 231100); The Pennsylvania State University.

Acknowledgments

The authors thank C. Butchkoski, C. Stihler, S. Darling, and the AMCC of the American Museum of Natural History for assistance with sample collection. In Pennsylvania, work with live bats was conducted by personnel of the Pennsylvania Game Commission in compliance with the National Environmental Policy Act of 1969 (42U.S.C. 4321–4347) and as authorized by Pennsylvania law Title 34: Sections 104 and 322.

Author notes

Data deposited at Dryad: http://dx.doi.org/doi:10.5061/dryad.2pn66

Corresponding Editor: Warren Johnson