Genomic Insights into the Atopic Eczema-Associated Skin Commensal Yeast Malassezia sympodialis

The draft high-coverage genome of M. sympodialis ATCC 42132 was assembled from a shotgun 454 data set and was extended and scaffolded using a 3-kb insert Illumina HiSeq mate-pair data set (Table 1) to a total of 65 scaffolds (L₅₀ = 511 kb), corresponding to a nuclear genome of 7.67 Mbp (Table 1). According to CEGMA (13) analysis, this assembly shows a high degree of completeness: 88.3 to 93.1%, comparable to 89.5 to 93.5% for the genome of the closely related species M. globosa (6). The estimated genome size of M. sympodialis is smaller than that of the M. globosa genome (8.96 Mbp), while both are in line with previously reported genome sizes from electrophoretic karyotyping experiments (14). The Malassezia genomes are among the smallest in the fungal kingdom, a feature probably related to their dependence on warm-blooded animals. Whole-genome alignments showed extensive synteny between the M. sympodialis and M. globosa scaffolds (Fig. 1A).

We predicted a total of 3,517 protein-coding genes in M. sympodialis using multiple lines of evidence (ab initio predictors, expressed sequence tags (ESTs) from M. globosa, and protein and nucleotide alignments) and manual annotations for specific genes. We applied mass spectrometry (MS) proteomics using peptide isoelectric focusing (IEF) on a broad and a narrow pH range to achieve high proteome coverage of M. sympodialis protein extracts and confirmed 3,176 (90%) of the predicted proteins. Of the MS-identified proteins, 98% have 2 or more peptide spectrum matches (PSMs) supporting their identification and 90% have 2 or more unique peptides identified (Fig. 1B) with sequence coverage of 13.5% or more for 75% of the proteins (Fig. 1C). The completeness of the annotation of protein-coding genes was estimated using unique peptides identified by mass spectrometry from both predicted protein coding genes and peptides generated by searching the theoretical tryptic peptidome of 6-reading-frame (6RF) translation of the M. sympodialis genome (Fig. 1D). A high proportion of identified peptides (30,559, i.e., 86% of the total number) overlapped between predicted protein-coding genes and 6RF translation of the genome. From the predicted proteome, tryptic peptides can be derived from sequence reaching over exon boundaries and other possible variants not present in the 6RF translation of the genome. In this study, 857 peptides were identified only in the predicted-protein database. In addition to the annotated proteome, 4,246 unique peptides were identified only in the 6RF translation. Thus, the comprehensive MS-based proteomics data confirm the majority of the predicted proteins and suggest that future analysis incorporating experimental peptide evidence has the potential to complement and refine the protein-coding-genome annotation.

When inspecting the synteny of orthologous proteins shared between M. globosa and M. sympodialis, we observed a tendency for genes which are distinct in the M. globosa genome to be fused into the same open reading frame (ORF) in the M. sympodialis genome. These discrepancies, as well as the difference between the two species regarding the total number of genes, mainly reflect differences between annotation platforms, the lack of RNA evidence for gene predictions in M. sympodialis, and, to a smaller extent, errors in sequencing. The annotation of the genomes of Malassezia species is particularly challenging, given the limited availability of transcript evidence (only 1,392 ESTs are available for M. globosa). Furthermore, the genomes of only a few species from this phylogenetic clade have been sequenced to date, such as U. maydis (15), Ustilago hordei (16), and Sporisorium reilianum (17). This renders the task of gene predictions via homology searches challenging, especially for fast-evolving and species- and genus-specific genes. For M. sympodialis, high-coverage sequencing of RNA extracted from distinct growth conditions is now under way to allow better annotation of the genome in a future update (A. Scheynius and T. L. Dawson, personal communication).

The mitochondrial genome was assembled from the 454 data and maps as a circular fragment of 38,622 bp with an estimated GC content of 32% (Table 1). Overall it is syntenic with the mitochondrial genome of M. globosa (Fig. 2). The genome contains all 15 expected protein-coding genes, 2 rRNA genes, and 25 tRNAs representing all 20 amino acids, with Ser, Leu, and Arg in two copies and Met in three copies. This coding content is the same as that of the related species U. maydis (accession number NC_008368). Although there is no synteny between Malassezia and Ustilago regarding gene order, genes are present on both strands in each species. In contrast to the mitochondrial genome of M. globosa, that of M. sympodialis has several group I introns (8 in total, including 3 that encode putative homing endonucleases).

A conspicuous feature of the M. sympodialis mitochondrial DNA (mtDNA) is a large (5.9-kb) inverted repeat containing the ATP9 gene and tRNA genes for Met, Leu, and Arg (Fig. 2). The inverted repeat is also present in the M. globosa mitochondrial genome, although it is shorter and poorly conserved between the two species, apart from regions corresponding to ATP9 and tRNAs. Large inverted repeats are common in chloroplast genomes (18), yet they occur infrequently in fungal mtDNAs and have been reported for only a few genera. In Candida species, homologous recombination between large inverted repeats in mtDNAs of Candida species has been proposed to play a role in replication and is associated with genome rearrangements (19, 20). Strand invasion structures in the inverted repeats of C. albicans mtDNA support a recombination-driven mechanism of DNA initiation (20), while comparative studies of the mitochondrial genomes of Candida species indicate that large inverted repeats are involved in conversions between circular and linear forms of the genome, as well as the formation of multipartite linear forms (19). Comparative analysis of mtDNAs of additional Malassezia species and clinical isolates is necessary to assess the functional significance of the inverted repeats and to determine whether mitochondria play a role in virulence, as has been reported with other fungal pathogens (21–23).

