Comparative host transcriptome in response to pathogenic fungi identifies common and species-specific transcriptional antifungal host response pathways

Graphical abstract


Introduction
Opportunistic fungal infections are worldwide causes of significant morbidity and mortality that equals either bacterial, viral, or parasitic infections [19,108]. Among fungal diseases, aspergillosis, candidiasis, and mucormycosis account for the majority of opportunistic infections in immunocompromised patients, with the notable exception of HIV infected patients who primarily suffer from cryptococcosis in addition to oral candidiasis. caused by the three major species Aspergillus fumigatus, Candida albicans, and Rhizopus oryzae, significantly increased in incidence over the past years from approximately 610.000 to 960.000 lifethreatening invasive infections annually, with mortality rates ranging between 30 and 95% collectively [14,19]. The number of immunocompromised patients is steadily increasing [19,35] due to evolving medical care with more invasive procedures, increased use of immunosuppressive drugs [26,75], and growing numbers of ICU admissions. These factors mutually contribute to the increasing incidence of aspergillosis, candidiasis, and mucormycosis.
Although C. albicans is a widespread commensal in the human gut [46,52], it can cause invasive candidiasis when the host immune system is impaired [55,78,106,111] and the microbiome is disrupted, for example, by the use of broad-spectrum antibiotics [11,117]. In contrast to C. albicans, Aspergillus species are saprotrophic fungi, which usually do not colonize the human host. However, all humans are exposed to the conidia of these environmental fungi on a daily basis, which can result in severe infections in immunocompromised hosts [57,86,107]. Compared to C. albicans and A. fumigatus, fungi of the genus Mucorales cause invasive disease in a more heterogeneous group of immunocompromised individuals, including diabetes and trauma patients [86]. Mucormycosis is considered one of the most aggressive forms of opportunistic fungal infection [19,72] and Rhizopus oryzae is the most commonly isolated species [85].
The outcome of opportunistic invasive fungal infections remains poor, and therapeutic approaches targeting the immune system have been recently recognized as a potentially lifesaving strategy for patients with invasive fungal infections [6,103]. To improve antifungal immunotherapy, a global insight into the protective antifungal host response is needed. A common denominator of these three phylogenetically distinct fungal species is that they can all cause infections in immunocompromised patients with similar predisposing factors, including neutropenia, myeloablative and immunosuppressive therapy [15,38,51,60,64].
We hypothesized that there is a core host response that governs resistance to these three opportunistic pathogens, which is debilitated in immunocompromised patients. Nevertheless, there are also intrinsic differences in the pathogenesis of candidiasis, aspergillosis, and mucormycosis among these patients. Speciesspecific host defense pathways in terms of pathogen recognition, cytokine release, and metabolism may be activated depending on the particular fungal species encountered by the immune system. The knowledge of these mechanisms will provide valuable insights for the selection of immunotherapies adapted against each species. To characterize the core host response and determine speciesspecific host responses against A. fumigatus, C. albicans, and R. oryzae, an unbiased transcriptional profiling approach was used to establish differentially expressed RNA transcripts in human peripheral blood mononuclear cells (PBMCs) in response to these 3 fungal species. Subsequently, the most relevant pathways core and species-specific pathways have been functionally validated.

Results
2.1. Transcriptional profiling of host response to A. fumigatus, C. albicans, and R. oryzae The genome-wide transcriptional response of human PBMCs in response to A. fumigatus, C. albicans, and R. oryzae was assessed using RNA sequencing (RNA-Seq). Unstimulated PBMCs served as reference, and the transcriptional profile was assessed after 4 h (4 h) and 24 h (24 h) of stimulation with inactivated yeast and conidia of the three fungal species. Transcripts that exhibited a Log 2 fold change >1 or <À1 compared to unstimulated PBMCs and an adjusted p-value of < 0.05 were considered differentially regulated. A. fumigatus differentially regulated 472 transcripts (290 up and 182 down) at 4 h, and 506 transcripts (467 up and 39 down) at 24 h (Fig. 1A, B). C. albicans differentially regulated 1055 transcripts (821 up and 234 down) at 4 h, and 2448 transcripts (1733 up and 715 down) at 24 h (Fig. 1A, B). R. oryzae differentially regulated a staggering amount of 10,104 transcripts (6043 up and 4061 down) at 4 h, and 9646 transcripts (4974 up and 4672 down) at 24 h (Fig. 1A, B). Unsupervised hierarchical clustering revealed that the transcriptional response of PBMCs from various donors show stimulus dependent clustering at both the 4 h and 24 h time points (Fig. 1C, D). Of note, the response induced by R. oryzae of a single donor at 24 h was observed to form a separate cluster (Fig. 1D). Following the observation that different expression profiles were induced by the three fungi in the unsupervised hierarchical cluster analysis, the inter-fungal heterogeneity in the response of PBMCs was investigated. Principle component analysis (PCA) revealed that the majority of the variance in transcription was stimulus-dependent, as demonstrated by the three distinct clusters corresponding to the three different fungal species. Of note, the response to A. fumigatus stimulation clustered closely to the unstimulated cells and the response to R. oryzae showed high inter-individual variation (Fig. 1E).

