Antifungal Potential of the Skin Microbiota of Hibernating Big Brown Bats (Eptesicus fuscus) Infected With the Causal Agent of White-Nose Syndrome

Little is known about skin microbiota in the context of the disease white-nose syndrome (WNS), caused by the fungus Pseudogymnoascus destructans (Pd), that has caused enormous declines of hibernating North American bats over the past decade. Interestingly, some hibernating species, such as the big brown bat (Eptesicus fuscus), appear resistant to the disease and their skin microbiota could play a role. However, a comprehensive analysis of the skin microbiota of E. fuscus in the context of Pd has not been done. In January 2017, we captured hibernating E. fuscus, sampled their skin microbiota, and inoculated them with Pd or sham inoculum. We allowed the bats to hibernate in the lab under controlled conditions for 11 weeks and then sampled their skin microbiota to test the following hypotheses: (1) Pd infection would not disrupt the skin microbiota of Pd-resistant E. fuscus; and (2) microbial taxa with antifungal properties would be abundant both before and after inoculation with Pd. Using high-throughput 16S rRNA gene sequencing, we discovered that beta diversity of Pd-inoculated bats changed more over time than that of sham-inoculated bats. Still, the most abundant taxa in the community were stable throughout the experiment. Among the most abundant taxa, Pseudomonas and Rhodococcus are known for antifungal potential against Pd and other fungi. Thus, in contrast to hypothesis 1, Pd infection destabilized the skin microbiota but consistent with hypothesis 2, bacteria with known antifungal properties remained abundant and stable on the skin. This study is the first to provide a comprehensive survey of skin microbiota of E. fuscus, suggesting potential associations between the bat skin microbiota and resistance to the Pd infection and WNS. These results set the stage for future studies to characterize microbiota gene expression, better understand mechanisms of resistance to WNS, and help develop conservation strategies.


