Pseudogymnoascus destructans transcriptome changes during white-nose syndrome infections

ABSTRACT White nose syndrome (WNS) is caused by the psychrophilic fungus Pseudogymnoascus destructans that can grow in the environment saprotrophically or parasitically by infecting hibernating bats. Infections are pathological in many species of North American bats, disrupting hibernation and causing mortality. To determine what fungal pathways are involved in infection of living tissue, we examined fungal gene expression using RNA-Seq. We compared P. destructans gene expression when grown in culture to that during infection of a North American bat species, Myotis lucifugus, that shows high WNS mortality. Cultured P. destructans was grown at 10 to 14 C and P. destructans growing in vivo was presumably exposed to temperatures ranging from 4 to 8 C during torpor and up to 37 C during periodic arousals. We found that when P. destructans is causing WNS, the most significant differentially expressed genes were involved in heat shock responses, cell wall remodeling, and micronutrient acquisition. These results indicate that this fungal pathogen responds to host-pathogen interactions by regulating gene expression in ways that may contribute to evasion of host responses. Alterations in fungal cell wall structures could allow P. destructans to avoid detection by host pattern recognition receptors and antibody responses. This study has also identified several fungal pathways upregulated during WNS infection that may be candidates for mitigating infection pathology. By identifying host-specific pathogen responses, these observations have important implications for host-pathogen evolutionary relationships in WNS and other fungal diseases.


