Voriconazole Treatment Induces a Conserved Sterol/Pleiotropic Drug Resistance Regulatory Network, including an Alternative Ergosterol Biosynthesis Pathway, in the Clinically Important FSSC Species, Fusarium keratoplasticum

Fusarium keratoplasticum is the Fusarium species most commonly associated with human infections (fusariosis). Antifungal treatment of fusariosis is often hampered by limited treatment options due to resistance towards azole antifungals. The mechanisms of antifungal resistance and sterol biosynthesis in fusaria are poorly understood. Therefore, in this study we assessed the transcriptional response of F. keratoplasticum when exposed to voriconazole. Our results revealed a group of dramatically upregulated ergosterol biosynthesis gene duplicates, most notably erg6A (912-fold), cyp51A (52-fold) and ebp1 (20-fold), which are likely part of an alternative ergosterol biosynthesis salvage pathway. The presence of human cholesterol biosynthesis gene homologs in F. keratoplasticum (ebp1, dhcr7 and dhcr24_1, dhcr24_2 and dhcr24_3) suggests that additional sterol biosynthesis pathways may be induced in fusaria under other growth conditions or during host invasion. Voriconazole also induced the expression of a number of ABC efflux pumps. Further investigations suggested that the highly conserved master regulator of ergosterol biosynthesis, FkSR, and the pleiotropic drug resistance network that induces zinc-cluster transcription factor FkAtrR coordinate the response of FSSC species to azole antifungal exposure. In-depth genome mining also helped clarify the ergosterol biosynthesis pathways of moulds and provided a better understanding of antifungal drug resistance mechanisms in fusaria.


Introduction
Fusarium keratoplasticum belongs to the Fusarium solani species complex (FSSC), a genetically diverse mould genus within the fungal order Hypocreales and Sordariomycetes class. Fusaria can be found in various ecological niches such as soil, water, plants, and man-made habitats. Fusaria are plant pathogens that cause significant diseases in a variety of agriculturally important crops [1][2][3]. Fusaria also frequently cause invasive mould infections (IMIs) in humans [4], and invasive fusariosis (IF) has a high mortality rate (56%) [5], especially in neutropenic patients (up to 100%) [6]. Infections caused by FSSC are problematic as most are innately resistant to commonly prescribed azole antifungals such as voriconazole (VRC) [7,8] which is one of the recommended treatment options for IF in addition to liposomal amphotericin B formulations or a combination of both [9][10][11]. Excellent guidelines for the diagnosis and management of rare mould infections, including IF, have recently been provided by the One World-One Guideline initiative [11]. Most FSSC species are susceptible to amphotericin B [7,8], but the use of this antifungal can of whether they were clinical or environmental isolates. Most intriguingly, however, this association held true across species boundaries with all VRC-resistant FSSC isolates (i.e., 3 F. keratoplasticum and 6 F. suttonianum isolates) possessing the 23 bp cyp51A promoter deletion [25]. We also recently reported a 372-fold induction of the mRNA levels of the pleiotropic drug resistance (PDR) transporter abc1 in response to VRC exposure [46]. When abc1 was expressed in S. cerevisiae it conferred up to 1024-fold increased resistance towards a range of xenobiotics including azoles [46].
The present study describes the transcriptional response of the clinically relevant FSSC species F. keratoplasticum to VRC exposure. The upregulation of a group of duplicated ergosterol biosynthesis genes and the annotation of all possible sterol biosynthesis homologs/orthologs in the sequenced genomes of related FSSC species enabled the proposal of an alternative ergosterol biosynthesis salvage pathway(s) that involve human cholesterol biosynthesis orthologs. Key players of the pleiotropic drug resistance network involved in azole antifungal drug resistance that seem conserved throughout the entire FSSC were also upregulated. The phylogenetic relationship of the manually curated ergosterol biosynthesis genes of the FSSC to their A. fumigatus and F. graminearum counterparts revealed possible alternative ergosterol biosynthesis pathways for these two model species as well.

Fungal Isolate and Culture Conditions
F. keratoplasticum strain 2781 is a clinical isolate obtained from a human nail as part of routine diagnostic procedures at the Hospital Canselor Tuanku Muhriz [25]. F. keratoplasticum 2781 showed high MICs to itraconazole (ITC; MIC > 32 mg/L) and VRC (MIC > 32 mg/L) [25]. Although antifungal clinical breakpoints have not yet been established for fusaria, epidemiological cut-off values have been defined using the CLSI method [7]. We have previously shown that exposure to high but sub-growth inhibitory concentrations of VRC (16 mg/L) caused increased expression of cyp51A [25] and the PDR transporter abc1 [46]. F. keratoplasticum strain 2781 was initially grown on potato dextrose agar, PDA (Merck & Co., Kenilworth, NJ, USA) at 28 • C for a week. The preparation of conidia suspensions, the growth and VRC induction, cell harvest and total RNA extraction were conducted as previously described [25]. In short, 50 mL potato dextrose broth, PDB (Merck & Co., Kenilworth, NJ, USA) was inoculated with 10 µL of F. keratoplasticum conidia suspension (∼1 × 10 8 cfu/mL) and cells were grown to mid-logarithmic phase for 21 h [47]. Then, cells were incubated for an additional 80 min in the presence or absence of VRC (16 mg/L) followed by cell harvest. VRC was added to the cell suspensions as a 100× stock (1.6 mg/mL) dissolved in DMSO and the control cells were incubated with an equal volume of DMSO only. The experiment included three biological replicates.