2.2.
Common transcriptional response to A. fumigatus, C. albicans, and R. oryzae Subsequently, the common transcriptional responses to the three different fungal species were investigated at both timepoints ( Fig. 2A). Despite major differences in the number of differentially regulated genes ( Fig. 1A-C), 82 genes (58 up-regulated and 24 down-regulated) were commonly regulated by all three fungal species after 4 h of stimulation, and 91 (73 up-regulated and 18 down-regulated) after 24 h of stimulation ( Fig. 2A* and B). Genes involved in host defense such as cytokine and chemokine signaling, phagosome maturation, as well as coagulation, and lncRNAs were overrepresented at both 4 h and 24 h. Pathway enrichment analysis using the commonly up-regulated protein-coding genes confirmed these findings (Fig. 2C). A significant enrichment of interleukin-10 signaling (p = 4.01 Â 10 À13 REACTOME), Chemokine receptors bind chemokines (p = 4.68 Â 10 -7 REACTOME), cytokine receptor binding (p = 2.43 Â 10 -9 GO Molecular Function), and Regulation of blood coagulation (p = 7.19 Â 10 -5 GO Biological Process) was observed at 4 h. Enriched up-regulated pathways at 24 h included Receptor CXCR2 binds ligands CXCL1 to 7 (p = 2.39 Â 10 -5 REACTOME), Iron uptake and transport (p = 6.68 Â 10 -4 REACTOME), sphingolipid catabolic process (p = 3.22 Â 10 -4 GO biological process) and PERK-mediated unfolded protein response (p = 1.27 Â 10 -4 GO biological process). Key coagulation genes (F3, PLEK, PLAU, and PLAUR) were up-regulated after 4 h of stimulation with all the three fungi, posing the coagulation cascade as an early core host response. Previous studies have demonstrated associations between coagulation and antifungal host defiance [31,37,102]. When exaggerated, this ''protective" host response can become detrimental, by causing fungal sepsis-induced coagulopathy, which can be thrombotic or hemorrhagic [56,74,82]. As the coagulation pathway is identified as a common response to phylogenetically diverse fungi, it provides a rationale for diagnostic assessment of coagulation in patients with invasive fungal infection, which may aid the allocation of thromboembolic protective therapeutic measures.
No significant enriched pathways were observed in the common down-regulated genes. Only 2 genes, namely CCR2 and CD101, were down-regulated both after 4 h and 24 h stimulation with all the three fungal species (Fig. 2B). Interestingly, using a whole blood infection model a recent study showed only a com-mon up-regulation of genes in response to four tested Candida species, while no common gene down-regulation was observed [49]. An overview of the common upregulated and downregulated genes can be found in the related Data in Brief article [20].
2.3. Functional validation of the core host response to A. fumigatus, C. albicans, and R. oryzae Despite not enriched as a pathway per se, numerous genes involved in phagocytosis, phagosome maturation, endocytosis, and actin remodeling were up-regulated in response to all three fungi. Plasma membrane phosphoinositides [109] and RAB GTPases [114] mediate engulfment and vesicle fusion dynamics with the phagosome to form the mature phagolysosome. Interestingly, RAB42 and RAB20, known to be involved in vesicular remodeling [114], were commonly up-regulated by the three fungi. To validate the involvement of these differentially regulated phagocytosis genes in antifungal host defense; genetic variation in the loci was associated with the risk of invasive candidiasis. To this end, a previously published genome wide association study (GWAS) on candidemia susceptibility [70] was used to test whether genetic variants in a 1 MB window around the genes of interest were asso-  ciated to candidemia susceptibility. Six intronic and intergenic genetic variants demonstrated significant association (p < 5 Â 10 -4 ) in proximity to the genes: ASGR2, LAMP1, RAB20, RAB42, GEM, and ATP6V1B2 (Table 1). The minor allele of rs4773114 increases susceptibility to candidemia (OR > 1), and the rest of the SNPs are protective against the disease (OR < 1) ( Fig. 3A-F, Table 1). Genetic variants within non-coding regions can regulate genes through different mechanisms. Intronic variants can affect mRNA splicing or could also be relevant as enhancers that could regulate expression [27], also known as expression quantitative trait loci (eQTL). To address this point, the QTLbase repository [119] has been consulted: the eQTL query from publicly available databases of the six discovered candidemia-associated SNPs have shown that rs4773114 is a known eQTL for RAB20 gene expression in blood (p = 0.031) and rs9994 influences the gene expression of LAMP1 in monocytes after LPS and IFNa exposure (p < 0.0001) (Supplementary Table 1).
To functionally compare the host response between the three fungal species, differences in the capacity of the three fungi to induce oxidative burst were evaluated in term of reactive oxygen species (ROS) production by the phagocytic NADPH oxidase complex. The Pentose Phosphate Pathway (PPP) is crucial for both fueling NADPH oxidase with NADPH to produce superoxide and for activating the cellular antioxidant systems [8,50]. On a transcriptional level, a significant upregulation of NCF1, encoding the p47 phox subunit of the NADPH oxidase complex was observed after 4 h stimulation with R. oryzae and at both time points for C. albicans. R. oryzae significantly up-regulated some PPP enzymes. After 24 h stimulation, both C. albicans and A. fumigatus significantly upregulated the NCF2, coding for the p67 phox subunit, while R. oryzae not only significantly down-regulated NCF2, but also all the other NADPH oxidase 2 subunits (Fig. 4A), as well as key enzymes of the PPP (PGLS, TKT, and RPIA). On a functional level, C. albicans most potently induced oxidative burst followed by A. fumigatus, while R. oryzae failed to induce any ROS release (Fig. 4B).
As cytokine and chemokine signaling were the key pathways induced as a core transcriptional response against A. fumigatus, C. albicans, and R. oryzae, the induction of cytokine and chemokines was further validated by systematically profiling the common features and differences in the cytokine response to these diverse fungi. PBMCs were stimulated with the three fungal species. In addition, swollen conidia were generated from A. fumigatus and R. oryzae spores, and hyphae from C. albicans yeast. Additional strains of the three fungal species were also assessed (Supplementary Fig. 1A).
In response to all the fungal species PBMCs released the proinflammatory cytokines MIP1a, TNF, IL-1b, GM-CSF, and IL-6 and the anti-inflammatory cytokine IL-1Ra. This is in line with the upregulation of the correspondent genes, with the exception of IL1RN for R. oryzae stimulation (Fig. 4D). High level of the abovementioned cytokines found also in supernatants from stimulation with diverse morphologies and strains of A. fumigatus, C. albicans and R. oryzae ( Supplementary Fig. 1A). Apart from MIP1a, which was released in response to all the three fungi, production of the abovementioned cytokines was higher in response to C. albicans and R. oryzae than A. fumigatus (Fig. 4C). IFNc was also released in response to all fungi, albeit at low level in response to A. fumigatus (Fig. 4C), which may relate to its poor transcriptional induction (Fig. 4D).
Interestingly, the chemokine IFN-gamma-inducible protein 10 (IP-10, also known as CXCL10) which is abundantly produced in the course of pulmonary viral infections [4] was released by unstimulated cells, but protein levels markedly decreased in PBMCs exposed to the three fungi (Fig. 4C, Supplementary Fig. 1A).
More detailed pathway visualization and descriptions of the pathways mentioned above can be found in the related Data in Brief article [20].

