The White-Nose Syndrome Transcriptome: Activation of Anti-fungal Host Responses in Wing Tissue of Hibernating Little Brown Myotis

Summary

White-nose syndrome is the most devastating epizootic wildlife disease of mammals in history, having killed millions of hibernating bats in North America since 2007. We have used next-generation RNA sequencing to provide a survey of the gene expression changes that accompany this disease in the skin of bats infected with the causative fungus. We identified possible new mechanisms that may either provide protection or contribute to mortality, including inflammatory immune responses. Contrary to expectations that hibernation represents a period of dormancy, we found that gene expression pathways were responsive to the environment. We also examined which genes were expressed in the pathogen and identified several classes of genes that could contribute to the virulence of this disease. Gene expression changes in the host were associated with local inflammation despite the fact that the bats were hibernating. However, we found that hibernating bats with white-nose syndrome lack some of the responses known to defend other mammals from fungal infection. We propose that bats could be protected from white-nose syndrome if these responses could be established prior to hibernation or if treatments could block the virulence factors expressed by the pathogen.

Introduction

White-nose syndrome (WNS) is an epizootic disease that has killed millions of bats in North America [1, 2]. WNS is caused by the psychrophile Pseudogymnoascus destructans (Pd) (formerly Geomyces destructans), an ascomycete fungal pathogen [3–5] that affects bats during hibernation. Pd grows at temperatures between 2 and 18°C and can infect bats while they hibernate [4, 6]. Pd is invasive and damages the cutaneous tissues of bats, including the wing [7], forming characteristic cupping erosions that are diagnostic of Pd infection [8]. Mortality rates due to WNS vary by species. In the little brown myotis, Myotis lucifugus, the mortality rate is up to 91% in affected caves [9, 10] whereas WNS resistance has been reported in the big brown bat, Eptesicus fuscus [11]. Bats in Europe are exposed to endemic Pd, but do not exhibit WNS mortality and appear to be resistant to the disease [12], despite cutaneous invasion by Pd [13].

Cutaneous infection by Pd causes some species of bats to arouse more frequently from torpor [5, 14, 15]. Although hibernating mammals spend less than 1% of their time euthermic [16], they use up to 90% of their stored energy during these periods [17, 18]. Because each arousal in little brown myotis utilizes an estimated 108 mg of stored fat [18], the increase in arousal frequency caused by WNS explains 58% of the morbidity rate associated with Pd infection [14]. Other factors that are also associated WNS pathology include effects of Pd infection on the integrity of wing tissue [7, 19], electrolyte balance and hydration [7, 20, 21], chronic respiratory acidosis [22], oxidative stress [23], and immune function [24]. The relative importance of each of these mechanisms in causing death in WNS is not clear, and the most likely model that has emerged is a multi-stage progression of WNS with contributions of several of these factors [22]. Differences in susceptibility to WNS between species in North America may be explained in part by different responses to Pd infection such as changes in thermoregulatory behavior. Understanding host responses to Pd infection may provide insight that could be useful for improving survival of affected species.

Cutaneous fungal infections in mammals are first recognized by components of the innate immune system, including C-type lectin receptors and Toll-like receptors [25]. Conserved components of the fungal cell wall activate pattern recognition receptors on phagocytes such as neutrophils, macrophages, and dendritic cells, and on epithelial cells [26]. Activation of these cells can lead to induction of the inflammasome, the production of inflammatory cytokines, and generation of reactive oxygen species that can mediate fungal cell killing [25]. The importance of the innate immune response to the initial recognition of fungal infections is demonstrated by the observation that deficiencies in these signaling pathways can lead to chronic fungal infections in humans [27, 28]. In the absence of invasion, colonization by commensal fungi can be maintained through tolerance mechanisms mediated by interactions with dendritic cells and epithelial cells in the skin [29]. Local activation of innate immune pathways can slow the growth of invasive pathogenic fungi and promote tolerance, possibly leading to a commensal relationship with the fungus [30], but is not usually sufficient to clear infections. Clearance of infections typically requires T helper (Th) cells, as demonstrated by the susceptibility of patients with acquired immune deficiency syndrome, immunosuppressant therapy, or chemotherapy to fungal infections [31]. These T cell responses can be mediated by Th17 cells [32, 33] or, in some cases, Th1 cells [34], with Th2 responses typically associated with greater susceptibility [35]. Th17 responses can contribute to clearance of invasive fungal infections through the actions of IL-17A and IL-22 [36] and the further recruitment and activation of neutrophils [37]. These T cell subsets have not been well characterized in bats, but those T-cell mediated immune mechanisms that have been studied appear to be conserved between bats and other mammals [38–41].

