Azole Resistance Mechanisms in Pathogenic Malassezia furfur

Malassezia spp. are emerging fungal pathogens causing opportunistic skin and severe systemic infection. Nosocomial outbreaks are associated with azole resistance, and understanding of the underlying mechanisms is limited to knowledge of other fungal species.

M alassezia, the lipid-dependent skin commensal, has also been shown to be an emerging pathogen not just in neonatal intensive care wards but also in Crohn's disease and pancreatic cancers (1)(2)(3)(4). Malassezia-associated superficial mycoses (dandruff, pityriasis versicolor, and folliculitis) affect up to 50% of the global population (5) and have the potential to cause severe systemic and invasive infection in immunocompromised individuals (6)(7)(8)(9). Outbreaks in neonatal and intensive care wards have been increasingly reported (6). Azole antifungals are the primary antifungal treatment for Malassezia-associated diseases due to their intrinsic resistance to some antifungals, particularly echinocandins (10). Because reduced in vitro antifungal susceptibility may be an indicator of clinical failure (11), there is a need to profile the susceptibility of Malassezia isolates derived from individuals with healthy or disease backgrounds to elucidate and design prevention strategies for the underlying mechanisms contributing to current and emerging trends of antifungal resistance.
Azole resistance may be primary (intrinsic) or secondary (acquired) (12). The former is found naturally without prior (known) antifungal exposure. The latter is a result of a previously susceptible strain being exposed to antifungals or other selective pressure and may be a result of altered gene expression, point mutations, or allelic variations The isolates further clustered into 2 groups, with disease isolates in cluster 2 ( Fig.  1A, green) having elevated MICs for almost all azole antifungals and disease isolates in cluster 1 (Fig. 1A, yellow) having high MICs specifically for TRB, CTZ, and MCZ (P , 0.001) ( Fig. 1B; also, see Table S1 in the supplemental material). ITZ had the lowest MICs of all antifungals tested (Fig. 1B), although disease isolates of M. furfur still had significantly higher MICs for KTZ (P , 0.05) (Fig. 1B). FLZ and VRZ MICs were higher for disease isolates in cluster 2 (Fig. 1A, green), although these values were not significant.
CYP51A1 Y67F mutants have elevated fluconazole and voriconazole MICs. The gene sequence of CYP51A1 in M. furfur was identified by a whole-genome-sequencing (WGS) BLAST search (BioProject ID 286710) of a known M. globosa CYP51A1 sequence (15) against M. furfur genomes and validated by the identification of the presence of the highly conserved heme-binding site (FGFGRHRCIG), EXXR motif (ERLR), and conserved threonine of the I helix involved in proton delivery (Fig. S1).
A tyrosine-to-phenylalanine mutation in residue 67 (Y67F) of CYP51A1 was identified in the blood isolate CBS 14141 ( Fig. 1C; Fig. S2). The mutation was also identified in four other isolates, i.e., the systemic-disease isolates Mal24, Mal26, and JLPK13 and the healthy-skin isolate 003 SC2. The three systemic-disease isolates showed elevated MICs ( Fig. 1B and C; Table S1) specifically for FLZ (P , 0.001), consistent with literature indicating the role of the tyrosine residue in drug binding (28). Healthy skin isolate 003 SC2 also showed slightly elevated MICs for MCZ and FLZ. This mutation is synonymous with the Y132F mutation in Candida albicans and is well known to affect binding to fluconazole due to the residue being 7 Å away from a key binding motif (29). It has been identified as a Y127F mutation in disease isolates of M. globosa (15).
Y67F knock-in and F67Y rescue mutants show no changes in azole susceptibility. The single CYP51 nucleotide change coding for the Y67F mutation was inserted into the wild-type, low-MIC strain CBS 7982 via Agrobacterium tumefaciens-mediated Skin catheter insertion PM315 Anal swab of neonate transformation (ATMT) using a nourseothricin (NAT) expression vector ( Fig. 2A to D). Similarly, a rescue F67Y CYP51 mutation was generated in the CYP51 mutant strain CBS 14141 ( Fig. 2A to D). Colonies were screened for the presence of the correct 1.5-kb flanking arms and insert (Fig. 2E). Sanger sequencing of the CYP51 mutation regions and internal transcribed spacer regions (ITS) was performed to validate the presence of the correct single nucleotide mutation and strain background of the respective transformants ( Fig. 2F and G). Insertion of the Y67F CYP51 mutation into the CBS 7982 strain did not yield any change in azole susceptibility ( Fig. 2H; Table 2). Similarly, an F67Y rescue mutant generated from the CBS 14141 CYP51 mutant strain did not show any reduction in azole MICs when tested against our panel of antifungals. The F67Y rescue mutation restored amphotericin B MICs to wild-type levels in CBS 14141, although introduction into wildtype CBS 7982 increased AMB MICs only slightly.
These observations suggest that reduced azole susceptibility may be mediated by factors other than CYP51 in M. furfur and that the Y67F mutation mainly impacts AMB  Antimicrobial Agents and Chemotherapy susceptibility in our selected strains. This is supported by the observation that the healthy-skin isolate 003 SC2 contains the Y67F CYP51 mutation but does not show the level of reduced azole susceptibility observed in disease isolate cluster 2.
Genes involved in metabolism and secondary metabolite production are upregulated in disease strains of M. furfur. To further elucidate the pathways potentially involved in M. furfur azole susceptibility, we performed RNA-seq in regular modified Dixon (mDixon) medium for representative candidate strains from each cluster (healthy isolates, disease isolate cluster 1, and disease isolate cluster 2) based on their antifungal susceptibility testing (AFST) profiles ( Fig. 1). They were defined as the wild type, susceptible strain CBS 7982 (healthy isolate cluster), the disease CYP51 mutant strain CBS 14141 (disease isolate cluster 2), and the intermediate disease strain CBS 7019 (disease isolate cluster 1).
A list of the top 20 upregulated genes commonly expressed in disease isolates CBS 7019 and CBS 14141 is presented in Table S2. Genes involved in metabolism, biological process, catalytic activity, and extracellular activity were the key genes upregulated in the disease strains (Fig. 3b). Statistically enriched KEGG pathways include secondary metabolite production (Fig. 3c). Upregulation of some exons coding for ABC and MFS transporter proteins, such as YBT1, ITR1, and OPT1, was also observed.
Healthy and disease isolates have different gene expression profiles following long-term exposure to clotrimazole in vitro. Long-term (up to 4 weeks) treatment of the wild-type healthy isolate CBS 7982 with clotrimazole (8 mg/ml) ( Fig. 4a; Table S3) induced up to an 8-fold increase in clotrimazole MICs (Fig. 4b). After treatment removal at 4 weeks, MICs were observed to fall back to close to starting values. This suggests that elevated MICs may be induced in vitro via a transient mechanism. For the strains CBS 7019 and CBS 14141, long-term clotrimazole treatment in vitro did not produce any meaningful change in MIC readings, as their starting MICs were already high.
Differential gene expression analysis was performed on RNA-seq of isolates exposed to the above-mentioned MICs of clotrimazole for the three M. furfur strains (CBS 7982, CBS 7019 and CBS 14141). Of 1,238 differentially expressed genes between 4-week clotrimazole-treated and nontreated CBS 7982, 689 were upregulated ( Fig. 4c and d).
Examples of genes commonly upregulated by clotrimazole treatment include DNM1, encoding the dynamin protein, and RBP1P, encoding the yeast RNA-binding protein, (Table S4), which have been reported to affect antifungal susceptibility and calcineurin-dependent azole tolerance in C. albicans, respectively (30,31). Metabolic pathways were commonly upregulated in 4-week clotrimazole-treated and nontreated groups for CBS 7982 and CBS 7019 ( Fig. 4e1 and e2). KEGG pathways upregulated in  clotrimazole-treated CBS 14141 were distinct from those in CBS 7982 and CBS 7019 and included the ubiquitin-mediated proteasome pathway and cell cycle (Fig. 4e3). A reverse trend was observed for KEGG pathways downregulated in the three strains, with metabolic pathways being downregulated in CBS7019 and CBS 14141 and the ubiquitin-mediated proteasome pathway being downregulated in CBS7982 ( Fig. 4f1 to f3). The ABC transporter encoded by PDR10 is upregulated after long-term exposure to clotrimazole in vitro. Narrowing RNA-seq gene identifications (IDs) to transporter pump proteins identified multiple reads matching exons mapping to "Similar to S. cerevisiae protein, SNQ2." These reads were consistently the single most highly upregulated transporter gene in all three M. furfur strains (Fig. 5). Further, multiple-sequence alignment matched the exon sequence to that of the multidrug transporter encoded by PDR10, as described by Ianiri et al. (27). Further analysis of PDR10 gene expression via reverse transcription-quantitative PCR (RT-qPCR) at baseline and at 3-, 4-, and 6-week time points showed that PDR10 expression increased over the 4week clotrimazole treatment period and dropped at 6 weeks, 2 weeks after treatment was withdrawn (Fig. 5b).
RNA-seq analysis also revealed that nine other transporter genes were upregulated in the clotrimazole-treated disease isolates CBS 7019 and CBS 14141 but not CBS 7982, including the previously reported mitochondrial ABC transporter, encoded by ATM1 ( Fig. 5a; Table S5). Other highly upregulated transporter genes in clotrimazole-treated CBS 7982 and CBS 7019 include the oligopeptide protein transporter gene OPT1, the multidrug transporter gene FLR1, and the caffeine resistance protein gene CRP5 (Fig.  5c), which have not yet been sequence validated in all isolates of the divergent M. furfur strains. These findings suggest that differential expression of transporter genes between different isolates may also affect antifungal susceptibility.
Deletion of PDR10 abrogates elevated MICs in disease isolate CBS 14141. The pleiotropic drug transporter gene PDR10 was identified to be significantly upregulated in clotrimazole-treated M. furfur strain CBS 7982 as described above. It was also found to play a role in Malassezia fluconazole resistance by Ianiri et al. (27), for which the CBS 14141 insertional mutants 7D9 and 2H11 and the CRISPR (clustered regularly interspaced palindromic repeats)/Cas9 pdr10D deletion mutant were constructed using the high-MIC mutant strain CBS 14141. 7D9 has a T-DNA insertion in the putative promoter region of PDR10, whereas 2H11 was previously described to have a chromosomal rearrangement mutant involving the gene ERG5 (27). We confirmed that the PDR10 gene was absent in the pdr10D mutant genome and that the pdr10D mutant showed null PDR10 expression (Fig. S3).
AFST was performed on the insertional mutants 7D9, 2H11, and the CRISPR/Cas9 pdr10D mutant alongside wild-type CBS 7982 and mutant CBS 14141 M. furfur strains ( Fig. 5d; Table 3). Deletion of PDR10 in CBS 14141 completely reversed the elevated MIC phenotype. The pdr10D strain showed an AFST profile close to that of the wild type (i.e., CBS 7982) (Fig. 5d). The 7D9 mutant exhibited intermediate susceptibility to  Wild-type and mutant strains show differences in rhodamine 6G efflux. The fluorescent small-molecule dye rhodamine 6G is commonly used to assess the efflux activity of ABC transporters (32). The pdr10D mutant and 7D9 were observed to have a higher uptake of rhodamine 6G (Fig. 6a) than CBS 7982 and CBS 14141 (Fig. 6a). This is consistent with literature which has reported that rhodamine 6G efflux is correlated with CDR1 expression (33). CBS 14141 was observed to have the highest rhodamine 6G efflux of the four strains, and CBS 7982 had the lowest (Fig. 6b). Surprisingly, the pdr10D mutant had rhodamine 6G efflux levels close to those of CBS 7982, suggesting that PDR10 is a key effector of rhodamine 6G efflux. The insertional mutant 7D9 showed intermediate rhodamine 6G efflux levels between those of CBS 7982 and CBS 14141.