Antifungal response pathways
After identification of the genes and pathways that dictate the core transcriptional response to A. fumigatus, C. albicans, and R. oryzae, genes exclusively regulated by each of the fungi were investigated in more detail ( Fig. 2A) to identify fungal specific transcriptional responses in cytokine signaling, pattern recognition, and metabolic pathways. A detailed overview of the species-specific induced genes and enriched pathways can be found in the related Data in Brief article [20].
Pathway enrichment analysis of the protein-coding genes specifically up-regulated after C. albicans stimulation confirmed type I interferon signaling (4 h: p = 3.21 Â 10 -30 , GO; 24 h: p = 1.3 3 Â 10 -19 Reactome) as unique for the response to C. albicans. It has been previously shown that the type I IFN pathway is a characteristic transcriptional signature of anti-Candida host defense when compared to bacteria [95]. The comparison of the transcriptional response to A. fumigatus and R. oryzae underscores that the type I IFN pathway is also a characteristic cytokine signaling signature of anti-Candida host defense even when compared to other fungi (Fig. 4E).
Overall a limited transcriptional (Fig. 1A, B, E) and cytokine (Fig. 4C, D) in response to A. fumigatus was observed, which can be explained by the fungal properties. Structural components of the conidia, particularly hydrophobins [1,79] and melanin [83,98], mask pathogen associated molecular patterns (PAMPs) and thereby mediate immune evasion.
The host transcriptional response to R. oryzae is divergent from the other two fungi, in terms of the number of differentially regulated genes, and of the pathways modulated. The sheer amount of differentially expressed genes by R. oryzae was 21-fold (4 h) and 19-fold (24 h) higher compared to A. fumigatus, and 9-fold (4 h) and 3-fold (24 h) higher compared to C. albicans respectively. Pathway enrichment analysis of the genes up-regulated at both time points (4 h and 24 h) revealed significant enrichment of pathways relevant to ribosomal function, endoplasmic reticulum (ER) stress and RNA processing.
Stimulation with the three different fungal species did not impact cell viability (Fig. 4F, Supplementary Fig. 1).  Differentially expressed genes in response to the three fungi are found to be located in candidemia-associated loci (p < 0.0006); related to Table 1. Genetic association study in candidemia patient to validate the importance of phagocytosis-associated genes that were differentially regulated in PBMCs of healthy volunteers. (A) Regional association plot of candidemia-associated rs1876730, in close proximity to AGSR2 gene on chromosome 17. (B) Regional association plot for candidemia-associated SNP rs73210879, which was mapped in a window of 1 MB around ATP6V1B2 gene on chromosome 8. (C) Regional association plot for candidemia-associated SNP rs12542176, in close proximity to GEM gene on chromosome 8. (D) Regional association plot for candidemia-associated SNP rs9994, in close proximity to LAMP1 gene on chromosome 13. (E) Regional association plot for candidemia-associated SNP rs12562688, in close proximity to RAB42 gene on chromosome 1. SNPs are plotted as the -Log of the p-value. (F) Regional association plot for candidemia-associated SNP rs4773114, which was mapped in a window of 1 MB around RAB20 gene on chromosome 13. Local linkage disequilibrium (LD) structure is reflected by the plotted estimated recombination rates (from HapMap) in the region around the associated SNP (purple diamond) and its correlation proxies. The correlation of the lead SNP to other SNPs at the locus is indicated by color. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Pathways involving pattern recognition receptors
The initiation of the antifungal host response is mediated through the recognition of PAMPs by pattern recognition receptors (PRRs) of the host cells [100]. Comparative transcriptional analysis of PRR expression in PBMCs stimulated with the three fungal species revealed profound differences and peculiarities. CLEC7A, encoding the dectin-1 receptor that plays a key role in antifungal immunity, was only significantly up-regulated by R. oryzae after 4 h, but not by other fungi at both time points. C. albicans induced expression of various C-type lectin receptors (CLRs; CLEC2B, CLEC2D, CLEC4E, CLEC4G, CLEC5A, and CLEC6A). Interestingly, the nucleotide oligomerization domain-like receptors (NLRs) NOD2 and NLRP3 stand out as uniquely up-regulated by C. albicans and not by the other fungi, although they also play a key role in modulating host responses against A. fumigatus [22,43]. Likewise, TLR2 and TLR4 were selectively up-regulated by C. albicans. The fact that expression of the intracellular nucleic acid-recognizing Toll-like receptors (TLRs; TLR3, TLR7, and TLR8), NLR (NOD2), together with the RIGi/MDA5 helicases IFIH1 and DDX58 are specifically induced by C. albicans further underscores that recognition of this fungus induces PRR pathways that drive type I IFN signaling (Fig. 5A).
Noteworthy, R. oryzae stimulation induced an early selective upregulation of some CLRs (CLEC4A, CLEC7A, CLEC10A, CLEC11A, CLEC12A), which are thereafter down-regulated after 24 h. In addition, the NLRs NLRP6 and NLRP9 were specifically up-regulated by R. oryzae. The receptor that mediates invasion of R. oryzae in endothelial host cells, GRP78 [65], also a regulator of ER stress and of the unfolded protein response (UPR) [61], was significantly up-regulated at both time points in response to R. oryzae and, less potently, to C. albicans (Fig. 5A). GRP78 was expressed in several of the pathways that were specifically up-regulated in response to R. oryzae. On the contrary, another receptor involved in host tissue invasion for Rhizopus species, epidermal growth factor receptor (EGFR) [3], was not significantly modulated by R. oryzae (Fig. 5A).
Compared to the other two fungi, A. fumigatus stimulation did not induce expression of PRRs in PBMCs; except for TLR5 (at 4 h), CLEC1B and CLEC12A (at 24 h). While CLEC1B and CLEC12A were not specific for A. fumigatus, TLR5 was the only PRR gene specifically up-regulated by A. fumigatus stimulation, while downregulated in response to 24 h stimulation with C. albicans and R. oryzae (Fig. 5A). Possible links between TLR5 and host defense against A. fumigatus have been demonstrated [88]. To validate the increased TLR5 expression in response to A. fumigatus, PBMCs were pre-incubated with heat-killed A. fumigatus conidia to induce TLR5 expression prior to stimulation with the TLR5 ligand flagellin. The Flagellin/TLR5-induced TNF response was significantly enhanced upon previous exposure to A. fumigatus, which does not induce abundant TNF responses itself (Fig. 5B). To determine whether TLR5 also plays a functional role in immune responses to A. fumigatus, common genetic variations in TLR5 were analyzed for their impact on A. fumigatus-induced cytokine responses in a functional genomics cohort [43,76]. The non-synonymous SNP rs5744174 (Phe616Leu) in the TLR5 ectodomain, which abrogates flagellin-induced signaling [92], has been associated with urinary tract infections [47] and hepatitis B infection [113], but it never been studied in the context of A. fumigatus infection. A reduced capacity to release IL-6 and an increased IL-1Ra release was observed in A. fumigatus-stimulated PBMCs of individuals carrying the minor allele (C) (Fig. 5C). Cytokine levels in bronchoalveolar lavage (BAL) samples from hematopoietic stem cell transplant patients (HSCT) with invasive aspergillosis (IA) were stratified based on the TLR5 (rs5744174) donor genotypes. Patients carrying the C allele showed a trend towards increased IL-1b and decreased IL-1Ra BAL levels compared to carriers of the homozygous reference TT genotype (Fig. 5D) suggesting that this polymorphism (rs5744174) influences pulmonary inflammation during aspergillosis. These data, collectively, suggest that TLR5 and genetic variation in the TLR5 gene differentially influences cytokine response to A. fumigatus and pulmonary inflammation in IA patients. In line with our data, TLR5 is up-regulated in THP1 cells exposed to A. fumigatus, and TLR5 neutralization was observed to impair conidial killing [89]. Other SNPs in the TLR5 coding region have been associated with an increased risk of aspergillosis, for example, insertion of a STOP codon (Arg392Ter; rs5744168) increased the risk of aspergillosis in HSCT patients [45]. This supports the notion that, TLR5 should be considered as a candidate gene for susceptibility to IA in at-risk patients. Future studies are also warranted to elucidate the A. fumigatus ligand that TLR5 recognizes.
Strikingly, many PRRs, in particular TLR genes were significantly down regulated upon R. oryzae stimulation (Fig. 5A). Concordantly, also the Toll-like receptor-signaling pathway was down-regulated (p = 7.32 Â 10 -4 , WikiPathways). In line with the reduced expression, PBMCs exposed to R. oryzae demonstrate decreased responsiveness to subsequent stimulation with the TLR4 ligand LPS or the TLR5 ligand flagellin (Fig. 5E). The mechanisms driving the downregulation of TLR, signaling pathway, and cytokine genes warrants further in-depth investigation. Endoplasmic reticulum (ER) stress responses, with subsequent ER-specific mRNA degradation and reduced influx of nascent proteins into the ER lumen can result in massive gene downregulation [23,39,61], leading to a reduced expression of crucial host defense pathways. Stimulation of PBMCs with R. oryzae significantly induced genes related to the unfolded protein response and cellular response to stress, including the key player GRP78 (Fig. 5A). Similar to R. oryzae, PBMCs stimulated with the ER-stress inducer tunicamycin demonstrated a downregulation of TLR gene expression, with the exception of TLR2, and strongly up-regulated GRP78 (Fig. 5F). These data underscore a potential connection between ER-stress  responses and downregulation of key immune receptors such as TLRs.