Fungal infections in animals are typically life-threatening only upon suppression of adaptive immune responses in the host, such as when chytrid fungus (Batrachochytrium dendrobatidis) blocks lymphocyte-mediated inflammatory responses [42]. Hibernation produces a natural suppression of some immune responses in mammal species where it has been studied. During hibernation, when the conservation of energy is critical, certain immunological mechanisms are downregulated while others remain unaffected [43–51]. Changes during hibernation can include depressed antibody responses [44, 52], decreased ability of T and B lymphocytes to proliferate in response to challenge [53, 54], and reduced complement activity [47]. Hibernation does not affect all immune responses equally, as shown in thirteen-lined ground squirrels (Ictidomys tridecemlineatus) that have a suppressed T-independent antibody response but are capable of mounting a T cell-dependent response during hibernation [44]. Studies of transcriptome-wide changes during hibernation in squirrels [55–59] have shown expression changes in genes involved in metabolism, oxidative stress, protein folding, ischemia/hypoxia, and other processes, but these studies were not examining an active immune response. Hibernation is also known to affect the distribution of leukocytes [45, 60] and platelets [61]. However, we have an incomplete understanding of how hibernation affects the suppression, or subsequent recovery, of immune responses [43], or how immune physiology in bats during hibernation may differ from that of rodents.

The cost of immune suppression during torpor is presumably outweighed by the benefits of energy conservation because most pathogens are not capable of proliferating at the low body temperatures of hibernating animals. However, the psychrophilic nature of Pd allows it to infect bats within hibernacula [2, 4]. The brief euthermic bouts of hibernating bats are shorter than most other hibernating mammalian species [14, 62] and it may not be possible for a bat naïve to a pathogen to mount a primary immune response in the few hours that it is euthermic throughout the hibernation season. We have observed antibody responses to Pd in bats, but these responses are strongest in active bats exposed to Pd after emergence from hibernation [63]. Therefore, hibernating bats may keep pathogens in check by relying on hypothermia, innate immune responses, and/or memory immune responses. The psychrophilic nature of Pd overcomes the first of these barriers to infection and the difficulty in fighting fungal pathogens with innate mechanisms alone may allow Pd to proliferate and invade the cutaneous tissues of bats.

The WNS panzootic has created an urgent need to understand if North American bat populations can persist in the presence of the fungal pathogen [1, 10]. Understanding the complete array of host responses mounted by bats afflicted with WNS may help illuminate sources of variation in survival within and among bat species. To determine which host responses are activated by Pd infection, we measured transcriptome-wide gene expression levels in bat wing tissue from hibernating bats affected by WNS. Gene expression was compared to bats that were hibernated in captivity in the absence of Pd exposure. We hypothesized that Pd infection would cause changes in gene expression that would reveal physiological responses during WNS that might be either protective or pathological. By using next-generation RNA sequencing to examine transcriptome-wide gene expression changes we expected to discover consistent patterns of host responses that occur in Pd-infected tissues. Combined with changes in gene expression within the Pd pathogen, these results have provided a survey of the host and pathogen interactions occurring during WNS.

Results

Gene Expression Changes Revealed by Next Generation RNA Sequencing

To determine the host response mounted by little brown myotis to Pd during hibernation, we measured changes in gene expression at the whole transcriptome level. Wing tissue samples were obtained from hibernating little brown myotis with no known exposure to Pd and bats exhibiting physical signs of WNS, as shown in Table 1. Histopathology [8] and quantitative PCR (qPCR) for Pd [64] were used to confirm the WNS status of each bat (Table 1). Cupping erosions diagnostic of WNS were found on all 6 bats captured in Kentucky, but on none of the 5 bats from states negative for WNS at the time of capture. Low levels of neutrophilic inflammation were found in all 11 wing samples (Table 1; Infl), although this inflammation was not associated with sites of Pd infection. All 6 WNS-affected bats tested positive for Pd by qPCR, although the fungal load measured on wing swabs (Table 1; qPCR) did not correlate with the number of cupping erosions found by histology (Table 1; WNS). As previously shown [5, 14, 15], WNS-affected bats had significantly lower body condition (Table 1; SMI; p = 0.017, t = 2.9255, df = 9).

Next generation RNA sequencing (RNA-Seq) was performed using poly-A selected RNA isolated from each RNAlater-preserved wing tissue sample (S1 Table). Using expression levels of Pd-derived transcripts, we confirmed that all 6 WNS-affected bats had abundant expression of Pd genes. The Pd-derived transcripts were not present at significant levels in any of the 5 samples from unaffected bats (S2 Table; p = 2.2x10^-6, t = 21.5, df = 5.33), including the MN090 sample that had tested positive for Pd by qPCR in one of the two replicates (Table 1). Because high levels of differential expression of Pd transcripts would make it more difficult to detect significant changes in host gene expression, the assembly was filtered [65] to remove Pd-derived sequences. Comparison of the filtered assembly with the original revealed that removing the Pd sequences did not significantly decrease the completeness of the assembly (S3 Table) as determined by BUSCO [66]. This filtered assembly (S1 Dataset) was used to calculate differential expression in host genes between the unaffected and WNS-affected samples.