We investigated the metabolic pathways for potential changes that could explain the in vitro nutritional requirements of Malassezia species. Comparing the genomes of M. globosa and M. sympodialis for genes involved in lipid metabolism, we found that both genomes had a similar complement of genes. Neither genome encodes a recognizable fungal-type fatty acid synthase, while each genome encodes a plethora of lipid-hydrolyzing enzymes, such as lipases, phospholipases C, and acid sphingomyelinases (Table 2). M. sympodialis has a slightly reduced number of lipases and phospholipases C compared to M. globosa. It is not clear whether this contributes to any physiological differences or if this pattern is simply due to the potentially incomplete predicted proteome of M. sympodialis. No Δ2,3-enoyl coenzyme A isomerase gene was found in either genome, suggesting defects in utilizing unsaturated fatty acids. Notably, a Δ9 desaturase gene, absent in both the M. globosa and Malassezia restricta genomes, was found in the M. sympodialis genome (MSY001_2159). This suggests that M. sympodialis, in contrast to M. globosa, may be able to add a double bond to generate oleic acid if provided with stearic acid in the culture medium. However, since oleic acid is included in the mDixon medium used to grow Malassezia, the difference between the two species regarding the Δ9 desaturase is unlikely to explain the reported better in vitro growth of M. sympodialis (24), unless oleic acid uptake is limiting in M. globosa and M. restricta.

A second explanation for the differences in growth between the two species could relate to genes involved in sugar assimilation, such as those encoding β-glucosidases. In contrast to M. globosa, M. sympodialis is positive in the β-glucosidase enzymatic assay (25). In fungi these enzymes belong to two families, glycoside hydrolase 1 and 3 (GH1 and 3); the number of genes belonging to each family varies in basidiomycetes, with 0 to 2 copies for GH1 and 3 to 7 for GH3 (15). Similar to the U. maydis and Cryptococcus neoformans genomes, Malassezia genomes do not code for a GH1-type fungal β-glucosidase but, notably, show even greater compactness, with only one gene coding for a putative GH3-type enzyme (see Fig. S1 in the supplemental material). As the GH3-type gene is present in both Malassezia genomes, differences in expression may explain why M. globosa does not show a β-glucosidase enzymatic activity.

A unique characteristic of Malassezia is the cell wall, which is very thick compared to the cell walls of other yeasts (2). The M. sympodialis genome contains representatives of the major polysaccharide biosynthesis genes (Table 3). These include genes encoding six chitin synthases, proteins associated with β-1,6-glucan synthesis, and only one β-1,3-glucan synthase catalytic subunit. Other basidiomycetes also have one Fks-like β-1,3-glucan synthase catalytic subunit, whereas ascomycetes such as S. cerevisiae generally harbor a number of alternative Fks subunits (Table 3). There are six predicted chitin deacetylases, which is similar to the number found in other basidiomycetes but larger than the number generally seen in ascomycetes. This suggests that a large proportion of cell wall chitin may be converted to chitosan, the deacetylated form of chitin. The chromosomal location of this gene family suggests that there have been three separate gene duplication events leading to expansion of the chitin deacetylase family. Studies in C. neoformans have shown that chitosan helps to maintain cell wall integrity (26), suggesting that chitin deacetylases and the chitosan made by them may prove to be excellent antifungal targets. There is also an abundance of putative exoglucanases with similarity to S. cerevisiae Exg1 (Table 3). The classical cell wall integrity or protein kinase C (PKC) pathway seems to be highly conserved in M. sympodialis (see Table S1 in the supplemental material). Putative enzymes involved in O-glycosylation are represented in the genome such as the O protein mannosyltransferase (Pmt) family, but there is little evidence of orthologs of N-glycosylation enzymes, in particular those that add sugars to the outer chains to N-glycan (see Table S1).

In yeasts such as S. cerevisiae and Candida albicans, the major class of cell wall-localized proteins comprises proteins that are modified by the addition of a glycosylphosphatidylinositol (GPI) anchor (27). The posttranslational addition of a GPI anchor targets proteins to the plasma membrane. Through a poorly understood mechanism, the GPI anchor of a subset of GPI proteins is cleaved and the proteins are translocated to the wall and covalently attached to β-1,6-glucan. Only ten M. sympodialis proteins (and 20 in M. globosa) were predicted to become GPI anchored (Table 3), and surprisingly, they are likely to be associated with the membrane rather than the cell wall. This is a very small number compared to that in other fungal species; pathogens of the genus Candida can have over ten times more predicted GPI-anchored proteins (28). In C. albicans, these GPI-modified proteins include important adhesins that are involved in adhesion to host epithelial and endothelial cells as well as to innate surfaces, such as indwelling catheters, and thereby also contribute to biofilm formation (27). Both adhesion and biofilm formation play a role in virulence in invading pathogens that cause systemic and bloodstream infections. GPI-modified proteins also include carbohydrate-active enzymes with important roles in cell wall construction and maintenance of cell wall integrity. These proteins (transglucosidases and chitin-glucan cross-linkers) act by modulating the cell wall polysaccharides chitin and β-1,3-glucan. Only one putative chitin-glucan cross-linker gene, a Utr2 homolog, was identified in the M. sympodialis genome. No Dfg5/Dcw1 family members were identified. These putative mannosidases are generally assumed to be involved in the cleavage of GPI anchors, consistent with the lack of GPI cell wall proteins. This prompted us to look for homologs of the proteins that synthesize the GPI anchor itself, and we found representatives of the enzymes that synthesize most steps in the pathway in S. cerevisiae (see Table S1 in the supplemental material).

From the above, we can surmise that there are β-1,3-glucan, β-1,6 glucan (probably β-1,3/β-1,6 glucan), chitin, and chitosan in the cell wall of M. sympodialis, and glycosylation may be limited primarily to O-glycosylation. Analysis of the genome provides no evidence of α-glucan synthesis and no or very few “classical” fungal cell wall proteins. These in silico results are in agreement with a previous cell wall carbohydrate analysis that revealed that the M. sympodialis cell wall is composed primarily of β-1,6 glucan, with trace amounts of branched β-1,6/β-1,3-glucan and mannan (29). We further performed high-pressure freezing--transmission electron microscopy (HPF-TEM) to corroborate our observations on the absence of cell wall proteins and proteins for outer N-mannosylation, which participate in forming a fibrillar layer. M. sympodialis indeed lacks the extensive outer fibrillar layer (Fig. 3) that is evident on the wall of S. cerevisiae and C. albicans (30, 31).

