Landscape Genetic Connectivity and Evidence for Recombination in the North American Population of the White-Nose Syndrome Pathogen, Pseudogymnoascus destructans.

White-Nose Syndrome is an ongoing fungal epizootic caused by epidermal infections of the fungus, Pseudogymnoascus destructans (P. destructans), affecting hibernating bat species in North America. Emerging early in 2006 in New York State, infections of P. destructans have spread to 38 US States and seven Canadian Provinces. Since then, clonal isolates of P. destructans have accumulated genotypic and phenotypic variations in North America. Using microsatellite and single nucleotide polymorphism markers, we investigated the population structure and genetic relationships among P. destructans isolates from diverse regions in North America to understand its pattern of spread, and to test hypotheses about factors that contribute to transmission. We found limited support for genetic isolation of P. destructans populations by geographic distance, and instead identified evidence for gene flow among geographic regions. Interestingly, allelic association tests revealed evidence for recombination in the North American P. destructans population. Our landscape genetic analyses revealed that the population structure of P. destructans in North America was significantly influenced by anthropogenic impacts on the landscape. Our results have important implications for understanding the mechanism(s) of P. destructans spread.

1. Introduction

White-Nose Syndrome (WNS) is caused by the fungal pathogen Pseudogymnoascus destructans (P. destructans) (formerly Geomyces destructans [1]. Within hibernating North America (N. America)n bats, infections are characterized by the presence of white mycelial growth on the muzzle or wing tissues of bats, leading to the formation of ulcers and the erosion of epithelial tissues [2]. The N. American P. destructans population likely originated from a single European P. destructans migrant strain [3,4]. The European P. destructans population is known to infect local bat species but has no known associated mortality [5,6]. Since the first recorded case in N. American bats in 2006, P. destructans has caused the deaths of millions of bats involving multiple species [7]. Coupled with the increasing bat deaths is the expanding geographic range of WNS in N. America, which is now found in 38 US States and seven Canadian Provinces, including the west coast of N. America [8,9]. Despite the drastic declines in bat populations, relatively little is known about how P. destructans and WNS have spread among geographic regions. Because there was very little genetic variation in the initial population of Pd when an accidental transmission from Eurasia to North America occurred sometime around 2005/2006, tracking the spread and analyzing the population structure of the N. American P. destructans population has been difficult [10,11].

Early investigations into genetic structure of P. destructans populations showed results consistent with clonal spread of a single genotype in N. America [12,13,14]. If the pathogen were spread in a step-wise fashion from the center of outbreak origin in New York State, we should expect a relationship between the amount of accumulated mutation during the spread of P. destructans and the geographic distance from the site of initial infection. However, recent studies have failed to identify large-scale geographic or temporal structure in the N. American P. destructans population [10]. For example, the P. destructans isolate from the western-most Washington State was more closely related to the isolate from the epicenter in New York State than some other isolates from within New York State are to the earliest known isolate [11]. However, such a result could be due to chance events in sampling where the genetic relationships between specific isolates may not be representative of the whole populations, and consequently causing misinterpretations of how P. destructans has been transmitted between bats and their hibernacula. Alternatively, the spread of P. destructans may not follow a strictly step-wise model of outward expansion from the epicenter, but may instead be a more complex processes, influenced by multiple modes of transmission, and both forward and backward migrations between infected regions. In addition, secondary introductions of P. destructans genotypes between hibernacula could also influence genetic variation and complicate population structure, making the patterns of genetic variations inconsistent with the step-wise model of clonal expansion. Indeed, the introduction of additional foreign genotypes into N. America could further complicate the population structure and genetic variation. This is especially true if additional introduced strain(s) from Eurasia have a complementary mating type (MAT1-2) to that within the N. American P. destructans population (MAT1-1), which could significantly contribute to generating genotypic and phenotypic diversity [15,16].

Much about when and how P. destructans spread in N. America remains unknown. Broadly speaking, P. destructans could be spread by three methods: primarily via host bats, and secondarily through alternative hosts such as human activities or through passive mechanisms such as wind currents via the dispersal of airborne spores. The majority of WNS spread occurs via the dispersal of P. destructans propagules, which are transferred through contact events between hosts. Additional opportunities for transmission may incorporate the movement of bat ectoparasites [17,18], or the spread of propagules by predators of bat species [19], but the extent of their contribution to the spread of WNS may be limited to relatively small geographical distances. In comparison, the distances covered by N. American bats can sometimes extend upwards of 500 km within the migration season [20]. Indeed, the early estimated rate for the expansion of WNS was about 300 km each year, within the range of bats’ movement range within each season [21]. However, the high bat mortality rates (upwards of 90%) associated with P. destructans infection in some species result in few surviving individuals within most hibernacula impacted by WNS. Over particularly long distances, the activity of secondary hosts may also contribute to the spread of P. destructans. Indeed, human-mediated transmission of P. destructans was the most likely cause for P. destructans spread from Europe to N. America [22], and it could have continued to contribute to the dispersal of P. destructans between caves in N. America [23]. Lastly, the spores of some fungi are known to easily disperse over thousands of kilometers by wind and air currents [24,25]. While short-distance dispersal events of P. destructans through the air are unlikely [26], aerosols collected from both inside and outside of hibernacula have been found to contain P. destructans spores [27]. Ultimately, the extreme sensitivity of P. destructans to temperature and radiation [28,29,30] are likely to further dampen any chance of wind-aided dispersal over extensive distances [31].

The host species impacted by WNS are distributed across a diverse landscape with highly variable climatic conditions. Understanding how landscape features influence the spread of P. destructans is an important component of assessing the risk to the spread of further infections. In this study, we employed microsatellite markers and genomic variants to assess the genetic relationships among P. destructans isolates and populations, and to identify the potential patterns of dispersal that may exist among the locations sampled across the WNS distribution range in North America. Specifically, we compare the potential influences of landscape features with the null model that geographic distance was the main contributor to genetic distance between isolates. Here, the landscape features that we analyzed include geographical, environmental, or anthropogenic factors. We consider patterns of wind conductance as an additional abiotic factor that may influence connectivity, not as a viable mechanism for the transmission of P. destructans in N. America. The relationships of these environmental variables with genotype distributions of P. destructans strains across the landscape can also help infer the presence/absence of barriers to P. destructans spread. Furthermore, any one of these mechanisms are not mutually exclusive, and may contribute to dispersal depending on local conditions of the hibernacula and/or the surrounding climate. Identifying the dominant influences on the connectivity between bat hibernacula could help formulate targeted approaches to WNS control and prevention within different regions.

2. Materials and Methods

2.1. Strain Information

The strains analyzed in this study came from two sources: (i) those from published literature that contained genome sequences and/or microsatellite genotype information and (ii) our own strains described in an earlier study [14]. Briefly, our own samples were collected from hibernacula between 2008–2013; the majority were isolated from live bats, cave-associated arthropods, or hibernaculum walls. A few were collected from deceased bats confirmed to have died of WNS. All samples were revived from glycerol freezer stocks held at −80 °C and cultured on nutrient agar media at 14 °C and purified through sub-culturing at the hyphal tip of mycelial growth. These purified cultures matured for one month on Sabouraud Dextrose Agar (SDA) before genomic DNA was extracted using a standard cetyl-trimethylammonium bromide (CTAB) and phenol-chloroform protocol typical for filamentous fungal cultures [32]. Complete strain information is presented in Table S1.

2.2. Multilocus Genotyping and Bioinformatic Analyses

For microsatellite genotyping of our own strains, we chose the nine most polymorphic microsatellite loci from the panel developed by [10]. We genotyped 108 isolates using a multiplex PCR method, consisting of three reactions with each reaction containing three sets of primers. In each reaction, each of the three primer sets contained either a HEX, TET, or FAM fluorophore to allow scoring for all three fragments in each reaction. The multiplexed PCR recipe included: 6 μL of 2× GoTaq Green MasterMix (Promega: Madison, Wisconsin), 3 μL nuclease-free water, 0.5 μL of 10 μmol forward primer (total of 1.5 μL for a reaction amplifying three loci), 1 μL of 10 μmol reverse primer (total of 3 μL for a reaction amplifying three loci), and 1 μL of template DNA. PCR products were diluted 10× prior to fragment analysis. All samples were run by Mobix Lab at McMaster University using a GeneAmp PCR SYSTEM 9700 machine (ThermoFisher: Waltham, Massachusetts). The resulting fragment sizes were binned using the Microsatellite Plugin (V1.4, Biomatters: Auckland, New Zealand) for Geneious (4 November 2018), based on the known fragment sizes of isolates in [10] as references.

To further understand the potential geographic patterns of genetic variation in P. destructans populations, we analyzed single nucleotide polymorphisms from 41 published genome sequences from the Short Read Archive on NCBI, plus the genomes of three newly sequenced strains [29]. The three new genomes were obtained from strains originating from hibernacula in the Canadian Maritimes and were sequenced on an Illumina MiSeq platform with paired-end libraries. All genomic reads were aligned to the NHWC 20631-21 reference genome [33] using the Burrows-Wheeler Alignment (BWA) bwa mem algorithm (V0.7.15; [34]. We sorted BAM files, removed PCR duplicates, and added read group identifiers using PicardTools (V1.131; [35]. Single Nucleotide Polymorphism (SNP) variants were called from our sample cohort using HaplotypeCaller in gVCF mode (GATK: V4.1.2) at a minimum calling threshold of 20. Variants within repeat-rich regions were discarded after aligning back to the reference genome using NUCmer (V3.23; [36]. Missing calls in our multi-sample VCF were resolved using alignment information in each BAM using FixVcfMissingGenotypes in jvarkit [37]. We excluded variants within our re-sequenced sample of the reference genome and variants with low quality (QUAL > 1000) and depth (DP > 10). The final set of variants included 131 SNPs among 45 N. American P. destructans strains (see results below). We added annotations and prediction of effects of SNPs using SnpEff [38].

2.3. Population Genetic Analyses

We calculated several population genetic statistics using our microsatellite and SNP datasets. All population genetic statistics were generated using the R package poppr [39]. Our microsatellite dataset consists of 139 isolates, including 108 P. destructans isolates we genotyped in this study, plus 31 N. American and 5 European isolates (those with geographical coordinates available) initially genotyped by [10]. For the SNP dataset, we excluded SNP loci if greater than 5% of isolates had no ’REF/ALT’ call. If we identified multiple isolates from the same hibernaculum having the same multilocus genotype, only one of the clonal isolates was used in subsequent population and landscape analyses.

The patterns of allelic relationships among loci were measured using the index of association tests standardized by the number of loci used (). The null hypothesis for index of association measures is that there is random association among alleles at different loci. Statistical significance was derived by comparing the observed index of association to a null distribution (assuming random recombination) over 999 permutations. After a randomization test on a shuffled dataset, simulated p-values below 0.05 indicates that the null hypothesis of random recombination should be rejected. To further reveal which pairs of loci might show evidence of recombination, we also measured linkage disequilibrium (LD) between all pairs of loci.

The proportion of pairwise loci that were phylogenetically incompatible, also known as the four-gamete test (FGT) [40], was determined using Multilocus (V1.3; [41]). In short, the FGT considers all possible combinations of alleles between two loci. Loci are assumed to be phylogenetically compatible if there is no evidence of homoplasy or recombination. We considered the results of a FGT against a null distribution of 999 permutations of a shuffled dataset. We used SNP loci to determine genomic regions where potential breakpoints from recombination exist using FGTpartitioner [42]. Lastly, we conducted tests for recombination in both datasets using SplitsTree4 (V4.13.1) [43]. Population genetic structure was estimated using the STRUCTURE (V2.3.4) and fastSTRUCTURE (V2.7) [44] for microsatellite and SNP datasets, respectively, obtained from N. American isolates and European isolates.

2.4. Landscape Genetics

To analyze how landscape features might be associated with genetic relationships among isolates, we first generated genetic distance metrics from our SNP and microsatellite datasets using the first two principal components based on individual-level allele usage [45]. We followed the best practices for resistance surface optimization using the all_comb function in ResistanceGA [46]. Resistance surfaces were created using the commuteDistance algorithm in gdistance [47]. Optimized resistance surfaces were used in a series of Maximum Likihood Population Effects mixed-effects models (MLPE), using the function MLPE.lmm [46]. We considered multiple hypotheses to test the patterns of genetic distance between our sampled sites: (i) Isolation by Distance (IBD) that proposes that gene flow is a function of the Euclidean distance among populations; (ii) Isolation by Environment (IBE) that states that a higher level of gene flow should occur among locations with similar climate; (iii) Isolation by Resistance (IBR) that proposes that gene flow is a function of the resistance distance; and (iv) a null model that assumes the absence of any geographic structure. To evaluate the relative support of the competing models, we assessed each model’s fit using the Akaike information criterion (AIC), and conducted bootstrapping replications by subsampling response and distance matrices, then fitting them back to the MLPE model. The percentage of bootstrap replications in which each IBD, IBE, or IBR model was the best-fitting determines the level of support for each model.

To test the different models, we collected rasters of climate, elevation, anthropogenic impact, and wind conductance across regions overlapping the sampling locations for the isolates that were investigated in this study. The consensus of climatic differences were determined by the combining of 19 bio-climatic variables (worldclim.org, variables “bio1-19”, accessed 18 March 2020) into a single surface using raster principal component analysis (PCA), availible in the R package RStoolbox [48] in order to reduce the amount of multicollinearity between variables. Elevation data were obtained from R package elevatr [49]. We collected raster layers of anthropogenic impact from the Commission for Environmental Cooperation which estimates the impact of human activity across the landscape with an index of human influence; a measure which encapsulates land-use, light pollution, and human density, among others, to gauge the degree of disturbance across the landscape. We generated a raster surface of wind conductance using monthly averages of wind speed and direction with rWind [50]. All variables were aggregated to a resolution of ∼ 20 (∼26 km) prior to optimization.

Lastly, we estimated effective migration rates based on genetic distances using the Estimated Effective Migration Surfaces (EEMS) software package [51]. This method visualizes the influence of IBD across geographic space and creates estimates of effective migration rates by interpolation. Complementary to our MLPE approach using resistance surfaces, EEMS creates effective migration surfaces across space without incorporating environmental variables. We conducted three independent simulations of 200 demes with 1 million burn-in Markov chain Monte Carlo steps over 2 million iterations. Estimated migration surfaces were created from averaging the results of all simulation replicates using the rEEMSplots package [51].

3. Results

3.1. Genetic Relationships among North American P. destructans

We obtained multilocus microsatellite genotypes from 108 N. American P. destructans isolates in our collection using primers targeting nine microsatellite loci (Table S2). These isolates belonged to 42 multilocus genotypes (MLGs). We present these results in combination with the MLG data from [10] to provide a greater context for the genetic relationships among P. destructans isolates in N. America (after clonal correction, , MLGs ; Figure 1). Within this combined dataset, and after correcting for clonal genotypes among different locations, we found an overall high allelic diversity (Simpson’s Diversity; ), and well over half the MLGs were exclusive to a single isolate. However, 35 isolates distributed among caves in New Brunswick, Prince Edward Island, New York, Ontario, Québec, and Vermont belonged to a single microsatellite genotype (Table 1). Aside from this dominant genotype, 12 other genotypes were also shared by two or more isolates, often from multiple US States/Canadian Provinces.

Figure 1

Map of sampling locations from across North America, with each point representing one of more isolate(s) genotyped using whole genome SNPs (blue), multiple microsatellite loci (green), or a combination of both (red). Many of the isolates used here are the result of genotyping in previous studies (See Table S1).

Table 1

Shared microsatellite genotypes present in multiple caves across North America. The isolates reported here for each genotype have been clone-corrected to exclude clonal individuals from the same sampling site.

MLG Code	Number of Isolates	Province/State, Country	Site
MLG.9	3	New Brunswick, Canada	Glebe Mine
			Markhamville Mine
			White Cave
MLG.11	3	New Brunswick, Canada	Glebe Mine
			Howes Cave
MLG.13	2	New Brunswick, Canada	Glebe Mine
			White Cave
MLG.15	2	West Virginia, USA	Pendleton
			Tucker
MLG.16	2	Nova Scotia, Canada	Falmouth
		New York, USA	Williams Hotel Mine
MLG.25	4	Not Available	Not Available
		North Carolina, USA	Yancey
		Ontario, Canada	Not Available
MLG.33	2	Tennessee, USA	Not Available
		Tennessee, USA	Not Available
MLG.35	2	New Brunswick, Canada	Dorchester Mine
		Prince Edward Island, Canada	Uigg
MLG.40	8	Connecticut, USA	Not Available
		Delaware, USA	Not Available
		Massachusetts, USA	Not Available
		Missouri, USA	Not Available
		Tennessee, USA	Not Available
		Virginia, USA	Not Available
		Not Available	Not Available
MLG.47	7	New Brunswick, Canada	Markhamville Mine
			Not Available
		Not Available	Not Available
MLG.69	2	New Brunswick, Canada	Berryton Cave
			Harbells Cave
MLG.72	7	New Brunswick, Canada	Berryton Cave
			White Cave
		Prince Edward Island, USA	Rocky Point
			Vernon Bridge
MLG.73	3	New Brunswick, Canada	Berryton Cave
			Harbells Cave
			White Cave
MLG.75	6	New Brunswick, Canada	Berryton Cave
			Dorchester Mine
			Markhamville Mine
			White Cave
		Ontario, Canada	Not Available
MLG.80	2	Indiana, USA	Not Available
		West Virginia, USA	Not Available
MLG.82	4	New Brunswick, Canada	Berryton Cave
			Dorchester Mine
			Markhamville Mine
		Prince Edward Island, USA	Murray River
MLG.83	35	New Brunswick, Canada	Berryton Cave
			Dorchester Mine
			Glebe Mine
			Harbells Cave
			Markhamville Mine
			White Cave
		North Carolina, USA	Avery
		New York, USA	Williams Hotel Mine
		Ontario, Canada	Not Available
			Not Available
			Not Available
			Not Available
			Not Available
		Prince Edward Island, Canada	Rocky Point
		Québec, Canada	Not Available
		Virginia, USA	Greely Mine
MLG.84	3	Ontario, Canada	Not Available
		Not Available	Not Available
MLG.87	6	New Brunswick, Canada	Berryton Cave
			Glebe Mine
			Harbells Cave
			Markhamville Mine
			White Cave

Complementary to the nine microsatellite marker data, we also analyzed the genome sequences of 41 N. American P. destructans strains. A total of 54 SNPs were identified and these SNPs resolved the 41 strains into 40 SNP genotypes. Two isolates from Glebe and Markhamville Mine in New Brunswick shared an identical SNP genotype. Similar to what we observed within microsatellite genotypes, SNP diversity was high ( = 0.966, Table 2) compared to previous observations in the N. American population ( = 0.69) [10]. All SNPs reported here were unique to N. America and have not been found within European P. destructans genomes.

Table 2

Microsatellite genotypes and population genetic metrics from 139 isolates presented in this study. All metrics were clone-corrected at the Province/State level. Tests of standardized index of association () are paired with simulated p-values generated from a randomization tests, conducted over 1000 permutations of the test statistic.

Region	Country	State/Province	N	MLG	λ †	Hexp ‡	r¯D	r¯Dp-Value
North America	Canada	New Brunswick	23	23	0.957	0.38	0.065	0.65
		Nova Scotia	1	1	0	-	-	-
		Ontario	11	11	0.91	0.32	−0.013	0.87
		Prince Edward Island	6	6	0.83	0.24	−0.069	0.88
		Québec	2	2	0.5	0.22	-	-
		Total	97 (not clone-corrected)	36	0.97	0.38	0.02	0.89
	USA	Connecticut	1	1	0	-	-	-
		Delaware	1	1	0	-	-	-
		Indiana	2	2	0.5	0.22	-	-
		Massachusetts	2	2	0.5	0.33	-	-
		MD	1	1	0	-	-	-
		ME	1	1	0	-	-	-
		Missouri	1	1	0	-	-	-
		North Carolina	3	3	0.67	0.33	0.83	0.78
		NJ	1	1	0	-	-	-
		New York	5	5	0.8	0.41	0.12	0.12
		OH	1	1	0	-	-	-
		PA	2	2	0.5	0.67	-	-
		Tennessee	3	3	0.67	0.22	0.83	0.77
		Virginia	3	3	0.67	0.37	0.4	0.16
		VT	2	2	0.5	0.44	-	-
		West Virginia	5	5	0.8	0.4	0.12	0.1
		Total	37 (not clone-corrected)	26	0.96	0.46	0.04	0.03
	Total		134 (not clone-corrected)	59	0.98	0.46	0.03	0.16
Europe	Total		5	5	0.8	0.51	0.08	0.27

† = Simpson’s Diversity and ‡ = Expected Heterozygosity.

3.2. Linkage Disequilibrium and Tests for Recombination

Standardized index of association () values are often used in tests of LD and are compared to a null distribution of values, representing linkage equilibrium (LD) as with random recombination in the genome. If the observed is significantly different from the null distribution, the null hypothesis of LD would be rejected and the observed would be consistent with LD. After censoring clonal isolates, analyses of microsatellite ( = 0.02, p = 0.22) and SNP genotypes ( = 0.006, p = 0.23) failed to reject the null hypothesis. These tests of index of association were also performed in a hierarchical manner, to identify patterns of linkage between pairs of loci (Table 2). With further partitioning, showed patterns of association not significantly different from LD in most loci combinations of microsatellite (72%) and SNPs (94%). Only ten of the 36 microsatellite pairs were inconsistent with LD (Figure 2A); whereas 83 of the 1 431 SNP pairs (∼6% of all combinations) were not in LD (Figure 2B). These results demonstrate that most loci are in LD within the N. American P. destructans population.

Figure 2

Pairwise comparison of the results from standardized index of association tests (), conducted on Microsatellite (A) and SNP (B) datasets. Squares of the heatmap outlined in black highlight combinations of loci where p-value < 0.05, consistent with null-expectations under clonality. Samples missing greater than 5% of genotype calls were excluded.

The above index of association and LD tests used random recombination as the null hypothesis to generate the expected genotype counts. We also investigated whether loci pairs are incompatible with strict clonal reproduction in the N. American P. destructans population using the FGT. Overall, 50% of microsatellite loci were phylogenetically incompatible (Figure 3A, p < 0.01). However, phylogenetic incompatibility was not evenly distributed, as combinations of alleles between some loci (i.e., Pd1/Pd11/Pd19 and Pd7/Pd13) had higher proportions of phyllogenetically compatible alleles (50–12.5%, respectively) compared to others (p < 0.01; Figure 3A). Considering all combinations of SNP loci, 25% were inconsistent with phylogenetic compatibility (p > 0.05; Figure 3B). Consistent with results from linkage equilibrium and phylogenetic compatibility analyses, the neighbour-network analysis of the N. American P. destructans isolates also showed prevalent closed loops (Figure 4) and evidence of recombination in the population. As expected, greater evidence for recombination in European isolates were found. For example, the FGT suggested the presence of four recombination breakpoint events in the N. American P. destructans population while 402 breakpoints were estimated in the European population (Figure S1). Taken together, our results indicate unambiguous evidence for recombination in the North American P. destructans population.

Figure 3

Pairwise comparison of bi-allelic loci where the presence of more than three allele combinations results in phylogenetic incompatibility. Within our microsatellite dataset (A) the numerator shows the number of pairs that are compatible, while the denominator lists the total number of allelic combinations between the loci. (B) SNP loci with greater than 5% of samples missing calls were dropped.

Figure 4

Neighbor-network constructed using SplitsTree using (A) 9 microsatellite loci and (B) 53 SNP loci. Parallel edges in this network indicate incongruence with a perfect monophyletic tree.

A neighbour-network with N. American P. destructans isolates showed the prevalence of closed loops which indicates incongruence with a perfect monophyletic tree (Figure 4), consistent with recombination and the presence of recombinant lineages in the population. However, tests for recombination conducted on these collections are not consistent with frequent recombination (SNPs = 0.138, p = 0.9). Whereas, with the inclusion of European strains in these tests, is consistent with frequent recombination (SNPs = 0.122, p < 0.01). FGT loci suggest the presence of four breakpoint events within the N. American P. destructans population, compared to the 402 FGT breakpoints estimated from EU loci in the P. destructans genome (Figure S1). Taken together, our results indicate unambiguous evidence for recombination within the N. American P. destructans population.

3.3. Population Structure

To visualize the population structure of P. destructans in N. America, we created unrooted phylogenies using N. American strains genotyped by microsatellites and SNPs, respectively. As revealed by microsatellite genotype information, common genotypes (e.g., MLG 83) and genotype clusters often include isolates from multiple locations (Figure 4). While we detected a statistically significant positive correlation between geographic distance and genetic distance based on microsatellite genotypes, that correlation was very weak ( = 0.05, p > 0.01; Figure 5A). Furthermore, there was no correlation between geographic distance and genetic distance based on SNPs ( = −0.002, p = 0.9; Figure 5B). These results are consistent with frequent gene flows between geographic populations.

Figure 5

Correlation between genetic distance of (A) microsatellites and (B) SNPs and logarithm of geographical distance of all pairwise combinations of isolates.

Using STRUCTURE, we identified the optimal number of clusters () for the P. destructans genomes analysed in this study from N. America and Europe. Interestingly, the two genotyping methods used here present different patterns in ancestral population structure. Based on microsatellite genotyping, the majority of USA isolates and all European isolates were assigned to one cluster separated from the other cluster that contained mostly Canadian and some US isolates (Figure S2A). In comparison, based on SNP markers, all isolates from Canada and US strains formed a cluster that also contained two isolates from Europe (Germany and Switzerland) while the other cluster contained only European isolates (Figure S2B).

3.4. Landscape Genetics

We examined the potential influence of multiple landscape factors on the patterns of P. destructans genetic distance in the microsatellite and SNP genotype datasets. We measured the performance of several models to test hypotheses of IBD, IBE, and IBR. The summary results of these tests are presented in Table 3. Interestingly, the model selection results differed between microsatellite and SNP datasets (Table 3). Models incorporating estimates of anthropogenic influence on the landscape (AICc = 1221.26, bootstrap = 88.9%, Table 3) was the best fitting model compared to all others considered for the microsatellite dataset, explaining patterns of genetic distance between individuals better than all other models. Overall, the microsatellite dataset demonstrated that genetic distance increases with greater resistance across human impacted landscapes (slope = 0.16, t-value = 10.6). These results suggest that the anthropogenic factors, such as differences in land use and population density, may have prominent impacts on the connectivity between bat hibernacula, and in-turn the genetic similarity between P. destructans isolates.

Table 3

Assessment of environmental variables contributing to P. destructans genotype distributions across landscape. Using a maximum likelihood population effects mixed-model, we measured the effects of various environmental variables/resistance surfaces to predict changes in pairwise genetic differences (PCA scores) generated from microsatellite or SNP alleles. AICc: adjusted Akaike information criterion; weight; marginal values of the fitted MLPE model. LL: log likelihood; Percent Top: the percentage of times each model was the best-fit model in over 1000 bootstrap replications.

	Variable(s)	Mean AICc	Mean Weight	Mean Marginal R^2	Mean LL	Percent Top	k
MSAT	Human Influence	1217.11	0.56	0.08	−604.01	90.5	4
	Climate × Human Influence	1221.43	0.07	0.08	−605.07	0	7
	Wind × Human Influence	1221.53	0.07	0.08	−605.12	0	7
	Climate	1222.17	0.06	0.06	−606.55	0.1	4
	Elevation × Human Influence	1221.64	0.06	0.07	−605.17	0	7
	Wind	1223.34	0.12	0.06	−607.13	9.4	4
	Climate × Elevation	1224.69	0.02	0.06	−606.7	0	7
	Climate × Wind	1225.9	0.02	0.06	−607.3	0	7
	Climate × Elevation × Human Influence	1227.26	0	0.07	−606.08	0	10
	Distance	1227.83	0.01	0.05	−609.76	0	2
	Elevation × Wind	1227.97	0	0.05	−608.34	0	7
	Elevation	1228.45	0.01	0.05	−609.68	0	4
	Climate × Elevation × Wind	1228.86	0	0.06	−606.88	0	10
	Climate × Wind × Human Influence	1228.93	0	0.05	−606.92	0	10
	Elevation × Wind × Human Influence	1229.25	0	0.06	−607.08	0	10
	Climate × Elevation × Wind × Human Influence	1234.15	0	0.07	−606.57	0	13
SNP	Distance	1211.72	0.25	0	−601.52	46.3	2
	Climate	1212.3	0.21	0.02	−600.9	23.7	4
	Human Influence	1212.09	0.24	0.01	−600.8	25	4
	Elevation	1213.18	0.13	0.01	−601.34	4.9	4
	Wind	1213.42	0.11	0	−601.46	0.1	4
	Climate × Human Influence	1218.24	0.01	0.02	−600.81	0	7
	Climate × Wind	1218.57	0.01	0.01	−600.98	0	7
	Wind × Human Influence	1218.4	0.01	0.01	−600.89	0	7
	Climate × Elevation	1218.58	0.01	0.01	−600.98	0	7
	Elevation × Human Influence	1218.43	0.01	0.01	−600.91	0	7
	Elevation × Wind	1219.38	0.01	0	−601.38	0	7
	Climate × Elevation × Human Influence	1231.92	0	0.01	−600.96	0	10
	Climate × Elevation × Wind	1232.13	0	0.01	−601.07	0	10
	Elevation × Wind × Human Influence	1231.96	0	0.01	−600.98	0	10
	Climate × Wind × Human Influence	1232.18	0	0.01	−601.09	0	10
	Climate × Elevation × Wind × Human Influence	1262.07	0	0.01	−601.04	0	13

Based on the SNP dataset, the MLPE model with the highest bootstrap support was IBD (51.3%), followed by IBR—Human influence (22.7%), and IBE (21%; Table 3). These results suggested IBD as an important driver of the relationships among strains in N. American P. destructans. However, the relationship between genetic and geographic distance is relatively weak (slope = −0.002, t-value = −0.01). Furthermore, models suggesting isolation via human influence on the landscape consistently had the lowest average AIC and highest average marginal scores compared to all other models tested (Table 3). Together, these analyses also indicate a positive relationship between genetic distance of P. destructans strains and resistance to human influence across the landscape (slope = 0.34, t-value = 1.26).

We estimated effective migration surfaces across the landscape for an area roughly   in size, equivalent to ∼1/5 of the total area of the United States. Figure 6 shows areas of high/low estimated migration rates between locations using both the SNP and microsatellite datasets. Consistent between these datasets is a large section of the distribution with low rates of migration, extending in a south-west/north-east direction (Figure 6A,B). As the estimated migration rates in this region are lower than what would be anticipated under IBD, this suggests reduced gene flow between these areas of the WNS distribution. Flanked by this region are parts of the distribution that have comparably higher estimated rates of migration, greater than expected under the null hypothesis of IBD (Figure 6A,B). Taken together, estimated migration surfaces for both the SNP and microsatellite datasets suggest a disconnect between the south west and north east corners of the distribution.

Figure 6

Estimated effective migration surfaces (EEMS) plot for P. destructans based on (A) SNP and (B) MSAT markers. Considering a maximum of 200 demes across the landscape, EEMS generates posterior mean migration rates (log10), indicating areas which have a higher (blue) or lower (orange) rate of migration (shown in blue) than expected under isolation by distance (IBD).

4. Discussion

This study investigated the patterns of genetic variation in the N. American population of P. destructans, including the potential landscape factors influencing strain relationships and population structure. Using microsatellite markers and whole-genome sequence data, we found substantial allelic diversity, including novel alleles not reported in prior studies with these markers [10,52]. In addition, we found that the N. American P. destructans population is inconsistent with the expectations of strict clonality. We built upon similar recent studies by including a collection of isolates from regions not previously investigated [10,11]. Our analyses revealed limited signals of population genetic structure based on a standard population genetic approach. Our analyses identified connectivity of P. destructans isolates strongly influenced by climate and human impacts on landscapes in N. America. Below we discuss the factors influencing the genetic diversity of P. destructans across the landscape.

4.1. Multilocus Genetic Variation

In clonal microbial pathogens, the rapid generation of novel MLGs via mutation can help overcome selective pressures and niche heterogeneity challenges [53,54,55]. For a recently introduced pathogen that primarily reproduces asexually, we expect the N. American P. destructans population to harbor low genotype diversity [56,57]. In comparison, an ancient population capable of both asexual and sexual reproduction would have high genotype diversity, with evidence of loci in LD and phylogenetically incompatible loci [58,59]. Our analyses based on both microsatellite and SNP genotype data are inconsistent with strict asexual reproduction of the N. American P. destructans population (Figure 2 and Figure 3). The occurrence of hyphal fusion between members of the same mating type is not uncommon within clonal populations, as parasexual recombination can be a major force of evolution among asexual fungi, causing increased genetic diversity and the emergence of new pathotypes [60]. Many other species of Ascomycota maintain cellular machinery necessary for the parasexual life cycle [61]. Although mitotic recombination has been described within other Pseudogymnoascus genera [62], more evidence is needed to distinguish the mechanisms responsible for the allelic patterns observed in this study.

Evidence of recombination has been reported in several seemingly asexual fungal populations where only a single mating type was present (e.g., [63]). A recent study by [64] also suggested evidence of recombination based on sequences from multiple loci in a collection of global P. destructans strains. However, they compared their data to only the null model of random recombination. With relatively limited sample size, it might be difficult to reject the null model. In contrast, our analyses based on a much larger sample size compared the data to both the strictly asexual reproduction model and the random recombination model. Thus, we believe that our results provide much more robust support for recombination in the N. American P. destructans population.

As different mating types have been found to co-occur within the same ecological niche [4], mating between complementary mating types and sexual reproduction might be common within the ancestral range of P. destructans in Europe and Asia. However, only one mating type (MAT1-1) has been reported in the N. American P. destructans population (Table S2) [4]. As a result, recombination in the N. American P. destructans population unlikely arose through traditional heterosexual mating. Instead, same-sex mating and/or parasexuality, as have been reported for other fungi [65], could have contributed to recombination in the N. American P. destructans population. Alternatively, there could be multiple introduction events of strains of the same mating type from Europe where recombination is likely common [66]. More extensive sampling is needed in order to critically evaluate these two possibilities.

The two types of genotype data provided complementary insights into the structure of P. destructans in N. America. However, there were several differences worth noting. For example, in this dataset, SNPs were useful for distinguishing European isolates of P. destructans and for identifying those from Europe that were closely related to N. American strains (Figure S2B). In addition, patterns of inferred migration differed slightly between the two genotyping methods (Figure 6). Some differences were likely due to the different sample sizes and distributions of the analyzed isolates by the two genotyping methods. Here, more locations were sampled at a higher depth using the microsatellite markers than the SNPs while very few isolates were genotyped using both methods. Furthermore, due to the different mutational processes for the two types of markers, models explaining population structure based on microsatellite markers may be more useful for revealing fine-scale genetic structure while SNP results are more suited for large-scale and long-term inferences.

Indeed, previous studies have suggested that the rate of mutation in microsatellite loci is significantly higher than nucleotide substitutions. Consequently, homoplasy will likely be less common for SNP markers than for microsatellite markers. Here, homoplasy refers to a microsatellite allele or a SNP shared by different strains due to parallel mutational events but not due to their common ancestor. Homoplasy could contribute to LD and phylogenetic incompatibility. While the true rate of homoplasy for either type of markers is unknown in the N. American P. destructans population, if present at a significant rate, we expect that the number of phylogenetically incompatible pairs of loci should be higher for loci with a higher number of alleles per locus. However, there is limited statistical support for this hypothesis (Figure 4 and Table S1), suggesting that homoplasy is unlikely to be the major cause of phylogenetic incompatibility and the observed evidence for recombination. Furthermore, we did not find any evidence for the assayed microsatellite markers or the SNPs being in regions of mutational hotspots. For example, we found that none of the SNPs were distributed close to the 1 kb region flanking microsatellite markers based on our whole-genome analyses.

4.2. Factors Contributing to Pseudogymnoascus destructans Population Structure

The spread of P. destructans in N. America was initially believed to be due to clonal expansion [13]. As such, we would expect P. destructans to accumulate genetic diversity as it spread from the epicenter of the outbreak. Overall, we found statistically significant support for a correlation between genetic distance based on microsatellite markers and geographic distance within N. America. However, the correlation coefficient was very weak ( = 0.05) and many isolates with identical genotypes were shared across distant sampling sites (Figure 5; Table 2). Indeed, almost half of the isolates genotyped using microsatellite markers belong to a single genotype (MLG 83) and this genotype was distributed across several Canadian Provinces/US States (Figure 4 and Table 2). Furthermore, there was no correlation between SNP-based genetic distance and geographical distance among isolates (Figure 5). Taken together, this suggests that gene flow is frequent between sites within regions impacted by WNS.

Our interpretation of P. destructans population structure suggests limited correlation with genetic distance (Figure 4 and Table 3). Even through the impact of P. destructans infections varies depending on species characteristics [67], the genetic structure of Myotis lucifugus has been found to be both consistent [68,69] and inconsistent [70] with the spread of infections. Indeed, the effect of geographic distance on P. destructans population structure may be limited [68,71] as long-distance dispersal events push the expansion front and back-transmission fills in behind [72,73].

Climatic patterns influence bat population structure via the seasonal regulation of migration behaviours, impacting when and where host bats congregate in social groups [74,75]. An understanding of the migration patterns of N. American Myotis between winter and summer roosts is crucial for understanding the patterns of WNS transmission [68,76,77,78,79,80,81,82]. Regional similarities in climate appear to be important in regulating the timing of seasonal migrations which may correspond to minimizing energy expenditure [83,84] or prey availability [85]. We report no substantial influence of seasonal variations in climate on the genetic connectivity between P. destructans strains (Table 3). Yet, estimated migration rates between sampled locations appear to be consistent with landscape resistance to climate (Figure 6B and Figure S3A); this result indicates that additional climatic measures (or other variables that vary with climate) need to be included to better explain population structure.

Variation in landscape topography or climate is an important driver of the phylogeographic patterns found in a wide variety of taxa [86]. For instance, in Eastern United States and Canada, the Appalachians present a major geological feature that overlaps with bat migratory routes [22]. The high plateau and extending mountain range of the Allegheny Front escarpment may act as a barrier to gene flow, as colonies of little brown bats located on the western side of the Appalachian high plateau were infected with WNS 1–2 years later than colonies in central or eastern Pennsylvania [8,68]. While this feature of the landscape definitely has a local impact on bat populations, we find no influence of altitude on the overall structure of P. destructans populations (Table 3). Similarly, while wind dispersal may be a common mechanism of passive dispersal for some fungal species [25,57], we find the genetic connectivity among P. destructans isolates is not consistent with prevailing currents (Table 3).

Variation in climate across the landscape could potentially impact the survival of P. destructans propagules [29,30] in environmental reservoirs during the summer months. Recent detection of P. destructans on free-flying bats during the summer months [23,87] suggests that cells of this fungus can persist on bats for some time. In comparison, positive detection of P. destructans appears to be much more common on the skin of European bats [88,89]. As such, either the persistence of P. destructans over the summer months or the seasonal shedding and recolonization of the skin microbiome from the environment could enable WNS transmission in N. America. As bat activity decreases at temperatures ≤ 10 °C [90], dormant P. destructans cells could then propagate when conditions are favorable [67,91,92].

Human activities could potentially influence the transmission of WNS, by directly facilitating the transmission of P. destructans between caves and mines [3,22]. Although commonly isolated from the sediments of caves impacted by WNS [93], P. destructans does not appear to be commonly picked up by cave visitors [23,94] and the frequency of human visits to rural caves is unknown. Instead, indirect influences on transmission may be impacted by the consequences of urbanization/land use [95,96].

Our results suggest that areas with high anthropogenic impacts have lower connectivity between locations with regards to P. destructans genotypes. Compared to forested areas, both species richness and relative abundance of bats are lower in urban areas [97,98]. At higher latitudes of the WNS distribution, the degree of human impact on the landscape drops off quickly (Figure S3C). Regional migration of bats occurs in all directions, not just latitudinally [99]. For instance, different migration strategies may impact the distance covered by certain species during seasonal migrations. Overwintering at northern latitudes is energetically expensive, and hibernating in social groups is one strategy that can reduce the amount of energy needed to survive the winter [100]. Instead, longer seasonal migrations may take place, in order to reach warmer hibernacula in the south. This has been recorded in Perimyotis subflavus, as individuals summering at northern latitudes have been found to migrate farther south [101]. Yet, some P. subflavus still overwinter in the hibernacula found at the northern edge of their distribution alongside other bat species [102]. These differing migration strategies may serve as an explanation for why these regions may remain well connected Figure 6 and the long distance dispersal of P. destructans.

Although some insectivorous bats do remain in large urban environments [103,104], some Vespertilionidae may have a higher preference for natural areas than urban areas in N. America [105]. With a higher degree of urbanization and fragmentation of suitable habitat, population sizes are decreased [106,107]. As such, high levels of urban development seem to restrict the frequency and length of migration events, reducing the connectivity between populations separated by anthropogenic impact on the landscape [105]. If bat migrations encounter resistance from anthropogenic change to the landscape, then by extension, they can impact the genetic structure of P. destructans.

The breadth of sampling for WNS may be limited, as the understanding of where particular environmental sources and sinks exist could be further developed. Much of the least populated areas are less likely to be surveyed as those near city centers. As our analyses suggest higher rates of migration that span areas with lesser human influence, where resistance to migration (and connectivity of P. destructans strains) is relaxed (Figure 6), more thorough sampling of these areas could reveal the impact of human activity on the spread of P. destructans.

5. Conclusions

The increasing prevalence of fungal epidemics emphasizes the importance of understanding how fungal pathogens spread and invade new ecological niches [108]. Rapid and effective genotyping methods are essential for monitoring wildlife diseases, especially WNS, as the high mobility of host species and airborne nature of pathogens themselves present significant challenges in tracking pathogen transmission. In recent years, high-throughput DNA sequencing technologies have improved the accuracy of genetic variation estimates, providing abundant SNPs often with greater discriminating power compared to a few microsatellite markers in most studies [109,110,111]. In this study, we have found that these genetic markers can provide complementary information for inferences about P. destructans population structure. Both types of genotype data showed deviations from the null expectations of strict clonality in N. America. In addition, we find that the N. American P. destructans population could not strictly be explained a stepping-stone model, involving IBD scenarios of dispersal. Instead, N. American population structure of P. destructans has been influenced by anthropogenic impacts across the landscape and IBD and IBE may play secondary roles.

More complex models incorporating seasonal changes across the landscape, including more robust estimates of P. destructans prevalence in environments, and measurement of fungal loads on free-flying bats are needed to better understand the population structure of P. destructans in N. America.