We compared host gene expression across all samples (S2 Dataset) using DESeq2 [67] to identify transcript clusters that were expressed at a minimum of 2-fold difference and significant at a false discovery rate (FDR) of 0.05 (S1 Fig). We found 1804 transcript clusters that were expressed at higher levels, and 1925 transcript clusters expressed at lower levels, in WNS-affected bat tissues (S4 Table). Hierarchical clustering (Fig 1) revealed that expression of these transcripts from all 5 bats without WNS was similar to each other. Gene expression in wing tissue from WNS-affected bats was different from unaffected bats and more similar to each other, as predicted. The normalized expression levels of the 3729 identified transcript clusters differentially expressed are listed in S4 Table.

Differential expression of individual gene isoforms was further analyzed using EBSeq [68], an empirical Bayesian approach to modeling gene expression. For each transcript cluster identified as differentially expressed by DESeq2, we used EBSeq to determine if any of the individual transcripts were differentially expressed at a posterior probability greater than 0.99 (S4 Table). Of the 3729 differentially expressed transcript clusters identified by DESeq2, EBSeq identified at least one differentially expressed transcript for 1427 (38% of total, 43% of upregulated genes and 33% of downregulated genes). These results indicate that differences in gene expression are likely due to alternative splicing or other isoform differences for many of the differentially expressed genes.

To annotate the functions of these genes and identify those likely to be involved with host responses to Pd infection, we used the Trinotate pipeline. BLAST was used to identify 1365 upregulated transcripts and 325 downregulated transcripts in WNS-affected tissues with significant homology to known genes from vertebrates in the Swissprot database. Of the 2295 remaining transcripts, 13 were mapped to genes from non-vertebrates in the Swissprot database, presumably due to environmental contamination or incomplete removal of Pd transcript sequences. Of the 2842 trinity transcript clusters without a BLASTx match in Swissprot, 2731 (96%) were found to align to sequences (e-value < 0.0001) in the little brown myotis genome. Of the aligned transcripts, 204 (7.4%) were found to correspond to previously identified non-coding RNA sequences. Of the 111 transcript clusters without a transcript that aligned to the little brown myotis genome or Swissprot, BLAST was used to align their transcripts to the UniRef90 database. We found that 7 genes aligned to vertebrate homologs, 9 aligned to fungal homologs, and 15 aligned to other metagenomic sequences. We were unable to identify homologous sequences for any transcripts from 80 (2.1%) of the transcript clusters that were differentially expressed.

Expression levels for the Swissprot-identified transcript clusters with the 100 lowest adjusted p values are shown in Fig 2 (see S4 Table for all results). Some of the differentially expressed genes with putative functions that were predicted to associate with host responses to a fungal pathogen are listed in Table 2. WNS caused dramatic changes in expression of genes involved in inflammation, immune responses, wound healing, metabolism, and oxidative stress, even though the bats were hibernating during the Pd infection. Most of these genes were upregulated in WNS-affected tissues, while a much smaller number of identified genes with putative functions in these categories were downregulated (Tables 2 and S4).

To determine if all 6 little brown myotis with WNS exhibit similar changes in gene expression, we performed clustering analysis of the differentially expressed transcripts (Fig 1). To confirm the significance of these patterns of gene expression, bootstrap analysis of clustering was performed [69]. The clustering of the unaffected samples together and the clustering of the WNS-affected samples together was verified with a confidence of 99% (Fig 3A). Principal component analysis was performed to better understand the relationships between the transcripts expressed in the 11 samples (Fig 3B). All 5 samples from unaffected bats were very similar based on the first three principal components identified, which account for 71% of the variance in these transcripts. The WNS-affected bat samples have more diverse gene expression (S5 Table) and PC1 (accounting for 44% of the variance) differentiates all 6 from the unaffected bat samples. The genes represented by PC1 include those that are more highly expressed in unaffected than WNS-affected wing tissue (Fig 2). PC2 (17% of the variance) and PC3 (10% of the variance) distinguish the KY19, KY23, and KY39 samples from the other two WNS-affected samples and from the unaffected samples. The rotation values of principal component analysis (S5 Table) reveal that inflammatory genes made the greatest contribution to PC2. Clustering analysis revealed diverse host responses among the bats infected with Pd.

Metabolic and Inflammatory Immune Pathways Associated with WNS

We next examined the functional pathways that were most affected in little brown myotis infected with Pd. For this gene ontology analysis, DESeq2 results on transcript isoforms were used with a higher FDR threshold of 0.1, as is typical for this type of analysis. From WNS-affected bat tissue, 3104 upregulated transcripts were aligned with BLAST to the human Uniprot database. Homologs for these transcripts were identified and a list of 1937 unique Ensembl IDs associated with upregulated genes was generated (S6 Table). GOrilla [70] was used to determine significantly upregulated gene ontology categories from the Uniprot GO ID database (Table 3 and S7 Table) and REVIGO [71] was used to visualize biological processes that were significantly overrepresented in the WNS-affected transcriptome (Fig 4). The functional analysis revealed that Pd infection increases expression of genes involved in metabolism, defense responses, and other pathways (Table 3).