INTRODUCTION
The skin is the first physical and immunological barrier against invading pathogens. It is also a complex and dynamic ecosystem inhabited by a rich community of microorganisms composed of bacteria, archaea, fungi, and viruses. This community contributes to host defenses by limiting colonization and persistence of pathogens that compete for resources, and by producing pathogen inhibitors (Grice and Segre, 2011;Walter et al., 2018;Woo et al., 2019). Here, we use the term 'microbiota' to refer to the taxonomic diversity of Bacteria and Archaea assessed using marker genes, rather than 'microbiome, ' which refers to both taxonomic and functional diversity of the complete community (Marchesi and Ravel, 2015). The microbiota also interact with the innate and adaptive immune system and contribute to maintenance of skin integrity and tissue repair (Lai et al., 2009;Curtis and Sperandio, 2011;Naik et al., 2012). Despite this importance, the complex interactions among the skin microbiota, skin invading pathogens and hosts remain poorly documented in animals, especially for newly emerged wildlife diseases causing population declines around the globe (Daszak, 2000;Belden and Harris, 2007;Jones et al., 2008;Fisher et al., 2012). Introduced fungal pathogens such as Batrachochytrium dendrobatidis and Batrachochytrium salamandrivorans have caused global declines and even extinctions in amphibians (Longcore et al., 1999;Martel et al., 2013). In recent decades, studies have highlighted the potential role of skin microbiota in patterns of resistance and susceptibility to these fungal pathogens and have pointed to potential management practices that might help conserve host populations (Harris et al., 2009;Bletz et al., 2013Bletz et al., , 2017Jani and Briggs, 2014;Woodhams et al., 2014;Bataille et al., 2016;Ange-Stark et al., 2019;Rebollar et al., 2019). However, despite this potential, management actions based on the skin microbiota have still not been widely applied in response to wildlife disease.
Similar to amphibians, hibernating North American bats have suffered dramatic impacts from a newly introduced fungal disease, white-nose syndrome (WNS) Frick et al., 2016). WNS is caused by the fungus Pseudogymnoascus destructans (Pd) (Gargas et al., 2009;Lorch et al., 2011) and has caused widespread declines of hibernating bats since being introduced from Eurasia sometime around 2006 (Warnecke et al., 2012;Frick et al., 2015Frick et al., , 2017Hoyt et al., 2016;Drees et al., 2017;Trivedi et al., 2017). In vulnerable bat species, Pd invades the skin creating lesions that alter fluid balance, thermoregulation, and gas exchange Cryan et al., 2010;Warnecke et al., 2012;Verant et al., 2014;McGuire et al., 2017). Hibernating bats survive the winter using stored body fat and prolonged energy-saving bouts of torpor characterized by dramatically reduced body temperature and metabolism (Jonasson and Willis, 2012;Czenze et al., 2013;Czenze and Willis, 2015). The fungal infection causes bats to warm up too frequently which, in turn, depletes their fat reserves (Reeder et al., 2012;Warnecke et al., 2012). Understanding how bat species respond to this disease is important for designing management strategies to protect bat populations both for biodiversity conservation, and preservation of the ecosystem services that bats provide Maine and Boyles, 2015;Wray et al., 2018). There is enormous variation in the impacts of WNS for different bat species, with some exhibiting little to no mortality, and others facing local extinction (Frick et al., 2010Langwig et al., 2015Langwig et al., , 2017. Mechanisms underlying this variation in susceptibility are not fully understood, despite potential benefits for disease management. Variation in susceptibility to Pd could reflect mechanisms that allow for resistance to infection (i.e., reduction or elimination of pathogen infection) and/or tolerance of infection (i.e., reduction of the harm caused by infection) (Råberg et al., 2007(Råberg et al., , 2009Svensson and Råberg, 2010) and potential mechanisms of resistance and tolerance are beginning to be addressed in the WNS context (Frick et al., 2016;Langwig et al., 2017;Cheng et al., 2019).
Big brown bats (Eptesicus fuscus) exhibit evidence of resistance with mild WNS symptoms compared to more susceptible species (Frank et al., 2014;Moore et al., 2018). E. fuscus are of particular interest as they often hibernate under similar environmental conditions as species that are highly vulnerable to WNS, including the little brown bat (Myotis lucifugus), the tricolored bat (Perimyotis subflavus), and the northern long-eared bat (Myotis septentrionalis), all three of which are now listed as endangered in Canada due to WNS (Canadian Wildlife Service and Committee on the Status of Endangered Wildlife in Canada, 2013). Physiological and behavioral factors could play a role in this resistance (Willis and Wilcox, 2014;Langwig et al., 2015;Frick et al., 2016) and understanding these mechanisms could help fulfill an important knowledge gap by identifying traits that contribute to WNS survival.
One underexplored factor that could affect resistance to WNS is the skin microbiota. Differences in lipid profiles affecting WNSresistance (Frank et al., 2016(Frank et al., , 2018) may contribute to different microbial profiles by providing different nutritional substrates. Hundreds of microorganisms isolated from wild bats and their corresponding habitats have been tested in controlled laboratory conditions for their inhibitory effects on Pd, by focusing on the actions of secreted compounds, contact inhibition or volatile molecules (Hamm et al., 2017;Micalizzi et al., 2017). Little brown bats persisting with Pd have proportionally more abundant Rhodococcus and Pseudomonas in their skin microbiota (Lemieux-Labonté et al., 2017), suggesting a possible role for these microorganisms in bat survival. Moreover, antifungal strains of Pseudomonas isolated from the skin of E. fuscus inhibit Pd growth in vitro Hamm et al., 2017), and improve survival of WNS-infected little brown bats when applied as a probiotic treatment in vivo Hoyt et al., 2019). All of these results are promising but they do not provide a complete picture of the complex skin community and its potential role in Pd resistance of some species.
We explored the skin microbiota of WNS-resistant E. fuscus inoculated with Pd. This species can hibernate under the same environmental conditions as more susceptible species like M. lucifugus but survives Pd infection with limited skin colonization (Walke et al., 2015;Moore et al., 2018). Therefore, E. fuscus represents a good model to study potential resistance mechanisms resulting from the skin microbiota in a controlled environment that accounts for variation in hibernation conditions. We tested two hypotheses about the skin microbiota of E. fuscus in the context of Pd infection: (1) that Pd inoculation is not a strong selective force on the skin microbiota of this species (Ange-Stark et al., 2019) and would therefore not disrupt the microbial community, and (2) that microbial taxa with antifungal properties, which are common in the hibernation sites of bats persisting after WNS (e.g., Rhodococcus and Pseudomonas, Lemieux-Labonté et al., 2017), would be abundant both before and after infection with Pd in this resistant species.

Collection of Bats
On January 18, 2017, we visited Richard Lake Mine, a hibernaculum approximately 100 km east of Kenora, ON, Canada (49 • 45 N, −94 • 28 W), which houses several hundred M. lucifugus and E. fuscus each winter. We first swabbed eight E. fuscus in the cave to test whether transport to the laboratory affected the skin microbiota (Supplementary Files 1, 2). We then collected 32 adult E. fuscus (16 males and 16 females), suspended them in cloth bags in a cooler lined with wet towels to maintain hibernation conditions, and transported them by car to a bio-secure animal facility at the University of Winnipeg. All handling of bats in the lab occurred in a biosafety cabinet. We swabbed the left wing of each torpid bat immediately after arrival in the lab to sample the "pre-captivity" hibernation microbiota. We swabbed the dorsal surface of the left wing (forearm and under forearm) in linear strokes for 20 s with a sterile Whatman Omniswab (Fisher Scientific) soaked in sterile 0.15 M NaCl (Lemieux-Labonté et al., 2017). Swabs tips were ejected into MoBio Powersoil DNA isolation Kit tubes (Mo Bio Laboratories), which were then transferred to 4 • C within 2 h and to −20 • C within 12 h of sampling. As a negative control, a humidified sterile swab was exposed to open air for 20 s, prior to ejecting its tip into a MoBio tube. We analyzed the skin microbiota between the inoculated and control groups at the beginning of the experiment to detect any pre-existing differences occurring by chance. After collecting the microbiota swab we then swabbed the right wing to determine Pd status.
All individuals were then randomly assigned to one of 4 cages with 8 individuals/cage in a 1:1 sex ratio (i.e., 4 males and 4 females/cage). We replicated the experiment in two incubators with one cage of Pd-inoculated bats, and one cage of shaminoculated controls per incubator. Incubators were maintained at 8 • C and 98% relative humidity. Pd-inoculated bats received a dose of Pd that has repeatedly been shown to cause symptoms of WNS in wild bats: 20 µl of inoculum, containing 500,000 conidia suspended in phosphate buffered saline (PBS) with Tween 20 to prevent clumping, pipetted onto the ventral side of their wings (Lorch et al., 2011;Warnecke et al., 2012;McGuire et al., 2016). Healthy controls were inoculated only with PBS-Tween 20 (Warnecke et al., 2012;McGuire et al., 2016).
After 11 weeks of captive hibernation, bats were swabbed again, as described above, to assess the "post-captivity" skin microbiota (left wing) and Pd load (right wing). We recorded body temperature immediately after removing bats from the incubator at the end of the study using a digital thermometer accurate to ±0.1 • C (SPER Scientific, Model 80008, Arizona, United States) and inserting a 1 mm diameter l-type thermocouple probe 4 mm into the rectum (Supplementary File 3). Bats were then swabbed within ∼1-2 min of body temperature measurement. Only a maximum of 2-3 min would have elapsed from the time bats started rewarming in response to our disturbance until the time we swabbed them. After obtaining additional measurements and samples for a range of complementary studies, bats were humanely euthanized by CO 2 inhalation under isoflurane anesthesia.
Pd swabs were processed and analyzed at the Pathogen and Microbiome Institute at the Northern Arizona University, whereas skin microbiota samples were shipped to Université de Montréal (Québec, Canada) for further processing. Areas of the wing that are infected and colonized by Pd hyphae fluorescence orange under UV lightning (Turner et al., 2014;Cheng et al., 2017). We took UV photographs of every bat's wings using a digital camera (Olympus © Tough TG-830) using a Canadian dime as a sizing reference, and then digitally quantified these areas as a measure of Pd load and infection intensity using the Image-J © software. All methods were approved by the University of Winnipeg Animal Care Committee (Protocol Number AEO08399).

DNA Extraction, Amplification and Sequencing
Bacterial genomic DNA was extracted from each swab using the MoBio Powersoil DNA isolation kit. We followed the manufacturer's protocol, modified following Castelino et al. (2017), by adding a 15-min incubation period at 70 • C after the addition of buffer C1 to increase the efficiency of microbial cell lysis. All procedures were conducted in a laminar flow hood to limit potential sample contamination, and extractions were randomized to avoid detecting false patterns (Salter et al., 2014). Four extraction blanks, two amplification blanks, and the HM-782D Human Microbiome Project mock community (BEI Resources) were also included to detect possible contamination and assess sequencing accuracy (Salter et al., 2014;Glassing et al., 2016). Because of contamination and low input DNA, 23 preand post-captivity E. fuscus skin microbiota samples, and five sampling negative controls were retained for subsequent analysis. Amplification and sequencing were performed as previously described (Preheim et al., 2013a;Lemieux-Labonté et al., 2017). Briefly, libraries were prepared using a two-step PCR. The first PCR amplified the hypervariable region V4 of the 16S small subunit ribosomal gene with forward primer U515 F and reverse primer E786 R (Caporaso et al., 2011). The amplifications were performed with a Mastercycler Nexus GSX1 (Eppendorf). Each sample was amplified in quadruplicate and pooled to limit possible PCR artifacts. The second PCR step consisted of adding primers containing a barcode (index) and Illumina adapter sequences to each DNA amplicon, and used forward primer PE-III-PCR-F and reverse primer PE-III-PCR-001-096 (Preheim et al., 2013b). This second amplification was performed in triplicate. After each PCR step, samples were pooled and purified with the PCR purification Agencourt AMPure XP (Beckman Coulter). Indexed samples concentration was measured with Qubit 2.0 Fluorometer (Invitrogen), and samples were pooled to obtain a final concentration range between 10 and 20 ng/µl. DNA was next diluted and denatured according to the manufacturer's protocol for paired-end sequencing using MiSeq Reagent Kit v2 (500 cycles) 2 × 250 bp on MiSeq (Illumina).

Data Analysis
For all sequences, quality filtering, trimming, dereplication, sample inference and merging of paired-end sequences, were performed in DADA2 version 1.14.1 (Callahan et al., 2016). The corresponding amplicon sequence variants (ASV) table providing ASV abundances was assigned in DADA2 using SILVA database release 132 (Quast et al., 2012).
We amplified 8,282,605 sequences classified into 43,436 ASV from the 67 samples sequenced. We filtered out mitochondrial and chloroplastic DNA sequences, as well as sequences from the genera Halomonas and Shewanella, the two most abundant taxa in negative controls. ASVs with abundance values smaller than 10 were filtered out leaving 5,176,773 sequences in 17,162 ASVs. After these filtration steps, only a small number of sequences (less than 11,000) were identified in the sampling control, extraction and PCR negative controls, and these samples were also excluded. One sampling negative control had 47,608 sequences, but it was excluded from our analysis because it diverged dramatically in composition from bat samples collected (see Supplementary File 4). A total of 5,098,811 sequences, classified in 17,162 ASVs, were amplified and analyzed from the 46 bat samples, with a mean of 110,843 sequences per sample (range: 34,416-300,474). We were able to match all expected genera in the mock positive control (see Supplementary File 5). The genera of the 20 expected mock taxa were the most abundant in the mock profile (see Supplementary File 5). We detected 733 false positives, but with very low relative abundance values ( <0.1%). All analyses were conducted in R version 3.5.0 (R Core Team, 2018).

Alpha Diversity
We first quantified diversity of the skin microbial community for all samples based on alpha diversity using the Shannon index (Shannon, 1948). The Shannon index, which includes both ASVs richness and evenness, was selected due to its reduced sensitivity to sample depth differences (Haegeman et al., 2013;Preheim et al., 2013a) (see Supplementary File 6). Normally distributed alpha diversity values were compared using linear models and linear mixed-effect models [lm() and lme() function] in R with Shannon index as the response variable, Pd-inoculation versus control (hereafter 'inoculation') and pre-versus post-captivity (hereafter 'captivity') as fixed effects and cage ID as a random effect. Significance was tested using ANOVA for linear model and a likelihood ratio test with a chi-square distribution for linear mixed-effect models (Pinheiro et al., 2017).

Beta Diversity
We used two distance measures to account for the phylogeny of microbiota in our samples (i.e., unweighted UniFrac and weighted UniFrac) (Lozupone and Knight, 2005;Lozupone et al., 2007). Distances were computed on rarefied data, as such measures could be sensitive to differences in sequencing depth Weiss et al., 2017). Computations were performed with the phyloseq package (McMurdie and Holmes, 2013) and results were visualized with principal coordinates analysis (PCoA) (Gower, 1966) using the ordinate() function in R. Distance matrices were checked with the function is.euclid() of the ade4 package (Dray and Dufour, 2007) prior to the ordination to ensure that all distances were Euclidian and properly representable by PCoA (Gower and Legendre, 1986). When required, square-root transformations were applied to obtain distance matrices satisfying the Euclidian condition (e.g., Weighted UniFrac). All phylogeny-based UniFrac distances were calculated using a phylogenetic tree constructed with FastTree 2.1.8 (Price et al., 2010).
To test for effects of inoculation and captivity we used distance-based redundancy analysis (db-RDA) (Legendre and Anderson, 1999). It is computed by first decomposing UniFrac distances (weighted or unweighted) into principal coordinates, and then applying RDA to the corresponding principal coordinates using the capscale() function of the R package vegan (Oksanen et al., 2019). We computed partial db-RDA to better understand the influence of Pd-inoculation and captivity on variation in microbial assemblages while controlling for possible confounding factors (i.e., incubator, sex, and cage ID) (Davies and Tso, 1982). Adjusted R-squared values (Ezekiel, 1930) were calculated to compare the explanatory power of different models. Significance of db-RDA and partial db-RDA was tested via 9999 permutations with the anova.cca() function of the R package vegan. We then performed an analysis of multivariate homogeneity (PERMDISP) (Anderson, 2006) with the betadisper() function to test whether groups differed in their dispersion. The null hypothesis of this test is that the average within-group dispersion is identical in all groups (Anderson, 2001). For these two tests, the number of permutations was set to 9999. To test for an effect of inoculation on any change in the microbiota from pre-to post-captivity, we calculated distances between paired pre-and post-captivity samples for each individual, and then used linear mixed effects models to test for the fixed effect of inoculation on this distance while controlling for cage ID as a random effect (Pinheiro et al., 2017). For all analyses, a p-value threshold of 0.05 was considered significant.