Genes coding for all 10 allergens previously cloned from M. sympodialis (Mala s 1 and s 5 to s 13) and for the three allergens from M. furfur (Mala f 2 to 4) (9) were identified in the M. sympodialis genome (Table 4). For some of the allergens, where the previously described sequence was incomplete (Mala s 13) or contained sequencing errors at the 5′ end (Mala s 11), the availability of the genome sequence and protein alignments allowed identification of the full coding sequence, including the start codon. Sequence identity between M. sympodialis and M. globosa at both the nucleotide and amino acid levels is generally high (Table 4). Despite the high degree of protein identity between putative orthologs, a molecular evolution analysis indicated high levels of nucleotide substitutions (dS > 1) for all of them (Table 4). A population of 56 clinical M. sympodialis isolates from healthy individuals and atopic eczema patients (see Table S2 in the supplemental material) was further analyzed for the presence of two major allergens, Mala s 1, an allergen of unknown function (32), and Mala s 12, showing sequence similarity to the GMC oxidoreductase family (33). The partial genes for Mala s 1 and Mala s 12 (encoding 81% and 45% of the mature proteins, respectively) were amplified (primers are listed in Table S3 in the supplemental material) in all clinical isolates and showed strong conservation. This finding suggests that the regions of the genes we examined are under high selective constraints, possibly reflecting the maintenance of essential roles in the interaction with the host.

The function of Mala s 1 is still an enigma despite the availability of its three-dimensional (3-D) structure. The 3-D structure indicates that Mala s 1 is a β-propeller-folded protein. This novel fold among allergens has structural similarity in the potential homologs Q4P4P8 and Tri 14, from the plant pathogens U. maydis and Gibberella zeae, respectively (34), suggesting that Mala s 1 and the plant pathogen proteins may have similar functions. Because gene deletion approaches have not been established for Malassezia, we investigated the role of the Mala s 1 ortholog in the related smut fungus U. maydis. Quantitative real-time PCR revealed that expression of the Mala s 1 ortholog in U. maydis (um04915), encoding a protein predicted to be secreted, is induced during colonization of maize seedlings (see Fig. S2A in the supplemental material). The gene um04915 was deleted from the genome of the solopathogenic strain SG200 (see the supplemental methods in the supplemental material). With respect to SG200, um04915 mutants were unaltered in growth sensitivity to various stressors, including H₂O₂, sorbitol, Calcofluor white, and Congo red, or in filamentation on charcoal-containing plates (data not shown). Virulence was determined in seedling infections by comparing four independent mutant strains and SG200. No significant differences in disease symptoms were noted in comparison to SG200 (see Fig. S2B in the supplemental material). Therefore, the um04915 gene is not directly related to U. maydis pathogenicity in seedling infections; however, we cannot exclude the possibility that the gene may be required for disease in different maize organs (35) or has evolved different roles in Malassezia, associated with a different host.

Enhanced release of Mala s allergens and particularly Mala s 12 has been observed when M. sympodialis is cultured at a higher pH, which reflects that of the skin of atopic eczema patients (12). Here we predicted four of the known allergens to be secreted proteins (Table 4). Combining this observation with proteomics experiments on the M. globosa orthologs (6) and results from a previous study that showed that Mala s 1 and 12 are expressed on the cell surface of Malassezia (36), we suggest that these allergens may be exported and/or loosely associated with the cell wall, for example, via disulfide bonds (27) or, for Mala s 1, via binding to phosphoinositides involved in membrane trafficking (34). Notably, the characterized ortholog of Mala s 1 from the wheat pathogen Fusarium graminearum (Tri 14) was proposed to be functionally associated (either as a regulator or as a transporter) with an adjacent gene cluster involved in biosynthesis of a mycotoxin (37). This observation is of interest, as in both the M. sympodialis and M. globosa genomes, the gene encoding Mala s 1 is located adjacent to the gene encoding Mala s 8. Furthermore, the first 10 genes located on the 5′ end of Mala s 8 code for proteins potentially involved in secondary metabolism, such as a putative monooxygenase, a permease, a taurine dioxygenase, and a cobalamin-independent methionine synthase. A compelling hypothesis that merits further investigation is that Mala s 1 and Mala s 8 are involved in cell wall or postsecretory modifications of an as-yet-unidentified secondary metabolite. A few pathways of secondary metabolism have been observed in Malassezia, with some evidence for contributions to pathogenesis (4). Another potentially secreted allergen is Mala s 7 (Table 4), an allergen of unknown function (38, 39), which we identified here as a member of a novel family. Indeed, most of the 13 allergens represent single-copy genes in M. sympodialis, with a few belonging to small gene families with 2 to 5 gene copies each. The M. sympodialis genome has four genes predicted to encode proteins similar to the protein sequence of Mala s 7, in contrast to three genes in M. globosa. The identified sequences of all Mala s 7-like proteins from M. sympodialis and M. globosa bear signal peptides and do not have transmembrane domains or GPI anchors. Two of the genes in M. sympodialis are highly similar: Mala s 7a, which codes for the published allergen (38, 39) and a second, termed here Mala s 7b; both are located at the ends of relatively short scaffolds, showing a nucleotide identity of ~90% over a region of 4 kb, including both the complete Mala s 7 genes as well as the surrounding intergenic sequences. We confirmed by PCR (primers are listed in Table S3 in the supplemental material) that the observed duplication does not represent an assembly artifact (data not shown). Additional sequencing and mapping to chromosomes is required to resolve whether this is a segmental duplication or whether the duplicated fragments lie on distinct chromosomes. A segmental duplication event could also be at the origin of the other two Mala s 7 copies in M. sympodialis and M. globosa, as the respective genes in both species are located only 1.5 kb apart. These genes are Mala s 7c and 7d in M. sympodialis and MGL_0968 and a gene missed in previous annotations (accession number JX857443) in M. globosa. Gene conversion often occurs between copies lying in close proximity in the genome, and the species-based clustering of genes for Mala s 7-like proteins in the gene tree (see Fig. S3 in the supplemental material) is in line with this scenario. Another interesting observation is that the Mala s 7-like proteins in M. globosa (MGL_0968, MGL_2673, and JX857443) are shorter and have low similarity with their M. sympodialis counterparts at both the C and N termini, which could indicate that Mala s 7 in M. sympodialis comes from a fusion of two smaller Mala s 7-like proteins. A molecular evolution analysis using branch models in PAML (40, 41) did not further elucidate this family’s complex history, due to inconclusively high dS values. Overall, the genome sequence revealed a gene family amplification for the Mala s 7 allergen in M. sympodialis, which merits further investigation. Gene duplication is a major force driving evolution of new traits, including virulence, and thus it will be of interest to determine the roles of the secreted Mala s 7-like proteins in the interaction with the host and their evolutionary history.