Metabolic rewiring of host cells: R. oryzae induces a dysfunctional glycolysis
Innate immune responses are energy-demanding processes and consequently involve reprogramming of cellular metabolism to meet energy demand [17]. The recognition of PAMPs induces distinct metabolic programmes in monocytes varying on the activated PRR pathway [54]. Recent works have shown rewiring of host metabolism and pointing out a crucial role for glucose metabolism in response to C. albicans [33,105] and A. fumigatus [41]. Even though dysregulated glucose homoeostasis (e.g. hyperglycemia or diabetic ketoacidosis) increases risk for mucormycosis [10], the host metabolic response to R. oryzae has not been studied in detail. Moreover, effects of different fungal species on host metabolism have not been compared so far. The significant upregulation of OXPHOS genes from each of the mitochondrial complexes after 4 h and enrichment of the respiratory electron pathway (OXPHOS) (p = 4.6 Â 10 -9 , WikiPathways) suggested that R. oryzae stimulation specifically induces oxidative metabolism. To validate increased OXPHOS activity, PBMCs were stimulated with A. fumigatus, C. albicans, or R. oryzae and subsequently the oxygen consumption rate (OCR) and extracellular acidification rate (ECAR) were analyzed. In contrast to the gene expression, OCR in host cells was similar in cells stimulated by the three different fungal species (Fig. 6A). Unexpectedly, when glycolytic capacity was measured, a significant reduction was observed only in R. oryzae-stimulated PBMCs (Fig. 6B). This correlated with a significantly decreased transcription of the ratelimiting glycolysis enzymes hexokinase-1 (HK1), hexokinase-2 (HK2), phosphofruttokinase-1 (PFKL, PFKM), and the terminal enzyme in anaerobic glycolysis lactate dehydrogenase (LDHA) in response to R. oryzae stimulation, but not by the other fungal species (Fig. 6C). Interestingly, the downregulation of HK2 was also found after 4 h of tunicamycin treatment in PBMCs (Fig. 5D), suggesting a role for ER-stress response induced by R. oryzae.
Glucose metabolism, in particular generation of the NAPDH via the PPP, is crucial for the generation of NADPH-oxidase-derived ROS in monocytes in response to microbial stimuli [44]. When comparing the direct fungal-induced ROS response, a defective Rhizopus-induced ROS-production was observed (Fig. 4B). To investigate the impact of the host metabolic changes induced by diverse fungal species on ROS production, PBMCs were exposed to A. fumigatus, C. albicans, or R. oryzae and subsequently zymosan-induced oxidative burst was quantified. Pre-incubation with either A. fumigatus or C. albicans slightly attenuated the capacity of PBMCs to induce an oxidative burst in response to Zymosan (Fig. 6E), which may correspond to the NADPH consumed to mount an oxidative burst in response to the pathogen by itself (Fig. 4B). Exposure to R. oryzae completely abolished capacity of PBMCs to induce any ROS in response to zymosan (Fig. 6E), even though it did not induce ROS in PBMCs by itself (Fig. 4B). This is in line with the significant downregulation of genes from the NADPH oxidase complex and key PPP enzymes after R. oryzae stimulation (Fig. 4A).

