# Pseudogymnoascus destructans invasion stage impacts the skin microbial functions of highly vulnerable Myotis lucifugus

Virginie Lemieux-Labonté, Jananan Pathmanathan, Yves Terrat, Nicolas

Tromas, Anouk Simard, Catherine Haase, Cori Lausen, Craig Willis,

François-Joseph Lapointe

► To cite this version:

Virginie Lemieux-Labonté, Jananan Pathmanathan, Yves Terrat, Nicolas Tromas, Anouk Simard, et al.. Pseudogymnoascus destructans invasion stage impacts the skin microbial functions of highly vulnerable Myotis lucifugus. FEMS Microbiology Ecology, 2024, 100 (11), 10.1093/femsec/fiae138. hal-05182126

HAL Id: hal-05182126

https://hal.science/hal-05182126v1

Submitted on 5 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Pseudogymnoascus destructans invasion stage impacts the skin microbial functions of highly vulnerable Myotis lucifugus

Virginie Lemieux-Labonte<fnref n="1" />*, Jananan S. Pathmanathan<fnref n="2" />, Yves Terrat<fnref n="1" />, Nicolas Tomas<fnref n="1" />, Anouk Simard<fnref n="3" />, Catherine G. Haase<fnref n="4" />,
Cori L. Lausen<fnref n="5" />, Craig K. R. Willis<fnref n="6" />, François-Joseph Lapointe<fnref n="1" />

## Abstract

The role of the skin microbiome in resistance and susceptibility of wildlife to fungal pathogens has been examined from a taxonomic perspective but skin microbial function, in the context of fungal infection, has yet to be studied. Our objective was to understand effects of a bat fungal pathogen site infection status and course of invasion on skin microbial function. We sampled seven hibernating colonies of *Myotis lucifugus* covering three-time points over the course of *Pseudogymnoascus destructans* (Pd) invasion and white nose syndrome (pre-invasion, epidemic, and established). Our results support three new hypotheses about Pd and skin functional microbiome: (1) there is an important effect of Pd invasion stage, especially at the epidemic stage; (2) disruption by the fungus at the epidemic stage could decrease anti-fungal functions with potential negative effects on the microbiome and bat health; (3) the collection site might have a larger influence on microbiomes at the pre-invasion stage rather than at epidemic and established stages. Future studies with larger sample sizes and using meta-omics approaches will help confirm these hypotheses, and determine the influence of the microbiome on wildlife survival to fungal disease.

**Keywords:** microbial function assemblage; shotgun metagenomic; skin microbiome; *Pseudogymnoascus destructans*; white-nose syndrome; *Myotis lucifugus*; little brown bats

## Introduction

Fungal diseases are on the rise worldwide, causing unprecedented declines in wildlife populations (Daszak 2000; Jones et al. 2008, Fisher et al. 2012, 2020). The emergence of fungal skin infections with physiological consequences is of particular concern as they raise new challenges in conservation biology, operating by mechanisms not yet fully characterized (Fischer et al. 2020). Chytridiomycosis, caused by the fungus *Batrachochytrium dendrobatidis* (Berger et al. 1998), has caused decline in at least 501 amphibian species (6.5% of described amphibian species), 18% of which are confirmed or presumed extinct and 25% of which have experienced declines in abundance exceeding 90% (Schelle et al. 2019). Similarly, *Ophidiomyces ophidiicola*, the fungal agent of ophidiomycosis, is a potential threat to the European snake population and has potentially caused declines in North America (Allender et al. 2011, 2015, Clark et al. 2011, Franklinos et al. 2017).

The fungus *Pseudogymnoascus destructans* (Pd) (Gargas et al. 2009, Lorch et al. 2011), pathogenic agent of white-nose syndrome (WNS), has impacted North American bats over the last decade (Blehert et al. 2009, Frick et al. 2010, Frick et al. 2016). Pd invades the skin during hibernation, forming lesions that alter fluid bal-

ance, thermoregulation, and gas exchange, resulting in increased hibernal arousal frequency (Meteyer et al. 2009, Cryan et al. 2010, Warnecke et al. 2013, Verant et al. 2014, McGuire et al. 2017). Three species, northern myotis bats (*Myotis septentrionalis*), little brown bats (*Myotis lucifugus*), and tricolored bats (*Perimyotis subflavus*) declined by >90% since WNS arrival in North America (Cheng et al. 2021). In Canada, these species are listed as federally endangered due to mortality rates of 75%-90% during the several-year invasion stage of the disease (Canadian Wildlife Service 2014). Despite dramatic declines during the invasion and epidemic stages of the disease, some hibernating colonies of *M. lucifugus* have persisted following fungal invasion (Dobony et al. 2011, Langwig et al. 2012; Canadian Wildlife Service 2014; Frick et al. 2015, Vanderwolf and McAlpine 2021, Hooton et al. 2023), with population counts stabilizing at ~5%-30% of their initial sizes (Langwig et al. 2012, Frick et al. 2015). Variation in susceptibility to Pd could reflect mechanisms that allow for resistance (i.e. reduction or elimination of pathogen infection) and/or tolerance (i.e. reduction of the damage caused by infection) (Råberg et al. 2007, Råberg et al. 2009, Svensson and Råberg 2010). WNS resistance and tolerance are increasingly being studied, with mechanisms ranging from immune

<footnote n="1">¹Département de sciences biologiques, Université de Montréal, Montréal, Québec, H2V 0B3, Canada</footnote>

<footnote n="2">²Institut de Systématique, Evolution, Biodiversité (ISYEB), Sorbonne Université, Paris, 75005, France</footnote>

<footnote n="3">³Ministère de l'Environnement, de la Lutte contre les changements climatiques, de la Faune et des Parcs, Québec, G1R 5V7, Canada</footnote>

<footnote n="4">⁴Department of Biology, Austin Peay State University, Clarksville, TN, 37044, United States</footnote>

<footnote n="5">⁵Wildlife Conservation Society Canada, Kaslo, British-Columbia, V0G 1M0, Canada</footnote>

<footnote n="6">⁶Department of Biology and Centre for Forest Interdisciplinary Research, University of Winnipeg, Winnipeg, Manitoba, R3B 2E9, Canada</footnote>

<footnote>*Corresponding author. Département de sciences biologiques, Université de Montréal, 1375 Avenue Thérèse-Lavoie-Roux, Montréal, Québec, H2V 0B3, Canada.</footnote>

<footnote>E-mail: virginie.lemieux-labonte@umontreal.ca</footnote>

<footnote>Editor: [Leluo Guan]</footnote>
system response, increase in fat accumulation, roost selection, torpor regulation, and skin microbiome response (Hoyt et al. 2015, Frick et al. 2016, Langwig et al. 2017, Lemieux-Labonté et al. 2017, Lilley et al. 2017, Cheng et al. 2019, Auteri and Knowles 2020).

In a fungal disease context, the microbiome could play an important role in species' resistance or vulnerability. Originally, the microbiome has been defined as a collection of microorganisms (bacteria, archaea, eukaryotes, and viruses), their genomes (i.e. genes), and surrounding environmental conditions (Marchesi and Ravel 2015). However, most microbiome studies only investigated bacteria and/or, to some extent, fungi. Archaea, other eukaryotes, and viruses have been poorly characterized in comparison to these two taxonomic groups, in particular due to limited methods and resources. Therefore, when referring to microbiome studies in this paper, we will mostly be referring to the bacterial and/or fungal component of the microbiome.

The microbiome could be a critical component of host health, influencing biochemical and physiological processes, including defense against pathogens (Cho and Blaser 2012). The microbiome contributes to host defense 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). It also interacts with the innate and adaptive immune systems, and contributes to the maintenance of skin integrity and tissue repair (Lai et al. 2009, Curtis and Sperandio 2011, Naik et al. 2012). Several studies have highlighted the potential role of the skin microbiome in the resistance and susceptibility to multiple fungal pathogens (Harris et al. 2009, Bletz et al. 2013, 2017, Jani and Briggs 2014, Woodhams et al. 2014, Bataille et al. 2016, Lemieux-Labonté et al. 2017, 2020; Rebollar et al. 2019, Ange-Stark et al. 2023). However, previous studies have largely focused on the taxonomic composition of bacterial communities to infer microbiome functions.