Analysis of Skin Microbiota Composition
We assessed the effect of inoculation and captivity on the composition of microbial taxa down to the genus level using the Analysis of Composition of Microbiomes (ANCOM 2.0) (Mandal et al., 2015). ANCOM is based on non-parametric tests (i.e., either Kruskal-Wallis test for independent samples, or Friedman test for dependent samples) and is appropriate for compositional data (Gloor et al., 2017). The test relies on point estimates of data transformed by an additive log ratio, where presumed invariant taxa are selected as the denominator. To identify differences between samples, the analysis was performed using an unrarified ASVs table only including taxa with a relative abundance larger than or equal to 0.1%. The relative composition of skin microbiota was assessed using the same ASVs table.

Pd Status and Bat Data
We selected Richard Lake Mine for our study, in part because, based on our surveillance in winter 2016 (swabs from bats and substrates) and fall 2016 (swabs from swarming bats), it was negative for Pd the year before this study. It was also >350 km from the nearest known WNS-positive hibernaculum at the time of our study (U.S. Geological Survey, 2019) and we observed no signs of WNS at the time of capture. We swabbed bats to confirm Pd-negative status as soon as we returned to the lab after capture but it was not possible to wait for result for these swabs before assigning bats to experimental groups. Unfortunately, after we began the experiment, qPCR results revealed that eight of the 23 study animals selected for microbiota analysis were just above the threshold (40 C t ) to be considered Pd positive (McGuire et al., 2016) at the time of collection (39.3 ± 0.74 C t , mean ± SD; Table 1 and Supplementary File 7). Half of these eight bats had been randomly assigned to the inoculated group and half to the control group. Among the bats assigned to the control group, three Pd-negative bats remained negative throughout the experiment but four bats that started in the control group as negative were positive by the end of the experiment presumably because of transfer of Pd from naturally infected control bats. One control bat that was Pd-positive at the time of capture was Pd-negative by the end of the experiment ( Table 1). We do not think this infection of captured bats, or contamination of the control group, influenced our results because loads at capture were extremely low, near the limit of detection and, on average, 18-fold lower than loads for the bats we experimentally infected (see Supplementary File 7). Moreover, proportions of the wings exhibiting orange UV fluorescence (see Supplementary File 7) of contaminated control bats at the end of the experiment (0.2 ± 0.4%) were 14-fold lower than those of inoculated bats (2.8 ± 3.3%) (Wilcox test, W = 6, p < 0.001) and Pd loads of inoculated bats (0.0003 ± 0.0003 ng) were six-fold that of the contaminated control bats at the end of the experiment (0.00005 ± 0.0001 ng) (Wilcox test, W = 24, p = 0.01). All bats that we inoculated with Pd in the laboratory, and used for subsequent analysis of the microbiota (n = 12), were Pd-positive at the end of the experiment.
After experiment, we detected orange fluorescence typical of Pd infection on E. fuscus skin of all bats we inoculated and on four control bats (Table 1). Interestingly, green and/or blue fluorescence was also detected. Green was always detected with co-occurrence of orange and/or blue fluorescence, whereas blue fluorescence was detected in the presence of orange, or alone.