Discussion
The primary objective of this study was to elucidate the core and species-specific transcriptional responses against the distinct opportunistic pathogenic fungi A. fumigatus, C. albicans, and R. oryzae. The identified core response consists of immune-related ''usual suspect" genes, involved in cytokine signaling, regulation of chemotaxis, signaling cascades, and phagocytosis/phagosome maturation. Additional genes expressed in the core response have been less thoroughly studied in the context of fungal immunology, such as vesicle remodeling, coagulation, and lncRNAs. Validation of cytokine and chemokine responses against the three fungi revealed homogeneity for several pro-inflammatory mediators, but the responses are dominated by species-specific peculiarities. Apart from the common core response, species-specific transcriptional responses were identified to the three different fungi in PRR pathways and immunometabolism. Collectively these data underscore that the host induces unique transcriptional responses to each of the three different fungi, but still exhibits a defined antifungal core host response.
The investigation of the core transcriptional response to A. fumigatus, C. albicans, and R. oryzae stimulation revealed that genes related to phagosome maturation and actin dynamics, which may contribute to phagocytosis and clearance, were commonly induced in response to all three fungi. The crucial role for phagocytosis in antifungal host defense is well established [5,36]. The upregulation of both RAB20 and RAB42 suggests a role for these RAB GTPases, which have not been previously associated with antifungal host defense. RAB42 can associate with the phagosome, but its function remains unclear [114]. RAB20 is an IFNc regulated RabGTPase involved in regulation of phagosome maturation [80,91]. Yet, RAB20 was demonstrated to be involved in killing of engulfed E. coli and macrophage phagolysosome acidification [118]. Interestingly, we showed that genetic variation in close vicinity of the phagocytosis/phagosome maturation-associated genes ASGR2, LAMP1, RAB42, GEM, ATP6V1B2, and RAB20 are associated with susceptibility to candidiasis. To explore the potential of SNP functionality in the context of fungal diseases further investigation is required. So far phagocytosis has been often studied on a species-specific basis: it has been shown that A. fumigatus [2,24], C. albicans [9,110], and R. oryzae [5,24] have unique features that impair phagosome maturation. Our data provides a rationale to study phagocytosis as a common denominator in antifungal host defense, which can be exploited to design immunotherapies potentially applicable to diverse fungal infections. This study highlights a role of IFNc, which has been shown to improve fungal killing of C. albicans [21,99] and A. fumigatus [7]. Interestingly, IFNc induces RAB20 and its endosomal association in macrophages [80], which points to RAB20 as a potential target via which IFNc mediates its positive effects on fungal killing. Assuming that RAB20 is a key player in phagocytosis, another immunotherapeutic strategy capable of improving phagocytosis is Vitamin D3, which can improve the antimicrobial capability of macrophages [25,63,77,94] and, interestingly, also induces RAB20 [104].
Functional validation of the common cytokine signaling signature revealed that MIP1a and IL-1Ra were released in response to all three fungi. These two cytokines could, therefore, potentially be exploited as a core antifungal immunotherapy. MIP1a is currently being tested in the cancer field [90]. IL-1Ra, by competitive binding of IL-1 type 1 receptor (IL-1R), can be an immunomodulatory strategy for fungal infections where disease and pathology are the result of inflammation-driven collateral damage (e.g. VVC, pulmonary aspergillosis, or influenza-associated pulmonary aspergillosis). Recombinant IL-1Ra is already in clinical use (Anakinra/KINERET) and has been proven to reduce inflammation in murine infection models mimicking candidiasis [16] and aspergillosis [42].
A comparative transcriptional approach coupled with functional validations, offers the possibility to investigate differences in key events in antifungal host defense in term of cytokine release, pathogen recognition, and metabolic remodeling between the three pathogenic fungal species. While investigating the fungalspecific cytokine signatures, our study confirmed a previously identified key player in the response to C. albicans. Type I IFN response related genes were specifically induced in response to C. albicans and not by the other fungi. While a previous study identified type I IFN signaling as C. albicans-specific in comparison to the response to bacteria [95], our research pinpoints that this pathway still is specific for the response to C. albicans, even when compared to other opportunistic pathogenic fungi. Strategies that modulate type I IFN pathway could be auspicious targets for candidiasis-directed immunotherapy. In support of this, experimental systemic candidiasis models showed protective neutrophil recruitment dependent on the type I IFN receptor IFNAR [30]. However, given that type I IFN also has been observed to promote immunopathology [69], careful evaluation of their use as an immunotherapeutic approach is warranted.
As for the metabolic remodeling of host cells, R. oryzae uniquely influences cell metabolism by causing a dysfunctional glycolysis: PBMCs exposed to R. oryzae exhibited a reduced glycolytic rate, which was associated with reduced expression of the ratelimiting glycolysis enzymes HK1, HK2, PKM, PFKM, PFKL, and LDHA. In contrast, an upregulation of the OXPHOS genes, a possible compensatory response, was observed. We hypothesized that the glycolytic shutdown in response to R. oryzae stimulation may additionally explain why poorly controlled diabetes mellitus and obesity are risk factors for mucormycosis [28,40,93,97]. Together with the high blood glucose levels and diabetic ketoacidosis (DKA) both of which induce GRP78 overexpression [5,65], host immune cells reduce glycolysis in response to R. oryzae by downregulating key glycolysis genes. This may worsen diabetic patients' glucose homeostasis, which is already impaired due to the inadequate regulation of glycolysis by insulin [18,112]. In light of the altered glucose metabolism in response to R. oryzae, the glucoselowering, insulin-sensitizing agent metformin, which promotes glycolysis [81] could be a promising therapeutic approach to change immunometabolism. In line with this, metformin treated obese flies prior to systemic Rhizopus infection show improved survival [93].
Despite the fact that we identified many pathways to be modulated in response to the fungal pathogens, only a fraction could be validated in the current study. We prioritized genes and pathways that may serve as potential targets for immune-directed therapy. Additionally, many interesting pathways emerged from this comparative transcriptomic analysis that warrant validation in future studies. Our study revealed many non-protein coding RNAs (lncRNAs) in the core fungal host response. However, whether and how these lncRNAs play a role in antifungal host defense remains to be elucidated. Many pathways induced by R. oryzae are associated to ER stress [34,116] whose induction in PBMCs lead to a downregulation of most TLR genes and an upregulation of GRP78, according to our results. GRP78 (HSPA5), the host receptor for R. oryzae, [65], is a key regulator of these mechanisms [61] and was among the genes enriched cellular stress response and UPR pathways. It is tempting to speculate that engagement of GRP78 [66] by R. oryzae may play a role for the activation of these pathways [58,115].
Another aspect that is of clinical interest is whether host characteristics (e.g. age, gender, diabetes score) influence the core and species-specific antifungal responses. Large cohort studies have revealed drastic influences of sex, age and other donor-specific parameters on circulating concentrations of anti-inflammatory mediators and monocyte-derived cytokine release [101]. Although the current study design does not facilitate such an analysis, detailed insights in the individual-specific characteristics of the core-antifungal host response could facilitate more personalized treatment strategies.
Collectively, this study identified genes and molecular pathways that comprise the core and unique antifungal host response against A. fumigatus, C. albicans, and R. oryzae. These insights can direct new studies on mechanisms that might eventually lead to novel host-directed treatment strategies.