For the transcripts that showed lower expression in WNS-affected tissue (Fig 2), the same analysis was performed. Of the 1152 identified transcripts that were downregulated at an FDR of 0.1, 694 homologous human genes were identified by BLAST and mapped to Ensembl gene IDs (S6 Table). GOrilla did not identify any Biological Process categories that were significantly downregulated in the WNS-affected bat tissue.

Host-Pathogen Interactions during WNS

To examine the gene expression of the Pd pathogen using a dual RNA-Seq approach [72], we separately generated a genome-guided Trinity assembly (S3 Dataset) with the Broad Institute G. destructans genome. The reads from each of the WNS-affected tissues were mapped onto this assembly with Bowtie and gene expression estimated using RSEM (S4 Dataset). Expression levels for the Pd genes with the greatest variance are shown in Fig 5. Hierarchical clustering (Fig 6) and principal component analysis (S2 Fig) of the differentially expressed transcripts indicated that Pd gene expression was most similar in the wing tissues from bats obtained from the same hibernaculum (Table 1). The expression patterns of Pd genes were more similar for KY06, KY07, and KY11, which corresponds to bats captured in Cave 1 in Kentucky, and for KY19 and KY23, which were captured from Cave 2.

The possible functions of the Pd genes expressed among the WNS-affected samples were analyzed by sequence homology. We first examined the expression levels of a family of secreted proteases that have been proposed to be involved in Pd virulence [73, 74] and found that these alkaline proteases were expressed by Pd in all 6 wing samples (Table 4). Destructin-2 was the most highly expressed isoform in all WNS-affected bat Pd samples.

We next examined the Pd transcript clusters for additional factors that could affect virulence. Alignment by BLAST to the Swissprot and Uniprot90 databases identified 12 056 transcripts with significant homology to known fungal genes (S8 Table). For the remaining 67 Pd transcript clusters, Trinotate was not able to identify known functional domains or signal peptides present in these previously uncharacterized Pd transcripts. The results from the BLAST alignment were examined for genes known to be involved in processes that could affect Pd virulence, such as secreted proteases [73–75], metal binding proteins [76], fungal cell wall remodeling [76, 77], and other virulence factors [75, 77, 78]. This analysis identified 46 Pd genes that could be involved in pathogenesis (Table 5), including additional secreted proteases that could be involved in tissue invasion.

Because the tissue samples were collected from bats from 6 different hibernacula for this study, it is possible that differences in host or pathogen gene expression reflect differences in the environmental conditions present in each location, including the microbiome. In addition, the housing of the unaffected bats in captivity for 13 weeks prior to analysis could also have affected the microbiome. To examine the differences in the skin microbiome between the bats, we used MG-RAST to identify the lowest common ancestor of metagenomic sequences present (S8 Table). Although there were some differences observed in the bacterial microbiomes present on the wings of the 11 bats, there were no significant changes between the WNS-affected and unaffected samples when bacteria were identified at the class level. Several strains of Pseudomonas fluorescens isolated from bat tissues have been identified with Pd growth inhibiting properties [79]. MG-RAST analysis showed that Pseudomonas species are present in all 11 samples (S8 Table). P. fluorescens transcripts represented 2.8±0.6% of transcripts identified from gammaproteobacteria and 0.40±0.05% of all bacteria on the wings of unaffected bats and 0.37±0.07% of all bacteria on WNS-affected bats. P. fluorescens was present on all little brown myotis sampled, but was rare and relative abundance was not statistically different between WNS-affected and unaffected bats (p = 0.49, t = -0.71, df = 9).

Discussion

The comparison of host gene expression between WNS-affected and unaffected little brown myotis clearly demonstrates that Pd infection causes physiological responses in wing tissue, where substantial fungal invasion of the skin occurs in WNS-affected bats [8]. The changes in transcript levels that we have observed indicate that host responses to fungal infection remain intact during hibernation and are similar to those observed during the initial stages of fungal infection in euthermic mammals [32]. These host responses include acute inflammation, wound healing, and metabolic changes. Pathogen gene expression varies among bats with WNS, suggesting host-pathogen interactions that mediate pathogenesis. Together, these results lay a foundation to determine which host and pathogen responses contribute to WNS resistance and susceptibility and identify targets to increase survival.

Host Response to Pd Infection