Introduction
Fungal pathogens have emerged as major threats to biodiversity 1 and human health. 2 The diversity of these infectious eukaryotes and their hosts present new challenges in characterizing the interactions between host, pathogen, and the environment that lead to pathogenesis. One successful approach is to use systems biology to compare whole-transcriptome changes in gene expression between the pathogen infecting the host, the host without the pathogen, and the pathogen without the host. [3][4][5] This dual RNA-Seq approach can be used to identify genetic factors from the pathogen that contribute to host colonization and manipulation of host-pathogen interactions.
Among fungal emerging infectious diseases, white-nose syndrome (WNS) in bats has spread from Eurasia, where it is endemic, to North America, [6][7][8] where it is decimating several species of hibernating bats. Susceptible species, such as the little brown myotis (Myotis lucifugus) have shown population declines up to 90% in affected hibernacula. [9][10][11] WNS is caused by Pseudogymnoascus destructans, a psychrophilic fungus that grows in cold hibernacula and causes cutaneous infections in bats while they hibernate. During WNS, P. destructans invades the skin tissue, forming subcutaneous lesions identified as cupping erosions by histopathology. 12 The infection disrupts the hibernation behavior of susceptible bats and leads to more frequent arousals from torpor, premature energy depletion, electrolyte imbalance, and death. [13][14][15][16] WNS does not affect all species of bats equally. Many, but notall,NorthAmericanspeciesarebeingseverelyaffected, 17,18 while most European bats can host P. destructans infections, but have low mortality from WNS. [19][20][21][22] Coevolution of P. destructans and Eurasian bats, such as Daubenton's myotis (M. daubentonii), appears to have adapted these populations to a commensal or parasitic relationship with lower pathology. 8 North American bats, on the other hand, have yet to benefit from such selection against extirpation of the host species 23 and some species face the possibility of regional extinctions. 10,18,24 The virulence of the P. destructans infection is controlledbyacombinationoftheenvironment(i.e.,temperature and humidity of the hibernaculum), the host (and the host'sresponsetoinfection),andthepathogen(andthepathogen'sresponsetothehost). 25 Inthisstudy,wefocusonthethird component of this epidemiological triangle by dissecting the genetic components that allow P. destructans to infect hosts andbecomeavirulentpathogen.
Whether P. destructans remains a commensal parasite or becomes pathogenic is determined by host-pathogen interactions. 8,26 We have previously examined the host response of the WNS-susceptible M. lucifugus to P. destructans infection in the wing membrane and found robust gene expression changes in the host during hibernation. 27 We now shift our focus to characterize previously hypothesized virulence attributes of the fungus that include immune evasion, nutrient acquisition, stress responses, and tissue invasion. 28,29 We measured P. destructans gene expression at the whole-transcriptome level, comparing expression patterns between the fungus when growing in culture and when infecting a North American bat species.

Results
Two different groups of samples were used to measure gene expression in P. destructans for this study (Table 1). Gene expression during infection of M. lucifugus was measured in the MyLu samples of wing tissue from 6 individual P. destructans-infected M. lucifugus collected 60-120 minutes after arousal from hibernation in caves in Kentucky, USA. Gene expression in P. destructans during infection was compared with 4 samples from the 20631-21 strain of P. destructans growing in culture at 10-14 C for 23 d on Sabouraud dextrose agar plates ( Table 1).

Comparison of infected and uninfected bats
Prior to comparing the expression of P. destructans genes during host infection to those in culture, we confirmed that infection levels in host tissues were sufficient to measure pathogen gene expression by quantifying the number of RNA-Seq reads that mapped to the P. destructans transcriptome (Table S1). Compared to a group of samples from M. lucifugus not infected with P. destructans ( Figure S1), the samples from the infected bats from Kentucky showed significantly higher levels of P. destructans transcripts (t D 8.84, p < 0.00001). In the wing samples from infected bats, we found that 5990 § 324 P. destructans genes were expressed at a minimum count of 1, representing 63% of all P. destructans genes (Table S1). These samples expressed 13 512 § 357 M. lucifugus genes, representing 52% of all bat genes. Using a minimum of 1 count in any sample, the cultured samples expressed 8825 genes and the wing samples expressed 7264 genes of 9575 known P. destructans genes (Table S1). These results indicate that sufficient read depth was obtained in this data set to measure P. destructans gene expression, at least for the majority of genes.

Comparison of P. destructans gene expression during WNS and culture
Using both hierarchical clustering (Fig. 1A) and principal component analysis (Fig. 1B), we found that the patterns of P. destructans gene expression were similar in each group of samples (cultured or WNS). We observed a small batch effect between the cultured samples that were grown at different times and sequenced differently (Table 1). We also found that samples from bats KY19 and KY23, which came from a different cave in the same county as bats KY06, KY07, and KY11, 27 clustered separately from these samples and from sample KY39, which came from a different county. These results suggest that some of the differences in gene expression that we observe within the 2 groups could be due to variations in the environmental conditions or genetic differences between the P. destructans isolated growing in different hibernacula. However, the largest differences appear to be due to the different growth conditions between culture and growth on bats. We then compared P. destructans gene expression during WNS infection of M. lucifugus to the 20631-21 strain of P. destructans grown in culture using both edgeR (Fig. 2, 3A, 3B) and DESeq2 (Fig. 3C). Because of the lower depth of sequencing for the WNS samples, we then filtered the results to exclude any P. destructans genes that were not expressed in at least 2 of the 6 MyLu samples. With a cutoff of 0.001 for FDR and a 2-fold minimum change, similar results were obtained using these 2 different analysis methods (Fig. 3D), with the majority of the genes identified as differentially expressed by edgeR also being identified by DESeq2.
Using the subset of genes identified by both edgeR and DESeq2, 94 P. destructans genes were identified as more highly expressed during WNS infection of M. lucifugus, and 117 genes were more highly expressed in P. destructans growing in culture (Fig. 3D, Table S2). Using our Trinotate annotation, we identified 39 genes that showed significant changes in expression during WNS whose putative functions could contribute to virulence by affecting tissue invasion, the heat shock response, nutrient acquisition, immune evasion, and other pathways ( Table 2).
We specifically examined the expression levels of secreted proteases, because they have been implicated in the pathogenesis of WNS. 30,31 Protease genes were identified by homology and by PFAM analysis 32 and the expression of these genes was compared in the 5 culture samples and 6 M. lucifugus WNS samples (Table S2). Table 3 lists selected protease genes and demonstrates that the genes for subtilase-family proteases are more highly expressed during culture than during tissue invasion. Other proteases are highly expressed during host infection, such as VC83_01361, the P. destructans homolog of the Aspergillus fumigatus major allergen Aspf2, show lower gene expression when P. destructans is growing in culture.
To further explore the functional pathways that regulate infection, gene ontology enrichment analysis was performed using the genes identified by edgeR at a maximum FDR of 0.05 and minimum fold-change of 2. We examined the annotated functions of P. destructans genes upregulated in either M. lucifugus infections or in culture (Table 4). This analysis determined that several pathways involved in peptide and nitrogen metabolism are significantly enriched in P. destructans during infection (FDR < 0.05). While growing in culture, P. destructans showed enrichment of oxidation-reduction and transport pathways (FDR < 0.001) and depletion of other metabolic pathways (FDR < 0.05).

Discussion
We determined how parasitism affects the expression patterns of P. destructans genes by comparing expression levels between the fungus in culture and during host infection. We used dual RNA-Seq data and an approach that simultaneously mapped the reads to both host and pathogen transcriptomes followed by the removal of reads that mapped to host transcripts. This approach allowed for the estimation of expression levels of P. destructans genes with high levels of confidence by using RSEM to control for the uncertainty of multi-mapped reads. We compared gene expression changes of the cultured 20631-21 North American strain of P. destructans to infection of a na€ ıve North American species. Although the data set had limited read depth for P. destructans genes in the M. lucifugus samples, we observed significant differential gene expression in 211 genes, or 2.2% of the 9575 known P. destructans genes. This initial study has validated this approach to identifying changing patterns of pathogen gene expression. Future studies will be needed to overcome some of the limitations of the currently available data sets by using greater read depth for the dual RNA-Seq data, better matching environmental conditions in vitro to those in hibernacula, and using the identical isolate of P. destructans for both data sets. Future work could also compare changes in P. destructans gene expression during infection of North American or European bat species that show more resistance to WNS mortality than M. lucifugus. 8,19,20,33,34 As expected, we found that the transition from abiotic to parasitic growth was accompanied by many changes in P. destructans gene expression. Differences in temperature and humidity could also contribute to the differences in gene expression that we observed. Some of the gene expression changes are also presumably due to alterations in nutrient availability, such as the increased expression of lipase (VC83_00616) in vivo due to the high lipid content of mammalian skin. Although the cultured P. destructans was not grown on the same substrate that it would find in the environment, many of the gene expression changes that we observed appear consistent with adaptation to the host environment, rather than changes due to nutrient sources. For example, the increased expression of heat shock genes is consistent with the response to arousal from torpor to euthermic The mean expression level (log 2 counts per million (CPM)) and the fold change (log 2 FC) are shown for each gene. Genes more highly expressed in culture are on the upper half of the graph and those more highly expressed in M. lucifugus tissue in the lower half. Blue points indicate differential expression (FDR 0.001 determined by edgeR) that are expressed in at least 2 MyLu samples. Red points indicate significant differential expression for genes that were not expressed in at least 2 MyLu samples. An interactive version of this graph is available as Data Set S2. After unzipping File S2 and opening the html file in a web browser, hover over each point to view the annotation metadata for that gene and the expression level (in log 2 CPM) for each sample. Individual genes can be found by searching, for example by entering VC83_01361 in the search box. (b) Expression of P. destructans genes is compared by edgeR between culture and M. lucifugus infection with a volcano plot. The fold change (log 2 ) and the FDR (log 10 ) are shown for each gene. Genes more highly expressed in culture are on the right half of the graph and those more highly expressed in M. lucifugus tissue in the left half. Blue points indicate differential expression (FDR 0.001 determined by edgeR), with colors as for (a). An interactive version of this graph is available as Data Set S3 and can be manipulated as described above. (c) Expression of P. destructans genes is compared by DESeq2 between culture and M. lucifugus infection with an MA plot. The mean expression level and the fold change (log 2 ) are shown for each gene. The red line indicates equal expression and the blue line indicate a 2-fold change. Genes more highly expressed in culture are on the upper half of the graph and those more highly expressed in M. lucifugus tissue in the lower half. Red points indicate differential expression (FDR 0.001 determined by DESeq2). (d) A Venn diagram compares the number of P. destructans genes identified as differentially expressed by edgeR and DESeq2. The number of genes shared by edgeR and DESeq2, or unique to each method, are shown using a maximum FDR of 0.001 and minimum fold change of 2 for genes upregulated in M. lucifugus infections or upregulated in culture after removing genes not expressed in at least 2 of the MyLu samples. Table S2 lists results for all P. destructans genes.
body temperatures that occurred 60 to 120 minutes before collecting the M. lucifugus samples. 27 Correspondingly, a single sample from a bat that was allowed to become euthermic only briefly did not show upregulation of P. destructans heat shock genes (unpublished results). Thermal stress caused by a febrile response in the human host has been shown to activate a heat shock response in Candida albicans, preventing deleterious protein unfolding and aggregation. 35 This heat shock response could be important for fungal survival in our system, as bats arouse to euthermic temperatures several times throughout hibernation (thus several times throughout P. destructans infection), and susceptible populations arouse from torpor more frequently during WNS. 13,14,36 Consistent with a response of the pathogen to evade host immune recognition, we also found large increases and decreases in expression of genes involved in fungal cell wall structures ( Table 2). The fungal cell wall is composed of an inner layer of chitin, a middle layer of b-glucans, and an outer layer of mannose. The cell wall provides rigidity and structure, however is also highly dynamic. The pattern recognition receptor Dectin-1 has been shown to be a receptor for fungal 1,3-b glucans and 1,6-b glucans, 37,38 thus cell wall components serve to alert the mammalian immune system of a fungal pathogen. We have observed that Dectin-1 and several other C-type lectin domain family members are significantly upregulated in bat tissues infected with P. destructans. 27 Consistent with this host observation, we detected significant alterations in P. destructans enzymes predicted to be involved in fungal cell wall remodeling (Table 2). VC83_00788 and VC83_04729, homologs of Endochitinase 1, an enzyme which randomly cleaves and breaks down chitin, are upregulated 11.6 and 6.3-fold, respectively, while VC83_05104, a homolog of Chitin synthase  Enrichment (e) or purification (p) detected. Enrichment indicates that the GO category is more highly represented than expected by chance and purification indicates that the category is underrepresented. 2 Number of differentially expressed genes in this category compared with total differentially expressed genes. 3 Adjusted p-value of enrichment or purification after Benjamini-Hochberg FDR correction. 4 is downregulated 3.4-fold in P. destructans during infection compared with culture. Two homologs to Glucan endo-1,3-b glucosidases were differentially regulated; VC83_07327 was upregulated in P. destructans during infection while VC83_09076 was upregulated during culture. These enzymes presumably regulate cell wall b-glycan turnover and catabolism of b-glycans 39 by removal of non-reducing terminal glucosyl residues from saccharides and glycosides. Additionally, 3 Mannan endo-1,6-a mannosidases that were differentially expressed between P. destructans actively infecting a host and growing in culture (Table 2). Two were upregulated in culture (VC83_00261 and VC83_01650), and one was upregulated during WNS (VC83_07145). Mannan endo-1,6-a mannosidases are required for normal synthesis of the cell wall and alkaline pH-induced hypha formation, as well as being responsible for random hydrolysis of a-mannosidic linkages in unbranched mannans. 40 It is likely that the changes in Glucan endo-1,3-b glucosidase and Mannan endo-1,6-a mannosidase gene expression that we observed upon the switch from abiotic growth to host colonization leads to substantial alterations in the cell wall structures. The resulting differences in saccharide and glycoside branching patterns in the cell wall could make the pathogen less recognizable to the mammalian immune system.
Alternatively, these changes in cell wall enzyme gene expression could be due to changes in metabolic pathways that accompany the shift from abiotic to infectious niches. Different carbon sources can modulate cell wall structure and virulence in C. albicans. 41,42 It is possible that changes in cell wall structures are caused by differences in metabolism when infecting bats, rather than direct adaptation to the host.
Alterations in cell wall structures also accompany shifts in the morphological growth type of fungi, such as a shift from yeast to hyphal phase in C. albicans. 37 However, P. destructans grows vegetatively as hyphae on both Sabouraud's dextrose agar medium in culture, 43,44 and when forming cupping erosions in the wing tissue of the host. 12,19 Thus there is no difference in morphotype between our cultured and WNS P. destructans samples that might explain the dramatic alterations in expression of cell wall remodeling enzymes that we observed. Consequently, we propose that changes in the b¡glucan landscape on the fungal surface via cell wall remodeling are a mechanism of immune evasion for P. destructans, similar to other fungal pathogens. 45 Alterations of the cell wall during infection could explain the ineffectiveness of antibodies that recognize the cell wall of cultured P. destructans in providing protection from WNS. 46,47 These results may also explain why immunization with either cultured P. destructans or a b-glucan vaccine 48 did not affect the susceptibility of M. lucifugus to WNS (J. Johnson, J. McMichael, D. Reeder, and K. Field, unpublished). The antigens provided by these immunizations may not be present on the surface of P. destructans during infection because of changes in the cell wall structure that accompany the transition from abiotic to parasitic growth.
Because tissue invasion is a hallmark characteristic of P. destructans infections during WNS, we expected that expression of genes involved in degradation of the extracellular matrix would be upregulated. Unexpectedly, we found that the previously characterized subtilase-family of secreted proteases 30,31 showed lower expression in P. destructans during infection than in culture. Instead, the homolog of the A. fumigatus vacuolar protease, major allergen Aspf2, showed high levels of expression during infection of M. lucifugus and was significantly upregulated compared with culture conditions. This suggests that other proteases may be better targets for preventing colonization than the subtilase-family proteases, although the possible role of Aspf2 in tissue invasion remains unknown. It is also plausible that subtilase-family proteases are regulated at a post-transcriptional level or are used by the fungus primarily during initial colonization. Therefore, further proteomic and expression time-course experiments may prove useful to further dissect the infection. Nevertheless, the abundant expression of Aspf2, known to be an A. fumigatus allergen in humans, 49 suggests that further investigation of IgE-mediated allergic reactions during WNS may be warranted.
Infection of hosts was also associated with changes in expression for several genes involved in the transport or homeostasis of metal ions, including zinc, iron, and copper. This fungal response may be due to limited availability of some of these micronutrients in the host, which is likely sequestering metal ions as a form of nutritional immunity. 27,50 Changes in micronutrient acquisition gene expression appear to be associated with host colonization, including increased expression of the zinc transporter Zrt1, the copper homeostasis factor ATX1, and a putative copper transporter, as well as the unexpected loss of siderophore import using MirB. 51 Homeostasis of these micronutrients is essential for normal fungal metabolism and for the ability of the pathogen to respond to the oxidative stress activated by the host immune response. 50 However, our gene ontology analysis (Table 4) indicates that genes involved in oxidation-reduction pathways are more highly expressed during growth in culture than host colonization. Enrichment of pathways involving peptide metabolism and translation in P. destructans infecting bats (Table 4) indicates that host colonization demands higher levels of protein expression than abiotic growth. Competition between the host and pathogen for micronutrients and the generation of oxidative stress likely varies over the course of infection 52 and further study is needed to dissect this time course.
Together, these results provide a model of gene expression changes in P. destructans that accompany the transitions from abiotic to parasitic growth ( Figure 4). This model provides a framework to understand how the pathogen responds with phenotypic plasticity to the environment and its host to adopt a virulent phenotype. Our results also suggest approaches to minimize virulence and/or colonization by targeting immune evasion, micronutrient acquisition, tissue invasion, or the heat shock response. Efforts to understand why some species are more susceptible to WNS than others will require further examination of host-pathogen interactions to determine if the pathogen responds differently in hosts that exhibit lower WNS susceptibility.

Sample collection
Two different data sets were used for this study ( Table 1). The samples for the first data set (MyLu) consisted of wing tissue from 6 individual P. destructans-infected M. lucifugus (little brown myotis) collected 60-120 minutes after arousal to euthermy from hibernation from caves in Kentucky, USA, as described previously. 27 Hibernacula temperatures were 4-6 C at the time of collection and, based on our previous experience, we estimate that skin temperature varied between 4 and 8 C during torpor and up to 37 C during periodic arousals. The second data set was obtained from the North American 20631-21 strain of P. destructans growing in culture by D. Akiyoshi and A. Robbins (Department of Infectious Disease and Global Health, Cummings School of Veterinary Medicine, Tufts University). The 20631-21 strain of P. destructans was obtained from D. Blehert (National Wildlife Health Center, US. Geological Survey, Madison, WI, USA). The fungus was grown in culture at 10-14 C for 23 d on Sabouraud dextrose agar plates (BD Diagnostics, #221180) ( Table 1). Sabouraud dextrose agar contains nutrient sources of dextrose, pancreatic digest of casein, and peptic digest of animal tissue. RNA was isolated using a Qiagen RNeasy Lipid Tissue Kit after disruption of the cells using Zymos BashingBead Lysis Tubes and a bead beater on maximum speed for 30 sec for 3 times and then 20 sec once, with cooling on ice between each.

RNA sequencing
RNA sequencing was performed using Illumina sequencing as summarized in Table 1. Prior to analysis all data sets were quality trimmed using Trimmomatic v.0.35 53 with the parameters SLIDINGWINDOW:4:5 LEAD-ING:5 TRAILING:5 MINLEN:25. For samples with paired-end sequencing, only reads with both pairs remaining after trimming were used for further analysis. Analysis of the reads using FastQC v0.11.5 54 and the results of STAR mapping indicate that there are no significant differences in the quality of the RNA in any of the cultured samples from the MyLu samples.

Differential expression analysis
The quality trimmed reads were aligned using STAR v.2.5.1b 55 to the concatenated genomes of M. lucifugus and P. destructans. For M. lucifugus, we used genome assembly Myoluc2.0 and gene models from Ensembl release 84. 56 For P. destructans, we used the genome assembly and gene models from Drees et al.. 57 RSEM v1.2.29 58 was then used to apply an expectation maximization algorithm to predict gene expression counts for each transcript. The expected count matrix for all samples is available in Data Set S1. To determine if the number of reads mapped to P. destructans transcripts provided sufficient statistical power to detect differential expression of these genes, we used Scotty 59 to analyze the expected counts generated by RSEM. We determined that 65% of P. destructans genes expressed at a minimum of 4-fold change could be detected with a p-value cutoff of 0.05. Transcripts per million (TPM) was calculated by normalizing read counts for the length of each transcript and adjusting for the library size of mapped reads for each sample. 58 The M. lucifugus transcripts were then removed from the analysis and differential expression was determined using only P. destructans transcripts.
Differential expression between conditions was determined using either DESeq2 v1.10.1 60 or edgeR v.3.12.1 61 after normalizing across samples using the trimmed mean of M-values (TMM) method 62 and a minimum expression level of 2 TPM combined across all samples. False discovery rate (FDR) was used to control for multiple comparisons using the Benjamini-Hochberg procedure. 63 Hierarchical clustering was performed using R stats package v3.3.1 with Pearson correlation complete-linkage clustering of Euclidean distances. Clustering was confirmed by bootstrap analysis using pvclust v2.0-0 64 at an a level of 99% and 100 000 iterations. Genes without expression (expected count < 1) in at least 2 MyLu samples were excluded from the final analysis. Annotations for each gene were determined by using Trinotate v3.0, NCBI BLAST v2.2.29C 65 with the UniProtKB/SwissProt database (E-value cutoff of 1 £ 10 ¡4 ), and InterProScan v.5.20-59.0. 66 Gene ontology annotations were extracted from the InterProScan results and gene ontology enrichment analysis was performed using GOATOOLS v0.6.9 67 with enrichment or purification measured by Fisher's exact test after FDR correction.

Disclosure of potential conflicts of interest
No potential conflicts of interest were disclosed.