The Mycobiome of Bats in the American Southwest Is Structured by Geography, Bat Species, and Behavior

Bats are widespread mammals that play key roles in ecosystems as pollinators and insectivores. However, there is a paucity of information about bat-associated microbes, in particular their fungal communities, despite the important role microbes play in host health and overall host function. The emerging fungal disease, white-nose syndrome, presents a potential challenge to the bat microbiome and understanding healthy bat-associated taxa will provide valuable information about potential microbiome-pathogen interactions. To address this knowledge gap, we collected 174 bat fur/skin swabs from 14 species of bats captured in five locations in New Mexico and Arizona and used high-throughput sequencing of the fungal internal transcribed (ITS) region to characterize bat-associated fungal communities. Our results revealed a highly heterogeneous bat mycobiome that was structured by geography and bat species. Furthermore, our data suggest that bat-associated fungal communities are affected by bat foraging, indicating the bat skin microbiota is dynamic on short time scales. Finally, despite the strong effects of site and species, we found widespread and abundant taxa from several taxonomic groups including the genera Alternaria and Metschnikowia that have the potential to be inhibitory towards fungal and bacterial pathogens.


Introduction
The microbiome, or the microorganisms that live on and in organisms, can affect host development, behavior, metabolism, and inflammation processes [13]. In addition, host-associated microbial communities can provide protection from infection and invasion by foreign microbes by outcompeting foreign taxa or through the production of antimicrobial compounds, conferring pathogen resistance in addition to the host's natural defense [30]. Much of what we know about host-microbe dynamics is derived from tightly controlled laboratory or captive animal experiments (e.g., [11], and studies examining microbiomes outside of the lab have revealed dynamic and diverse communities [78]. Complex factors act to shape host-microbiome interactions including geography, host-microbe evolutionary relationships, host disease state, and microbe-microbe interactions on/in hosts [19,78]. On a global scale, host microbiomes are structured by microhabitat [70], climate conditions [38], and immune complexity of hosts [79].
Investigations into host-associated microbes have primarily focused on gut communities [14,33,45], while skin-associated microbiomes have received less attention [35,36,41,49]. The surface microbiome is a spatially and temporally dynamic system [5] that lies at the interface of host and environment, providing the first line of protection against foreign microbes that are not behaviorally avoided [8]. For bats, ecoregion characteristics, the dynamics of temperature during hibernation, physical interactions with other bat species, variable skin secretions, non-temperature related roosting site/hibernacula factors, feeding sites, and behavior are all potential factors shaping bat microbiomes [3,41,47,78]. The emergence of fungal skin infections including bat white-nose syndrome [25] provides a pressing need to understand the spatial and temporal forces that shape skin microbiomes and the role that bat microbiomes play in disease.
White-nose syndrome (hereafter WNS) is an introduced, emerging fungal disease of bats caused by the psychrophilic 1 3 fungus Pseudogymnoascus destructans. WNS has led to the mortality of millions of bats in the eastern North America and in 2016, WNS was documented near Seattle, WA in a little brown bat (Myotis lucifugus), a major geographic jump [45]. The fungus P. destructans was detected in 2017 in Texas (https:// tpwd. texas. gov/ newsm edia/ relea ses/? req= 20170 323c). Bats are typically exposed to P. destructans during roosting in hibernacula and the fungus infects the skin, muzzle, ears, and wings of bats. WNS has a high mortality rate in some species of bats (e.g., Myotis lucifugus) and in areas P. pestructans spreads to, population declines of 70% have been reported [57]. [72,73] predict that WNS will impact 12 of the 13 species based on a variety of traits previously studied, while predicting that Myotis velifer bats will likely be resistant. Despite the devastating nature of WNS, several microbial derived compounds have demonstrated efficacy in inhibiting P. destructans [16,17] with recent data demonstrating the beneficial nature of bat associated Streptomyces species in vitro [27]. The American Southwest is one of the most bat species rich areas of the USA, possessing over half of the 45 American bat species [6,78]. Given the high mortality rate and expanding footprint of WNS in the USA and Canada (https:// www. white noses yndro me. org/ where-is-wns), understanding the bat-associated microbiome in unaffected areas may provide valuable information about potential host-pathogen interactions as WNS spreads.
Studies investigating bat-associated microbial communities have primarily focused on bacterial communities of the skin [3,43,82] and gut [31,81]. Knowledge of skinassociated fungal taxa is limited [32,35,36,79], although mycobiome traits may predict disease vulnerability [72,73]. To address this knowledge gap, we surveyed skin and furassociated fungal communities from 174 bats collected from five sites in southeastern, central, and northwestern New Mexico and northwestern Arizona. We hypothesized that species and site, key drivers of microbial community composition in host-associated and environmental communities [69,82], would strongly structure bat-associated fungal communities. Our results demonstrate the dominating effect of geography and species on fungal community composition. Furthermore, analysis of widespread and highly abundant taxa suggests that in healthy bats, potentially fungal-inhibiting taxa are important members of the bat microbiome.

Sample Collection and Processing
We swabbed bats captured during post-hibernation by netting at surface sites near known bat roosts, and by plucking from cave walls across New Mexico and Arizona at the following NPS (National Park Service) and BLM (Bureau of Land Management) sites: in New Mexico at El Malpais National Monument, El Malpais Conservation Area, Bureau of Land Management (BLM) Caves 45 and 55, Fort Stanton-Snowy River Cave National Conservation Area, and Carlsbad Caverns National Park; and in Arizona in Grand Canyon-Parashant National Monument (SI Fig. 1). Sampling was performed from 2011 to 2014 during spring and summer months (March to August) in caves and at surface sites outside of the hibernaculum. Surface bats were captured using mist nets [37], while cave captures were conducted by plucking bats from cave walls using nitrile glove covered hands and changing gloves after each bat was handled. All captures were conducted using approved protocols listed on the Scientific Collecting After capture, bats were placed in individual, sterilized bags until swabbing, which was conducted near the site of capture, either in the cave or near the nets on the surface sites. Skin and fur of the bats were systematically swabbed immediately after capture, wearing new nitrile gloves that were changed between bats, and using a single sterile polyester fiber-tipped application swabs (Falcon) moistened with sterile double-distilled water or Ringer's solution. Each swab was placed in a sterile 1.7 ml snap-cap microcentrifuge tube containing 100 ml of RNAlater®, and immediately frozen in a liquid nitrogen dry shipper or placed on dry ice. Samples were transported to the University of New Mexico and stored in a − 80 °C freezer until they were sent for sequencing with a few days or weeks. We sampled 186 bats belonging to 14 species (SI Table 1) as part of a 16S rRNA gene survey [78]. Bat species was determined by bat biologists, each with decades of experience with southwestern bats, who did the captures. Their determination of bat species was based on their knowledge of these species, and the measurements of key traits of the bats that they took after swabbing were completed. These samples came from the five study locations in the Southwest noted above (SI Table 1). Once the bats were swabbed, a bat biologist collected standard data for each animal including: species, sex, wing condition, and ectoparasite presence [61]. The wings, muzzle, ears, and uropatagium were assessed for any tissue damage (necrosis), lesions, scarring, or skin mottling that might suggest infection by P. destructans [18,62].
Sequencing samples were sent to MR DNA Molecular Research LP, Shallowater, Texas (http:// www. mrdna lab. com/) for genomic DNA extraction and 454 sequencing diversity assays of the fungal internal transcribed spacer region (ITS). The 186 samples were sequenced in nine runs. Barcoded amplicon sequencing processes were performed by MR DNA® under the trademark service (bTEFAP®). The ITS1 and ITS4 primer sets were used in a singlestep 30 cycle PCR using the HotStarTaq Plus Master Mix Kit (Qiagen, Valencia, CA, USA) under the following conditions: 94 °C for 3 min, followed by 28 cycles (five cycles used on PCR products) of 94 °C for 30 s, 53 °C for 40 s, and 72 °C for 1 min, after which a final elongation step at 72 °C for 5 min was performed. Sequencing with the 27F primer was performed at MR DNA on a Roche 454 FLX titanium following the manufacturer's guidelines. During processing (extractions and PCR), samples were checked for fungal contamination with fungal primers and gel electrophoresis highlighted above and no contamination was detected.

Sequence and Data Analysis
Reads were demultiplexed and denoised by the sequencing facility (www. mrdna lab. com). Briefly, barcodes and primers were removed and reads < 200 bp, with ambiguous bases, and homopolymers runs exceeding 6 bp were removed. Reads were processed using the UPARSE pipeline [22] at 97% identity using a previously published workflow [63]. Reads were first clustered into operational taxonomic units (OTUs) against the UNITE database (Nilsson et al. 2018) and reads that appeared only once (singletons) were removed. Reads that failed to hit the database were then clustered in de novo mode and checked for chimeras with UCHIME [23].
Due to potential misclassification of fungal reads [49] by using a single database, we assigned taxonomy to OTUs with the RDP classifier [15] using the UNITE database, Warcup Fungal ITS training set 2 [76], and the NCBI ITS Database as a reference. Due to many taxonomic assignments not going beyond the family level for all databases, we compared the relative abundance of fungal families for each sample using a Spearman's correlation. We did not observe any statistical evidence for differences between databases for the abundance of fungal families, so we will only present and use identifications from the UNITE database in this manuscript.
We calculated Bray-Curtis similarity values in R [60] with the vegan package [54] and visualized the results with a non-metric multidimensional scaling plot. We assessed differences in composition with a non-parametric multivariate analysis of variance. To assess changes in beta diversity over geographic distance, we used linear regression in R to examine pair-wise Bray-Curtis similarity as a function of distance (in km), calculating pairwise distance in R with the geosphere package. Furthermore, we compared Bray-Curtis similarity values within and between species at a single site and in between sites to understand the role of bat species in structuring fungal communities. Significant differences in pairwise beta diversity values were calculated with a Welch's t-test with a Benjamini-Hochberg correction.
We calculated alpha diversity (Shannon Diversity and OTUs observed) in R with the vegan package [54]. We analyzed factors that influence microbial richness with a generalized linear mixed effects model in R with the rstanarm package [50]. We used a Gaussian model with a weakly informed priors with 10,000 iterations similar to previous work on the bacterial communities of bats sampled in this study [78]. The model included latitude as well as site, bat species, bat feeding behavior, bat sex, bat diet, bat size, and the location of capture (inside a cave or surface netted) as random effects. Convergence of Monte Carlo Markov Chain (MCMC) chains was evaluated with the Rhat statistic. Finally, we evaluated potential collinearity in the model by calculating variance inflation values (VIF). All VIF values were near one for comparisons, thus indicating the absence of collinearity.

Identifying Widespread and Abundant Taxa and FunGuild Analysis
To identify taxa both widespread (> 55% of samples) and highly abundant taxa (> 10% relative abundance), we constructed a percent occupancy graph. Finally, we examined the functional potential of fungal taxa by assigning taxa to putative functional guilds using FunGuild [51]. We assessed significant differences in the abundance of fungal guilds between sites and species of bats with a Welch's t-test. However, we found no significant differences and thus display data as aggregated percentages of total reads.

Community Composition
An NMDS of Bray-Curtis similarity values highlighted the significant role of sampling site (non-parametric PERMANOVA; p < 0.001, F = 24.96, R 2 = 0.555) in structuring fungal communities (Fig. 1A, Tables 1 and 2). To further elucidate the role of site in structuring fungal communities, we analyzed pair-wise Bray-Curtis dissimilarity values as a function of distance (Fig. 1B). We found a significant relationship (y = 0.0028x + 0.529, p < 0.001 for slope and intercept) between distance and changes on community composition, providing further support for the NMDS results. In addition to site, we found significant effects of location where the bat was caught (inside or outside of a cave); however, this effect explained less variance than site (R2 = 0.03, Table 1), suggesting a lesser effect.
PERMANOVA results revealed a significant effect (F = 29.43, p < 0.001, R2 = 0.39) of bat species on community composition. We further investigated this pattern by comparing pairwise Bray-Curtis dissimilarity values (Fig. 2). We found, when comparing bat species within a single site and not between sites, there were significantly lower dissimilarity values than when comparing between species, suggesting a role for bat species in structuring fungal communities.

Taxonomic Composition and Diversity
We first compared taxonomic assignments of fungal OTUs between three commonly used databases (SI Fig. 2). We found significant correlations of fungal family abundance (Spearman's rho > 0.85, p < 0.05) between all three databases. As such, for the remainder of the manuscript, we will only discuss taxonomic identifications from the UNITE Database. Taxonomic composition of fungi strongly varied as a function of both location and bat species (SI Fig. 2). Bat associated fungal communities were dominated by fungi from the phylum Ascomycota (mean relative abundance (RA) = 90.1%), Mortierellomycota (4.4%), and Basidiomycota (2.5%). We found that a relatively small proportion (2.1%) of taxa could not be assigned to a fungal phylum. At the class level, taxa from the classes Dothideomycetes were highly abundant across nearly all sites and species, and when these taxa were in low abundance, members of the Saccharomycetes and taxa with no assigned phylum replaced them. Finally, it is worth noting that while taxa from the class Dothideomycetes were highly abundant in our data, no taxa from the family of Pseudogymnoascus destructans (Pseudeurotiaceae), the causative agent of whitenose syndrome, were found.  Fungal diversity, as measured by the Shannon Diversity Index, ranged from 0.78 to 3.8 across all bats and sites (SI Fig. 3). When we had enough replication to test, no significant differences between species were observed (Welch's t-test, Benjamini-Hochberg corrected p > 0.05) and sites did not significantly differ in diversity. Mixed effects modeling of fungal richness (OTUs observed) revealed contrasting patterns of site and bat characteristics on fungal diversity (SI Table 2

Abundant Taxa and FunGuild Analysis
To determine taxa that were important to the bat mycobiome, we analyzed the abundance and distribution of taxa across the dataset (Figs. 3 and 4). We identified 55 OTUs that were found in > 55% of samples or in abundance > 10% relative abundance. The taxa were primarily from the phylum Ascomycota and of the 55 taxa, 26 could not be assigned to a functional group or taxonomy past the family level. Of these taxa, most possessed a saprotrophic or partially saprotrophic lifestyle and a colony morphology of a microfungus. Fungi from the genus Alternaria and Metschnikowia were common as well. To complement analysis of core taxa, FunGuild analysis of the whole dataset revealed that when a functional guild could be assigned, taxa were predominantly saprotrophic (decomposer), pathotrophic (pathogenic), or symbiotrophic (mutualistic).

Discussion
In recent years, host-associated microbial communities have received extensive attention due to their important role in host wellness [13]. Host-associated microbial taxa can modulate the host immune response through direct effects on the host (i.e., inflammatory response, [9,48], as well as through the production of antimicrobial compounds. Bacteria are commonly investigated components of host microbiomes; however, information about fungal-host interactions is lacking. In this study, we investigated the composition and diversity of fungi associated with the skin of 186 bats in the American Southwest. Our results revealed strong effects of geography and bat species, as well as a significant effect of location where the bats were caught (outside/inside a cave) on community composition.
This suggests that the mycobiome of bats is dynamic over the course of a day, despite no significant effect of seasonality (spring-fall) or year sampled. Host-associated microbial communities can be highly dynamic on short (day) or long (months to years) time scales due to natural variations including host activity, diet, or health [10]. While we have documented the variability of the bat mycobiome, more information is needed on the dynamics and immune interactions [4]. While FunGuild analysis suggests that most taxa employ a saprotrophic lifestyle, the functional capacity of these communities is not well understood, in particular with respect to interactions with white-nose syndrome. The strongest drivers of bat skin/fur mycobiome were geography and bat species. Host microbiomes commonly exhibit patterns of morphological, functional, and genetic differentiation concurrent with the distribution of their hosts [80]. Furthermore, environmental microbes display similar variations over spatial scales and, as with host-associated taxa, are governed by local environmental conditions that select for composition, diversity, and abundance. Despite a similar roosting environment (cave), the skin microbiome of bats was most strongly driven by the local abiotic conditions of each site, suggesting local cues and reservoirs of fungi drive fungal community composition of bat-associated fungi. A similar result was noted in recent work on bat skin fungal and bacterial communities [2,3,47,76,77,82]. In addition to geography, the species of bat also had a significant effect on fungal composition. Host microbiomes display host specificity at broad (family or greater, [35,69], narrow (e.g. genus, [34,59], and subspecies [7,67] taxonomic resolutions. Within a given site, bat species was a significant driver of fungal composition suggesting, despite no effect of bat size or sex, each species of bat harbors a distinct community. Physiological differences between bat species such as skin secretions [56], Pannkuk et al. 2021;[64], the physical properties of the skin [68], and host-microbe immune interactions [4,43] likely play a role in differentiating fungal communities among bat species. However, more information is needed to support this notion. Fig. 3 Plot of the mean log 10 abundance of OTUs and the percent occupancy (proportion of samples, the OTUs are found in). Points colored in black (n = 55) are found in > than 55% of samples or with greater mean abundance of 100 reads per sample (10% relative abundance) Fig. 4 Bar plot of the results from the FunGuild analysis of fungal trophic feeding guild. The numbers above each bar indicate the percentage of total reads belonging to each category. The "unknown" category is comprised of taxa that could not be assigned to a feeding guild due to poor taxonomic resolution. No significant difference for any category was found between sites or species, so aggregated values are presented A key service that host-associated microbiomes provide to their host is the protection against invasion by foreign and potentially pathogenic microorganisms. While these interactions can often be mediated by competition for niche space [58], host-associated taxa can produce antimicrobial compounds to inhibit or kill off potential invaders [16,17,30]. Analysis of widely distributed fungi revealed that several taxa, including fungal taxa from the genera Alternaria and Metschnikowia, have been demonstrated to have antimicrobial properties [20,39]. For example, members of the genus Metschnikowia display antagonism towards other fungi (Sisti et al. 2014) and Alternaria produce numerous toxic compounds [39]. The commonality of these taxa across a wide range of species, sites, and individuals suggests that they likely plat an important role in the bat skin mycobiome. Furthermore, a parallel study of the bacterial microbiome of the same group of bats [78] revealed Actinobacteria to be a dominant component of the microbiome. Actinobacteria have long been known to produce antimicrobial compounds (i.e., Streptomycin) including antifungal compounds [12,30,81]. In comparison to recent work on bat skin fungal communities, Vanerwolf et al. (2021) documented similar collection of fungi from the classes Dothideomycetes and Eurotiomycetes. Similar results were documented by Ange-Stark et al. [2] and Zhang et al. [82]. The widespread nature of these taxa across multiple regions and species of bats suggests a potentially important role for these fungal groups on bat skin.
Many of the taxa identified (e.g., fungi from the order Chaetothyriales) are closely related to recently identified fungi that display significant inhibitory activity against the amphibian killing fungus, Batrachochytrium dendrobatidis [34]. Furthermore, recent work by [72,73] revealed a handful of fungal taxa from the genus Aureobasidium (Class: Dothideomycetes, Family: Dothioraceae), a common group on bat skin, were capable of producing anti-fungal metabolites. While the read-length of ITS fragments from this study did not allow for high-levels of taxonomic classification, many of the bats sampled in this study had > 75% of their fungal communities from the class Dothideomycetes (SI Fig. 3). In addition, it is important to note that no taxa from the family of Pseudogymnoascus destructans (Pseudeurotiaceae), the causative agent of white-nose syndrome, were found. Furthermore, recent work in these cave systems detected ideal habitat (temp/humidity) for P. destructans in 50% of surveyed caves and species from the genus Pseudogymnoascus were detected with PCR/culturing in 70% of caves, but not in sampled soils. Overall, this suggests that while the conditions are right for the proliferation of P. destructans in the American Southwest, it simply has not taken off. Our lack of knowledge of where bats hibernate in the southwest and the wide distribution of bats across a variety of habitats may have prevented widespread detection of WNS. A recent discovery of white-nose syndrome symptoms in southwestern New Mexico in live bats indicates that the disease is present, but confirmation awaits a dead specimen. While more work is necessary to elucidate the presence of potential WNS inhibiting fungal taxa on wild bats, overall, our results suggest that bats, like other animals, may potentially recruit possible pathogen-inhibiting taxa.
In conclusion, our study has demonstrated the important role of geography and species in structuring fungal communities. Bat associated fungi were highly variable; however, we have identified cosmopolitan Ascomycetes that are likely important to the healthy bat mycobiome. Furthermore, our results demonstrate co-occurring bacterial and fungal taxa that display anti-microbial capabilities to both bacterial and fungal pathogens. Our results suggest that despite the heterogeneity of bat-associated fungal communities, bats harbor anti-microbial taxa that likely play key roles in host defense. As WNS progresses across the Western United States, it will be of great importance to understand the host-pathogenmicrobiome dynamics and the role that potential probiotic taxa play in the alleviation of WNS.
Author Contribution ASW and DEN collected and processed all samples in coordination with other local scientists and bat biologists. PJK and DCW performed all sequences/statistical analyses with input from ASW and DEN. All authors contributed to authorship of the manuscript.