We next addressed the function of Mala s 6, which is a member of the cyclophilin panallergen family (42). The sequence of the M. sympodialis Mala s 6 protein exhibited highest similarity to the cytoplasmic form of cyclophilin A, which is in agreement with the predicted absence of a secretion signal peptide (Table 4). Mala s 6 was found to be the most conserved protein in pairwise comparisons between M. globosa and M. sympodialis (Table 4), in line with its reported conservation across the tree of life. Using the C. neoformans model system, we investigated the functions of Mala s 6 using immunological, enzymatic, and drug inhibition assays. Recombinant Mala s 6 (43) reacted with a Cpa1-specific antisera (Fig. 4A) that successfully recognizes both the Cpa1 and Cpa2 cyclophilin A proteins in C. neoformans (44). Using a chymotrypsin-coupled assay that measures the cis-to-trans isomerization of a synthetic peptide (45), we found that the Mala s 6 recombinant protein shows robust cis-trans peptidyl-prolyl isomerase activity (Fig. 4B). Furthermore, the recombinant Mala s 6 protein is sensitive to inhibition by cyclosporine A (Fig. 4B), an effective immunosuppressive natural product that targets cyclophilin A and calcineurin (44). Cyclosporine A shows beneficial effects in atopic eczema, but due to side effects, its use is limited to patients with severe refractory disease (11). In conclusion, we demonstrated that Mala s 6 is a bona fide cyclophilin A and targeted by cyclosporine A.

We identified in the M. sympodialis genome assembly two scaffolds corresponding to the A (or pheromone/receptor [P/R]) and B (or homeodomain [HD]) mating type (MAT) loci based on sequence similarity shared with MAT genes of the closely related species M. globosa (6). The A and B mating-type loci in M. sympodialis were not linked in our assemblies; however, aligning them to the M. globosa MAT locus, where the two alleles are linked (6), identified three additional scaffolds which may lie between these loci and which could define a contiguous mating-type locus in M. sympodialis (Fig. 5A). Alignments of raw Illumina reads onto these scaffolds and PCR experiments (see the supplemental methods) provided further evidence for this contiguous assembly. Thus, similar to M. globosa (6), the M. sympodialis A and B loci appear to be physically linked and lie about ~141 kb apart.

In basidiomycetes, linkage of the A and B loci, along with strict biallelism, indicative of absence of recombination between these alleles, commonly defines bipolar mating systems, while tetrapolar mating systems contain unlinked and multiallelic A and B loci (46). However, in contrast to expectations for bipolar species, both the A and B MAT locus alleles of M. sympodialis showed extended flanking synteny with the corresponding M. globosa MAT locus regions (Fig. 5A). In comparison, in the bipolar species U. hordei, where the ancestral A and B MAT regions lie 430 to 500 kb apart, separated by a region that is highly rearranged between the two mating types, there is synteny on the 5′ end flanking the B locus and on the 3′ end flanking the A locus but not on their other flanks (47). Sequencing of the M. sympodialis A and B MAT alleles in a population sample of isolates (see Table S2 and the supplemental methods in the supplemental material) further indicated that the mating system of this species might not fit either of the traditional bipolar and tetrapolar systems, similar to what was previously reported for the pseudobipolar species Sporidiobolus salmonicolor (48). Below we present evidence suggesting that M. sympodialis has an intermediate mating system, termed the pseudobipolar mating system, where multiallelism might occur despite physical linkage (for a comparison of mating systems, see Fig. S4 in the supplemental material).

The A locus present in the genome of ATCC 42132 encodes a pheromone and a pheromone receptor arranged as two adjacent divergently oriented genes (Fig. 5A) and was designated a1. A MAT A locus allele that was distinct in terms of both sequence conservation and orientation of the genes was identified by sequencing the same region in a collection of isolates (see Table S2 and supplemental methods in the supplemental material) and designated a2. Both a1 and a2 are embedded in highly conserved, syntenic regions (Fig. 5B). In contrast to the a locus in the closely related species U. maydis, which is biallelic, the a1 and a2 MAT alleles from M. sympodialis do not contain a pheromone pseudogene, which in U. maydis is thought to descend from a tri-allelic P/R mating type system shared with the last common ancestor with S. reilianum, which is also tri-allelic (49–51). Moreover, the M. sympodialis A MAT locus alleles both lack the lga2 and rga2 genes that are present in the U. maydis a2 and S. reilianum a2 alleles, and which are involved in governing uniparental mitochondrial inheritance in those species (52, 53).