The gene expression changes we observed in the wing tissue of WNS-affected bats are similar to those observed in other cutaneous fungal infections [80]. Cutaneous Candida albicans infections in humans and mice typically initiate an immune response by activating pattern recognition receptors of the C-type lectin family [81–83] and the toll-like receptor family, both of which we found upregulated in WNS-affected bat wing tissue (S4 Table). These included C-type lectin domain (CLEC) family members CLEC4D (MCL), CLEC4E (MINCLE), CLEC7A (Dectin-1), CLEC6A (Dectin-2), and Toll-like receptor 9. In mice and humans, protective host responses to C. albicans are usually characterized by many of the same cytokines and chemokines [29] that we have found upgregulated in WNS-affected wing tissue, including the cytokines IL-1β, IL-6, G-CSF, IL-23A, and IL-17C. Little brown myotis infected with Pd are increasing transcription of the key genes necessary for initiating a host response that provides protection from fungal infection. This clearly demonstrates that hibernation does not prevent innate immune responses in bats infected with Pd and that, although they are not closely related to rodents and primates [41], bats respond to fungal infections similarly to these other mammals.

The responses to Pd infection within bat wing tissue may be mediated by keratinocytes in the epithelial tissue. Activation of pattern recognition receptors by fungal ligands is expected to induce keratinocytes to produce many of the cytokines that we have found upregulated at the transcript level in WNS-affected bat wing tissue [84]. In addition to the cytokines typically involved in C. albicans responses described above, keratinocytes are also known to express the chemokine Ccl2 and the cytokines IL-20 and IL-24 in response to pattern recognition receptor activation [85]. Keratinocytes and fibroblasts are also known to exhibit a paracrine loop of IL-1 and IL-6 activation [86] that enhances wound healing and host defense to microbial infection and we found evidence of IL-1 and IL-6 receptor activation in the increased RNA levels for transcription factor p65, NFκB, and P-selectin glycoprotein ligand 1 (S4 Table). Another important cytokine produced by epithelial cells in response to infection is IL-17C [87]. This is an atypical IL-17 family member that is expressed by epithelial cells and causes autocrine responses in the epithelial cells that also express the IL-17RA and IL-17RE heterodimeric IL-17 receptor [87]. The wing tissue transcriptomes from WNS-affected and unaffected bats show similar expression levels of both IL-17RA and IL-17RE (S2 Dataset) and would, therefore, be expected to be responsive to IL-17C. The gene ontology analysis also found evidence for functional enrichment of genes involved in keratinocyte differentiation, presumably due to wound healing responses. Keratinocytes or other epithelial cells in bat wing tissue appear to have responded to the invasion of the epidermis by fungal hyphae.

Genes for pro-inflammatory mediators characterized the innate immune response that we observed in the wing tissue of Pd infected bats. Under euthermic conditions this would be expected to provide protection by the recruitment of monocytes and neutrophils, mediated by G-CSF, IL-23A, Ccl2, IL-17C and IL-6 [88], and the initiation of an adaptive Th17 or Th1 response. However, under the constraints of hibernation, responses that require leukocyte migration do not appear to occur in Pd-infected bats. We do not find strong evidence of increased expression for genes characteristic of either innate or adaptive leukocytes, except for L-selectin, which is expressed on T cells, and CD177, which is expressed on neutrophils. Lower than expected levels of monocyte, neutrophil, Th1, and Th17 cell recruitment may be related to the sequestration of leukocytes during hibernation [45]. However, we have observed neutrophil recruitment in hibernating little brown myotis in response to another fungal infection (Table 1). In the histological examination of the current samples, we found neutrophilic inflammation in both WNS-affected and unaffected wing tissue (Table 1). However, this inflammation did not occur at the sites of Pd infection. Curiously, we found a significant increase in WNS-affected tissue for transcripts for CD3γ and CD45 that could be expressed by gamma-delta T cells or other innate lymphocytes that reside in the skin [89]. It is possible that Pd is specifically suppressing neutrophil and/or T cell recruitment by interfering with chemotactic signals, similar to the suppression of inflammatory immune responses during chytridiomycosis in amphibians [42]. However, analysis of tissue levels of the cytokines and chemokines is necessary to confirm the secretion of these proteins. Because neutrophils and T cells do not appear to be recruited to sites of Pd infection during hibernation, only local inflammatory mediators may be available and they appear to be unable to control the infection in little brown myotis.

In addition to immune responses, hibernating bats also respond to Pd infection in other ways. We found transcripts for proteins from many pathways involved in metabolism, signaling, gene expression, transport, migration, and differentiation that were altered in WNS-affected bats (Fig 4). We cannot exclude the possibility that some of these differences were due to the different hibernation conditions of the two groups of bats. However, the differential expression of the genes in these pathways demonstrates that they are subject to regulation during hibernation and can respond to infection, tissue damage, and/or environmental changes.

Host responses to fungal infection can be influenced by changes in the pathogen, including gene expression changes in the colonizing fungus, such as C. albicans [29]. We found significant variability in the gene expression by Pd, which is particularly interesting because all Pd in North America is presumed to be a clone of the same mating type [90]. The pathogen has adopted different gene expression profiles in the 6 bat tissues (Fig 5), perhaps in response to differences in the host environments. Correspondingly, host gene expression patterns also show differences between the WNS-affected tissue samples. Of particular interest is the observation that the cytokine and chemokine genes found in principal component 2 of our PCA analysis (Fig 3B and S5 Table) are expressed at very different levels in the 6 Pd-infected samples. From this study we cannot determine whether the differences in pathogen gene expression are driven by differences in the host environment or vice versa. Although all 6 WNS-affected bats had visible signs of WNS, had similar Pd burdens, and similar histopathology, it is possible that the differences in host or pathogen gene expression that we observed may have affected progression of WNS and survival.