Alpha Diversity
After assuring no difference in Shannon diversity between inoculation groups prior to the experiment (F = 0.62, p = 0.44), no effect of inoculation was detected on Shannon diversity at the end of the experiment while controlling for cage ID (χ2 = 0.05, p = 0.82). There was also no effect of inoculation on the change in alpha diversity from capture to the end of the experiment (F = 0.74, p = 0.40). Although there was no effect of inoculation, we did detect a significant difference in alpha diversity between the start and the end of the experiment (F = 93.15, p < 0.0001). Mean Shannon diversity values decreased by 24% (1.59, ± = 0.16) in post-inoculated samples indicating a strong effect of either hibernation or captivity (or both) on alpha diversity of the skin microbiota.

Beta Diversity
By chance, the inoculation group factor explained redundant variation in weighted UniFrac distances when we controlled for sex pre-captivity ( Table 2). However, this redundant variation was not persistent throughout the experiment. Indeed, there was no redundant variation between inoculation and UniFrac distances by the end of the experiment (i.e., for post-captivity samples) after controlling for incubator, sex and cage ID ( Table 2). This indicates that beta diversity changed over the course of the experiment, although it was not caused by inoculation. Yet, cage ID was the factor that explained most of the variation in beta diversity when included in the model post-captivity ( To account for temporal effects and better assess whether inoculation led to more or less change in the microbiota over the course of the experiment, we calculated distances between preand post-captivity samples paired within each individual. After controlling for cage ID, we found significantly higher dispersion of distances for Pd-inoculated bats based on both unweighted and weighted UniFrac distances (χ 2 = 6.52, p = 0.01, χ 2 = 8.92, p = 0.003) (Figure 1). Post-captivity samples from inoculated bats were 11% (0.09 ± 0.03) more distant than samples from controls based on unweighted UniFrac distances, and 22% (0.10 ± 0.03) more distant based on weighted UniFrac (Figure 1). In other words, as for pooled results above, within individuals the skin microbiota of Pd-inoculated bats changed more over time than that of sham-inoculated controls. Captivity was an important factor shaping community structure of the skin microbiota (Figure 2). About 26% of the variation in community structure was explained by captivity based on unweighted UniFrac distances, and about 42% of this variation was explained by captivity based on weighted UniFrac ( Table 2). Community structure of pre-captivity samples was more dispersed than that in post-inoculated samples based on unweighted UniFrac distances (PERMDISP, F = 12.34, p = 0.0008) and weighted UniFrac (PERMDISP, F = 30.33, p = 0.0001) (Figure 2).

Analysis of Skin Microbiota Composition
Only one taxon (i.e., Pseudonocardia) changed in abundance over the course of captivity at the genus level and it was more abundant in pre-captive individuals ( Table 3). No taxa were significantly more abundant in Pd-inoculated individuals at the end of the experiment ( Table 3).

DISCUSSION
Research on the skin microbiota of wildlife facing emergent skin diseases is a growing area of research (Daszak, 2000;Belden and Harris, 2007;Jones et al., 2008;Fisher et al., 2012), and in the specific context of North American bats, understanding responses of the skin microbiota could contribute to management of WNS Hoyt et al., 2019). We inoculated, or sham inoculated E. fuscus with Pd and maintained them under controlled conditions for 3 months to test whether the skin microbiota of a WNS-resistant species might evolve in response to fungal infection, and help understand the potential contribution of the microbiota to WNS resistance.
We found mixed support for our first hypothesis that Pd inoculation would have little impact on the skin microbiota of this WNS-resistant species. No taxa were significantly more abundant on Pd-inoculated bats compared to controls, at the end of the experiment, and no difference was detected in diversity as observed in the wild for this species (Ange-Stark et al., 2019). Interestingly, the observed differences in beta diversity (weighted UniFrac distances) between bats assigned to our two experimental groups disappeared by the end of the experiment. This suggests that inoculation has no effect on the skin microbiota diversity. However, our results indicate that Pd inoculation reduced the capacity of the host (or its microbiota) to regulate the skin microbiota, leading to more dispersion within and between individuals as observed in studies of Bd, the fungal pathogen affecting amphibians (Zaneveld et al., 2017). Jani and Briggs (2014) found that outbreaks of Bd increased temporal changes in the skin microbiota of Rana sierra and suggested that, even frogs that could tolerate Bd infection could be sensitive to disruption of the microbiota by the pathogen, although whether this disruption was protective or harmful was not clear. Similarly, our results reflect a transformation of the skin microbiota due to Pd infection, although we do not know if this disruption would be advantageous, detrimental, or neutral for infected E. fuscus. In Lithobates catesbieanus, an amphibian species with low susceptibility to Bd (Eskew et al., 2015), the fungus was a selective force on the microbial community, but the microbiota also affected Bd infectivity and, consequently, probably exhibited a negative effect on host fitness (Walke et al., 2015). Our results could not establish whether the skin microbiota affects the fungal infection, or if the fungal infection could in turn affect the microbiota, and this still needs to be investigated.
Consistent with our second hypothesis, the composition of the microbiota included several taxa known to inhibit Pd growth in vitro, and which have been found on bats persisting after WNS [e.g., Rhodococcus (Cornelison et al., 2014;Hamm et al., 2017), Pseudomonas ]. Pseudonocardia was also the most abundant microbial taxa in pre-captivity samples and, although it declined in abundance by the end of the experiment, it is known to have antifungal activity (Sen et al., 2009). The proportional abundances of Rhodococcus and Pseudomonas were stable from pre-to post-captivity and they were among the most abundant taxa at both time points. These bacteria may play a role in the persistence of some individuals of highly susceptible bat species (e.g., M. lucifugus) after Pd invasion (Lemieux-Labonté et al., 2017) and the fact that they represented a large proportion of the microbiota in the resistant species we studied also suggests their importance in the response of bats to WNS. Interestingly, the blue and green fluorescence we detected on E. fuscus skin, is typical of pyocyanine and pyoverdine pigments, which have been observed in some Pseudomonas such as Pseudomonas aeruginosa and Pseudomonas fluorescens (Reyes et al., 1981;Albesa et al., 1985). Pyoverdine is a molecule that scavenges iron from the environment and could compete with pathogens such as Pd that may be limited by iron availability (Mascuch et al., 2015;Reeder et al., 2017;Sass et al., 2017). Pyocyanine is an exotoxin with antimicrobial potential and could also be detrimental to Pd (Hassan and Fridovich, 1980;Baron and Rowe, 1981). Green fluorescence (typical of pyoverdine) was always observed co-occurring with orange fluorescence (typical of Pd) and/or blue fluorescence (typical of pyocyanine), whereas blue fluorescence was detected both alone or alongside orange fluorescence. These patterns are consistent with antifungal activity by Pseudomonas on the skin of E. fuscus. However, it is important to note that strains of the potential genera discussed above are not obligatorily antifungal, and that the metabarcoding method used in our study did not allow us to assess the proportion of antifungal metabolites produced as suggested by a previous study (Antony-Babu et al., 2017). The proportional overview of the community is also biased by the different copy numbers of ribosomal 16S rRNA gene present in different microorganisms (Klappenbach, 2001). Further studies should focus on the functional influence of these bacteria and other microorganisms (e.g., fungi and viruses) on the hostpathogen interaction between bats and Pd using approaches such as whole metagenomics sequencing.
The fact that some E. fuscus were, unexpectedly, already Pdpositive when we collected them before the experiment could have influenced our results. Since all collected bats were roosting in the same site, it is possible that all bats sampled were either carrying Pd at levels below our detectability threshold, or had previously been exposed. However, the facts that the mine was Pd-negative the year before our study, that we observed no clinical signs of WNS during our capture session, and that intensity of infection was very low for the few positive bats we did collect, suggest the cave only recently became infected with Pd and that impacts on our results were likely minimal. We also found that the skin microbiota shifted during captivity. This finding supports previous results in which samples from the same bat collected over time were more different than samples from multiple bats collected at the same time and place (Kolodny et al., 2019). The decrease in microbial diversity we observed in post-captivity samples likely reflects the controlled conditions in the laboratory setting. Captivity has been observed to negatively affect community diversity in the microbiota of wild amphibians (Kueneman et al., 2016;Sabino-Pinto et al., 2016) and it has been linked to sterile conditions in captivity (Loudon et al., 2014).
Our results suggest that the skin microbiota profile of E. fuscus may contribute to the apparent resistance of this species to Pd. FIGURE 3 | Relative abundance of different genera in the microbiota on bat skin from before (pre-captivity) and after (post-captivity) experiment. Analysis was performed on unrarefied ASVs table of taxa with relative abundance higher or equal to 0.1%.
Antifungal taxa predominated the microbiota prior to infection and remained stable in the skin microbiota after inoculation and 3 months of hibernation. However, despite apparent resistance to WNS and the abundance of anti-fungal bacteria, Pd could still affect host health in this species as it was shown to create skin lesion and inflammation (Moore et al., 2018) that may impact skin barrier integrity and lead to subsequent vulnerability to other pathogens. Although certain bacterial taxa remained abundant, there was destabilization of the microbiota that could be harmful for bats and this requires further study to truly understand impacts of infection on this species. We recommend that future studies try to better characterize post-Pd colonization patterns and, especially, characterize functional changes resulting from these patterns. WNS resistance and tolerance may reflect a range of mechanisms, from storage of extra body fat in fall , endogenous fat in epidermis composition (Frank et al., 2016(Frank et al., , 2018, variation in individual behavior (Frick et al., 2016) to variation in immune response (Maslo and Fefferman, 2015). Our results suggest that the skin microbiota could be another mechanism helping some bats survive infection with Pd. We recommend further research on gene expression by the microbiome of resistant and susceptible bat species to help understand disease mechanisms in WNS and contribute to conservation strategies for North American bats.

ETHICS STATEMENT
The animal study was reviewed and approved by The University of Winnipeg Animal Care Committee (Protocol Number AEO08399).

AUTHOR CONTRIBUTIONS
VL-L carried out experiments, analyzed the data, and wrote the manuscript. CW and ND designed the experiments, performed sampling, and wrote the manuscript. F-JL, VL-L, CW, and ND wrote the manuscript. All authors read and approved the final manuscript.