Based on taxonomical analysis obtained by amplicon sequencing, the bat microbiome, contrary to other mammals, is mostly influenced by ecological factors, and not by phylogenetic relationships (Lutz et al. 2019). Microbial community assemblage of gregarious bats changes over time at the colony level (Kolodny et al. 2019), and is predominantly influenced by the environment (Avena et al. 2016, Lemieux-Labonté et al. 2016, Winter et al. 2017, Grisnik et al. 2020). In accordance with a protective effect of the skin microbiome, hundreds of microorganisms isolated from wild bats and their corresponding habitats have shown inhibitory effects on Pd by secreted compounds, contact inhibition, or volatile molecules (Hamm et al. 2017, Micalizzi et al. 2017). Additionally, compared with the microbiome of the cave environment, the bat microbiome is enriched with antifungal taxa (Grisnik et al. 2020). Indeed, antifungal strains of *Pseudomonas* isolated from the skin of *E. fuscus* inhibit Pd growth in vitro (Hoyt et al. 2015, Hamm et al. 2017), and improve survival of WNS-affected little brown bats when applied to the skin as a probiotic treatment in vivo (Cheng et al. 2017, Hoyt et al. 2019). The survival of *M. lucifugus* may be due to commensal bacteria, as persisting colonies with Pd exhibit a less diverse microbiome, and proportionally more antifungal bacterial taxa on their skin (Lemieux-Labonté et al. 2017). Accordingly, with the hypothesis that the microbiome protects against Pd, Grisnik et al. (2020) found that Pd-negative bats exhibited microbial assemblages with more antifungal taxa than Pd-positive bats in *P. subflavus*. Moreover, one study has observed a significant reduction of bacterial community richness and evenness after invasion by Pd in *M. lucifugus* (Ange-Stark et al. 2023). These changes could be

indicative of bacterial selection, but they could also be indicative of altered microbiome ability to perform functions that inhibit or prevent the growth of pathogens and/or opportunistic bacteria, and may in turn potentially amplify disease severity. Taken together, these findings show that the response of host microbial communities with respect to fungal pathogens is still unclear.

While amplicon sequencing can provide information about the effects of Pd on the taxonomic composition of the microbiome, we still know little about its effects on microbiome functions. In *P. subflavus*, Grisnik et al. (2022) reported different patterns for taxonomical versus functional assemblages. While taxon-based analysis revealed differences between Pd-positive and Pd-negative bats, functional assemblages remained unchanged by Pd-infection (Grisnik et al. 2022). Other recent work suggests that functional characterization of the microbial community is more sensitive than taxonomic assessment, and that genes, rather than species or genera, may be the appropriate parameter for understanding patterns of diversity in many microbial communities (Green et al. 2008, Burke et al. 2011, Louca et al. 2016, Ma et al. 2019). Furthermore, metagenomics yield greater resolution and accuracy in detection of putative functional genes compared to amplicon sequencing (Poretsky et al. 2014, Louca et al. 2016, Ranjan et al. 2016, Brumfield et al. 2020, Durazzi et al. 2021). Metagenomics should thus be used to study microbiomes in relation to taxonomy and/or functions (genes) (Langille 2018), and functional investigation is required to unravel the role of microbiome in vulnerable bat species facing WNS.

Previous WNS studies have investigated bat microbiomes based only on positive or negative individual Pd status. For example, Grisnik et al. (2022) found no difference in microbiome functions between Pd-positive and Pd-negative *P. subflavus*. However, responses of the microbiome to Pd could also change over the course of disease invasion. Following initial invasion when Pd prevalence is low but increasing, WNS, like many pathogen outbreaks, is characterized by an epidemic phase during which the host population declines and pathogen prevalence is moderate to high (Langwig et al. 2015a). Following the epidemic phase, Pd appears to become established with, for *M. lucifugus*, high and stable prevalence and stable (though smaller than original) host population size (Langwig et al. 2012, 2015a, Frick et al. 2015). Understanding how Pd infection status and disease stage interact with skin microbial function assemblage is important for understanding mechanisms underlying persistence of bats.

We assessed the functional potential of the skin microbiome of wild hibernating *M. lucifugus* colonies using shotgun metagenomics. We compared bats from colonies at three stages of the WNS outbreak (Pd-negative/pre-invasion; Pd-positive/epidemic; Pd-positive/established) (Langwig et al. 2015b) to (1) disentangle effects of site infection status (i.e. Pd-positive vs. Pd-negative) from the stage of pathogen invasion on the gene functions in the microbiome and (2) identify differences in microbial functions that could affect, or be affected, by Pd. Based on previous results in *M. lucifugus* (Lemieux-Labonté et al. 2017, Ange-Stark et al. 2023), we hypothesized that Pd would affect this species microbiome functional assemblage, with a decrease in functional diversity during the epidemic phase of rapid population decline. This reduced diversity will persist into the established phase with stable populations of persisting, potentially resistant or tolerant bats. We also hypothesized that genes associated with antifungal properties, or which act in competition with fungi, would be more prevalent in persisting bats following Pd establishment.
<table><caption>**Table 1.** Pooled bat skin microbiome samples information collected at different sites for different Pd invasion stages.</caption><thead><tr><th>Site Pd status</th><th>Invasion stage</th><th>Site</th><th>Province /State</th><th>Hibernaculum</th><th>Rock</th><th>Number of bats per site</th><th>Number of pools</th><th>Months of collection</th><th>Collection years<fnref n="1" /></th><th>Collection method<fnref n="2" /></th></tr></thead><tbody><tr><td rowspan="3">Negative</td><td rowspan="3">Pre-invasion</td><td>Alberta</td><td>Alberta</td><td>Cave</td><td>Limestone</td><td>10</td><td>2</td><td>January</td><td>Early 2018</td><td rowspan="3">M1</td></tr><tr><td>Montana</td><td>Montana</td><td>Cave</td><td>Limestone</td><td>10</td><td>2</td><td>September</td><td>Late 2019</td></tr><tr><td>Patate</td><td>Québec</td><td>Cave</td><td>Limestone</td><td>10</td><td>2</td><td>February</td><td>Late 2018</td></tr><tr><td rowspan="4">Positive</td><td rowspan="2">Epidemic</td><td>St-George (SG)</td><td>Ontario</td><td>Cave</td><td>Limestone</td><td>5</td><td>1</td><td>November</td><td>Early 2015</td><td rowspan="4">M2</td></tr><tr><td>Richard Lake</td><td>Ontario</td><td>Mine</td><td>Granite</td><td>15</td><td>3</td><td>February</td><td>Late 2016</td></tr><tr><td rowspan="2">Established</td><td>St-George (SG)</td><td>Ontario</td><td>Cave</td><td>Limestone</td><td>10</td><td>2</td><td>March</td><td>Late 2018</td></tr><tr><td>Lafleche</td><td>Québec</td><td>Cave</td><td>Calcite</td><td>10</td><td>2</td><td>April</td><td>Late 2019</td></tr><tr><td></td><td>Lames</td><td>Québec</td><td>Cave</td><td>Calcite</td><td>10</td><td>2</td><td>February</td><td>Late 2019</td><td></td></tr></tbody></table>

# Materials and methods

## Skin microbiome collection