Responses that May Contribute to WNS Mortality

Because the increased frequency of arousals from torpor appears to be a primary cause of WNS mortality [5, 14, 22], we considered possible mechanisms that could affect torpor bout length. The increased gene expression of IL-1, IL-6, and other pro-inflammatory cytokines mediates a local acute inflammatory response to Pd. These cytokines also have systemic effects that modify behavior and thermoregulation [91]. In addition to cytokine and chemokine transcript increases, we also found increased transcripts for the enzyme cyclooxygenase-2 (prostaglandin G/H synthase 2) and both secreted and cytosolic phospholipase A2 that form critical inflammatory lipid mediators such as prostaglandin H₂. The eicosanoids generated by these enzymes, along with the actions of the upregulated genes kallikrein-6 and cathepsin S, are expected to generate pain and itching by locally activating neuronal nociceptors [92, 93]. This, in turn, could affect torpor bout length and/or behavior during periodic arousals. Indeed, we have documented significantly more grooming in WNS-affected bats infected in the wild [94], although a different study on laboratory-infected bats did not find similar behavior changes [95]. Together, the upregulated genes will likely generate an inflammatory microenvironment within the wing that may contribute to the robust wound healing response that we observe in WNS-affected bats. However, inflammation can also play a detrimental role in some diseases [96]. Further tissue damage and subsequent wound healing occurs in surviving bats upon emergence from hibernation [19]. These local affects of inflammation (pain and itching) as well as systemic effects are likely to play a key role in WNS pathology.

In addition to the gene expression changes that may contribute to acute inflammation locally within the epithelial tissues invaded by Pd, the systemic release of febrile cytokines such as IL-6 could affect the signals that control hibernation arousal. However, an exogenous pyrogen, lipopolysaccharide, is not able to provoke arousals in hibernating golden-mantled ground squirrels [97], so it may be unlikely that inflammation or febrile cytokines can directly trigger arousal in WNS-affected bats. Intracerebroventricular injection of prostaglandin E₂ in golden-mantled ground squirrels induces arousal from torpor and a febrile response during an extended periodic arousal [97]. Our observation of increased expression of the enzyme that generates prostaglandin H₂ may provide a mechanism that explains the shortened torpor bouts in WNS-affected bats, if it can be shown that this enzyme is active in the tissue and produces enough prostaglandin H₂ to act systemically.

In addition to the changes in expression of genes involved in immune responses and wound healing, we also found significant changes in metabolic genes. We found evidence of gene expression changes consistent with increased fat metabolism, including changes in transcripts for apolipoproteins, lipid transport proteins, protein metabolism, and carbohydrate metabolism. Of particular interest, we found increases in the expression of hydroxycarboxylic acid receptors 2 and 3 that are known to mediate adiponectin secretion [98]. This suggests that infection with Pd may directly trigger changes in lipid and carbohydrate metabolism that contribute to WNS pathology. These changes stand in contrast to the changes that have been seen in the brain transcriptome of hibernating horseshoe bats (Rhinolophus ferrumequinum) [99] and the brain proteome of hibernating Rickett’s big-footed bats (Myotis pilosus) [100], which show decreased fat metabolism during torpor. By leading to premature depletion of fat stores, these gene expression changes could contribute to WNS mortality.

The other changes in host gene expression that we observed are consistent with a multi-stage progression model of WNS [22]. We also found support for changes in genes involved in oxidative stress [23] and body fluid levels, which may contribute to WNS progression. Together, the pattern of gene expression changes that we find in little brown myotis with WNS suggests that a combination of maladaptive responses may contribute to mortality. However, the number of upregulated genes involved in the acute inflammatory response suggests that excessive inflammation may also be a factor contributing to pathology even prior to emergence from hibernation when it is suspected to contribute to wing damage [101].

Implications for Future Studies

The changes in host transcript levels that we have found are presumably caused by physiological responses of the host to infection. However, caution must be used when extending these transcriptional responses to functional mechanisms because the current study does not measure protein or metabolite levels directly. Future studies will be necessary to determine which of the gene expression changes observed affect which host response mechanisms.