Fungal strains
Several well-described fungal strains were selected for characterization of the fungal specific transcriptional response, A. fumigatus (Ku80 / Af1163 / CBS144.89), C. albicans (ATCC MYA-3573 / UC 820) [59], and R. oryzae (delemar) (RA 99-880 / ATCC MYA-4621) [68]. C. albicans was grown from glycerol stock on Sabouraud dextrose (SD) plates. For all experiments, a single colony was grown in Sabouraud broth at 37°C for 12 h to log phase while shaking at 120 rpm. The harvested C. albicans yeast were washed twice in sterile phosphate buffered saline (PBS) and counted using a Bürker hemocytometer. A. fumigatus and R. oryzae were grown from glycerol stock on potato dextrose agar for 7 days at 37°C. Conidia were harvested in sterile PBS-tween-20 (0.1%), subsequently washed twice in sterile PBS, and counted using a Bürker hemocytometer.
C. albicans hyphae were generated by inoculating RPMI-1640 culture medium adjusted to pH 6.4 with hydrochloric acid with a density of 1 Â 10 6 /mL Candida yeast and culturing overnight at 37°C while shaking at 120 rpm. Subsequently the hyphae were washed in PBS and brought to a concentration originating from 1 Â 10 5 /mL yeast for the experiments. A. fumigatus and R. oryzae swollen conidia were generated by inoculating RPMI-1640 culture medium with a density of 1 Â 10 6 /mL conidia, and culturing 4 h at 37°C while shaking at 120 rpm. Subsequently the swollen conidia were washed in PBS and brought to a concentration of 1 Â 10 7 /mL for the experiments. Fungi were inactivated for 30 min in 4% paraformaldehyde (A. fumigatus and R. oryzae) or heat-inactivated (C. albicans). For all experiments a concentration of 1 Â 10 6 /mL C. albicans yeast was used or 1 Â 10 7 /mL A. fumigatus or R. oryzae conidia.