DISCUSSION
In this study, we profiled the antifungal susceptibility of multiple M. furfur isolates collected from healthy and disease states and elucidated mechanisms underlying their antifungal susceptibility. Strains isolated from individuals with skin or systemic disease had elevated antifungal MICs relative to strains isolated from healthy individuals. Isolates of M. furfur from healthy and disease individuals cluster into groups which reflect their health/disease origin based on their antifungal susceptibility profiles. Isolates from individuals with systemic disease (e.g., blood, urine, and catheter isolates) were the most likely to have elevated MICs. However, exceptions include CBS 9374, which is from a healthy individual but still shows elevated MICs of terbinafine and miconazole. It is likely that the patterns of antifungal susceptibility observed in our disease isolates are stable, given that many of the disease strains have been maintained in culture and passaged repeatedly in fungal banks without antifungal exposure. However, a lack of information on the clinical background of many disease isolates (i.e., antifungal treatment prior to culture isolation and treatment dosage/duration)  hampers further interpretation regarding how reduced azole susceptibility could have arisen and whether it is intrinsic or acquired.
Although mutations in CYP51 are well documented to result in reduced azole susceptibility in fungi (17) (15), although the extent to which they confer azole resistance was not reported. The largescale construction of transformants from different strain backgrounds is currently limited by the need for strain-specific sequence data (preferably whole-genome sequencing data) for the design of 1.5-kb flanking arms for homologous recombination required in Agrobacterium tumefaciens-mediated transformation.
The presence of the Y67F mutation in the healthy-scalp isolate 003 SC2 suggests that the mutation alone is insufficient to confer reduced azole susceptibility. A synonymous Y134F mutation in the plant rust fungus Puccinia triticina was documented to have limited impact on its susceptibility to epoxiconazole (34). Our F67Y mutation in CBS 14141 was observed to result in increased amphotericin B susceptibility. This is consistent with observations in Leishmania mexicana, in which a single nucleotide mutation in CYP51 was also reported to induce amphotericin B resistance (35), and suggests a role for this gene in amphotericin B susceptibility.
Rapid transient increases in the mRNA expression of CDR efflux pump genes in C. albicans after exposure to fluconazole have been reported (36,37). Transporter pump genes such as PDR10 were observed to be upregulated after exposure to clotrimazole based on RNA-seq data. These genes were validated based on known homology to similar proteins in Saccharomyces cerevisiae or based on existing literature. While the transporter genes PDR10 (27), FLR1 (38), and ITR1 (39) have been reported to be associated with azole resistance, there is little information on the role of other transporter genes, such as OPT1 and CRP5, in multidrug resistance. A putative paralog of PDR10 was identified immediately adjacent and 39 to PDR10 and was identified as PDR10_2 by Ianiri et al. (27). This warranted further gene validation to ensure that the pdr10D mutant did not include a knockout of PDR10_2 (Fig. S3). However, we confirmed that a single knockout of the PDR10 gene (i.e., pdr10D) was sufficient to reduce MICs to wild-type levels.
While we demonstrated an increase in clotrimazole MICs after 4 weeks of successive treatment, the effect was dose dependent, and the phenomenon is best observed at antifungal concentrations two to four steps higher than the known MIC for the specific strain. While higher treatment doses may promote higher fold changes in MICs, they are less well tolerated, and the subsequent low viable inoculum impedes susceptibility testing.
At present, PDR10 expression in CBS 7982 appears to be strongly inducible by the presence of clotrimazole and likely azoles and other stress factors such as chemicals, UV exposure, elevated temperature and nutrient-limiting conditions (27). As described previously, a deletion in PDR10 for CBS 14141 resulted in increased sensitivity to fluconazole and benomyl, which are antifungal drugs that have different mechanisms of action (27). This suggests a critical role of PDR10 in nonspecific (or pleiotropic) cellular detoxification in M. furfur CBS 14141, most likely through active efflux of xenobiotics, including antifungal drugs, in line with the known function of ABC transporters. Further transcriptional analysis and deletion models in more strains and species of Malassezia are required to validate the role of ABC transporters in Malassezia azole resistance. This current progress is limited by the poor practical efficiency of generating ATMT plasmid constructs and the need for strain-specific knowledge of genome regions upstream and downstream of the gene of interest. While this has been improved by the development of a CRISPR/Cas9 system (27), challenges remain in the targeting of specific genes (versus unknown homologs) and the role of random insertions as confounding factors.
While rhodamine 6G assays have been useful in giving us an approximation of the relative ABC transporter efflux activity in strains from the same background (i.e., mainly CBS 14141), many transporter pumps are pleiotropic, and more specific binding assays are likely required to narrow down the role of specific pumps across multiple M. furfur strains. For the same reasons, the use of generic ABC transport inhibitors, such as verapamil and carbamazepine, was not useful for comparison across different strains.
In summary, we identified differences in azole antifungal susceptibility patterns of 26 isolates of M. furfur (13 from healthy subjects and 13 from those with disease). While a Y67F mutation was identified as a contributor to intrinsic resistance in some systemic-disease isolates, its impact on antifungal susceptibility must be interpreted on a strain-by-strain basis in context with other physiological factors. In CBS 7982 and CBS 14141, acquired resistance to clotrimazole and likely other azoles appear to be largely associated with differential activity of ABC transporters, particularly that encoded by PDR10. This has been validated by gene expression and functional deletion studies. There are likely additional factors contributing to interstrain differences, which are a result of multiple environmental stresses and other selection pressures. Understanding the functional mechanisms underlying azole resistance in Malassezia across a spectrum of strains, sources, and genetic backgrounds will be useful for the identification of new therapeutic targets to prepare for the emergence of new resistant strains.

