We surveyed five sites in Virginia and New York during the fall of 2011 to collect samples from bats and cave substrates at the start of the hibernation period (Table 1). Sites in both locations in the eastern U.S. were within the white-nose syndrome (WNS) epidemic area at the time of sampling, however no bats showed active infection at the time of sampling and all were tested for the presence of P. destructans. In addition, during the summer and fall of 2012 and 2013, bats and cave substrates were sampled from 14 sites across the state of Colorado. Collecting permits were granted for each state (CO, NY, VA) and work was conducted with an approved IACUC protocol (#1305.06 from the University of Colorado). Sampling was conducted opportunistically at different sites (including caves, mines, and buildings) where bat populations or cave habitats were being monitored by different state and federal agencies. We sampled four species of hibernating bats in the eastern United States: the little brown bat (Myotis lucifugus), tricolored bat (Perimyotis subflavus), northern long-eared bat (Myotis septentrionalis), and the Indiana bat (Myotis sodalis). In Colorado, we sampled eight species across a variety of habitat types: long-legged myotis (Myotis volans), little brown bat (Myotis lucifugus), western small-footed bat (Myotis ciliolabrum), Yuma myotis (Myotis yumanensis), big brown bat (Eptesicus fuscus), western long-eared bat (Myotis evotis), hoary bat (Lasiurus cinereus), and Townsend's big-eared bat (Corynorhinus townsendii). One species, M. lucifugus, commonly inhabited sites in all three regions (CO, VA, and NY), allowing us to assess geographical and environmental effects in greater depth.

To target the bacterial communities of bat skin as well as the surrounding environment, a sterile swab protocol was established and used across sampling teams. Each bat was individually sampled using a sterile rayon or cotton swab that was moistened with autoclaved (sterile, DNA free) water, and care was taken to use sterile technique to prevent cross-contamination. Bats were captured using either harp or mist nets that were sterilized prior to use at each site. Captured animals were removed from traps and handled with clean, sterile gloves that were changed after every individual. The skin of the bat along the forearm and muzzle were swabbed firmly with the same swab tip five times along each targeted area, for a total of ten passes (Langwig et al., 2015). After skin swabs were taken, the weight, sex, species, and hibernaculum site were recorded for each individual. If access to the roost site was available, environmental samples of roost walls (rock, mine shaft, house) were taken by passing the swab ten times over the surface in areas where bats were roosting. If soil was available in the roost site, small samples were taken directly beneath roosts using a sterilized scoopula and placed in small, sterilized plastic bags and were sampled using swabs for final analysis. In downstream analyses, environmental sample refer to those sites where it was possible to collect any of the above types of samples. Samples were stored on ice until they could be shipped or transported to the University of Colorado Boulder, where they were kept in a −20°C freezer prior to DNA analysis and sequencing.

DNA extraction, library preparation, and sequencing of the 16S rRNA marker gene was conducted as in Fierer et al. (2012). Briefly, DNA was extracted using the MoBio PowerSoil extraction kit following the manufacturer's protocol. PCR amplification was conducted using the primers 515F/806R to target the V4 region of 16S rRNA and contained 12 bp barcodes as well as Illumina sequencing adapters following Caporaso et al. (2011). The PCR reactions contained 11 μL PCR water, 10 μL 5 Prime Master Mix, 1.0 μL each of the forward and reverse primers, 1.0 μL MgCl₂, and 1.0 μL genomic DNA. The thermocycling conditions for PCR consisted of an initial denaturation step of 94°C for 3 min, followed by 35 cycles at 94°C for 45 s, 50°C for 60 s, 72°C for 90 s; and final extension of 10 min at 72°C. Each sample was amplified in triplicate and combined. Amplicons were quantified using the Quant-IT Picogreen dsDNA reagent in 1X TE buffer. A composite sample for sequencing was created by combining equimolar ratios of amplicons from the individual samples and was cleaned using the MoBio UltraClean PCR clean up DNA purification kit. The final sample including aliquots of the sequencing primers were sequenced at the Biofrontiers Next-Gen Sequencing Facility (University of Colorado Boulder). Samples collected in the Eastern US (n = 55) were sequenced using the Illumina HiSeq platform in 2011. Samples collected in Colorado (n = 147) were sequenced using the Illumina MiSeq platform in 2013.

We tested for P. destructans DNA using real-time quantitative PCR following the protocol of Muller et al. (2013). We quantified P. destructans based on the cycle threshold (C_t) value to estimate the fungal load on each bat, with a C_t cut-off of 40 cycles. The standard curve for quantification was generated using genomic DNA from P. destructans ATCC MYA-4855 quantified with the Quant-IT PicoGreen double-stranded DNA assay kit (Life Technologies, Carlsbad, CA) in conjunction with a DynaQuant 300 fluorometer (Harvard Bioscience, Inc., Holliston, MA). Serial dilutions of the DNA from 10 ng to 1000 fg were prepared and analyzed with IGS qPCR, resulting in a significant curve from 17.33 to 30.74 C_t (C_t = −3.348^*(Log10 P. destructans [ng]) + 22.049).

Both sets of HiSeq (from the eastern US) and MiSeq (from Colorado) sequences were de-multiplexed using a custom python script “prep_fastq_for_uparse.py,” available at: https://github.com/leffj/helper-code-for-uparse). We trimmed the MiSeq forward reads to 100 bp in length (Prober et al., 2015) in order to combine them with the HiSeq data. Both HiSeq and MiSeq sequences were then combined into a single data file. The combined dataset was processed using the UPARSE pipeline (Edgar, 2013). After de-multiplexing, sequences were filtered against a maximum per-sequence expected error frequency value of 0.5 and singleton sequences were removed to filter out low quality sequences. The resulting sequences were clustered at the ≥97% similarity threshold to generate a de novo OTU database. The raw, combined sequences were then mapped to the de novo database to generate sequence counts per OTU and sample. OTUs were assigned taxonomic classifications using the RDP classifier (Wang et al., 2007) trained on the Greengenes database (August 2013 version; McDonald et al., 2012). All OTUs classified as mitochondria or chloroplasts were removed from analyses. Samples were then rarefied to 9800 sequences per sample, resulting in a total of 23,782 OTUs in the final data set.

Unless otherwise noted, bacterial sequence data were analyzed using R (Version 3.2.2, R Core Team, 2015). The combined data set was then used to generate three distinct data subsets: All Samples (202 bat skin and environment samples collected from a site where bats were sampled and soil samples that are not paired with bat samples), Bats Only (164 only bat skin samples, no environmental samples or soil samples), and Paired Samples (93 samples where bats and their environments (cave wall) could be sampled together). Within the Bats Only subset, a smaller subset of the 70 M. lucifugus samples was generated in order to better describe site and state variation within a single species.

We examined alpha diversity in the Bats Only and Paired Only data sets, respectively, for significant differences across sample type, site, species, and region. We used a custom script in the mctoolsr (https://github.com/leffj/mctoolsr), biom, and ggplot, and vegan libraries to calculate Shannon diversity, richness, and Kruskal-Wallis tests of significance. To determine which bacterial OTUs were shared across bat species, the Paired Only samples were used with the core_microbiome.py script in QIIME (http://qiime.org). In order to assess differences in community composition, we used a Bray-Curtis dissimilarity matrix for each of the 4 data subsets (All Samples, Bats Only, and Paired Only, and M. lucifugus Only). We used permutational multivariate analysis of variance (PERMANOVA) via the adonis function within the vegan package to examine the differences between communities using host species, site, and region as explanatory factors. To test for differences in dispersion amongst groups, we used betadispr in the vegan package. In addition, we used pairwise PERMANOVAs to examine which bat species, sites or regions were driving significant differences. We applied false discovery rate (FDR) corrected p-values to the pairwise PERMANOVA tests. The Bray-Curtis dissimilarity matrix was used to create ordinations using non-metric multidimensional scaling (NMDS) for each data set in order to visualize the effects of state, site, and host species on bacterial communities. We used a multiple linear effects model to determine which bacterial classes were driving differences in community composition among bat species, with site as a random effect. Relative abundances of bacterial classes were rank transformed to meet the assumptions of the models. All p-values were FDR corrected based on the number of taxa tested (all bacterial classes with mean relative abundances ≥1% in any species).

To compare the skin bacteria observed on bats to other vertebrates, we performed a literature search and identified 9 studies representing 4 taxa with the following criteria: (1) targeted 16S rRNA gene for bacteria, (2) comparable high-throughput sequencing techniques within the last 7 years, and (3) a measure of the bacterial abundances on a specific area of skin. We mined the 5–10 most abundant bacterial taxa per host species, and calculated percent relative abundance using available supplemental materials or graphs as provided by published papers. Remaining taxa, excluding the 5–10 most abundant bacterial taxa, were labeled as “Other Classes” or “Other Phyla,” respectively.

The clustering step (>97%) to a de novo database ultimately produced 28,487 OTUs from 235 samples. The average number of sequences per sample was 52,250. Samples were rarefied to 9000 sequences per sample, which removed 29 samples from the dataset that had fewer sequences than the threshold and four samples which were negative controls. In total, 202 bat and environment samples achieved high quality sequencing and were included in downstream analyses. (Supplementary Table 1). These samples were used in three data sets: Bats Only (166), Paired Only (93 samples, 19 environmental and 74 bats), M. lucifugus only (70) and All Samples (202).

The top classes represented in the data set for all samples were: Gammaproteobacteria, Alphaproteobacteria, Actinobacteria, Betaproteobacteria, Bacilli, Flavobacteria, Cytophagia, and Thermoleophilia. For both bat and environmental samples, the top 15 OTUs by relative abundance comprise 89% of the sequences. The top OTUs for just bats from all states fall within the classes: Gammaproteobacteria, Alphaproteobacteria, Actinobacteria, Betaproteobacteria, Bacilli, Flavobacteria, Saprospirae, Thermoleophilia, and Clostridia (Figure 1). For bats, there were more than 604 OTU genera represented (with an abundance of more than 0.0001%). Two common genera found on bats are Pseudomonas (9%), and Acinetobacter (5%), and there are many genus-level that remain unclassified at that level such as Xanthamonadaceae: Other (4%), Bacteria: Other (2%), Gammaproteobacteria: Other (2%), Bacilli: Other (2%), Actinomycetales: Other (1.8%), Enterobacteraceae: Other (1.8%), Sphingomonadaceae: Other (1.7%) and Pasteurellaceae: Other (1.5%).

Across multiple vertebrate taxa, bats share the most classes with other mammals, particularly dogs (Table 2). At the phylum level, Proteobacteria and Bacteroidetes are represented across all taxa, with particularly high abundances (11–61%) of Proteobacteria in all skin microbiomes. Bats and salamanders have skin microbiomes that are largely composed of Proteobacteria, at 65 and 60% relative abundance respectively. At the class level, Gammaproteobacteria are broadly represented across host taxa. Bats and dogs showed a stark difference in their relative abundance of class Gammaproteobacteria, which dominated the bat skin community at more than 40% but was 12.5% in domestic dogs (Canis familiaris). Alphaproteobacteria, Betaproteobacteria, and Actinobacteria are also shared across multiple taxa, however in lower abundances overall. One class, Thermoleophilia, was only found on furred mammals (dogs and bats) and Cytophagia was the only class found exclusively on bats. Overall, much of the skin microbial community remains unclassified, even at the class level (between 7 and 59% in the species compared in this data set).

In order to examine how diversity varied across species, sites, and regions, Shannon diversity indices were calculated for paired environment and bat samples as well as within bat samples only. For the paired samples across all sites, Shannon diversity of environmental samples was marginally different from the bat skin samples (p-value = 0.06), and there are observed differences in alpha diversity metrics between sites and regions. When comparing differences between bats at different sites, there were significant differences in alpha diversity observed among all sites and species in both Shannon diversity and richness (Shannon Diversity: Kruskal-Wallis chi-squared = 20.15, df = 10, p-value = 0.02, Richness: Kruskal-Wallis chi-squared = 24.70, df = 10, p-value = 0.005).

We measured patterns in community composition across sample regions, sample sites, and host species (beta diversity) in the Bats Only data set of 164 individuals found across 11 different sites in all 3 regions (NY, VA, and CO). Among the bat samples, seven OTUs were found across 85% of all samples, suggesting that these could be considered as core taxa. They were classified as belonging to the class and family, Gammaproteobacteria: Pasteurellaceae, Gammaproteobacteria: Enterobacteriaceae, Alphaproteobacteria: Sphingomonadaceae, and Betaproteobacteria: Burkholderiales. We used PERMANOVA analyses to assess differences in the bacterial community and find that region, site, and species all significantly explained variability (p-value < 0.001 for each variable, R² 0.12, 0.27, and 0.15). These differences are apparent when the NMDS plot is colored by species (Figure 2). Betadispr analysis of each of the factors (site, sample, and region) within the Bats Only data set showed that there are significantly different dispersion levels for each (each p-value < 0.001). To further explore these patterns, we used pairwise PERMANOVAs. Many of the comparisons between community composition were significant across data sets (i.e., Bats Only) and sample groups (i.e., site). (Supplementary Table 2). Using a multiple liner effects model, Betaproteobacteria were found to be the class that had the strongest effect on differences in community composition on bat skin since it was the only one that varied significantly across species (FDR corrected p-value = 0.009). The significant (p-value < 0.001) differences in beta diversity between sites for the Bats Only data set are shown in an NMDS plot in Figures 3A,B, with an accompanying dendrogram showing the geographic pattern of beta diversity across sites. Lastly, the significant (p-value < 0.001) difference in beta-diversity between regions is graphed using an NMDS of the paired samples analysis, with each region highlighted in a different color and each sample type (bat skin vs. environment) represented by a shape (circle or triangle, Figure 4). These results demonstrate a high amount of exchange from the environment (at both the site and regional level) as well as the bat host species in the bacterial communities of bat skin. However, when considering the overall dissimilarities between skin bacterial communities between and within sample site beta diversity and species beta diversity, there is more variation between sites than between species (Figure 6). This suggests that while environmental and host factors are both important in determining the composition of the skin microbiome, differences between sites more strongly predict differences in bacterial communities than differences across host species.

To better understand the differences between the effects of species on the host bat community and the effect of environment (on both the regional and site-specific level), we used a subset of the data set that included a single species. Samples from 70 individuals of M. lucifugus from all three states and multiple sites in each state were used for further analysis. Using a Bray-Curtis distance matrix, patterns by site and region were visualized using an NMDS ordination and the dissimilarity patterns between host bacterial communities were analyzed using PERMANOVA (Figure 5A). Site and region were both significant factors for the M. lucifugus subset, suggesting that there may be site- and region-specific microbes that influence the development of the microbiome shared within a species (R²: 0.33, p-value < 0.001). The NMDS plot, colored by site and state, shows clear differences in the bacterial communities between different sites where bats were sampled. A dendrogram of the beta diversity dissimilarity shows that the most similar communities are grouped by state (Figures 5A,B).

A subsample of the data set that included paired environment (cave wall or soil) samples as well as bat samples (skin swabs) for sample site were analyzed to see which taxa were shared amongst the two groups. Soil and bats did not differ significantly in their richness or Shannon diversity (Kruskal-Wallis test of significance) (Figures 7B,C). Based on these results, a heatmap was generated using relative abundances to compare the top 10 classes found on both bats and environment, which shows that Gammaproteobacteria, Alphaproteobacteria, and Actinobacteria are the most abundant and are also shared within both soil and bats (Figure 8). However, a PERMANOVA analysis of the difference in community composition of the microbial communities between bat skin and the environment using a Bray-Curtis dissimilarity matrix found a significant difference (Figure 7A) (R² = 0.033, p-value < 0.001), suggesting that there are differences below the class level that distinguish soil and bat skin.

WNS testing yielded five positive individuals of the 68 bats sampled from the eastern US portion of the data set, although only four of these individuals were used in the final analysis. WNS status was not a significant factor in the community composition of bat skin community when comparing either the entire sample set or when subsampling only from within sites with M. lucifugus (ANOSIM using QIIME, p > 0.05). Supplementary Figure 1 shows PCOA plots of infected and uninfected individuals across all species and within just M. lucifugus. There are no strong patterns by infection status.

The bat skin microbiome is influenced by both species and the habitat in which the host is found. Many of the same bacterial genera are found in similar abundances across bat species (Figure 1). Many of the bacterial taxa found on bat skin are also dominant taxa in soil and plant material, suggesting that these groups are shared readily between the host and its environment (Figure 8). The community of bacteria found on bat skin contains high diversity, with even the most abundant bacteria comprising less than 10% of the overall sample. Pseudomonas (9% of total abundance in the paired bat and environmental data set) is a commonly found genus in the environment, along with Acinetobacter (5%) and the plant-associated family Xanthamonadaceae (4%). As bats are the most vagile mammal species, it is possible that during activities such as foraging and migration they are acquiring a diverse microbial community from the environment. Two other groups of interest in this data set that are known to be host-associated. Enterobacteriaceae (1.8%) is an environmental as well as gut-associated bacterial family and suggests that there may be contamination of the skin with fecal material, which is not surprising given the close association with the roost environment. Previous culture-based studies have also detected Pseudomonas on bat skin, and found multiple strains in the Pseudomonas fluorescens complex which exhibited anti-fungal properties against P. destructans (Hoyt et al., 2015). Pasteurellaceae (1.5%) is usually associated with the oral microbiome, suggesting that these microbes may be deposited on the skin during grooming bouts. Both Enterobacteriaceae and Pasteurellaceae were found to be shared between 85% of all the bat species as part of the core microbiome analysis. Overall, the knowledge of how many of these taxa are specialized to the bat host remains limited by the availability of data for host-associated bacterial genera (Table 2). However, the general pattern emerges that much of the diversity of bacterial taxa is shared between both bats and the environment, and many of those taxa are known to be either common symbionts of other mammalian hosts or found in the environment.

In a review of the available literature on the skin microbial communities of four mammals (humans, dogs, whales and bats), two birds (turkey vultures and black vultures), two amphibians (wild frogs and red-backed salamanders) and marine fishes, bats share some of the key skin associations with other terrestrial mammals (dogs) but have few bacterial taxa in common across these groups. The important classes shared between dogs and bats include Alphaproteobacteria, Betaproteobacteria, Gammaproteobacteria, Actinobacteria, Bacilli, and Thermoleophilia (Table 2). However, bats and dogs showed a difference in their relative composition of class Gammaproteobacteria, which dominated the bat skin community at more than 40% but was only 12.25% in dogs. Gammaproteobacteria was found in moderately high abundance in whales (33.9%) and 50% in salamanders, and was shared in varying amounts amongst all the wild vertebrates sampled. One class, Thermoleophilia, is only found on bats and dogs, suggesting a relationship to terrestrial mammals, however it is a known soil organism (Suzuki and Whitman, 2012; Crits-Christoph et al., 2013). For all of the vertebrates examined for this study, 7.6%–49.4% of bacteria were not in the most abundant taxa. Although broad comparisons could be made, synthesis across studies was complicated by disparate sampling schemes, analysis types, and OTU reporting methods. For comparative analyses in the future and the ability to understand shared microbial taxa between study groups, it is essential that future studies provide taxa tables to aid in comparative analyses to allow an understanding of broader patterns across taxa.

Many factors were found to be significant drivers of the host-associated microbial communities of bats. There were significant differences in the community composition of the bacterial communities between the three regions (Colorado, New York, and Virginia) (Figure 4). While there are shared taxa among these three areas, New York and Virginia bats appear to share more in common with each other than Colorado, suggesting a regional pattern of soil- and host- associated microbial communities, though broader regional sampling would be needed to confirm this conclusion. Using the data set containing only bat samples, there are significant differences amongst bacterial OTU abundances (alpha diversity) and between these microbial communities (beta diversity) by both species and site (Figures 1, 3). Patterns of host-associated bacterial communities within a single species (specifically M. lucifugus) also show differences between regions and sites. Taken together, this represents a complicated picture of how bats acquire and maintain their skin microbial community. Along with the results of the paired analysis, which suggest that alpha diversity is not different between bats and their environment, this suggests a shared relationship both among bats (inter- and intra-species) and between the bat hosts and their surrounding environment.

The results from the bat microbial community suggest that many of the bacteria that are dominant in the environment of hibernacula are also shared with bats, which may provide one mechanism for why we see general regionalized patterns as well as site-specific differences as the strongest factors explaining variation in the bat skin microbiome (Figure 8). However, a significant difference in the bacterial community composition amongst the 14 bat species sampled does exist, though less pronounced than the site effects. In amphibians, species tends to be a better predictor of the host-associated skin community rather than site (McKenzie et al., 2012; Kueneman et al., 2014). For a single whale species (Megaptera novaeangliae), sampling location is known to be an important factor shaping the host skin community (Apprill et al., 2014). In humans and domestic dogs, we know that the composition of the skin microbial communities can be affected by many factors, such as disease status, temporal variation, and body site location (Costello et al., 2009; Grice et al., 2009; Grice and Segre, 2011; Hoffmann et al., 2014). In homes where domestic dogs and humans share the same environment, many of the bacterial taxa will swap hosts to form a mixed microbial community that is significantly different from other humans or domestic dogs not within the same household (Song et al., 2013). In a recent study of carrion-eating birds, many of the taxa found on the exposed dermis on the heads of two species of vultures closely resembled prey items, as well as taxa that were specific to bird skin (Roggenbuck et al., 2014). The results from the bat microbial community suggest that many of the bacteria that are dominant in the environment of hibernacula are also shared with bats, which may provide one mechanism for why we see general regionalized patterns as well as site-specific differences as the strongest factors explaining variation in the bat skin microbiome (Figure 8). However, host species effects are also present, suggesting that the bat skin microbiome is influenced by complex and interacting factors. Broader studies that examine whether there are generalizable patterns that drive vertebrate skin microbiome communities are needed.

For this study, all of the bats sampled were captured at the end of the summer, during the fall swarming period, or at the very start of hibernation. The pathogen only begins to invade and colonize the dermis of the host during hibernation (Langwig et al., 2015) when the temperature of the skin is lowered for sustained periods. None of the bats included in this study were visibly infected, and P. destructans loads were low in positive individuals (Supplementary Figure 1). Our dataset does not, therefore, permit a robust investigation of the relationship between P. destructans and the bat skin bacterial community.

Bats share much of their abundant microbial taxa with the environment, and this may drive the patterns observed in the alpha- and beta-diversity differences between sampling sites. There is more variation observed between site beta-diversity than between species (Figure 6). To a lesser extent, but still significant, species play an important role in determining the composition of the bat skin microbiome. Overall, the complex factors regulating the bat microbiome suggest that the environment and host factors are important, but more research is necessary to understand the relative contributions and functions of these host-associated microbial communities. We are just beginning to understand the patterns of microbial diversity on bat hosts, and understanding the ecology of these associations will build a foundation for future work on the influences on host health and interactions with the environment.

Overall, our understanding of host-associated skin microbial communities is still in its infancy. Compared to the amount of research available on the human microbiome, the microbiome of animals, and particularly skin, remains relatively underexplored, and many of the patterns that we observe need to be rigorously evaluated. The skin microbial community is not as abundant nor as diverse as the gut community, and therefore requires particular care when collecting and processing samples from the field. The emergence of new fungal diseases, several of which invade the skin of their host and cause morbidity and mortality, highlights the importance of fungi for understanding microbiomes (Fisher et al., 2012). Understanding the patterns by which these communities assemble can shed light on how to approach these infectious diseases from a microbiome perspective and possibly enhance treatment of these diseases.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.