Healthy volunteers and patients
For RNA-Seq experiments buffy coats were obtained from anonymized healthy donors after written consent (Sanquin Blood Bank, Nijmegen, the Netherlands). For validation experiments blood was similarly obtained from buffycoats or collected from healthy volunteers by venous blood puncture after informed consent was obtained. All experiments were performed and conducted in accordance to Good Clinical practice, the Declaration of Helsinki, and the approval of the Arnhem-Nijmegen Ethical Committee (nr.2010/104). For the validation of the findings related to several genes, healthy individuals that were recruited in the 200FG cohort study (www.humanfunctionalgenomics.org) [43,76]

RNA-Seq
PBMCs (n = 8 donors) were stimulated in 6 well plates (1 Â 10 7 cells/well) with 1 Â 10 7 /mL A. fumigatus (Ku80, PFA fixed) 1 Â 10 6 / mL C. albicans (UC820, heat killed), 1 Â 10 7 /mL R. oryzae (RA 99-880, PFA fixed), or were left unstimulated for 4 h or 24 h at 37°C with 5% CO 2 . After incubation, the supernatant was removed and the cells were lysed and RNA was isolated using the mirVANA RNA isolation kit (Applied Biosystems) according to the protocol supplied by the manufacturer, before RNA-Seq was performed. All samples were analyzed and profiled in one batch. Before data analysis, batch effect was carefully checked. From the PCA analysis, all the samples were well separated by time (1st PCA) and by stimulation (2nd PCA). No batch effect was detected. The RNA-Seq analysis of this dataset was performed as previously described [62]. Briefly, sequencing reads were mapped to the human genome using STAR (version 2.3.0) [32]. The aligner was provided with a file containing junctions from Ensembl GRCh37.71. HTSeq-count of the Python package HTSeq (version 0.5.4p3) was used (the HTSeq package, (http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html) to quantify the read counts per gene based on annotation version GRCh37.71, using the default union-counting mode. After quality control 2 donors stimulated for 4 h and three donors stimulated for 24 h with R. oryzae needed to be excluded from the analysis. Differentially expressed genes were identified by statistics analysis using DESeq2 package from bioconductor [67]. The statistically significant threshold (FDR p 0.05 and Fold Change ! 2) was applied. The data sets generated from this study have been deposited in GEO Series accession number GSE162746 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE162746).