MATERIALS AND METHODS
Strains and culture conditions. No active primary culture isolation was performed in this study. All 26 isolates used in this study were either obtained from the Westerdijik Fungal Diversity Institute or obtained from previously published studies (40)(41)(42). Based on data available in fungal bank records and from previous publications, (13 of the strains were isolated from healthy individuals and 13 from patients with preexisting skin or medical conditions) ( Table 1). Strains CBS 14141, JLPK13, Mal18, Mal24, Mal25, Mal26, Mal32, and PM315 were a kind contribution from Bart Theelen and Claudia Cafarchia. All Malassezia furfur strains were maintained on modified Dixon agar or broth at 32°C as described previously (10,40).
Antifungal susceptibility testing. Antifungal susceptibility testing was performed using a broth microdilution method as described by Leong et al. (10). Briefly, 200Â drug stock dilutions were prepared at a 2Â concentration in fresh OptiMAL medium. Amphotericin B, terbinafine, clotrimazole, miconazole, itraconazole, fluconazole, voriconazole, and ketoconazole were purchased from Sigma-Aldrich, Singapore. Stock and drug plate dilutions were prepared in accordance with CLSI and EUCAST guidelines. Yeast inocula were obtained from 4-to 7-day-old strains of Malassezia spp. A 50-ml yeast inoculum was added to 50 ml of 2Â concentrated antifungals to achieve a final cell density of 5 Â 10 3 to 5 Â 10 4 CFU/ml. A 10-ml portion of 2Â yeast inoculum that had been diluted 10 times was also plated onto a modified Leeming-Notman agar plate and incubated for 4 to 7 days at 35°C for postverification of the CFU inoculum (10 to 100 colonies per 10 ml). Each assay was performed in triplicate plates for a single culture at every individual time point or reading.
Long-term in vitro treatment with clotrimazole. CBS 7982 was maintained in triplicate cultures of mDixon broth containing 8 mg/ml of clotrimazole for 4 weeks at 32°C, with fresh medium and antifungal supplied every alternate day. After 4 weeks, medium was replaced with fresh medium only (i.e., no antifungal) every alternate day for another 2 weeks. A 100-ml portion of culture was removed every 7 days for antifungal susceptibility testing as described above. For CBS 7019 and CBS 14141, clotrimazole concentrations were both 256 mg/ml. RNA extraction was performed on aliquots of each triplicate samples at log phase in week 0 (before addition of clotrimazole) and week 4 (after 4-week clotrimazole treatment) for RNA-seq analysis.
Gene analysis. All primer sequences used in this study are listed in Table S5. Unannotated gene sequences for CYP51 and PDR10 were derived from a nucleotide BLAST of the whole-genome sequencing (WGS) data (NCBI BioProject database no. PRJNA286710 [43]) or by sequencing the specific gene using the self-designed primers (Table S5). PCR was performed using Platinum Taq DNA polymerase (Thermo Fisher, Singapore), TaKaRa HS DNA polymerase (TaKaRa Bio USA Inc.) for 1.5-kb flanking arms, and LongAmp (New England Biolabs [NEB], Singapore). Gel extraction of the PCR product was performed with a QIAquick gel extraction kit per the manufacturer's instructions (Qiagen, Singapore). Sanger sequencing was performed on the purified PCR product using BigDye (Thermo Fisher, Singapore) per the manufacturer's instructions. RT-qPCR was performed using the GoTaq one-step RT-qPCR system (Promega, Singapore) or a Superscript transcriptase III kit (Thermo Fisher, Singapore) with the respective gene-specific primers using the Malassezia actin gene ACT1 as the housekeeping gene. PCR with reverse-transcribed cDNA for actin was used to detect genomic-DNA contamination. Relative gene expression was analyzed using the Q-Gene module as described by Muller et al. (44).
Agrobacterium tumefaciens-mediated transformation. For ATMT, 1.5-kb upstream and downstream flanking arms of CYP51 were derived from PCR (TaKaRa Bio USA Inc.) of CBS 7982 or CBS 14141 genomic DNA and assembled with CBS 7982 CYP51 gene blocks (TTC/TAC) (Integrated DNA Technologies, Singapore) in the binary vector, pGI3, and the nourseothricin (NAT) and neomycin (NEO) plasmid cassettes pAIM2 and pAIM6 under the control of the M. sympodialis ACT1 promoter and terminator. Insertional mutagenesis was performed using the A. tumefaciens strain EHA105 as described previously (46).
RNA-seq. Total RNA was extracted and purified from log-phase cultures of CBS 7982, CBS 7019, and CBS 14141. Three biological replicates of each sample were analyzed. Briefly, cell pellets were washed three times in sterile phosphate-buffered saline (PBS) before resuspension in TRIzol and frozen at 280°C overnight. On thawing, samples were subjected to bead beating followed by use of the Direct-zol RNA miniprep kit (Zymo Research, USA), following the manufacturer's instructions. An additional DNase step was performed with the Turbo DNase kit (Thermo Fisher, Singapore) per the manufacturer's instructions. RNA sequencing and differential analysis services were provided by Novogene AIT, Singapore. One microgram of RNA was used for library preparation using a NEBNext Ultra RNA library prep kit for Illumina (NEB, USA) per the manufacturer's instructions. Next, 125-bp/150-bp paired-end reads were generated from an Illumina-based sequencing platform and processed with the appropriate quality controls. Differential expression analysis (untreated versus clotrimazole treated, healthy versus disease isolate) was performed using the DESeq R package (1.18.0).
Data and statistical analysis. All assays were performed as three independent experiments, each with triplicate readings unless otherwise stated. For pairwise comparisons, a paired two-tailed Student's t test was performed with Microsoft Excel. For grouped comparisons, a one-way analysis of variance (ANOVA) with Dunnett's test was performed with GraphPad Prism 8 (GraphPad Software, CA, USA). For heat map plotting of MICs, all MICs were normalized across the different antifungal concentrations such that the highest values were set at 1 and the lowest value were calculated as 1/total number of concentrations. Heat maps and Venn diagrams were plotted using the gplots package in R (version 3.03).
Rhodamine 6G assay. The rhodamine 6G efflux assay was performed as described previously (30). Briefly, equal volumes of cells were normalized to an optical density at 600 nm (OD 600 ) of 0.2 and incubated with 10 mM rhodamine 6G (Tee Hai Chem, Singapore) in PBS for 30 min at 32°C. Next, cells were transferred to 4°C and spun down to collect the supernatant. Pellets were washed twice with PBS and resuspended in equal volumes of PBS containing 2% glucose and incubated at 32°C. An aliquot of solution was removed at 0-, 5-, 10-, 15-, and 20-min intervals and spun down to collect the supernatant. Each sample was plated in triplicate.