The little brown myotis chosen for the WNS-affected samples were exhibiting WNS pathology and appeared unlikely to survive at the time of sample collection. For this reason, it is presently uncertain which of the gene expression changes that we have observed are contributing to protection and which are pathological. Another factor that likely contributes to the variation in gene expression that we observed among the samples collected from free-ranging bats is the time since the most recent arousal from torpor. Prior to collection of each wing tissue sample, bats were artificially aroused for 30 to 120 minutes. This period of arousal is similar in duration to the natural arousals during hibernation for little brown myotis [14], and presumably of sufficient duration for some innate immune responses to occur and for transcript levels to be altered. One reason for this procedure was to avoid disparities between the elapsed time from the most recent arousal bout until tissue collection. For the WNS-affected bats we could not determine when the most recent natural arousal would have occurred, but it would have likely been more recently than in unaffected animals, as affected animals arouse from torpor more frequently [14]. In the current study we cannot resolve whether the changes in gene expression that we observed occurred during the most recent arousal, during previous periodic arousals, or during torpor. Future studies will be needed to determine which of the changes in gene expression that we observed during WNS in bats in the wild also vary in controlled captive hibernation conditions when prior arousal patterns are known. Further studies are also needed to compare the physiological responses in bats exhibiting WNS morbidity to responses in less susceptible bats, such as European species, North American species that are less susceptible like the big brown bat [11], and the remnant populations of little brown myotis that appear to have developed tolerance or resistance to Pd [1]. Such studies should point to a path forward for bats in North America to persist in a landscape where Pd is endemic.

Conclusions

Little brown myotis mount a host response to Pd infection during hibernation. Which components of this response are protective or contribute to WNS pathology remains to be resolved. The innate immune response we have observed would be expected to promote a Th17-directed adaptive immune response that could clear the infection. However, the energetic constraints of hibernation may prevent little brown myotis from execution of the Th17- and neutrophil-mediated phases of the immune response. This may lead to excessive inflammatory responses, either during hibernation or upon emergence. The changes in host gene expression that we observed demonstrate that during Pd infection, little brown myotis also alter other defense responses, metabolic pathways, and transcription. Numerous Pd genes that may contribute to virulence were identified and these represent potential pathogen responses to host defense. Hibernation does not prevent a host response to infection and a better understanding of the differences between host and pathogen responses in bats susceptible to WNS and those resistant may lead to ways for increasing survival.

Materials and Methods

Ethics Statement

This study was carried out on bats from non-endangered species in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All methods were approved by the Institutional Animal Care and Use Committee at Bucknell University (protocol DMR-016). Animals were humanely euthanized by isoflurane anesthesia overdose followed by decapitation. In Illinois, animal collection was conducted by state wildlife officials and a numbered permit was not required. Scientific collector’s permits were obtained in Michigan (SC1448), Minnesota (201174), and Kentucky (SC1411147).

Samples