Identification of core and specific host response using Venn diagram analysis
Genes were sorted based on a significance threshold of a corrected p-values < 0.05 and a Log 2 FC > 1 for up-regulated genes or a Log 2 FC of < À1 for down-regulated genes ( Fig. 1A-C red genes). Within this list of differentially regulated genes by the 3 fungi a VENN-analysis, using the online VENN analysis tool Venny (https://bioinfogp.cnb.csic.es/tools/venny/), genes commonly upregulated or down-regulated (significant for all 3 fungi, same direction and a Log 2 FC of > 0.9 or < À0.9) by all three fungi were identified as the core host response. Subsequently genes were identified that were induced or down-regulated for each fungus uniquely (significant for the specific fungus and a Log 2 FC of > 0.9 or < À0.9), last the sets of genes that were commonly induced or down-regulated by at least two of the fungi were identified (significant for 2 fungi, same direction and a Log 2 FC of > 0.9 or < -0.9), subsequently proportional Venn diagrams were drawn with the eulerAPE application v2.0.3 [73,87]. The less strict threshold to group the genes for the core and species-specific responses was selected to increase the likelihood of identifying enriched pathways in the common host response, and reducing potentially false-positive enriched pathways in the species-unique responses.

Pathway enrichment analysis
Pathway enrichment analysis was performed using Cytoscape software with the ClueGO v.2.5.5 and CluePedia v.1.5.5 plugin [13]. To interpret and visualize the functionally group terms in the form of gene networks and pathways we have used Reactome, Wikipathways, and Gene Ontology (GO) categories, by excluding the annotations with the IEA code (Inferred from Electronic Annotation, which are assigned automatically computationally inferred based in sequence similarity comparisons. Degree of functional enrichment was determined by sorting enriched terms based on a p-value threshold of < 0.05. To reduce the redundancy of GO terms we applied the GO term fusion of related terms with similar associated genes. We used GO tree intervals between levels 3 and 8 and a kappa score of 0.4. The statistical test used for the enrichment was based on a right-sided hypergeometric option. The hypergeometric test p-values are further corrected for multiple testing using the Benjamini Hochberg multiple testing correction [12]. For R. oryzae modulated genes due to the high number of modulated genes, we performed pathway analysis on the list of overlapping genes between 4 h and 24 h. The complete pathways analysis results and settings are available in the related Data in Brief article.

Cytokine profiling with Luminex
To validate homogeneity in cytokine responses as a common pathway, PBMCs were stimulated for 24 h with A. fumigatus, C. albi-cans and R. oryzae as described above, supernatants were collected after 24 h of stimulation, and stored at À20°C till analysis. Cytokine protein levels were detected with a 25-plex human cytokine panel (LHC0009M; Thermo Fischer Scientific) according to the protocol supplied by the manufacturer. Bead counts were read on a Luminex TM MAGPIX TM instrument (Thermo Fischer Scientific).

In vitro PBMC viability assay
To determine the effect of fungal stimuli on viability of PBMCs, apoptosis and necrosis were assessed using Annexin-V and propidium iodide (PI) staining after 24 h of stimulation with the various fungal species. The cells were transferred to a v-bottom plate, washed with PBA (PBS pH 7.4, 1% w/v bovine serum albumin), and CD45 stained was stained using anti-CD45-BV510 (clone HI30; Biolegend) for 25 min at 4°C in the dark. After washing twice with PBA, cells were stained with FITC-conjugated Annexin-V (Biovision) for 15 min at 4°C in the dark in the presence of 5 mM CaCl 2 . Immediately afterwards, PI (Invitrogen Molecular Probes) was added and the cells were incubated an additional 5 min at 4°C in the dark. The cells were measured using a CytoFLEX cytometer (Beckman Coulter). Data analyses were performed in FlowJo (vX.0.7). First, CD45 + cells were gated to exclude fungal particles from the analysis. In the CD45 + population single cells were selected by subsequent SSC/FSC and FSC-H/FSC-A gating. The percentage of cells positive for only FITC-Annexin-V were designated as apoptotic and cells additionally positive for PI were designated as necrotic. Etoposide-treated (50 mM; Sigma-Aldrich) and heatkilled cells were used as positive controls for apoptosis and necrosis, respectively.

SNP genotyping of IA patients
Genomic DNA was isolated from whole blood of patients enrolled in the IFIGEN study using the QIAcube automated system (Qiagen). Genotyping of rs5744174 in TLR5 was performed using KASPar assays (LGC Genomics) in an Applied Biosystems 7500 Fast Real-Time PCR system (Thermo Fisher Scientific), according to the manufacturer's instructions.

GWAS on candidemia susceptibility
The GWAS on candidemia susceptibility was previously described [48]. In short, upon quality control per SNP and sample, this GWAS was performed in a cohort of 161 candidemia cases and 152 disease-matched controls of European ancestry whose demographic and clinical characteristics have been previously described [53]. DNA was genotyped using Illumina HumanCoreExome-12 v1.0 and HumanCoreExome-24 v1.0 BeadChip SNP chips. Geno-types were imputed using the human reference consortium (HRC) panel [71] using the Michigan imputation server [29]. In total, 5,426,313 SNPs were tested for disease association using Fisher's exact test with PLINK v1.9 [84]. Detailed results and statistics can be found in Table 1.

BAL fluid collection
BAL specimens were collected from patients using a flexible fiberoptic bronchoscope following local anesthesia with 2% lidocaine (Xylocaine), when infection was clinically suspected. Samples were obtained by instillation of a 0.9% sterile saline solution (20 mL twice). The sampling area was determined based on the localization of lesion on chest imaging. BAL specimens with comparable recovery rates were used.

ELISA for cytokine in PBMCs supernatant and patient BAL
Cytokine levels were measured in the cell culture supernatants from PBMCs using commercially available ELISA assays (R&D systems) according to the protocol supplied by the manufacturer. For human cytokines, IL-1b, TNF, IL-6, and IL-1Ra were measured after 24 h of stimulation.
Cytokines were quantified in BAL samples using customized Human Premixed Multi-Analyte Kits (R&D Systems, MN, USA). All cytokine determinations were performed in duplicates, and concentrations were reported in pg/mL.

ER-stress induction in PBMCs and qPCR
PBMCs were incubated with 2 ug/mL tunicamycin in RPMI 1640 for 4 h. RNA was isolated using RNeasy Mini Kit (Qiagen), according to the manufacturer's instructions. For reverse transcriptionquantitative PCR (RT-qPCR), isolated RNA (500 ng) was treated with DNase I (Fermentas) and subsequently transcribed into cDNA using 0.5 mg Oligo(dT)12-18 Primer, 200 U Superscript TM III Reverse Transcriptase and 40 U RNaseOUT TM Recombinant RNase Inhibitor (Thermo Fischer Scientific). RT-qPCR was performed with GoTaq Ò qPCR Master Mix (Promega) in a CFX96 thermocycler (Bio-Rad) and the expression levels were normalized against b-actin. All the primers used are listed in Supplementary Table 2. 4.15. Seahorse metabolic analysis 2 Â 10 6 PBMCs were seeded in a polypropylene tube (BD Falcon) in 2 mL of RPMI 10% serum and stimulated as previously described for RNA-Seq for 12 h at 37°C, 5% CO 2 . After that, 4 Â 10 5 cells were plated to overnight-calibrated cartridges in assay medium (either DMEM supplemented with 2 mM glutamine, 11 mM glucose and 1 mM pyruvate, or DMEM supplemented with 1 mM glutamine depending on the assay [pH adjusted to 7.4]) and incubated for 1 h in a non-CO 2 -corrected incubator at 37°C. Oxygen consumption rate (OCR) and extracellular acidification rate (ECAR) were measured using Cell Mito Stress test (for OCR) kit and the Glycolysis Stress test (for ECAR) in an XFe96 Analyzer (Seahorse Bioscience), with final concentrations of 1 mM oligomycin, 1 mM FCCP, and 1.25/2.5 mM rotenone/antimycin A for the Mito Stress test, and 11 mM glucose, 1 mM oligomycin, and 22 mM 2-DG for the Glycolysis Stress test.

Statistical analysis
For the RNA-Seq data and the pathways analysis statistical details can be found in the dedicated method sections. When comparing two differentially stimulated groups, or groups with different genotypes, we used Student's t tests or the Mann-Whitney U