Similarly, the M. sympodialis B MAT locus of the reference strain ATCC 42132 (Fig. 5A) was designated allele b1 and contains divergently oriented genes encoding the homeodomain transcription factors bE and bW (by analogy with U. maydis). Again, PCR with flanking primers and sequence analysis revealed the presence of two additional alleles, b2 and b3, distinguished from b1 by a series of substitutions, many of which lie in the N-terminal region or in the protein-protein interaction domains of the two homeodomain factors and are nonsynonymous (Fig. 5B; also, see Fig. S5 in the supplemental material). The differences observed in the M. sympodialis B locus alleles may be sufficient to represent different mating types based on two lines of reasoning (54, 55). First, comparison to similar amino acid changes that are naturally occurring in U. maydis B locus alleles of different mating types reveals a similar degree of substitution and in similar regions of the proteins. Second, amino acid substitutions found in U. maydis mutants that exhibit compatibility with different-partner homeodomain proteins also exhibit a similar pattern to the changes observed between the b1, b2, and b3 mating type alleles of M. sympodialis. The finding of biallelism for the A locus and triallelism for the B MAT locus suggests that recombination could occur between these regions. Further evidence for this is based on an allele compatibility test that revealed five of the six possible allelic configurations predicted for two unlinked MAT loci (a1b1, a1b2, a2b1, a2b2, and a2b3) in a collection of M. sympodialis isolates (see Table S2 in the supplemental material). In contrast, if the three B alleles observed were simply the result of drift, we would have expected linkage between the A and B alleles (a1b1 and a2b2), but not the recombinant a1b2 or a2b1 combinations. These results suggest that the ~141 kb region separating the A and B MAT loci does not suppress recombination as in the biallelic MAT locus of U. hordei (56).

Overall, our comparative genomic and polymorphism analyses are consistent with a pseudobipolar mating system for M. sympodialis, distinct from tetrapolar mating systems such as the one observed in U. maydis and from strict biallelic bipolar systems, such as in the species U. hordei (47, 48, 54, 56–59). The pseudobipolar mating system of M. sympodialis may have arisen recently from a tetrapolar ancestor, so that large-scale rearrangements between the two linked loci that erase flanking synteny have not yet occurred, in contrast to other derived bipolar species, such as U. hordei. Multiple independent transitions from an ancestral tetrapolar state to a bipolar state and possibly also a pseudobipolar derived state in the Basidiomycota have frequently been reported (47, 60). In this view, M. globosa could also represent a pseudobipolar state derived from a tetrapolar ancestor; however, evidence for multiallelism in isolates of this species is needed to confirm this hypothesis.

Pheromone response during mating in fungi is regulated through a MAP kinase cascade, coupled to a heterotrimeric G protein consisting of α, β, and γ subunits. The genes of the MAP kinase module (e.g., FUS3, STE7, STE11) and the receptor for a-factor pheromone (STE3) are generally conserved in Malassezia (Table 5; for a full list, see Table S4 in the supplemental material). In contrast, the G protein subunits show more diverged profiles: of the three (or four in U. maydis [61] and S. reilianum) Gα subunits present in filamentous ascomycetes and many basidiomycetes (GPA1-4) (62), only one is present in Malassezia species (Table 5). This protein is most closely related to Gpa3 from U. maydis, required for mating (61). The Gβ subunit, associated with mating in many fungi (for examples, see reference 63), is an ortholog of Ste4 of S. cerevisiae (64); homologs of this protein are present in U. maydis, S. reilianum, and M. globosa but not in the M. sympodialis genome (Table 5). However, mating is regulated in U. maydis by a WD-40 protein related to Gβ subunits, Rak1 (65); Rak1 is conserved in all the Ustilagomycotina, including Malassezia species (see Table S4). It is possible that this protein interacts with the single Gα (Gpa3) in Malassezia as part of the mating signaling response. We could not identify a Gγ subunit (Ste18) ortholog in either of the Malassezia species. However, Gγ proteins are small and poorly conserved and can be difficult to find in genomic sequences.

Of the 29 genes defined as “core” for meiosis in eukaryotes (66, 67), only 19 are unequivocally present in the Malassezia genomes (Table 5). However, many losses are not specific to these species. For example, in most organisms, strand invasion is promoted following formation of double-strand breaks (DSB) by the activity of two proteins that arose from a gene duplication event that preceded the evolution of eukaryotes: Rad51, required during recombination in meiosis and for repairing DSBs in somatic cells, and Dmc1, which functions only in meiosis (68). Similarly to the genomes of other Ustilagomycotina species (M. globosa, U. maydis, and S. reilianum), the yeasts Candida guilliermondii and Candida lusitaniae (28, 69), the microsporidian species Encephalitozoon cuniculi, Caenorhabditis elegans, and Drosophila melanogaster, the M. sympodialis genome contains only a Rad51 homolog, no Dmc1. The loss of the DMC1 gene is correlated with the absence of genes (see Table S4 in the supplemental material) coding for the assembly factors Sae3 and Mei5 (70, 71). It is therefore possible that Rad51 alone is required for recombination in Malassezia, as previously proposed for C. elegans and D. melanogaster (72). Initiation of meiotic recombination in eukaryotes requires the formation of double-strand breaks in DNA by a complex containing Spo11. Spo11 orthologs are conserved in almost all eukaryotes, apart from one lineage of protists (73). It was initially difficult to identify the SPO11 gene in M. sympodialis (MSY001_2221). However, preliminary RNA-Seq data (T. Holm and A. Scheynius, unpublished data) provided evidence for a gene structure that contains six introns, one of which has a noncanonical splice site. SPO11 is expressed at a low level in two-day cultures of M. sympodialis, suggesting that the organism has retained the capacity to undergo meiosis. A putative multi-intron SPO11 candidate is also found in the M. globosa genome (J. Xu and C. W. Saunders, unpublished data).

Most eukaryotic genomes contain two paralogs of the essential cohesion complex, namely, REC8 and RAD21; the corresponding proteins are required for cohesion of sister chromatids during mitosis and meiosis (66). Rec8 is a meiosis-specific component (74). Surprisingly, the Malassezia genomes contain only one paralog gene, and this is more closely related to RAD21 than to REC8 (see Fig. S6 in the supplemental material). Both paralogs are present in the other Ustilagomycotina (U. maydis and S. reilianum) and in other basidiomycetes (e.g., C. neoformans) (see Table S4 in the supplemental material). The loss of Rec8 is unusual in fungi; however, this protein is also missing from the protists (66). It is possible that Rad21, which plays a role in meiosis in mammals (75, 76), may substitute for Rec8. Other meiotic subunits of the cohesion complex (Smc1, Smc3, Scc3, and Pds5) are present in Malassezia and related basidiomycetes (see Table S4 in the supplemental material).