We collected hibernating little brown myotis from cave or mine walls at the locations listed in Table 1. Bats collected from all locations are expected to be from the same genetic population of eastern little brown myotis [102]. For bats unaffected by WNS, little brown myotis were first swabbed on the left forearm for quantitative PCR analysis. After measurements were taken, bats were individually placed in cloth bags and hung in constant temperature thermoelectric coolers (Koolatron PC-3) maintained at ~7°C. Water-saturated sponges were placed in the bottom of each cooler to maintain humidity during transportation to Bucknell University. Bats were housed for 13 weeks in a Percival (model I36VLC8) environmental chamber with conditions set to 4°C and 95% relative humidity. Bats were provided water throughout hibernation. Bats were aroused from hibernation for 30–120 minutes prior to euthanasia. For WNS-affected bats, little brown myotis were collected in the field, measured, swabbed for quantitative PCR, and humanely euthanized after being aroused from hibernation for 60–120 minutes. Scaled mass index (SMI) was calculated using the formula (mass(in g))*(38.01/(forearm length(in mm))^1.406 [103]. Wing tissue was placed in formalin for histology and placed in RNAlater (Sigma-Aldrich) for gene expression analysis. RNAlater samples were stored at ambient temperature for up to 24 hours before long-term storage at -80°C. RNA was purified from 50 mg of wing tissue using a QIAGEN RNeasy Mini Kit. All samples used for RNA sequencing had RNA integrity values greater than 7.0 using an Agilent Bioanalyzer.

Verification of WNS Status

Wing skin tissue was removed from the bones of the arm and digits and rolled onto 2 cm paraffin wax logs. The logs were then fixed in 10% neutral buffered formalin for at least 24 hours. Each log was cut into 3 pieces that were processed into paraffin blocks overnight in a Tissue-Tek VIP processor (Sakura Finetek). The pieces were embedded in paraffin blocks, sectioned at 3 microns, and stained with periodic acid Schiff with a hematoxylin counterstain [8]. WNS lesions (Table 1; WNS) were identified as cupping erosions with fungal hyphae and conida present. Inflammatory foci (Table 1; Infl) were identified as clusters of infiltrating neutrophils and were not associated with the asymmetrical curved conidia of Pd.

To determine presence or absence of Pd on bats unaffected by WNS, each swab was tested twice by quantitative PCR [64] by Jeffrey T. Foster at University of New Hampshire. A cycle-threshold less than 40 was used as a positive result. One of the 5 unaffected bats had one positive and one negative test (Table 1), but histology (Table 1) and subsequent RNA sequencing determined this to most likely be a false positive (S2 Table; p = 2.2x10^-6). For bats affected by WNS, we performed quantitative PCR to measure the Pd load, in genomic equivalents normalized to swabs spiked with 10 000 Pd conidia, that were detected on each bat [15].

Next Generation RNA Sequencing

The Genome Sequencing and Analysis Facility at the University of Texas at Austin performed all library preparation and quality control procedures. Directional RNA libraries were prepared with poly-A mRNA enrichment, dUTP/UDG strand-specific labeling, fragmentation, and 200 base pair size selection. RNA-Seq was performed in two lanes of an Illumina HiSeq 2500 with 101 base pair length reads obtained.

Transcriptome Assemblies

The paired reads from all samples were preprocessed by removing adapters and using trimmomatic PE [104] with settings of Illumina clip:2:30:10, seed mismatches:2, palindrome threshold:30, clip threshold:10, leading:5, trailing:5, minlength:36. The remaining paired reads were then combined and Trinity (v2.0.4) was used in strand-specific mode (RF) to construct a de novo assembly [105]. K-mer in silico read normalization with maximum coverage of 50 resulted in 22 482 456 read pairs that were used for assembly out of 177 755 004 total. The assembly was then filtered to remove Pd sequences using the program Deconseq [65] with the Broad Institute Geomyces destructans genome 20631–21 used to identify pathogen sequences and with the little brown myotis genome (Myoluc2.0) used to retain host sequences. Bowtie 1.0.1 [106] was used to determine the number of reads that mapped to each transcript in the assembly.

Differential Expression

The script align_and_estimate_abundance.pl included in the Trinity v2.0.6 distribution [105] was used to estimate expression levels for each transcript. Bowtie 1.0.1 [106] was used to map reads (including unpaired reads after quality trimming) from each sample onto the assembly. RSEM v1.2.20 [107] was used to apply an expectation maximization algorithm to predict gene expression counts for each transcript. Expression levels are presented after trimmed mean of M-values (TMM) normalization in fragments per kilobase of transcript per million mapped reads (FPKM). DESeq2 v1.8.1 [67] was used to determine the probability of differential expression for each Trinity transcript cluster that had a minimum RSEM-estimated count, before normalization, of 5 across all samples. For DESeq2 analysis, the default values for removing outliers and filtering lowly expressed transcripts were used. An alpha value of 0.05 was used instead of the default of 0.1 to decrease the number of differentially expressed genes identified. Posterior probabilities of differential expression for individual transcript isoforms were estimated using a Bayesian approach with EBSeq v1.8.0 [68]. False discovery rate [108] was used to control for multiple comparisons. NCBI BLAST v2.2.29+ [109] was used to identify the highest-ranking match for each isoform in the UniProt Swissprot database (downloaded on Sep 17, 2014) with an e-value cutoff of 1x10^-5.

Hierarchical clustering of samples and genes was performed within R 3.1.2 using the hclust function with the complete linkage method. Bootstrap analysis of clustering was performed using the pvclust 1.3–2 package and 1000 replications [69]. Principal component analysis was performed using the prcomp function and visualized with the rgl 0.93.1098 package.

Gene Ontology

NCBI BLAST v2.2.29+ [109] was used with an e-value cutoff of 1x10^-5 to identify homologs in the Uniprot Swissprot human protein database (downloaded on Nov 25, 2014) for transcripts significantly upregulated in WNS-affected bat wing tissue with an FDR of less than 0.1 (in order to increase the number of genes prior to subsequent analysis with higher stringency FDR). Unique Ensembl gene IDs were identified for 1144 of the 1922 upregulated transcripts and 481 of the 1356 downregulated transcripts. GOrilla [70] was used with a p value cutoff of 0.001 to identify upregulated or downregulated biological processes by comparison to the background list of 12 828 human genes identified by BLAST in the Trinity assembly. Multiple testing correction [108] was used with an FDR cutoff of 0.01. Results were visualized as a treemap with REVIGO [71].

Pd Gene Analysis

Trinity v2.0.4 was used to generate a Pd assembly in genome-guided mode with jaccard clipping and using the Broad Institute G. destructans genome 20631–21. This assembly was used to assess pathogen gene expression in the samples from WNS-affected bats using RSEM v1.2.20 [107]. Trinotate v2 was used to annotate the Pd transcripts by using NCBI BLAST v2.2.29+ [109] and both the Swissprot and Uniref90 databases (downloaded on Sep 17, 2014).

Metagenome Analysis

Reads for each sample were analyzed using MG-RAST v.3.5 [110] to identify metagenomic sequences after filtering against the B. taurus genome (the taxonomically closest genome available for filtering). For assignment of organism abundance, the best hit classification was used with the M5NR database, maximum e-value cutoff of 1x10^-5, minimum identity cutoff of 60%, and minimum alignment length cutoff of 15.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14

Attachment 15