We sampled skin microbiome from *M. lucifugus* in seven different hibernacula from 2015 to 2019 that were in either the pre-invasion, epidemic (one year after first Pd detection and observed population decline) or established phases (9 years after first Pd detection and noticeable WNS decline with actual stable population trend) of WNS. We sampled one of the seven sites longitudinally during both the pre-invasion and epidemic stages. We collected samples from 10 to 15 bats and one negative field control (humified swab open in the air) per hibernaculum. All methods were approved by the Animal Welfare and Ethics Committee at Université de Montréal (Protocol Number #16-015), Montana Fish, Wildlife and Park (Permit Number #2018-008-W), and the University of Winnipeg Animal Care Committee (Protocol Number AEO5639). We used two different protocols depending on the site, as the sites were part of different monitoring programs. We controlled for these procedural differences in our statistical analysis. For the first method (M1), bats were swabbed in the hand because samples also needed to be collected for other purposes. Bats were captured by gently pulling them off the wall and microbiome samples were collected while handling the bats using a clean pair of exam gloves for each bat. We swabbed the dorsal side of the right wing (forearm and wing skin) in overlapping linear strokes for 12 s with a sterile Puritan HydraFlock swab (Fisher Scientific) soaked in sterile 0.15 M NaCl 0.1% Tween-20 buffer. Swab tips were ejected into a 1.5-ml microtube containing the same buffer. Tubes were transferred to 4°C within 2 hr and to -20°C within 24 hr of sampling (see Table 1 for more information on the samples). All samples were shipped to Université de Montréal (Québec, Canada) for further processing. This method was used at one hibernaculum located 60 km south of Great Falls, Montana, USA (47°30'N, -111°1'W), and a second hibernaculum located in Alberta, Canada (53°5'N, 116°34'W). Both hibernacula were in the pre-invasion stage and Pd-negative at time of collection in 2018 (geographically isolated hibernacula from known positive sites and no decline observed at sampling time). Twenty bats (10 per site) were collected using this method. The Montana samples were heated to 56°C for 30 min to ensure eradication of potential rabies viruses prior to being shipped on dry ice, as requested by USDA animal products export certificates 0579-0256 and Canadian Food Inspection Agency importation permit A-2017-06316-4.

We used a slightly different sampling protocol (M2) for the other five hibernacula, as samples were collected as part of an-

nual visual survey where bats were not handled. For these sites, we randomly selected hibernating bats from those we could reach from the ground. Samples were collected without handling the bats by swabbing the back (fur) and wing forearm (skin) of each bat in overlapping linear strokes for 20 s with a sterile Whatman Omniswab (Fisher Scientific) soaked in sterile 0.15 M NaCl. We swabbed these bats longer as swabbing was harder to perform while not handling the bats to ensure quantities of biomass collected are similar to M1 methods. Swab tips were ejected into MoBio Powersoil DNA isolation Kit tubes (MoBio Laboratories), which were transferred to -20°C within 24 h of sampling. We used this method at one site (Patate, pre-invasion) on Anticosti island, Québec, Canada (49°29'N, -63°00'W), one site near Lake St-George (SG, 2016/pre-invasion and 2019/epidemic), Manitoba, Canada (51°22'N, 97°14'W), one site at Richard Lake (RL, epidemic) ~100 km east of Kenora, Ontario, Canada (49°45'N, 94°28'W), and two sites (Lames and Laflèche, established) located 60 km north of Gatineau, Québec, Canada (45°28'N, 75°42'W). Sixty bats were collected with this method and used for analysis. Pd site status was determined at epidemic sites by collecting separate swabs and performing qPCR (Table S0). For established sites, multiple qPCR assessment has been performed years before and after sampling, and they showed positive results for qPCR (Table S0 shows results for 2020, i.e. one year from microbiome collection). All samples with values crossing baseline below 40 cycles threshold (40 Ct) were considered positive. Mass bat mortalities were recorded in newly Pd-positive sites indicating that colonies were in the Pd epidemic invasion stage. In Pd-established sites, bat mortality rates were stable and similar to pre-invasion rates. The collection years indicate whether the samples were collected in early or late winter of the corresponding year.

## DNA extraction, purification and metagenomic sequencing

Bacterial genomic DNA was extracted from each swab using the ChargeSwitch gDNA Mini Bacteria Kit (Invitrogen™). To increase DNA yield, we modified the manufacturer's protocol by adding a 15-min heating period and eluting in 50 µl of elution buffer. All procedures were conducted with gloves in a laminar flow hood to limit sample contamination, and extractions were randomized to avoid detecting false patterns (Salter et al. 2014). Because of low input DNA and the fact that bat microbiomes show low inter-individual variation within colonies (Kolodny et al. 2019), samples were pooled by sites and collection years. Samples were

<footnote n="1">¹Late indicates samples collected in late winter (January to April) and early indicates samples collected in early winter (September to November).</footnote>

<footnote n="2">²Method M1: bats were handled and swabbed 12 s on the right wing. Method M2: bats were swabbed 20 s on back and wing while hanging on roost wall. Each pool is a mixed of five individuals collected at the same site.</footnote>
randomly pooled within site. Therefore, each line in Table 1 repre-
sents a pooled sample from five individuals occupying the same
site. Pooled samples were purified using a purification kit spin col-
umn (Zymo Research™), and sample concentrations were mea-
sured with a Qubit 2.0 Fluorometer (Invitrogen). Samples with
DNA concentration below 0.01 ng/µl were discarded (including
field negative control), except for extraction controls that were
used to keep track of contamination. To improve quality, we added
a mix of 20 µl nuclease-free water, 51 µl NEB3 buffer, and 3 µl
RNase to each sample prior to heating at 37°C for 60 min, followed
by 70°C for 20 min. Samples were then sent to Génome Québec
facilities at McGill University (Québec, Canada) for shotgun DNA
library preparation with NEB Ultra II kit (New England BioLabs
Inc.). Ten samples failing the quality control test were discarded,
including all three pooled extraction controls. Metagenomic DNA
sequencing was performed at Génome Québec facilities on No-
vaSeq6000 S4 flow cells of type 150 bp paired-end (Illumina).

Sequence processing

We obtained a total of 320 713 730 shotgun sequence reads and 96 855 546 460 bases for the 16 bat sample pools under study (Table S1 for details on the raw sequences). No sequences were obtained from extraction negative controls. Sequence reads were mapped to the *M. lucifugus* genome, and only unmapped reads representing bat skin microbial genes (i.e. genes from other eukaryotes, prokaryotes, and viruses) were selected for further analysis. First, quality trimming of raw reads was performed using the SolexaQA v3.1.7.1 program (Cox et al. 2010) with default parameters. Trimmed reads shorter than 75 nt were removed for further analysis. Reads with identical leading 20 bp were removed as artificial duplicates. From the trimmed high-quality reads, gene fragments were predicted using FragGeneScan-Plus v3.0 (Kim et al. 2015). Predicted protein fragments were clustered at a 90% similarity level using cd-hit v4.8.1 (Fu et al. 2012). One representative of each cluster was further used for a similarity search on the M5nr database (https://github.com/MG-RAST/myM5NR) using the DIAMOND sequence aligner (Buchfink et al. 2015). To assess function of protein fragments, we retrieved predicted functions of best hits through Clusters of Orthologous Groups (COG) of proteins databases (Tatusov et al. 2000). There are three COG levels ranging from general to more precise function: COG1 has 4 classification levels, COG2 level has 22 classes, and COG3 level has greater than 20 324 classes (Galperin et al. 2015). The COG database is updated periodically as new genomes become available. A matrix of function abundance by samples was produced to analyze the diversity and composition profiles in R Studio v1.1.447 (R Core Team 2018).

Alpha diversity

We first quantified alpha diversity of gene functions for all sam-
ples with the Shannon index (Shannon 1948), using the phyloseq
package (McMurdie and Holmes 2013) at the COG2 and COG3 lev-
els. The Shannon index, which includes both richness and even-
ness, was selected due to its reduced sensitivity to sample depth
differences (Haegeman et al. 2013, Preheim et al. 2013). To ensure
normality, we first transformed data with Tukey's ladder transfor-
mation (Tukey 1977). Alpha diversity values were compared using
linear mixed-effect models (lme functions) in R, with alpha di-
versity as the response variable, site Pd status and invasion stage
at collection time as fixed effects, and site, collection years, heat-
ing, and collection methods as random effects. Significance was
tested by performing 999 permutations of the data, and the re-

sulting models were computed to generate a distribution of AIC
statistics. The reference AIC (not randomized) was considered sig-
nificant using a P-value threshold of 0.05.

Beta diversity

Microbiome data are compositional because they are relative. The total number of counts per sample is highly variable and constrained by the maximum number of sequences reads. This induces strong dependencies among feature abundances, as an increase of one feature implies the decrease of counts for another feature so that the total number of counts does not exceed the specified sequencing depth (Calle 2019). Ignoring the compositional nature of microbiome data may lead to erroneous results. Consequently, we used the coherent Aitchison distance that respects the compositional properties of microbiome data (Aitchison 1982, Aitchison et al. 2000, Gloor et al. 2017, Calle 2019). The Aitchison distance is defined as the Euclidian distance after a centered-log transformation (CLR). However, since microbiome data contain a high number of zeros, they should be corrected prior to the centered-log transformation to ensure validity of results. We used the function transform from the microbiome package to perform the CLR transformation (Lahti and Shetty 2017). This function applies a pseudocount before taking log results. Beta diversity was then visualized with principal component analysis (PCA) using the ordinate and plot ordinate functions of the phyloseq package (McMurdie and Holmes 2013).

To test for effects of site Pd status and invasion stage, we used a
distance-based redundancy analysis (db-RDA) (Legendre and An-
derson 1999). This approach first decomposes distances into prin-
cipal coordinates, and then applies RDA to the corresponding prin-
cipal coordinates using the capscale function of the R package
vegan (Oksanen et al. 2019). We also computed a partial db-RDA
to better understand the influence of site Pd status and invasion
stage, while controlling for potential confounding factors (Davies
and Tso 1982). To select the constrained variables and to max-
imize degrees of freedom, we performed first a dbRDA not con-
trolled on potential confounding factors (i.e. sites, year of collec-
tion, collection method, heating) (Table S3) and only kept signifi-
cant factors for further analysis. By doing so, sites and collection
years were used to constrain the partial db-RDA at COG3, whereas
only collection years were used at COG2 level. Significance was
tested via 999 permutations with the anova.cca function of the R
package vegan. For all analyses, a P-value threshold of 0.05 was
considered significant.

Skin microbial functions composition

We assessed the effect of invasion stage on the composition of microbial functions down to the COG3 level using the Analysis of Composition of Microbiomes (ANCOM 2.0) (Mandal et al. 2015). ANCOM is based on the analysis of difference in pairwise log-ratios of function abundance/relative abundance, between comparison groups of interest. It is appropriate for compositional data (Gloor et al. 2017), under the assumption that a minority of features (<25%) are differentially abundant (Mandal et al. 2015). This analysis does not output a P-value. Instead for each function, we computed a W statistic indicating the number of significantly different pairwise log-ratios, while controlling for false discoveries. Significance level was set to 0.05 and multiple corrections were performed with the Benjamini-Hochberg procedure (Benjamini and Hochberg 1995). ANCOM analysis was conducted on unrarefied COG3-level table of the more abundant functions (>1%).
## Bipartite graph analysis

We used a bipartite graph analysis to better explore the microbial functions without a priori group definition. The aim is to determine how samples group and by what discriminant function. A bipartite network is a collection of nodes, grouped in two separate sets, with no interactions among nodes within a set. In the present case, the bat samples were pooled by sites and year of collection to explore relationships with functional traits at the COG3 level. We used the R package bipartite and bipartiteD3 to plot bipartite graphs and to compute modules that are aggregated sets of interacting components (Dormann et al. 2008, Terry 2018). Their defining feature is that interaction within modules is more prevalent than between modules (Newman 2003, Newman and Girvan 2004, Fortunato 2010). In other words, modules are link-rich clusters of genes in a group. In addition, we chose to compute the species specificity index, the coefficient of variation of interactions, normalized to values between 0 (low specificity) and 1 (high specificity) to target the most discriminant functions between modules (Julliard et al. 2006, Poisot et al. 2012). We performed the modules and species specificity index computation on the most abundant (>1%) functions of the eight samples pooled by year and site (177 functions).

## Results

### Alpha diversity functional profile

The site Pd status model at the COG2 level (22 functions overall) was not significant (Table 2 and Fig. 1A), whereas the invasion stage model was significant (Table 2 and Fig. 1B). Mean-transformed Shannon diversity values differed by 6% (–855 ± 566) between pre-invasion and epidemic stage (Fig. 1B), but this was not significant (Table 2). Shannon diversity values from Pd-established site were significantly different from pre-invasion and epidemic stages. Alpha diversity of pre-invasion groups was significantly lower by 6% (920 ± 577) from established stage (Table 2 and Fig. 1B). Between epidemic and established stages, mean-transformed Shannon diversity was higher by 13% (1208 ± 492). The top model (based on ΔAIC) included the invasion stages: epidemic vs. established. Models at the COG3 level (2324 functions overall) were not significant (Table S2).

### Beta diversity functional profile

Clustering according to invasion stage is not clearly defined for the COG2 level (Fig. 2A), but is much clearer for the COG3 (Fig. 2B). No models were significant at the COG2 level (Table 3). Invasion stage clusters are observable in the PCA at the COG3 level, with a clear separation on the first principal axis that explained 27.4% of variation between bats in the established stage vs. bats in the epidemic stage (Fig. 2B). However, according to the db-RDA, the confounding effect of sites was predominant and explained up to 9% of the functional variation (Table 4). These results indicate that invasion stage is nested into the site factor effect, since without site, invasion stage explained up to 8% of variation (Table 4). Even when unconstrained by site, site Pd status explained no redundant variation in the functional variation (Table 4).

### ANCOM on functional profiles

We conducted the ANCOM analysis on the 332 more abundant COG3 functions, representing >1% of total abundance contained within the 16 pooled samples. Eleven functions differed according to invasion stage (Fig. 3 and Table S4). The epidemic samples appeared to be the most distinct. Three functions differed

between the epidemic group and the pre-invasion group. Ten functions differed between the epidemic group and established group (Fig. 3). The epidemic bats exhibited the lowest abundance for the outer membrane protein porin, the ABC type transport system involved in lipoprotein release permease component, ribulose 5 phosphate 4 epimerase and related epimerases and aldolases, phospholipase C, putative theorine efflux protein, high affinity Fe2 Pe2 permease, Mn2 and Fe2 transporter of the NRAMP family, and transcriptional regulators compared to bats from established stage. On the contrary, the functions asparagine synthase glutamine hydrolyzing and type II restriction enzyme methylase subunits had higher abundances in the epidemic bats compared to the bats from established stage (Fig. 3). The pre-invasion and established groups appeared to have similar function abundances, except for A 4 amino 4 deoxy L arabinose transferase and related glycosyltransferases of PMT family, which exhibited larger functional abundances mean in the bats from the established invasion stage (Fig. 3).

## Bipartite analysis

We explored the microbial functions without a priori group definition using a bipartite graph analysis of the samples pooled by site and collection year at the COG3 level. We performed the modules and species specificity index computation on the most abundant (>1%) functions of the eight pooled samples (see Table S5 for modules and the ten functions with higher species specificity by module). Out of these 177 functions, the nine functions with higher species specificity index were considered as the most discriminating between modules (Fig. 4).

The first module clusters pre-invasion bats collected from Patate with the bats in epidemic stage from St-George collected in 2019 (SG19) (Fig. 4) and is characterized by 71 functions (see Table S5). Among the nine more discriminant functions, three belong to this module, including ATP-dependent Lon protease bacterial type, arylsulfatase A and related enzymes, and type I site-specific restriction modification system R restriction subunit and related helicases (Fig. 4). Pre-invasion bats collected in 2016 (SG16) are part of a second module associated with 33 functions (Table S5), including two among the nine more discriminant functions: catalase and outer membrane receptor proteins mostly Fe transport (Fig. 4). The epidemic bats from the Richard Lake site are forming a third module characterized by 25 functions (Table S5 and Fig. 4). The last module clusters pre-invasion bats collected from sites in Montana and Alberta with colonies in the established invasion stage collected at sites Lames and Laflèche. Collectively, these four sites are characterized by 42 functions (Table S5), including four among the nine most discriminating functions: high affinity Fe₂ Pb₂ permease, transposase and inactivated derivatives TnpA family, outer membrane protein porin, and ABC type transport system involved in lipoprotein release permease component (Fig. 4).

## Discussion

In this study, we used shotgun metagenomics to identify the functional potential of the skin microbiome of hibernating *M. lucifugus* colonies not affected by Pd and at different stages of pathogen invasion. Our objectives were to investigate the effect of site infection status and invasion stage on the microbial functional assemblage and to understand how the microbiome could potentially influence or be influenced by fungal infection. Overall, microbial functional assemblages of bat colonies were not impacted by site Pd-positive (epidemic and established) or Pd-negative status (pre-invasion) alone. Our metagenomic results on *M. lucifugus* are con-
<table><caption>**Table 2.** Linear mixed models of functional alpha diversity at the COG2 levels.</caption><thead><tr><th>Model</th><th>Model AIC</th><th>Δ AIC</th><th>Likelihood ratio</th><th>P-value</th></tr></thead><tbody><tr><td>Epidemic vs established</td><td>127.940</td><td>0.000</td><td>-56.970</td><td>0.006<sup>1**</sup></td></tr><tr><td>Pre-invasion vs established</td><td>150.296</td><td>22.356</td><td>-68.148</td><td>0.024<sup>1*</sup></td></tr><tr><td>Pre-invasion vs epidemic</td><td>181.806</td><td>53.866</td><td>-83.903</td><td>1<sup>1</sup></td></tr><tr><td>Invasion stage</td><td>225.065</td><td>97.125</td><td>-104.533</td><td>0.004<sup>**</sup></td></tr><tr><td>Site Pd status</td><td>248.963</td><td>121.023</td><td>-117.482</td><td>0.925</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Functional Shannon diversity of skin microbiome samples transformed at the COG2 level according to (A) site Pd status and (B) invasion stage. Pairwise comparisons among invasion stages are significantly different for pre-invasion vs. established, and epidemic vs. established (Table 2).</figcaption></figure>

<table><caption>**Table 3.** Db-RDA model tests of Aitchison distance of *M. lucifugus* skin microbial functions at the COG2 level.</caption><thead><tr><th>Model formula</th><th>F-statistic</th><th>P-value</th><th>R<sup>2</sup> adjusted</th></tr></thead><tbody><tr><td>~site Pd status</td><td>0.42</td><td>0.81</td><td>NA</td></tr><tr><td>~Invasion stage</td><td>1.12</td><td>0.37</td><td>NA</td></tr><tr><td>~site Pd status + invasion stage</td><td>1.25</td><td>0.37</td><td>NA</td></tr></tbody></table>

sistent with Grisnik et al. (2020, 2022), who reported considerable functional and taxonomic redundancy among Pd-positive and Pd-negative individuals but in *P. subflavus* bats from epidemic and established sites. However, our results appear to differ from what was found on *M. lucifugus* taxonomic studies, which found differences in bacterial taxa among Pd-status (Lemieux-Labonté et al. 2017, Ange-Stark et al. 2023). We must stress that our results are reliable at a site/colony level and may differ from previous studies as samples from established and epidemic sites were not individually accessed for Pd. Consequently, pools in these sites were

probably formed of negative and positive individuals. We suggest future studies to assess individual Pd load and metagenomes to obtain a finer-scale resolution of fungal infection intensity effect on skin microbiome functions.

Still, when we considered the disease invasion stage and not just the site positive or negative status, we did find an effect on the microbiome functional assemblages. We observed significantly lower diversity in pre-invasion and epidemic stage bats compared to the established disease stage but we did not observe a significant lower diversity in epidemic stage bats compared to pre-

<footnote n="1">¹P-value adjusted with Bonferroni method.</footnote>

<footnote>Asterisks indicate significant results.</footnote>

<footnote>*** ≤0.001.</footnote>

<footnote>** ≤0.01.</footnote>

<footnote>* ≤0.05. Site, collection year, heating, and collection method were random effects. Significance tested by 999 permutations.</footnote>

<footnote>Asterisks indicate significant results.</footnote>

<footnote>*** ≤0.001.</footnote>

<footnote>** ≤0.01.</footnote>

<footnote>* ≤0.05.</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Principal component analysis of Aitchison skin microbiome functional profiles according to invasion stage and sites for (A) the COG2 level and (B) the COG3 level. Each point represents 5 bats.</figcaption></figure>
<table><caption>**Table 4.** Db-RDA model tests of Aitchison distance of *M. lucifugus* skin microbial functions at the COG3 level.</caption><thead><tr><th>Model formula</th><th>F-statistic</th><th>P-value</th><th>R<sup>2</sup> adjusted</th></tr></thead><tbody><tr><td>~ site Pd status | site, collection method, collection years</td><td>0</td><td>NA</td><td>NA</td></tr><tr><td>~Invasion stage | site, collection method, collection years</td><td>0</td><td>NA</td><td>NA</td></tr><tr><td>~ site Pd status + invasion stage | site, collection method, collection years</td><td>0</td><td>NA</td><td>NA</td></tr><tr><td>~site | collection method, collection years</td><td>1.50</td><td>0.05*</td><td>0.09</td></tr><tr><td>~ site Pd status | collection method, collection years</td><td>0</td><td>NA</td><td>NA</td></tr><tr><td>~Invasion stage | collection method, collection years</td><td>2.31</td><td>0.007**</td><td>0.08</td></tr><tr><td>~ site Pd status + Invasion stage | collection method, collection years</td><td>2.31</td><td>0.005**</td><td>0.08</td></tr><tr><td>~Invasion stage + site | collection method, collection years</td><td>1.50</td><td>0.04*</td><td>0.09</td></tr></tbody></table>

Asterisks indicate significant results.

***≤0.001,

**≤0.01,

*≤0.05.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Count abundance of functions detected as significantly divergent by ANCOM according to invasion stages at 0.8 thresholds (Table S4 for W statistics). Asterisks indicate significant differences between corresponding groups.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Bipartite graph associating samples pooled by sites and collection years with the more abundant functions at the COG3 level (>1%). Only the nine functions with the higher species specificity index are presented (for other functions see Table S5). The four different modules identified in the complete bipartite graph are highlighted to illustrate the relationships with the corresponding functions. Wider links indicate functions with higher relative abundance in a sample.</figcaption></figure>

invasion bats. Furthermore, our results suggest that bat colonies persisting after invasion may regain Shannon diversity at higher levels than pre-invasion bats (Lemieux-Labonté et al. 2017). This could indicate an enrichment in new functions in response to the disease and/or a better evenness after invasion. Yet, this trend was only observed at higher functional levels (COG2), which suggests change in more general function groups. Consistent with these findings, our first hypothesis is not supported, as we did not observe a decrease in diversity between pre-invasion bats and epidemic stage bats, and that an increase in diversity was observed in bats at the established stage of invasion. The functional assemblage diversity differed the most between bats collected at epidemic and established invasion stages. This pattern could represent selection for a particular functional profile that could help bats survive persistent disease during the established stage or a return to a normal state similar to that occurring pre-invasion.

By analyzing functional abundance, we observed that predominant genes detected are parts of various biological processes, such as molecular functions or cellular components (UniProt 2019). We also discovered changes in predominant functions between invasion stages. Overall, the established group had a community composition similar to that of pre-invasion bats and the epidemic

stage appears to be the most different. Many functions with anti-fungal potential are less abundant in epidemic stage bats. Some of these functions are involved in iron (Fe2+) and manganese (Mn2+) acquisition and transport [Mn2 and Fe2 transporter of the NRAMP family-like protein, high affinity Fe2 Pb2 permease, and outer membrane receptor proteins mostly Fe] (UniProt 2019). Iron and manganese are essential micronutrients that play a central role in infection processes because they serve as cofactors in many reactions, including many with direct and indirect roles in pathogens virulence (Gerwien et al. 2018). Hosts can also fight pathogens by deploying toxic levels of certain metals (Hood and Skaar 2012). Upregulation or downregulation of function influencing Fe2+ and Mn2+ bioavailability could influence Pd invasion (Reeder et al. 2017, Gerwien et al. 2018). We also detected less phospholipase C, which can be toxic and cause lysis in eukaryotic cells (Titball 1993), at the epidemic stage. However, as these functional abundances were similar between pre-invasion and Pd-established colonies, our hypothesis on the selection of anti-fungal functions following Pd establishment is not supported. The low abundance of these functions in newly infected epidemic colonies is consistent with a new hypothesis that disruption by the fungus could cause a decrease in these functions with potentially negative effects on the micro-
biome and bat health. Still, more work is needed to understand
the implications of changes in microbiome function on disease
severity and survival of bats with WNS.

Moreover, we also observed an interesting trend in pre-invasion
bats. This group appears to have a wider dispersion in their func-
tional composition abundance and diversity. Indeed, functional
abundance appears more variable between pooled samples in the
pre-invasion group than in the epidemic and established groups.
This is also true for beta diversity analysis, where one can ob-
serve a wider dispersion for this group, with pooled samples being
more different among sites. The pre-invasion group has the largest
dispersion, followed by the epidemic and the established groups,
supported by the significant interaction between invasion stage
and site. Accordingly, sites might have a larger influence on pre-
invasion microbiomes compared to the microbiome collected at
epidemic and established sites. Unfortunately, this could be an ar-
tifact of the relatively few sites sampled within each category and
large geographic dispersion of pre-invasion group, and we suggest
further investigation to verify this trend.

Although our results are consistent with variation in effects of
WNS on microbiome function depending on the stage of pathogen
invasion, we must stress that shotgun metagenomics only eluci-
dates the functional potential of putative COGs, and does not re-
flect what is truly expressed by the corresponding genes (Jansson
and Hofmockel 2018). Transcriptomic and other meta-omics ap-
proaches will be the next step to elucidate the true skin microbial
potential in the WNS context. It is also important to note that pat-
terns we observed in this study could be site- (i.e. colony-)specific
and might not reflect invasion stage, owing to microbiome dif-
ferences associated with geographical/environmental variation.
This site effect on function assemblage is consistent with what
was found in previous taxon-based microbiota studies (Avena et
al. 2016, Lemieux-Labonté et al. 2017). Thus, we could expect
site-specific selection on microbiome function due to WNS. Bats
might also show mechanisms of resistance or tolerance to WNS
independent of the microbiome due to physiological or behav-
ioral mechanisms such as environmental selection and physio-
logical adaptation (Cheng et al. 2019, Langwig et al. 2017). Thus,
there is a possibility that microbiome may act synergistically with
other mechanisms to help mitigate Pd infection. Another poten-
tial limitation of our study is that time (i.e. year) and WNS-status
(positive or negative) are confounded. However, we controlled for
confounding effects such as site, collection year, and collection
method in our analyses and found support for the hypothesis of
an effect of invasion stage on skin functions assemblage.

This study provides evidence of changes in microbiome func-
tional profiles that may be attributed to WNS. It represents an
important step towards understanding the role of skin micro-
biome in wild bat colonies facing an emergent fungal pathogen.
Overall, our results support that microbiome functional profiles
vary among Pd invasion stages but that changes are site-specific.
Such diversified responses could help explain why some colonies
may be affected differently by Pd infection during the initial epi-
demic phase of WNS (Frick et al. 2017). Interestingly, we found
that skin microbial functions of persisting bats in the established
phase of WNS were similar to sites that were pre-infection, pos-
sibly indicating reestablishment of microbiomes. It is still un-
clear whether the skin microbiome could protect against Pd, or
whether the changes we observed reflect neutral or maladaptive
effects of the disease. Future studies design should use concur-
rent sampling from multiple positive and negative sites in dif-
ferent phases of WNS invasion and assess individual change in
Pd load in relation to functional microbiome. Being able to have

more samples and perform sequencing of individual bat samples
would increase statistical power and help support or reject new
hypotheses brought by this research. However, given how far the
disease has spread and the range of bat species and communi-
ties now involved, this may no longer be possible, at least in the
eastern half of North America. We thus suggest local authori-
ties or researchers doing annual winter survey in western sites
not affected by the disease to collect additional swabs from bats
skin for microbiome study, anticipating infection, and document-
ing wing microbiomes throughout the infection stages. As we ex-
pect WNS to spread further, these samples will allow monitor-
ing the microbiome throughout disease stages. This kind of sam-
pling regime could provide better insight to understand interac-
tions between Pd and the microbiome. Understanding the roles of
the microbiome and developing new tools to better manage the
disease in the future, may hinge on use of combined meta-omics
approaches, and consideration of other resistance and tolerance
mechanisms.

Acknowledgments

We would like to thank all members of the Willis lab at the Univer-
sity of Winnipeg for their help with this project. We thank Olivier
Cameron-Trudel, Dave Hobson, Cory Olson, Nicole Dorville, and
Kaleigh Norquay for sampling assistance. We also thank the Mi-
crobial Evolutionary Genomics Laboratory of Jesse Shapiro at Uni-
versité de Montréal for technical support.

Author contributions

Virginia Lemieux-Labonté (Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Visualization, Writing – original draft, Writing – review & editing), Jananan S. Pathmanathan (Formal analysis, Resources, Software, Writing – review & editing), Yves Terrat (Formal analysis, Resources, Software, Validation, Writing – review & editing), Nicolas Tromas (Conceptualization, Methodology, Resources, Writing – review & editing), Anouk Simard (Data curation, Resources, Validation, Writing – review & editing), Catherine G. Haase (Data curation, Resources, Validation, Writing – review & editing), Cori L. Lausen (Data curation, Resources, Validation, Writing – review & editing), Craig K. R. Willis (Data curation, Resources, Validation, Writing – review & editing), and François-Joseph Lapointe (Conceptualization, Funding acquisition, Supervision, Validation, Writing – review & editing)

**Supplementary data**

Supplementary data is available at *FEMSEC Journal* online.

Conflict of interest: The authors declare that they have no competing interests.

Funding

This work was supported by the Discovery Grants from the Natu-
ral Sciences and Engineering Research Council (NSERC Canada)
[RGPIN-2015-05219] to François-Joseph Lapointe, as well as the
*Fightwns* Micro Grants for Microbats Program [2017] to Virginie
Lemieux-Labonté. The funders had no role in study design, data
collection and analysis, decision to publish, or preparation of the
manuscript.
## Data availability

The raw sequencing data supporting the conclusions of this article are available here: https://figshare.com/s/8c742fca1ec5e317 324a. The scripts are available here: https://figshare.com/s/5354 5eafd27ff4a20d53.

## References

Aitchison J, Barceló-Vidal C, Martín-Fernández JA et al. Logratio analysis and compositional distance. *Math Geol* 2000;**32**:271-5.

Aitchison J. The statistical analysis of compositional data. *J R Stat Soc Ser B Methodol* 1982;**44**:139-60.

Allender MC, Dreslik M, Wylie S et al. *Chrysosporium* sp. infection in eastern massasauga rattlesnakes. *Emerg Infect Dis* 2011;**17**:2383.

Allender MC, Raudabaugh DB, Gleason FH et al. The natural history, ecology, and epidemiology of *Ophidiomycyse ophiodiicola* and its potential impact on free-ranging snake populations. *Fungal Ecol* 2015;**17**:187-96.

Ange-Stark M, Parise KL, Cheng TL et al. White-nose syndrome restructures bat skin microbiomes. *Microbiol Spectr* 2023;**11**:e02715-23.

Auteri GG, Knowles LL. Decimated little brown bats show potential for adaptive change. *Sci Rep* 2020;**10**:1-10.

Avena CV, Parfrey LW, Leff JW et al. Deconstructing the bat skin microbiome: influences of the host and the environment. *Front Microbiol* 2016;**7**:1753. https://doi.org/10.3389/fmicb.2016.01753.

Bataille A, Lee-Cruz L, Tripathi B et al. Microbiome variation across amphibian skin regions: implications for chytridiomycosis mitigation efforts. *Microb Ecol* 2016;**71**:221-32.

Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. *J R Stat Soc Ser B Methodol* 1995;**57**:289-300.

Berger L, Speare R, Daszak P et al. Chytridiomycosis causes amphibian mortality associated with population declines in the rain forests of Australia and Central America. *Proc Natl Acad Sci USA* 1998;**95**:9031-6.

Blehert DS, Hicks AC, Behr M et al. Bat white-nose syndrome: an emerging fungal pathogen? *Science* 2009;**323**:227.

Bletz MC, Loudon AH, Becker MH et al. Mitigating amphibian chytridiomycosis with bioaugmentation: characteristics of effective probiotics and strategies for their selection and use. *Ecol Lett* 2013;**16**:807-20.

Bletz MC, Perl RGB, Bobowski BTC et al. Amphibian skin microbiota exhibits temporal variation in community structure but stability of predicted Bd-inhibitory function. *ISME J* 2017;**11**:1521-34.

Brumfield KD, Huq A, Colwell RR et al. Microbial resolution of whole genome shotgun and 16S amplicon metagenomic sequencing using publicly available NEON data. *PLoS One* 2020;**15**:e0228899.

Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. *Nat Methods* 2015;**12**:59-60.

Burke C, Steinberg P, Rusch D et al. Bacterial community assembly based on functional genes rather than species. *Proc Natl Acad Sci USA* 2011;**108**:14288-93.

Calle ML. Statistical analysis of metagenomics data. *Genomics Inform* 2019;**17**:e6.

Canadian Wildlife Service. Committee on the Status of Endangered Wildlife in Canada. COSEWIC Assessment and Status Report on the Little Brown Myotis (Myotis lucifugus), Northern Myotis (Myotis septentrionalis), Tri-Colored Bat (Perimyotis subflavus) in Canada. 2014. https://publications.gc.ca/collections/collection_2014/ec/C W69-14-688-2014-eng.pdf (21 October 2024, date last accessed).

Cheng TL, Gerson A, Moore MS et al. Higher fat stores contribute to persistence of little brown bat populations with white-nose syndrome. *J Anim Ecol* 2019;**88**:591-600.

Cheng TL, Mayberry H, McGuire LP et al. Efficacy of a probiotic bacterium to treat bats affected by the disease white-nose syndrome. *J Appl Ecol* 2017;**54**:701-8.

Cheng TL, Reichard JD, Coleman JTH et al. The scope and severity of white-nose syndrome on hibernating bats in North America. *Conserv Biol* 2021;**35**:1586-97.

Cho I, Blaser MJ. The human microbiome: at the interface of health and disease. *Nat Rev Genet* 2012;**13**:260-70.

Clark RW, Marchand MN, Clifford BJ et al. Decline of an isolated timber rattlesnake (*Crotalus horridus*) population: interactions between climate change, disease, and loss of genetic diversity. *Biol Conserv* 2011;**144**:886-91.

Cox MP, Peterson DA, Biggs PJ. SolexaQA: at-a-glance quality assessment of Illumina second-generation sequencing data. *BMC Bioinf* 2010;**11**:1-6.

Cryan PM, Meteyer C, Boyles JG et al. Wing pathology of white-nose syndrome in bats suggests life-threatening disruption of physiology. *BMC Biol* 2010;**8**:135.

Curtis MM, Sperandio V. A complex relationship: the interaction among symbiotic microbes, invading pathogens, and their mammalian host. *Mucosal Immunol* 2011;**4**:133-8.

Daszak P. Emerging infectious diseases of wildlife–threats to biodiversity and human health. *Science* 2000;**287**:443-9.

Davies PT, Tso MK-S. Procedures for reduced-rank regression. *Appl Stat* 1982;**31**:244.

Dobony CA, Hicks AC, Langwig KE et al. Little brown myotis persist despite exposure to white-nose syndrome. *J Fish Wildl Manag* 2011;**2**:190-5.

Dormann CF, Gruber B, Fründ J. Introducing the bipartite package: analysing ecological networks. *R News* 2008;**8**:8-11.

Durazzi F, Sala C, Castellani G et al. Comparison between 16S rRNA and shotgun sequencing data for the taxonomic characterization of the gut microbiota. *Sci Rep* 2021;**11**:3030.

Fischer NM, Dool SE, Puechmaille SJ. Seasonal patterns of *Pseudogymnoascus destructans* germination indicate host-pathogen coevolution. *Biol Lett* 2020;**16**:20200177.

Fisher MC, Gurr SJ, Cuomo CA et al. Threats posed by the fungal kingdom to humans, wildlife, and agriculture. *mBio* 2020;**11**:e00449-20.

Fisher MC, Henk DA, Briggs CJ et al. Emerging fungal threats to animal, plant and ecosystem health. *Nature* 2012;**484**:186-94.

Fortunato S. Community detection in graphs. *Phys Rep* 2010;**486**:75-174.

Franklinos LHV, Lorch JM, Bohuski E et al. Emerging fungal pathogen *Ophidiomycyse ophiodiicola* in wild European snakes. *Sci Rep* 2017;**7**:1-7.

Frick WF, Cheng TL, Langwig KE et al. Pathogen dynamics during invasion and establishment of white-nose syndrome explain mechanisms of host persistence. *Ecology* 2017;**98**:624-31.

Frick WF, Pollock JF, Hicks AC et al. An emerging disease causes regional population collapse of a common North American bat species. *Science* 2010;**329**:679-82.

Frick WF, Puechmaille SJ, Hoyt JR et al. Disease alters macroecological patterns of North American bats: Disease alters macroecology of bats. *Global Ecol Biogeogr* 2015;**24**:741-9.

Frick WF, Puechmaille SJ, Willis CKR. White-nose syndrome in bats. In: *Bats in the Anthropocene: Conservation of Bats in a Changing World*. Cham: Springer, 2016, 245-62.

Fu L, Niu B, Zhu Z et al. CD-HIT: accelerated for clustering the next-generation sequencing data. *Bioinformatics* 2012;**28**:3150-2.
Galperin MY, Makarova KS, Wolf YI et al. Expanded microbial genome coverage and improved protein family annotation in the COG database. *Nucleic Acids Res* 2015; **43**:D261–9.

Gargas A, Trest MT, Christensen M et al. *Geomyces destructans* sp. nov. associated with bat white-nose syndrome. *Mycotaxon* 2009; **108**:147–54.

Gerwien F, Skrahina V, Kasper L et al. Metals in fungal virulence. *FEMS Microbiol Rev* 2018; **42**:fux050. https://doi.org/10.1093/femsre/fux050.

Gloor GB, Macklaim JM, Pawlowsky-Glahn V et al. Microbiome datasets are compositional: and this is not optional. *Front Microbiol* 2017; **8**:2224.

Green JL, Bohannan BJM, Whitaker RJ. Microbial biogeography: from taxonomy to traits. *Science* 2008; **320**:1039–43.

Grice EA, Segre JA. The skin microbiome. *Nat Rev Micro* 2011; **9**:244–53.

Grisnik M, Bowers O, Moore AJ et al. The cutaneous microbiota of bats has in vitro antifungal activity against the white nose pathogen. *FEMS Microbiol Ecol* 2020; **96**:fiz193.

Grisnik M, Grinath JB, Munafo JP et al. Functional redundancy in bat microbial assemblage in the presence of the white nose pathogen. *Microb Ecol* 2023; **86**:713–26.

Haegeman B, Hamelin J, Moriarty J et al. Robust estimation of microbial diversity in theory and in practice. *ISME J* 2013; **7**:1092–101.

Hamm PS, Caimi NA, Northup DE et al. Western bats as a reservoir of novel *Streptomyces* species with antifungal activity. *Appl Environ Microb* 2017; **83**:e03057–16. /aem/83/5/e03057-16.atom.

Harris RN, Brucker RM, Walke JB et al. Skin microbes on frogs prevent morbidity and mortality caused by a lethal skin fungus. *ISME J* 2009; **3**:818–24.

Hood MI, Skaar EP. Nutritional immunity: transition metals at the pathogen-host interface. *Nat Rev Micro* 2012; **10**:525–37.

Hooton LA, Adams AA, Cameron A et al. Effects of bat white-nose syndrome on hibernation and swarming aggregations of bats in Ontario. *Can J Zool* 2023; **101**:886–95.

Hoyt JR, Cheng TL, Langwig KE et al. Bacteria isolated from bats inhibit the growth of *Pseudogymnoascus destructans*, the causative agent of white-nose syndrome. *PLoS One* 2015; **10**:e0121329.

Hoyt JR, Langwig KE, White JP et al. Field trial of a probiotic bacteria to protect bats from white-nose syndrome. *Sci Rep* 2019; **9**:9158.

Jani AJ, Briggs CJ. The pathogen *Batrachochytrium dendrobatidis* disturbs the frog skin microbiome during a natural epidemic and experimental infection. *Proc Natl Acad Sci USA* 2014; **111**:E5049–58.

Jansson JK, Hofmockel KS. The soil microbiome—from metagenomics to metaphenomics. *Curr Opin Microbiol* 2018; **43**:162–8.

Jones KE, Patel NG, Levy MA et al. Global trends in emerging infectious diseases. *Nature* 2008; **451**:990–3.

Julliard R, Clavel J, Devictor V et al. Spatial segregation of specialists and generalists in bird communities. *Ecol Lett* 2006; **9**:1237–44.

Kim D, Hahn AS, Wu S-J et al. FragGeneScan-Plus for scalable high-throughput short-read open reading frame prediction. 2015 IEEE Conference on Computational Intelligence in Bioinformatics and Computational Biology (CIBCB). New-York, USA: IEEE, 2015, 1–8.

Kolodny O, Weinberg M, Reshef L et al. Coordinated change at the colony level in fruit bat fur microbiomes through time. *Nat Ecol Evol* 2019; **3**:116–24.

Lahti L, Shetty S et al. Tools for microbiome analysis in R. 2017. ht tp://microbiome.github.io/tutorials/ (21 October 2024, date last accessed).

Lai Y, Di Nardo A, Nakatsuji T et al. Commensal bacteria regulate Toll-like receptor 3-dependent inflammation after skin injury. *Nat Med* 2009; **15**:1377–82.

Langille MGI. Exploring linkages between taxonomic and functional profiles of the human microbiome. *MSystems* 2018; **3**:e00163–17.

Langwig KE, Frick WF, Bried JT et al. Sociality, density-dependence and microclimates determine the persistence of populations suffering from a novel fungal disease, white-nose syndrome. *Ecol Lett* 2012; **15**:1050–7.

Langwig KE, Frick WF, Reynolds R et al. Host and pathogen ecology drive the seasonal dynamics of a fungal disease, white-nose syndrome. *Proc R Soc B* 2015; **282**:20142335.

Langwig KE, Hoyt JR, Parise KL et al. Resistance in persisting bat populations after white-nose syndrome invasion. *Phil Trans R Soc B* 2017; **372**:20160044.

Langwig KE, Voyles J, Wilber MQ et al. Context-dependent conservation responses to emerging wildlife diseases. *Front Ecol Environ* 2015; **13**:195–202.

Legendre P, Anderson MJ. Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments. *Ecological Monographs* 1999; **69**:1–24.

Lemieux-Labonté V, Dorville NAS-Y, Willis CKR et al. Antifungal potential of the skin microbiota of hibernating big brown bats (*Eptesicus fuscus*) Infected with the causal agent of white-nose syndrome. *Front Microbiol* 2020; **11**:1776. https://doi.org/10.3389/fmicb.2020.01776/fmicb.2020.01776.

Lemieux-Labonté V, Simard A, Willis CKR et al. Enrichment of beneficial bacteria in the skin microbiota of bats persisting with white-nose syndrome. *Microbiome* 2017; **5**:115.

Lemieux-Labonté V, Tomas N, Shapiro BJ et al. Environment and host species shape the skin microbiome of captive neotropical bats. *PeerJ* 2016; **4**:e2430.

Lilley TM, Prokkola JM, Johnson JS et al. Immune responses in hibernating little brown myotis (*Myotis lucifugus*) with white-nose syndrome. *Proc R Soc B* 2017; **284**:20162232.

Lorch JM, Meteyer CU, Behr MJ et al. Experimental infection of bats with *Geomyces destructans* causes white-nose syndrome. *Nature* 2011; **480**:376–8.

Louca S, Parfrey LW, Doebeli M. Decoupling function and taxonomy in the global ocean microbiome. *Science* 2016; **353**:1272–7.

Lutz HL, Jackson EW, Webala PW et al. Ecology and host identity outweigh evolutionary history in shaping the bat microbiome. *mSystems* 2019; **4**:e00511–19. https://doi.org/10.1128/mSystems.00511-19.

Ma X, Zhang Q, Zheng M et al. Microbial functional traits are sensitive indicators of mild disturbance by lamb grazing. *ISME J* 2019; **13**:1370–3.

Mandal S, Van Treuren W, White RA et al. Analysis of composition of microbiomes: a novel method for studying microbial composition. *Microb Ecol Health Dis* 2015; **26**:27663. https://doi.org/10.3402/mehd.v26.27663.

Marchesi JR, Ravel J. The vocabulary of microbiome research: a proposal. *Microbiome* 2015; **3**:31.

McGuire LP, Mayberry HW, Willis CKR. White-nose syndrome increases torpid metabolic rate and evaporative water loss in hibernating bats. *Am J Physiol Regul Integr Comp Physiol* 2017; **313**:R680–6.

McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. *PLoS One* 2013; **8**:e61217.

Meteyer CU, Buckles EL, Blehert DS et al. Histopathologic criteria to confirm white-nose syndrome in bats. *J Vet Diagn Invest* 2009; **21**:411–4.

Micalizzi EW, Mack JN, White GP et al. Microbial inhibitors of the fungus *Pseudogymnoascus destructans*, the causal agent of white-nose syndrome in bats. *PLoS One* 2017; **12**:e0179770.

Naik S, Bouladoux N, Wilhelm C et al. Compartmentalized control of skin immunity by resident commensals. *Science* 2012; **337**:1115–9.
Newman MEJ, Girvan M. Finding and evaluating community structure in networks. *Phys Rev E* 2004;**69**:026113.

Newman MEJ. The structure and function of complex networks. *SIAM Rev* 2003;**45**:167–256.

Oksanen J, Blanchet FG, Friendly M et al. vegan: community ecology package. 2019. https://cran.r-project.org/web/packages/vegan/vegan.pdf (21 October 2024, date last accessed).

Poisot T, Canard E, Mouquet N et al. A comparative study of ecological specialization estimators. *Methods Ecol Evol* 2012;**3**:537–44.

Poretsky R, Rodriguez-R LM, Luo C et al. Strengths and limitations of 16S rRNA gene amplicon sequencing in revealing temporal microbial community dynamics. *PLoS One* 2014;**9**:e93827.

Preheim SP, Perrotta AR, Friedman J et al. Computational methods for high-throughput comparative analyses of natural microbial communities. *Methods Enzymol* 2013;**531**:353–70.

R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. 2018. Vienna, Austria. https://www.R-project.org/.

Råberg L, Graham AL, Read AF. Decomposing health: tolerance and resistance to parasites in animals. *Phil Trans R Soc B* 2009;**364**:37–49.

Råberg L, Sim D, Read AF. Disentangling genetic variation for resistance and tolerance to infectious diseases in animals. *Science* 2007;**318**:812–4.

Ranjan R, Rani A, Metwally A et al. Analysis of the microbiome: advantages of whole genome shotgun versus 16S amplicon sequencing. *Biochem Biophys Res Commun* 2016;**469**:967–77.

Rebollar EA, Bridges T, Hughey MC et al. Integrating the role of antifungal bacteria into skin symbiotic communities of three Neotropical frog species. *ISME J* 2019;**13**:1763–75.

Reeder SM, Palmer JM, Prokkola JM et al. Pseudogymnoascus destructans transcriptome changes during white-nose syndrome infections. *Virulence* 2017;**8**:1695–707.

Salter SJ, Cox MJ, Turek EM et al. Reagent and laboratory contamination can critically impact sequence-based microbiome analyses. *BMC Biol* 2014;**12**:87.

Scheele BC, Pasmans F, Skerratt LF et al. Amphibian fungal zootic causes catastrophic and ongoing loss of biodiversity. *Science* 2019;**363**:1459–63.

Shannon CE. A mathematical theory of communication. *Bell Syst Tech J* 1948;**27**:379–423.

Svensson EI, Råberg L. Resistance and tolerance in animal enemy-victim coevolution. *Trends Ecol Evol* 2010;**25**:267–74.

Tatusov RL, Galperin MY, Natale DA et al. The COG database: a tool for genome-scale analysis of protein functions and evolution. *Nucleic Acids Res* 2000;**28**:33–6.

Terry C. bipartiteD3: interactive bipartite graphs. R package version 0.1.0. Website HttpsCRAN R-Proj Orgpackage BipartiteD3. 2018. https://cran.r-project.org/web/packages/bipartiteD3/bipartiteD3.pdf (21 October 2024, date last accessed).

Titball RW. Bacterial phospholipases C. *Microbiol Rev* 1993;**57**:347–66.

Tukey JW. Exploratory Data Analysis. Mass.: Addison-Wesley Publishing Company Reading, 1977.

UniProt. UniProt: A worldwide hub of protein knowledge. *Nucleic Acids Res* 2019;**47**:D506–15.

Vanderwolf KJ, McAlpine DF. Hibernacula microclimate and declines in overwintering bats during an outbreak of white-nose syndrome near the northern range limit of infection in North America. *Ecol Evol* 2021;**11**:2273–88.

Verant ML, Meteyer CU, Speakman JR et al. White-nose syndrome initiates a cascade of physiologic disturbances in the hibernating bat host. *BMC Physiol 2014*;**14**:10.

Walter J, Maldonado-Gómez MX, Martínez I. To engraft or not to engraft: an ecological framework for gut microbiome modulation with live microbes. *Curr Opin Biotechnol* 2018;**49**:129–39.

Warnecke L, Turner JM, Bollinger TK et al. Pathophysiology of white-nose syndrome in bats: a mechanistic model linking wing damage to mortality. *Biol Lett* 2013;**9**:20130177.

Winter AS, Hathaway JJM, Kimble JC et al. Skin and fur bacterial diversity and community structure on American southwestern bats: effects of habitat, geography and bat traits. *PeerJ* 2017;**5**:e3944.

Woo V, Eshleman EM, Rice T et al. Microbiota inhibit epithelial pathogen adherence by epigenetically regulating C-type lectin expression. *Front Immunol* 2019;**10**:928.

Woodhams DC, Brandt H, Baumgartner S et al. Interacting symbionts and immunity in the amphibian skin mucosome predict disease risk and probiotic effectiveness. *PLoS One* 2014;**9**:e96375.