Overall, the Malassezia species have lost a substantial number of genes associated with mating and meiosis (Table 5; also, see Table S4 in the supplemental material). However, few of the gene losses are unique; rather, most are shared with other species with an extant sexual cycle. Examples include genes encoding proteins involved in formation of the synaptonemal complex (SC) (HOP1, ZIP2, ZIP3, and RED1). These genes are absent not only in Malassezia spp. (see Table S4) but also in U. maydis (77, 78), a species in which sexual reproduction is well documented, which may suggest that the Ustilagomycotina do not form SCs at all, similar to Schizosaccharomyces pombe and some Candida species (28, 79). The genomic evidence is consistent with the possibility that there may be a sexual cycle in Malassezia resembling that of other Ustilagomycotina.

In this study, we identified independent lines of evidence for sexual reproduction in the Malassezia genus: (i) the presence of a MAT locus with apparently intact MAT alleles in the genomes of three Malassezia species (Fig. 5A) (6); (ii) evidence for recombination in a population of isolates of M. sympodialis, as shown by the discovery of segregating polymorphisms in the A and B regions of the MAT locus (Fig. 5B; also, see Fig. S3 in the supplemental material); and (iii) conservation of genes required for meiosis and signaling in the mating process in both M. globosa and M. sympodialis (Table 5; also, see Table S4 in the supplemental material). Although some of the core meiotic genes are absent, these are not consistently defined as “core” genes in the literature (for example, see reference 72), and their absence does not necessarily suggest an absence of sex (80). For example, genes absent from Malassezia but present in the other Ustilagomycotina include MSH4, MSH5, and MER3, required for the resolution of Holliday junctions during recombination in S. cerevisiae. All three are also missing from sexual fungi such as C. guilliermondii, C. lusitaniae (28), and S. pombe and from other eukaryotes with intact sexual cycles, including Plasmodium and Drosophila (66).

One way to confirm sexual reproduction and establish the role of MAT alleles in this process would be to conduct mating assays for fertility. In an attempt to detect an extant sexual cycle, a collection of M. sympodialis isolates (see the supplemental methods and Table S2 in the supplemental material) have been cocultured in pairwise and more complex mixtures under a variety of different media and environmental conditions (see the supplemental methods). However, despite these efforts, to date none of the morphological features associated with basidiomycete sexual reproduction (e.g., dikaryotic hyphae, clamp connections, basidia, and basidiospores) have been observed. Thus, the right combination of strains or conditions to detect an extant sexual cycle, if one exists, has not yet been established. The Malassezia sexual cycle might differ morphologically from that of other fungi and could require genetic approaches with marked strains to detect responses to pheromones, cell-cell fusion, or genetic recombination. No sexual cycle is known for any Malassezia species, but previous studies provide evidence for M. furfur hybrids that may result from mating between varieties or cryptic species (81; T. Boekhout, personal communication). Furthermore, recombination was observed in allozyme studies of the species M. pachydermatis (82).

In summary, the genome of M. sympodialis reported here, combined with previous studies of the M. restricta and M. globosa genomes, provides a rich foundation for future studies to elucidate the unique features of these ubiquitous commensals of human skin associated with multiple disease states. Moreover, apparent transitions in the mating type locus suggest that much remains to be learned about how, when, and where sexual reproduction might occur. Sexual reproduction may occur on human skin, with implications for antigens presented by the yeast, which might provoke immune reactions, leading to disease. This study provides insight into a number of hypotheses related to the life cycle of these fungi. The development of M. sympodialis transformation protocols for gene replacement studies will be a critical next step toward assessing the roles of genes potentially required for Malassezia species to become specialized to live on the skin as commensals but also to provoke disease.

DNA was isolated from M. sympodialis ATCC 42132 cultured on Dixon agar (24) modified to contain 1% (vol/vol) Tween 60, 1% (wt/vol) agar, and no oleic acid (mDixon) at 32°C for 4 days using the QIAamp DNA minikit (Qiagen GmbH, Hilden, Germany) according to the manufacturer’s instructions with small modifications. Briefly, glass beads were added to the cell suspension, which was vortexed for 4 min prior to the lysing incubation at 56°C for 3 h with additional vortexing during the incubation period.

Pyrosequencing was performed using 454 Titanium chemistries (Roche/454 Life Sciences, Branford, CT). Illumina libraries were made with the Illumina mate pair kit (3-kb insert) according to the manufacturer’s instructions, followed by 50-bp paired-end sequencing on one lane of an Illumina HiSeq instrument. Genome assembly of the 454 data was accomplished with the GS De Novo assembler, v. 2.3 (Roche Diagnostics, Basel, Switzerland). Contig extension and scaffolding were based on the Illumina data using SSPACE, v. 1.0 (83). We noted that a large fraction (~50%) of the Illumina mate-pair data in reality represented standard noncircularized paired ends, but no attempts to confirm or extend the scaffold structure by long-range PCRs were done at this point. The mitochondrial assembly was performed with Newbler 2.6, using a random subset of 40,000 reads from the 454 shotgun sequencing, representing ~25× coverage of the mitochondrial genome. Scaffolding and identification of the inverted repeat were aided by mapping of the Illumina 3-kb jumping library on the mitochondrial contigs. Finally, the genome was fully completed by manual identification and addition of 454 reads spanning the assembly gaps flanking the inverted repeat. CEGMA analysis was run using a set of 248 core eukaryotic genes (CEGs) as queries against the assemblies of M. sympodialis and M. globosa, with the completeness reported as a percentage reflecting the number of CEGs found as complete or partial genes, respectively.