Total RNA Extraction and Sequencing
Total RNA extraction (hot phenol method), DNase treatment (PureLink DNase kit; Invitrogen Inc., Carlsbad, CA, USA) and quality checks of RNA samples by agarose gel electrophoresis were conducted as previously described [25]. Total RNA samples (3 each of VRC treated and untreated control samples) were submitted to Genewiz Biotechnology Co. Ltd. (now Azenta Life Science), Suzhou, China, for transcriptome sequencing. Sample quality control checks were performed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and a NanoDrop spectrophotometer (Thermo Fisher Scientific Inc. Waltham, MA, USA). RNA sequencing (RNA-seq) of 150 bp paired-end reads with a total of~40 million reads was performed on NovaSeq 6000 (Illumina, San Diego, CA, USA).

De Novo Transcriptome Assembly and Analysis of Differentially Expressed Genes (DEGs)
Quality assessment of the sequencing data was evaluated using FastQC version 0.11.4 (Simon Andrew, Cambridge, UK, https://www.bioinformatics.babraham.ac.uk/projects/ fastqc/) [48] followed by the removal of adapter sequences using Cutadapt version 1.9.1 (Martin Marcel, https://cutadapt.readthedocs.io/en/v1.9.1/) [49]. The 5 and 3 end bases containing Ns or sequences with quality values < 20 were removed. Sequencing reads that were less than 75 bp long after trimming were also removed. De novo transcriptome assembly of the clean data was performed using Trinity version 2.2.0 (Brian J. Haas, https: //github.com/trinityrnaseq/trinityrnaseq/releases) [50] followed by clustering into long non-redundant unique sequences. The prediction of coding regions within the assembled transcripts was performed using the TransDecoder version 3.0.0 (Brian J. Haas, https: //github.com/TransDecoder/TransDecoder/releases?page=2) [51]. To determine their putative functions, transcripts were aligned with sequences in the following publicly available databases: NCBI non-redundant protein sequences (Nr), a manually annotated and reviewed protein sequence database (Swiss-Prot), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Cluster of Orthologous Groups of proteins (COG).
To generate read densities for gene expression analysis, sequencing reads were aligned to the assembled transcripts using Bowtie2 version 2.2.6 (Langmead et al., https://sourceforge.net/ projects/bowtie-bio/files/) [52] with default settings. The level of gene expression normalized to read density was calculated using RNA-seq by Expectation Maximization (RSEM) version 1.2.4 (Li et al., https://github.com/deweylab/RSEM/releases) [53] and quantified as fragments per kilobase of transcripts per million mapped reads (FPKM). The analysis of DEGs between the two groups of RNA samples was performed using the Bioconductor package DESeq2 version 1.6.3 (Love et al., https://bioconductor.org/packages/release/bioc/html/DESeq2.html) [54]. The results were further analysed to determine significant DEGs based on the criteria of fold change > 2 and a q-value (false discovery rate adjusted p-value) < 0.05.
In order to classify the biological functions of significant DEGs, GO functional enrichment analysis was performed using GOSeq [55]. Gene length and read count biases were included in the GO analysis using a filtering threshold for overrepresented sequences by applying a p-value ≤ 0.05.

Search for Homologs of F. keratoplasticum Genes in Other Fungi
We identified and manually curated all possible ergosterol and cholesterol biosynthesis homologs of A. fumigatus and F. graminearum and the three FSSC species F. vanettenii, F. solani FSSC5 and F. keratoplasticum 2781 by searching their entire genome or, in the case of F. keratoplasticum 2781, the entire transcriptome of cells grown either in the presence or absence of VRC using all known ergosterol or cholesterol biosynthesis protein sequences of A. fumigatus and S. cerevisiae or H. sapiens as query sequences (see Supplementary File S1).
Coding sequences (CDSs) within the F. keratoplasticum mRNA transcripts were used as query sequences with the online NCBI BLASTx to search the non-redundant protein sequence database (Nr) of Fusarium taxa (taxid: 5506). CDSs conserved across the Fusarium taxa were identified from the alignments and translated into protein sequences. The BLAST+ software version 2.13.0 [58] was used to conduct a local BLASTp search against the protein sequences of Fusarium solani FSSC 5 [59], F. vanettenii 77-13-4 [60], F. graminearum PH-1 [61], A. fumigatus Af293 [62] and S. cerevisiae S288C [63] downloaded from the publicly available MycoCosm [64] database (accessed on 29 June 2022) of the Joint Genome Institute (JGI). F. solani FSSC 5 and F. vanettenii 77-13-4 were selected because they are the closest relatives to F. keratoplasticum within the FSSC clade. The F. keratoplasticum genome has yet to be sequenced. F. graminearum is a more distant relative belonging to a neighbour clade Fusarium sambunicum species complex. A. fumigatus was included as a pathogenic mould for which azole resistance mechanisms have been well studied. S. cerevisiae was included as a well characterised yeast species.

Validation of RNA-Seq DEGs Using qPCR
First strand cDNA synthesis of 1 µg total RNA and the quantification of mRNA expression levels of individual genes was conducted by real-time qPCR as previously described [25,46]. The three cyp51 paralogs (cyp51A, cyp51B and cyp51C) and abc1 were selected for validation of the RNA-seq results [25,46]. The erg6A and erg6B paralogs encoding sterol C24-methyltransferases were also included. qPCR of three biological replicates with three technical replicates each was performed to confirm the RNA-seq results.

De Novo Transcriptome Assembly and Annotation of Transcripts
De novo transcriptome assembly of sequencing reads obtained from all six RNA samples generated a total of 9,011,979 contigs, of which 37,274 were identified as unique transcripts. The sequence length distribution of the 37,274 transcripts ( Figure 1a BLAST+ software version 2.13.0 [58] was used to conduct a local BLASTp search against the protein sequences of Fusarium solani FSSC 5 [59], F. vanettenii 77-13-4 [60], F. graminearum PH-1 [61], A. fumigatus Af293 [62] and S. cerevisiae S288C [63] downloaded from the publicly available MycoCosm [64] database (accessed on 29 June 2022) of the Joint Genome Institute (JGI). F. solani FSSC 5 and F. vanettenii 77-13-4 were selected because they are the closest relatives to F. keratoplasticum within the FSSC clade. The F. keratoplasticum genome has yet to be sequenced. F. graminearum is a more distant relative belonging to a neighbour clade Fusarium sambunicum species complex. A. fumigatus was included as a pathogenic mould for which azole resistance mechanisms have been well studied. S. cerevisiae was included as a well characterised yeast species.

Validation of RNA-Seq DEGs Using qPCR
First strand cDNA synthesis of 1 µg total RNA and the quantification of mRNA expression levels of individual genes was conducted by real-time qPCR as previously described [25,46]. The three cyp51 paralogs (cyp51A, cyp51B and cyp51C) and abc1 were selected for validation of the RNA-seq results [25,46]. The erg6A and erg6B paralogs encoding sterol C24-methyltransferases were also included. qPCR of three biological replicates with three technical replicates each was performed to confirm the RNA-seq results.

De Novo Transcriptome Assembly and Annotation of Transcripts
De novo transcriptome assembly of sequencing reads obtained from all six RNA samples generated a total of 9,011,979 contigs, of which 37,274 were identified as unique transcripts. The sequence length distribution of the 37,274 transcripts ( Figure 1a

Genes Differentially Expressed in the Presence of VRC
A total of 316 genes were significantly differently expressed (fold-change > 2 and q-value < 0.05) in the presence of VRC; 233 were upregulated and 83 were downregulated ( Figure 2a). To assess the reproducibility of the 316 significant DEGs in each sample from the two treatment groups (VRC-treated and untreated), we compared their FPKM profiles. Figure 2b shows that the triplicate repeats for cells grown either in the presence (T1, T2, T3) or absence (UT1, UT2, UT3) of VRC clearly grouped into two distinct clusters, confirming good reproducibility and the reliability of the results produced by the biological replicates.
Of the 316 DEGs, GO enrichment analysis revealed 117 genes were associated with various cellular components, 67 were involved in biological process and 22 were involved in molecular functions as shown in Figure 3. the two treatment groups (VRC-treated and untreated), we compared their FPKM profiles. Figure 2b shows that the triplicate repeats for cells grown either in the presence (T1, T2, T3) or absence (UT1, UT2, UT3) of VRC clearly grouped into two distinct clusters, confirming good reproducibility and the reliability of the results produced by the biological replicates.
Of the 316 DEGs, GO enrichment analysis revealed 117 genes were associated with various cellular components, 67 were involved in biological process and 22 were involved in molecular functions as shown in Figure 3.

Inventory of All F. keratoplasticum Sterol Biosynthesis Gene Transcripts and Their
Orthologs in F. solani, F. vanettenii, F. graminearum, A. fumigatus and S. cerevisiae The search for all possible fungal ergosterol biosynthesis gene orthologs in F. keratoplasticum mRNA transcripts revealed the presence of eight duplicated genes (erg10, erg13, erg7, erg6, erg24, erg25, erg3 and erg5) in addition to the azole drug target, cyp51 which is known to have three paralogs in Fusarium [25]. In order to 'correctly' assign 'A' and 'B' orthologs (e.g., cyp51A or cyp51B) for the paralogs of erg10, erg13, erg7, erg6, cyp51, erg24, erg25, erg3 and erg5, we compared F. keratoplasticum protein sequences with all homologous proteins from three different Fusarium species (F. solani FSSC 5, F. vanettenii, and F. graminearum), A. fumigatus and S. cerevisiae ( Figure 4). The two FSSC species F. solani and F. vanettenii had the same eight duplicates as F. keratoplasticum and three cyp51 genes. A. fumigatus had nine ergosterol biosynthesis gene duplicates (erg10, erg13, hmg, erg7, cyp51, erg24, erg25, erg26 and erg4) and it had three erg3 paralogs. F. graminearum had seven gene duplicates (erg10, erg6, cyp51, erg24, erg26, erg3 and erg5) and S. cerevisiae had just one (HMG). After careful annotation of all ORF sequences (some introns and ATG start codons were clearly misannotated; see Supplementary File S1 for further details) their phylogenetic relationships were determined to help identify the corresponding orthologs for all six fungal species. To remain consistent with the commonly accepted nomenclature for mould cyp51 orthologs (i.e., 'A' orthologs are usually expressed at low levels under normal growth conditions and are less conserved but they are upregulated in response to azole stress), 'B' was assigned to the more conserved ortholog across all species and 'A' was assigned to the less conserved gene duplicate. Thus, the F. keratoplasticum 'B' orthologs of all duplicated genes had the highest homology to their orthologs in F. graminearum, A. fumigatus and S. cerevisiae (Figure 4). For example, of the two Fk-erg13 paralogs, Fk-erg13B had a much higher homology to F. graminearum erg13 (95%) than did Fk-erg13A (69%). A. fumigatus also had two erg13 paralogs, Af-erg13A and Af-erg13B. As expected, Fk-erg13B was significantly more conserved (81%) than Fk-erg13A (71%) when compared to their respective A. fumigatus orthologs or to S. cerevisiae ERG13 ('A' 68% and 'B' 75%). The phylogenetic relationships between all F. keratoplasticum genes in the later stages of the ergosterol biosynthesis pathway, beginning with erg7 (lanosterol synthase), and their homologs in the remaining five fungal species are presented in Supplementary Figures S1a,b. Phylogenies of the human cholesterol gene (dhcr7 and dhcr24) homologs and the mevalonate pathway genes (erg10, erg13, hmg1) are available in Supplementary Figure S1c,d, respectively. The percentage identities and homologies for all genes across all species analysed in this study can be found in the 'Fker FSSC2 (protein)' sheet in Supplementary File S1.
To our surprise, we also found five genes in all three fusaria ( Figure 4, green font) that were homologs of human cholesterol biosynthesis genes [65]: the ∆8,∆7-sterol isomerase homolog of emopamil binding protein (ebp1) [66], the 7-dehydrocholesterol reductase (dhcr7) and three paralogs of 24-dehydrocholesterol reductase (dhcr24) which suggests that fusaria can also synthesize cholesterol or cholesterol-like compounds [67]. Further investigations revealed two ebp1 homologs and one dhcr24 homolog, ortholog of dchr24_1, in A. fumigatus and the same five epb1, dhcr7 and three dhcr24 orthologs in F. graminearum. Phylogenetic trees demonstrating the evolutionary relationship between these human cholesterol biosynthesis homologs are presented in Supplementary Figure S1b and their homologs in the remaining five fungal species are presented in Su Figures S1a and S1b. Phylogenies of the human cholesterol gene (dhcr7 and ologs and the mevalonate pathway genes (erg10, erg13, hmg1) are available in tary Figures S1c and S1d, respectively. The percentage identities and hom genes across all species analysed in this study can be found in the 'Fker FS sheet in Supplementary File S1.  For genes with no obvious orthologs in other fungal species the gene most homologous to the F. keratoplasticum gene is indicated inside the rectangle. NA indicates a homolog was not found. Genes (ebp1, dhcr7 and dhcr24) with characteristics similar to human cholesterol biosynthesis [65] genes are highlighted in green font. Asterisk (*) indicates only a partial gene transcript was assembled/sequenced.

Identification of F. keratoplasticum ABC Transporters and Transcription Factors Differentially Expressed in VRC Treated Cells
Among the 316 F. keratoplasticum significant DEGs, seven encoded proteins (DN33072, DN3265, DN5825, DN16847, DN36936, DN21770 and DN29261) that belonged to the ABC transporter superfamily. These proteins had homologs in other fusaria (Table 1) with percentage identities between 75-100%. F. keratoplasticum abc1 (DN16847) and abc3 (DN36936) were the closest homologs of A. fumigatus PDR transporter abcG1, also known as abcC or cdr1B, a major efflux pump involved in azole resistance in A. fumigatus [34,35] (JGI protein ID: 1300; Supplementary File S1), with percentages identity/homology of 60/73 and 63/76 (Table 1), respectively. F. keratoplasticum DEGs also included transcription factors. All six (DN25544, DN12705, DN14162, DN37099, DN16099 and DN13837) differentially expressed transcription factors were also present in all other fusaria. F. keratoplasticum atF (DN14162) orthologs were present in F. solani FSSC5 and F. vanettenii but the F. graminearum atF homolog was not as highly conserved as the other five TFs, suggesting that Fg-atF is possibly not an ortholog of Fk-atF and has evolved its own unique function in F. graminearum. Homologs of moc3 (DN12705), atF (DN14162), SR (DN37099) and SREBP1 (DN13837) were not found in A. fumigatus. No homologs of the F. keratoplasticum transcription factors were found in S. cerevisiae.

Transcriptional Response of F. keratoplasticum Cells to VRC Exposure Suggests the Presence of an Alternative Route in the Ergosterol Biosynthesis Pathway
Among the 37,274 transcripts (Figure 1a), 32 were genes from the ergosterol biosynthesis pathway of which 22 were significantly upregulated (Figure 5a, red shading in the rectangles on the right column of the heatmap; and Table 2) in response to VRC. In fact, erg6A, cyp51A and, most interestingly, the human cholesterol biosynthesis ortholog, ebp1 (Table 2), were the three most upregulated DEGs in response to VRC; 912-fold, 52-fold and 20-fold upregulated, respectively. To determine the transcription factors possibly involved in F. keratoplasticum 2781 ergosterol pathway gene regulation, we compared the log2 fold change values to those caused by transcription factor deletion in F. graminearum-strains PH-1 ∆FgSR [33] and ∆Fg-atrR [36] (Figure 5a). We found that 16 out of the 22 genes in F. keratoplasticum upregulated in response to VRC were downregulated in F. graminearum ∆FgSR suggesting that FkSR is also the main transcription factor involved in the upregulation of the ergosterol biosynthesis genes in F. keratoplasticum upon azole antifungal treatment. The unchanged FkSR expression level upon VRC exposure (Table 2) is not unexpected because the FgSR expression levels also remained unchanged upon tebuconazole exposure of F. graminearum PH-1 cells [33]. The F. keratoplasticum atrR ortholog of Fg-atrR, however, was the most significantly upregulated (4.5-fold) transcription factor upon VRC exposure ( Table 2). This too was expected as Fg-atrR transcript levels were also increased in response to tebuconazole treatment of F. graminearum cells [68]. The F. graminearum ∆Fg-atrR expression profile demonstrates that FgAtrR is also involved in the upregulation of some ergosterol biosynthesis genes, but it has a smaller target range as only six genes (erg1, cyp51A, erg26A, ebp1, erg6A and erg5A) were downregulated upon its deletion (Figure 5a).   3 DEGs with log2 fold changes above 2.00 are highlighted with a grey background. 4 Asterisk (*) indicates no significant change in expression (log2 fold changes with a q-value of > 0.05). 5 Genes are listed according to their predicted position in the ergosterol biosynthesis pathway (see text for further details). 6 Fan et al. [38] claimed that F. graminearum cyp51C does not encode sterol 14αdemethylase. 7 The cluster grouping of the ABCG family is based on previous reports by Lamping et al. [69] and James et al. [46]. 8 AtrR transcription factor orthologs have been reported to mediate azole resistance in A. fumigatus [34,35] and F. graminearum [36]. 9 The novel zinc-cluster transcription factor, FgSR, is the master regulator of ergosterol biosynthesis in F. graminearum [33].
The final clue for the proposed alternative ergosterol biosynthesis pathway in Figure 5b came with the discovery of the human ebp1 ortholog which was present in all fusaria species investigated (Supplementary File S1 and Supplementary Figure S1b). Human EBP orthologs perform the same ∆8,∆7-sterol isomerase reaction as fungal Erg2 orthologs [66]. Thus, like its human counterpart FkEbp1 most likely acts on sterol intermediates that are not methylated at their C24 position. Methylation at the C24 position, a characteristic feature differentiating ergosterol from cholesterol, is catalyzed by Erg6. In S. cerevisiae, Erg6 is thought to act after the Erg25/26/27/28 complex converting zymosterol to fecosterol which is the substrate for Erg2. However, mould Erg6B orthologs prefer lanosterol as a substrate [38,43] which places them before the Cyp51B orthologs in the ergosterol biosynthesis pathway (Figure 5b). These observations together with the facts that (i) Fk-erg6A, Fk-cyp51A and Fk-ebp1 were the most dramatically upregulated genes in response to VRC; (ii) they were the only three genes, apart from erg5A, that were significantly downregulated upon the deletion of both FgSR and Fg-atrR; and (iii) the likelihood that ebp1 acts on non-C24 methylated sterol intermediates led to the proposed alternative ergosterol biosynthesis pathway shown in Figure 5b. Fk-erg7A, Fk-erg25A and Fk-erg5A may not take part in either of these two pathways as neither is expressed under either of the growth conditions investigated ( Table 2).  To compare the expression levels of each paralog in the presence and absence of VRC, we examined their FPKM values (Table 2) in each treatment group. We found that in the absence of VRC, expression levels of both erg6A (3) and cyp51A (2) were~70 times lower than their 'B' paralogs (207 and 104, respectively). Conversely, when exposed to VRC, expression levels of erg6A (4186) and cyp51A (2982) were~3 times higher than their 'B' paralogs (1694 and 899, respectively). This trend was not observed in the other seven (erg10, erg13, erg7, erg24, erg25, erg3 and erg5) gene duplicates where their relative expression levels (i.e., expression level of 'B' paralog was higher than the 'A' paralog, except for erg3 which was the opposite) were similar in both the presence and absence of VRC.
Five transcription factors were differentially expressed in F. keratoplasticum when exposed to VRC of which three (DN25544, DN12705, DN14162) were up-regulated and two (DN16099, DN13837) were down-regulated (Table 2). Interestingly, the highest (4.4-fold) upregulated transcription factor DN25544 is clearly the ortholog of atrR, an autoregulatory Zn 2 -Cys 6 zinc cluster transcription factor and master regulator of cyp51 orthologs and/or ABC transporters involved in azole resistance in F. graminearum [36] and A. fumigatus [34,35]. The transcription factor DN13837, which was downregulated 2.2-fold, is a sterol regulatory element binding protein 1 (SREBP1) homolog; these proteins are the main regulators of sterol biosynthesis in many fungal species (e.g., SrbA in A. fumigatus [74]) and in higher eukaryotes including humans [75]). DN16099 (2.1-fold downregulated) is an AP-1 like transcription factor (yap5), homologs of which are involved in stress response and yeast iron regulatory mechanisms [76].

Correlation between Gene Expression Levels Measured by qPCR and RNA-Seq
The expression levels for key genes identified from the RNA-seq analysis were measured by qPCR. The Cq values were normalized with GAPDH (gpd1, Figure 6a), and their fold up-or down-regulation in response to VRC (Figure 6b) was determined using the 2 −∆∆Cq method [77]. Clear differences in expression levels ( Figure 6a) were observed between the cells grown in the presence and absence of VRC. The qPCR analysis revealed a 2100-fold upregulation of erg6A, a 9-fold upregulation of erg6B and a 2400-, 10-and 4-fold upregulation of cyp51A, cyp51B and cyp51C, respectively, and a 120-fold upregulation of abc1 in response to VRC (Figure 6b). These VRC responses showed similar trends to the expression levels determined by RNA-seq (R 2 = 0.712; Figure 6c). In particular, erg6A and cyp51A were both expressed far less (i.e., 2142-fold and 2431-fold less), at almost undetectable levels, than their 'B' paralogs under normal growth conditions, but they were induced much more than their 'B' paralogs in response to VRC so that they reached comparable ( Figure 6) or even 2-to 3-fold higher ( Table 2) expression levels than their 'B' paralogs in the qPCR and RNA-seq experiments, respectively.

Discussion
A recent landmark study by Liu et al. [33] highlighted a novel transcription factor in F. graminearum, FgSR, which binds to a 16 bp cis-element present in the promoter of 119 genes, 20 of which were involved in ergosterol biosynthesis including cyp51A, cyp51B, cyp51C, erg6A, and erg6B. Fourteen of these 20 genes were also significantly downregulated in the ∆FgSR strain [32]. Fg-erg5A and Fg-erg4 were two notable exceptions. Although FgSR did not bind to either promoter, both genes were also significantly downregulated in the ∆FgSR strain [32]. Deletion of FgSR in F. graminearum had an inverse effect to that of VRC exposure on all these orthologs in F. keratoplasticum (Figure 5a). Again, the only exceptions were erg5A and erg4 that were either not expressed under either growth condition (erg5A) or fell just below (0.98; erg4) the log2-fold change threshold of 1.0 (Supplementary File S1). These striking similarities in the regulation of each of these ergosterol biosynthesis orthologs indicates that FkSR is also likely the key transcription factor that regulates ergosterol biosynthesis in F. keratoplasticum, and it too possibly functions via

Discussion
A recent landmark study by Liu et al. [33] highlighted a novel transcription factor in F. graminearum, FgSR, which binds to a 16 bp cis-element present in the promoter of 119 genes, 20 of which were involved in ergosterol biosynthesis including cyp51A, cyp51B, cyp51C, erg6A, and erg6B. Fourteen of these 20 genes were also significantly downregulated in the ∆FgSR strain [32]. Fg-erg5A and Fg-erg4 were two notable exceptions. Although FgSR did not bind to either promoter, both genes were also significantly downregulated in the ∆FgSR strain [32]. Deletion of FgSR in F. graminearum had an inverse effect to that of VRC exposure on all these orthologs in F. keratoplasticum (Figure 5a). Again, the only exceptions were erg5A and erg4 that were either not expressed under either growth condition (erg5A) or fell just below (0.98; erg4) the log2-fold change threshold of 1.0 (Supplementary File S1). These striking similarities in the regulation of each of these ergosterol biosynthesis orthologs indicates that FkSR is also likely the key transcription factor that regulates ergosterol biosynthesis in F. keratoplasticum, and it too possibly functions via phosphorylation mediated by the kinases of the high osmolarity glycerol response (HOG) pathway [33]. A similar 16 bp cis-element binding site is present in the cyp51A and cyp51B but not the cyp51C promoters of F. keratoplasticum [25]. We also found this motif in all cyp51A and cyp51B promoters of two other FSSC species [25] as well as in F. graminearum, F. verticillioides, and F. oxysporum [33].
As expected, FkSR, the ortholog of FgSR (Table 2), was not differentially expressed by VRC exposure which is consistent with the findings in Liu et al. [33] where F. graminearum cells exposed to tebuconazole or depleted in ergosterol in a mutant deleted of erg4 did not alter FgSR protein levels. This confirms that the novel sterol regulator SR is indeed likely to be a characteristic feature of all Leotiomycetes and Sordariomycetes species [33].
Fk-atrR was the transcription factor most significantly upregulated (4.4-fold; Table 2) in response to VRC. FkAtrR is clearly the ortholog (63% protein identity) of AfAtrR (ABC transporter regulator; their phylogenetic relationship is shown in Figure 1 of [36]) that mediates resistance to azoles in A. fumigatus by possibly autoregulating its own expression [33] and co-regulating srbA, cyp51A and abcG1/cdr1B expression [34,35,45]. AtrR orthologs act similarly in fusaria species. Deletion of F. graminearum atrR (FGSG_06810; 84% protein identity) [36] caused the downregulation of cyp51A (FGSG_04092; 19-fold), erg5A (FGSG_03686; 29-fold) and erg6A (FGSG_05740; 46-fold) which are clearly the orthologs of Fk-cyp51A (81% protein identity), Fk-erg5A (80%) and Fk-erg6A (86%), respectively (Supplementary Figure  S1a,b). Unfortunately, FGSG_03686 and FGSG_05740 were named Fg-erg5B and Fg-erg6B and not Fg-erg5A and Fg-erg6A, respectively, in that publication, even though, like for the Fg-cyp51A and Fg-cyp51B orthologs, they are less well conserved than Fg-erg5B (FGSG_01959) and Fg-erg6B (FGSG_02783) and expressed at much lower levels than Fg-erg5B and Fg-erg6B under normal growth conditions and two of the most downregulated genes in the ∆FgSR strain ( Figure 5a and 'Fgra (gene)' sheet in Supplementary File S1). Thus, for consistency and to avoid confusion, we modified the nomenclature for all F. graminearum ergosterol biosynthesis gene duplicates so that all fusaria orthologs are labelled with the same 'A' (salvage pathway) or 'B' (default pathway) extensions (Figures 4 and 5b). Aspergillus species have only one erg6 ortholog which appears to be the common ancestor of all fusaria erg6B orthologs (Supplementary Figure S1b). As expected for the FgAtrR ortholog of AfAtrR which possibly autoregulated its own expression and co-regulated srbA expression in A. fumigatus [34], Fg-atrR expression was not affected in the ∆FgSR knock-out strain [33], but it was upregulated (3.3-fold) in response to tebuconazole [68], which is very similar to our findings with Fk-atrR being 4.4-fold upregulated in response to VRC.
To date, detailed mechanisms of ergosterol biosynthesis remain underexplored in mould pathogens like fusaria. The discovery of the human ∆8,∆7-sterol isomerase ortholog ebp1 in F. keratoplasticum and all related FSSC species, and even two paralogs ebp1A and ebp1B in A. fumigatus, was rather surprising given that fungi have evolved their own characteristic ∆8,∆7-sterol isomerase erg2. Plants and animals, on the other hand, use EBP orthologs for the ∆8,∆7-sterol isomerase reaction [78,79]. Even though murine EBP could complement a yeast mutant deleted in ERG2 which exhibited a very similar sterol composition to wildtype S. cerevisiae, murine EBP and ScErg2 share no similarities and are phylogenetically completely different proteins [66]. erg2 encodes a typical bitopic membrane protein with an N-terminal transmembrane span while Fk-ebp1 and all related fungal ebp1 homologs encode proteins with five predicted transmembrane spans, respectively (Supplementary File S1). Perhaps this is why these ebp1 homologs have been overlooked for so long in the model moulds A. fumigatus and F. graminearum.
Because Fk-ebp1 was one of the sterol biosynthesis genes most highly upregulated in response to VRC and because its presence forced us to ascertain how it might possibly fit into the commonly assumed ergosterol biosynthesis pathway in F. keratoplasticum, we were interested to find out whether this enzyme had simply been missed in previous investigations, by searching their data. That was indeed the case. The F. graminearum ebp1 ortholog FGSG_13888 was actually the most dramatically down-regulated sterol biosynthesis gene (919-fold) in the ∆FgSR strain (we found it among the list of DEGs in the Supplementary Data 3 file of [33]) and also one of the most down-regulated (8.6fold) sterol biosynthesis genes in the ∆FgAtrR strain (we found it among the list of DEGs in the Supplementary Table S2 of [36]). The cyp51A and ebp1 orthologs were also the highest upregulated genes (163-fold each) in response to exposing F. graminearum cells to tebuconazole [68]. In addition, one of the two A. fumigatus orthologs, ebp1B, was one of only four sterol biosynthesis genes (i.e., erg24A, erg25A, erg3B and ebp1B) that was dramatically downregulated in a ∆SrbA strain when the cells were grown under hypoxic growth conditions [74]. Despite its obvious significance, ebp1 was not recognized as the 'missing link' that it appears to be for an alternative ergosterol biosynthesis salvage pathway that is induced by azole antifungal treatment of fusaria and possibly also A. fumigatus. However, the rather different set of erg paralogs in A. fumigatus suggests a different salvage pathway is induced by azole exposure in this Eurotiomycetes model species [42,80].
We propose two major changes to the ergosterol biosynthesis pathway in response to the inhibition of cyp51B by VRC. Erg6A acts after Ebp1 because it possibly prefers non-C24 methylated substrates, like its mammalian homologs (C24 methylation is a typical feature of fungal and plant sterols), and Ebp1 acts after the Erg25/26/27/28 complex right where it is placed in the cholesterol biosynthesis pathway. Cyp51A possibly prefers lanosterol as a substrate. erg7A, erg25A and erg5A (in orange fonts in Figure 5b) possibly play no part in either pathway because their expression levels are simply too low under either growth condition (i.e., with and without VRC). Erg24A and Erg24B are both upregulated to metabolise both types of sterol intermediates, those that are C24 methylated by Erg6B and processed by the only partially inhibited Cyp51B (F. keratoplasticum 2781 had MIC VRC > 32 mg/L [25]) and those that are not methylated. These two intermediates are further processed by the Erg25B/26/27/28 complex which can process both types of intermediates to generate fecosterol and zymosterol, respectively. Fecosterol and zymosterol are then processed by the two ∆8,∆7-sterol isomerases Erg2 and Ebp1, respectively. Further progression in the two pathways is not clear. However, because there are two erg3 paralogs, both of which are expressed at similar levels under normal growth conditions (Table 2) and erg3A is significantly more induced by VRC (5.9-fold) than erg3B (2.1-fold), there are possibly two separate pathways leading all the way to the production of ergosterol. The final four enzymatic reactions necessary for ergosterol biosynthesis may involve a different series of intermediate steps (hence the question marks for the placements of erg6A and erg3A and an X for different intermediates in the proposed salvage pathway in Figure 5b). Three alternative pathways were proposed for the conversion of fecosterol to ergosterol in A. fumigatus. Each pathway involved a different series of enzymatic reactions carried out by Erg2, three Erg3, Erg5 and two Erg4 paralogs, each preferring different intermediates as their respective substrates (see Figure 3 in [42]). In all three of these proposed pathways the first step involved the conversion of fecosterol to episterol by Erg2. However, our discovery of two ebp1 paralogs, one of which, ebp1B, and erg3B were required for adaption to hypoxic growth [74], suggests one or both of these ebp1 paralogs may play a part in these alternative pathways.
The possible contribution of the third Cyp51C paralog to ergosterol biosynthesis remains to be investigated even though it was claimed to have no sterol 14-α demethylase activity [38]. This assumption was based on the fact that Fg-cyp51A and Fg-cyp51B, but not Fg-cyp51C, could complement a S. cerevisiae mutant deleted in ERG11. However, the authors did not explore whether Fg-cyp51C was actually expressed nor did they consider the possibility that FgCyp51C may not accept lanosterol as its preferred substrate. The same authors also created null-mutants for all three Fg-cyp51 paralogs and the double deletion mutants for the A,C and B,C paralog combinations. However, the authors did not create a double deletion mutant for the A,B paralog combination because they "considered it to be lethal" [38]. Yet, all indications are that cyp51C is part of the ergosterol biosynthesis pathway. It too is upregulated in response to azole stress, deletion of Fg-cyp51C also affected the sterol composition of F. graminearum cells and it was required for virulence on host wheat ears [38]. It is therefore conceivable that a Fg-cyp51A/Fg-cyp51B double deletion mutant is viable and possibly uses another alternative ergosterol biosynthesis pathway to produce enough ergosterol or sterol intermediates for cell survival. F. keratoplasticum erg7A, erg25A and erg5A were not expressed under either growth condition tested, yet they are conserved in all FSSC species investigated. This suggests the possibility of an additional ergosterol biosynthesis salvage pathway that may be induced by hypoxia or other growth conditions not tested in the current study.
Dhcr7 and Dhcr24 catalyse the final two steps of cholesterol biosynthesis by reducing the double bonds at C7 and C24, respectively [67]. The presence of one dhcr7 and three dhcr24 orthologs in fusaria ( Figure 4) indicates that fusaria may also be able to produce cholesterol or cholesterol-like compounds under certain stress conditions. A viable 'cholesterol yeast' was successfully engineered by deletion of S. cerevisiae ERG5 and ERG6 and replacing them with fish DHCR7 and DHCR24, respectively [81]. Chytridiomycota are examples of fungi that naturally synthesise cholesterol and other ∆5 sterols as their predominant membrane sterols [82]. In addition, an entire cholesterol biosynthesis pathway has also recently been reported for tomato plants that can produce significant quantities of cholesterol [79]. Fk-dhcr7 and Fk-dhcr24_3 were not expressed to any significant amount under either of the two growth conditions investigated. However, Fk-dhcr24_2, and to a lesser extent Fk-dchr24_1, were expressed at reasonable levels ( Table 2) indicating they are functionally relevant under normal growth conditions. Perhaps F. keratoplasticum produces cholesterol or cholesterol-like compounds at higher temperatures or during host invasion by upregulating Fk-dhcr7 and Fk-dhcr24_3 and downregulating Fk-erg5B and Fk-erg6B, respectively. Cholesterol may provide a more rigid membrane environment enabling cells to grow at higher temperatures. No further possible plant, fungal or human sterol biosynthesis homologs were found in any of the fungal genomes of F. graminearum, F. solani FSSC5 or F. vanettenii.

Conclusions
The highlight of this study is the discovery of mould ebp homologs that are possibly an integral part of an alternative ergosterol biosynthesis salvage pathway that is upregulated in fusaria species upon azole antifungal exposure. VRC exposure of F. keratoplasticum also caused the upregulation of a number of ABC transporters including the major multidrug efflux pump genes abc1 and abc3. The antifungal drug resistance response of F. keratoplasticum appears rather well conserved in all fusaria species with orthologs of the sterol biosynthesis master regulator, SR, and the pleiotropic drug resistance regulator, AtrR, working together to protect fusaria from the effects of azole antifungals. We also discovered the elements of a conserved cholesterol biosynthesis pathway that may be expressed at higher temperature or under different stress conditions. Future biochemical studies are required to confirm these discoveries.