The genome of M. sympodialis was annotated using the program MAKER versions 2.10 and 2.25 (84, 85). As evidence for gene annotations we used (i) protein alignments to a set of 67,086 publicly available proteins derived from the species M. globosa, U. maydis, C. neoformans, Fusarium graminearum, Magnaporthe grisea, Neurospora discreta, Neurospora tetrasperma, and Sordaria macrospora and clustered using Cd-hit (86) version 4.02 and a protein identity threshold of 90%; (ii) nucleotide alignments to 1,392 EST sequences coming from an M. globosa sequenced library (accession number IBEST_028020) previously used for gene predictions in M. globosa (6); and (iii) the ab initio predictors Augustus (87), using a U. maydis model, Genemark-ES (88), trained with the M. sympodialis scaffolds, and SNAP (89), trained within MAKER as follows. MAKER was run four consecutive times, and each annotation output file from MAKER (genome gff file) was converted into a model using the instructions from the SNAP documentation and provided as an input model in the next run. We complemented the predicted genes from the fourth MAKER run with (i) a set of 542 models, identified by protein alignments of the M. globosa proteome against ab initio models not retained by MAKER using the “pass-through” method implemented in MAKER version 2.25, and (ii) a few genes manually retrieved (e.g., SPO11, genes for Mala s 7-like proteins, and a pheromone gene) using tBLASTn and manual curation. The gene models from the above analyses have not been further curated. Mitochondrial assembly was performed with Newbler 2.6, using a random subset of 40,000 reads from the 454 shotgun sequencing, representing ~25× coverage of the mitochondrial genome. Scaffolding and identification of the inverted repeat were aided by mapping of the Illumina 3-kb jumping library on the mitochondrial contigs. The genome was completed by manual identification and addition of 454 reads spanning the assembly gaps flanking the inverted repeat. The mitochondrial DNA was annotated following steps described in reference 90; the mitochondrial map was generated using Geneious Pro, v.5.6.4 (Biomatters, Auckland, New Zealand). Details on identification of specific genes presented in tables are found in the supplemental methods.

Four replicates of M. sympodialis (ATCC 42132) were cultured on mDixon agar (see above) at 32°C for 2 and 15 days. Cells were harvested and washed twice with PBS (phosphate-buffered saline) by centrifugation at 1,200 × g for 5 min. Pellets were frozen at −80°C. To extract proteins for mass spectrometry analyses, 30 mg of every pellet was dissolved in 200 µl PBS and transferred to tubes containing 200 µl 425 to 600 µm acid-washed glass beads (Sigma-Aldrich, Sweden). The cells were disrupted by homogenization in a Precellys 24 tissue homogenizer (Bertin Technologies, France). Approximately 100 µl cell suspension from every sample was removed and transferred to a new tube and a 1:1 volume of lysis buffer was added to obtain a final concentration of 4% (wt/vol) SDS, 1 mM DTT (dithiothreitol), 25 mM HEPES (pH 7.6). The samples were heated at 95°C for 5 min followed by sonication twice for 30 s each time. Protein concentration was determined with a DC protein assay (Bio-Rad, Sweden). Samples were subsequently reduced by dithiothreitol and alkylated by iodoacetamide followed by overnight trypsination (Promega). The four biological replicates of each sample from the 2- and 15-day cultures were pooled and separated by immobilized pH gradient-isoelectric focusing (IPG-IEF) on a narrow-range pH 3.7 to 4.9 strip and a 3- to 10-gel strip as described previously (91). Extracted fractions from the IPG-IEF were separated using an Agilent 1200 Nano-LC system coupled to a Thermo Scientific LTQ Orbitrap Velos. Proteome discoverer 1.3 with Sequest-Percolator (Thermo Scientific) was used to search predicted proteins or a 6-reading frame translation of the M. sympodialis genome merged with the Bos taurus database (UniProt canonical sequences, 120,727) for protein identification, limited to a false discovery rate of <1%. Peptide matches to the B. taurus database were considered remnants from the culture medium and removed.

HPF of Malassezia isolates was carried out as described previously (31). Briefly, samples were prepared by high-pressure freezing with an EMPACT2 high-pressure freezer and rapid transport system (Leica Microsystems Ltd., Milton Keynes, United Kingdom). After freezing, cells were freeze-substituted in substitution reagent (1% [wt/vol] OsO₄ in acetone) with a Leica EMAFS2. Samples were then embedded in Spurr resin and additional infiltration was provided under a vacuum at 60°C before embedding in Leica FSP specimen containers and polymerizing at 60°C for 48 h. Semithin survey sections, 0.5 µm thick, were stained with 1% toluidine blue to identify areas containing cells. Ultrathin sections (60 nm) were prepared with a Diatome diamond knife on a Leica UC6 ultramicrotome and stained with uranyl acetate and lead citrate for examination with a Philips CM10 transmission microscope (FEI UK Ltd., Cambridge, United Kingdom) and imaging with a Gatan Bioscan 792 (Gatan United Kingdom, Abingdon, United Kingdom).

Orthologous genes were aligned using ClustalW 2.1 (92) or Muscle (93). Nonsynonymous (dN) and synonymous (dS) substitution frequencies were calculated using the method described in reference 94, as implemented in yn00 in the PAML package (94), except for the Mala s 7 family, where branch models were tested in PAML (40, 41). Gene trees for this analysis (see Fig. S3 in the supplemental material) and for the β-glucosidase (see Fig. S1 in the supplemental material) were constructed with PhyML software (95) using the LG model for amino acid equilibrium frequencies, allowing estimation of invariable sites and setting the rate categories number to 4. An optimized tree topology search was performed with a starting BioNJ tree and using the best-of-NNI and SPR method. Bootstrapping analysis was performed with 1,000 datasets.

The native clinical M. sympodialis isolates utilized for amplification of allergens, mating-type genes, and mating assays (see Table S2 in the supplemental material) were obtained from healthy individuals and from patients with moderate to severe atopic eczema at the Dermatology Unit, Karolinska University Hospital, Stockholm, Sweden, and the protocol was approved by the local ethics committee. The participants were instructed not to wash their upper back on the day of isolation. Samples were taken by holding a contact plate containing modified Leeming and Notman agar medium (96) against the skin of the upper back for 15 s. The contact plates were incubated at 32°C for 6 days. One colony from each plate was transferred to mDixon agar plates (see above) and cultured for 4 days at 32°C. The isolates were identified as M. sympodialis using the primers listed in Table S3 in the supplemental material.

For 56 clinical isolates (see Table S2 in the supplemental material), DNA corresponding to the partial gene sequences of Mala s 1 (915 bp) and Mala s 12 (913 bp) was amplified by PCR using primers (see Table S3 in the supplemental material) designed according to published sequences (32, 33). The PCR amplifications were carried out with Phusion high-fidelity DNA polymerase (New England Biolabs, Ipswich, MA) under the following cycling conditions: an initial denaturation at 98°C for 2 min followed by 30 cycles of 10 s at 98°C, 20 s at 64°C, and 15 s at 72°C and a final elongation step at 72°C for 7 min. The PCR products were purified using the QIAquick PCR purification kit (Qiagen) and sequenced at a verified core facility (KIGene, Karolinska Institutet, Stockholm, Sweden) using the same primers as the PCR amplification. Sequencing reactions were performed on an ABI 3730 Prism DNA analyzer (Applied Biosystems, Foster City, CA) using the BigDye terminator, v.3.1 (Applied Biosystems). The retrieved forward and reverse sequences were aligned using Geneious software, version 5.5.7 (Biomatters Ltd.), and the resulting consensus sequences were used for further analysis. The copies of Mala s 7a and 7b were confirmed in M. sympodialis ATCC 42132 by PCR amplification using the polymerase noted above and primers listed in Table S3 in the supplemental material. The following cycling conditions were used: an initial denaturation at 98°C for 1 min followed by 35 cycles of 15 s at 98°C, 15 s at 64°C, and 2.5 min at 72°C and a final elongation step at 72°C for 10 min. The product was analyzed by electrophoresis in a 1% (wt/vol) agarose gel (Invitrogen, Groningen, Netherlands) with SYBRSafe DNA gel staining (Invitrogen) in 1× TAE (Tris-acetate-EDTA).

Total proteins (50 µg) from C. neoformans strains, including wild-type H99 and cpa1 and cpa1 cpa2 mutants (44), were fractionated by 18% (wt/vol) SDS-PAGE in parallel with a protein extract from M. sympodialis ATCC 42132 (97) and recombinant Mala s 6 (43). Proteins were transferred to PVDF (polyvinylidene difluoride) membranes (Bio-Rad) and incubated with an anti-Cpa1 rabbit antiserum diluted 1:2,000 (44). Membranes were developed using the enhanced chemiluminescence (ECL) advanced detection kit (Amersham). Peptidylprolyl isomerization activity of Mala s 6 was assayed by an improved method described previously (45). Briefly, 1 ml reaction buffer containing 0.5 mg/ml chymotrypsin (Sigma) and 500 ng of the recombinant Mala s 6 protein (43) or 500 ng of the C. albicans Cyp1 cyclophilin A (Y. L. Chen and M. E. Cardenas-Corona, unpublished data) was pre-equilibrated to 10°C and then rapidly mixed in chilled cuvettes containing 10 µl of the substrate peptide (from a stock solution of 0.5 mM N-succinyl-Ala-Ala-Pro-Phe-p-nitroanilide [Sigma] dissolved in trifluoroethanol containing 470 mM LiCl). The cuvettes were immediately placed in the spectrophotometer, and the release of p-nitroanilide was monitored at 395 nm and 10°C with a Beckman DU-600 spectrophotometer. Cyclosporine A (LC Laboratories) from a 100 mM stock solution in methanol was added to 1-ml reaction mixtures at a final concentration of 1 µM prior to the addition of substrate peptide to test for inhibition of the enzyme activity.

The genomic scaffolds containing the MAT locus were retrieved using tBLASTn using as queries the M. globosa genes pra1 (MGL_0964), a putative pheromone gene (MGL_0963), bW1 (MGL_0883), and bE1 (MGL_0884). Whole-genome alignments between M. globosa and M. sympodialis genomes with Mummer (98) allowed identification of the three additional scaffolds located between the A (P/R) and B (HD) loci in M. globosa, which were not linked in the original M. sympodialis assembly. To further characterize the M. sympodialis MAT locus, we designed primers (see Table S3 in the supplemental material) to amplify the A and B regions from M. sympodialis isolates (see Table S2) cultured on modified Dixon agar (24) to containing 1% (wt/vol) peptone, 1% (wt/vol) desiccated ox bile, 1% (wt/vol) Tween 60, 2% (wt/vol) agar, and no oleic acid at 30°C for 4 days. Primers were designed based on alignments with M. globosa using Primer3 (99). The PCR products (~3 kb for the A locus and ~4 kb for the B locus) were purified on a 1% (wt/vol) agarose gel using a Qiaquick gel DNA extraction kit and used for direct DNA sequencing by primer walking. Sequencing reactions were carried out at the Genome Sequencing & Analysis Core Facility at the Duke Institute for Genome Sciences & Policy (IGSP), at Eton BioScience (Research Triangle Park, NC) and at CBS Fungal Biodiversity Centre, Utrecht, The Netherlands. For details on linkage analysis of MAT loci, mating assays, and mating type identification, see the supplemental methods.

The nuclear and mitochondrial genomes of M. sympodialis have been deposited in the EMBL database and were assigned accession numbers HE999549-HE999613 and HF558646, respectively. The EST data from M. globosa are accessible under accession number LIBEST_028020. ITS (internal transcribed spacer) and MAT locus sequences of different M. sympodialis isolates were deposited in GenBank under accession numbers JX964840 to JX964847, JX964800 to JX964802, and JX964848 to JX964850 (see Table S2 in the supplemental material).