Exploring the genomic diversity of black yeasts and relatives (Chaetothyriales, Ascomycota).

The order Chaetothyriales (Pezizomycotina, Ascomycetes) harbours obligatorily melanised fungi and includes numerous etiologic agents of chromoblastomycosis, phaeohyphomycosis and other diseases of vertebrate hosts. Diseases range from mild cutaneous to fatal cerebral or disseminated infections and affect humans and cold-blooded animals globally. In addition, Chaetothyriales comprise species with aquatic, rock-inhabiting, ant-associated, and mycoparasitic life-styles, as well as species that tolerate toxic compounds, suggesting a high degree of versatile extremotolerance. To understand their biology and divergent niche occupation, we sequenced and annotated a set of 23 genomes of main the human opportunists within the Chaetothyriales as well as related environmental species. Our analyses included fungi with diverse life-styles, namely opportunistic pathogens and closely related saprobes, to identify genomic adaptations related to pathogenesis. Furthermore, ecological preferences of Chaetothyriales were analysed, in conjuncture with the order-level phylogeny based on conserved ribosomal genes. General characteristics, phylogenomic relationships, transposable elements, sex-related genes, protein family evolution, genes related to protein degradation (MEROPS), carbohydrate-active enzymes (CAZymes), melanin synthesis and secondary metabolism were investigated and compared between species. Genome assemblies varied from 25.81 Mb (Capronia coronata) to 43.03 Mb (Cladophialophora immunda). The bantiana-clade contained the highest number of predicted genes (12 817 on average) as well as larger genomes. We found a low content of mobile elements, with DNA transposons from Tc1/Mariner superfamily being the most abundant across analysed species. Additionally, we identified a reduction of carbohydrate degrading enzymes, specifically many of the Glycosyl Hydrolase (GH) class, while most of the Pectin Lyase (PL) genes were lost in etiological agents of chromoblastomycosis and phaeohyphomycosis. An expansion was found in protein degrading peptidase enzyme families S12 (serine-type D-Ala-D-Ala carboxypeptidases) and M38 (isoaspartyl dipeptidases). Based on genomic information, a wide range of abilities of melanin biosynthesis was revealed; genes related to metabolically distinct DHN, DOPA and pyomelanin pathways were identified. The MAT (MAting Type) locus and other sex-related genes were recognized in all 23 black fungi. Members of the asexual genera Fonsecaea and Cladophialophora appear to be heterothallic with a single copy of either MAT-1-1 or MAT-1-2 in each individual. All Capronia species are homothallic as both MAT1-1 and MAT1-2 genes were found in each single genome. The genomic synteny of the MAT-locus flanking genes (SLA2-APN2-COX13) is not conserved in black fungi as is commonly observed in Eurotiomycetes, indicating a unique genomic context for MAT in those species. The heterokaryon (het) genes expansion associated with the low selective pressure at the MAT-locus suggests that a parasexual cycle may play an important role in generating diversity among those fungi.

Introduction

The order Chaetothyriales (Pezizomycotina, Ascomycetes) harbours melanised, non-lichenised fungi with a large morphological diversity. The order is included in the subclass Chaetothyriomycetidae along with the lichenised orders Verrucariales, Pyrenulales, and Celotheliales. Within the Chaetothyriales, at least five families are recognized: Chaetothyriaceae, Cyphellophoraceae, Epibryaceae, Herpotrichiellaceae, and Trichomeriaceae (Batista and Ciferri, 1962, Réblová et al., 2013), while some clades are as yet unassigned. The members of Chaetothyriales exhibit a complex ecological variation, and species are found in habitats characterised by extreme and adverse conditions, e.g. on rock surfaces in hot, arid climates, in toxic niches with hydrocarbons and heavy metals, and remarkably often occur in vertebrates as opportunistic pathogens (de Hoog 2014). Some species cause mutilating or even fatal infectious diseases, often in apparently healthy individuals. Recent studies sequenced rDNA from a large number of undescribed melanised fungi from ant colonies that clustered in various families of Chaetothyriales (Voglmayr et al., 2011, Nepel et al., 2014). The asexual morphs of members of Chaetothyriales show large morphological diversity, whereas the sexual morph shows limited variation over the entire order. Some genera produce budding cells or are entirely yeast-like, and hence the order is often referred to as “black yeasts and relatives” (BY) (Fig. 1).

Fig. 1

Phylogenomic distribution of Chaetothyriales and related ascomycetes used for comparative genomics. The majority of species are placed in the families Cyphellophoraceae (C) and Herpotrichiellaceae (H). The main characteristics such as niche, isolation source (red boxes – anthropophilic pathogens, orange boxes, zoophilic pathogens and green boxes geophilic), anamorphs, teleomorphs and sexual locus organization are displayed for each compared species.

The family contains species that generally are epiphytes, growing on the surface of plant leaves, but it is still unclear whether those species are plant pathogens or symbionts. The mycelium resides on the surface of plant leaves without truly penetrating the host plant cuticle (Chomnunti ). Members of this family are mainly distributed in tropical regions and are characterised by producing a sooty melanised mycelium resembling a loose network of hyphae covering the substrate. Ascomata are formed below the mycelial web and are easily released from the plant cuticle. Asexual Chaetothyriaceae are only reported for genera Chaetothyrium (Merismella) and Ceramothyrium (Stanhughesia) (Hyde ).

The family harbours a vast diversity of polyphyletic asexual morphs, which include both saprobic species on plant debris and clinically important species (Fig. 1) (Untereiner & Naveau 1998). Among the latter are causative agents of chromoblastomycosis, phaeohyphomycosis, disseminated infections, and primary cerebritis (McGinnis, 1983, Garnica et al., 2009). Main asexual genera are Cladophialophora, Exophiala, Fonsecaea, Phialophora, and Rhinocladiella, which all include opportunistic pathogens that cause a wide array of clinical syndromes in cold- and warm-blooded vertebrates (Crous et al., 2007, Seyedmousavi et al., 2013). Most species reproduce asexually with conidia generated by a filamentous phase, while members of the genus Exophiala show yeast-like budding. Occasionally meristematic growth is observed (Fig. 1) (de Hoog ). Muriform cell segmentation is the unique invasive form inside host tissue in chromoblastomycosis (da Silva et al., 2002, da Silva et al., 2008). Capronia is the homothallic sexual genus covering all asexual members of Herpotrichiellaceae. Ascomata are setose containing 8–32-spored asci; ascospores are pale to dark brown and are generally transversally septate or muriform (Untereiner 1995). Species of the family are generally found in nutrient-poor habitats such as showers and sinks in bathrooms or washing machines and dishwashers (Hamada and Abe, 2010, Lian and de Hoog, 2010, Zalar et al., 2011, Zupancic et al., 2016), while some thrive in extreme environments such as on rocks or in toxic niches (Badali et al., 2011, Seyedmousavi et al., 2011, de Hoog, 2014). A significant number of ant-associated undescribed species from carton galleries is also affiliated with this family (Voglmayr et al., 2011, Nepel et al., 2014).

The family is a small monophyletic group of species which are known through their asexual morphs only (Réblová ). Conidia may be hyaline and one-celled, but several species have pale brown, curved conidia with thin cross walls. Conidiogenous cells are inconspicuously phialidic and are cylindrical and intercalary, or swollen and lateral. This family includes mild opportunists on human skin and nails in Cyphellophora and Phialophora (Fig. 1) (Feng et al., 2012, Gao et al., 2015).

The family is composed by epiphytic species (Chomnunti ) and a large clade of rock-inhabiting species recently added (Isola ). Remarkably also the genus Arthrocladium clusters in the family, known for a single strain causing a fatal disseminated human infection (Nascimento ). The single sexual morph in the family is the genus Trichomerium, which morphologically is very similar to Capronia above (Chomnunti ). Trichomerium was first placed within the Chaetothyriaceae on the basis of morphological similarities of sooty mould-like mycelium, but later a separate family was erected using improved phylogenetic analyses (Chomnunti ). Ascomata of the Trichomerium species are spherical, covered by long, scattered setae, and contain 8-spored asci with septate, often brownish ascospores. Recently, phylogenetic studies also added some paraphyletic taxa, which morphologically are very deviant, such as the asexual species Brycekendrickomyces acaciae (Crous ). Also, some simple morphology known in the Herpotrichiellaceae is recurrent in the Trichomeriaceae in Cladophialophora modesta and Cl. proteae (Badali ). Meristematic, non-sporulating species were classified in the genera Knufia and Lithophila, a group of largely rock-inhabiting species (with the exception of the lichenicolous species Knufia peltigerae) within the Trichomeriaceae (Isola ). Numerous undescribed species of ant-associated fungi characterised by sooty mould-like mycelium are also contained within this family (Voglmayr et al., 2011, Nepel et al., 2014).

A recently proposed family is (Gueidan ), covering the genus Epibryon and the asexual morph Leptomeliola ptilidii, as well as some more simply structured asexual morphs that morphologically are classified as Cladophialophora sylvestris, Cl. humicola and Cl. minutissima (Badali et al., 2008, de Hoog et al., 2011, Gueidan et al., 2014). Several species are bryophilous fungi, but some have a rock-inhabiting life style, or occur in soil or on vascular plants. Ascomata are located superficially on or penetrating leaf tissue. Straight or curved dark setae cover globose to ovoid or pyriform, ostiolate, pale to dark brown to black ascomata, and the 8-spored asci are ovoid, ellipsoidal or subcylindrical, without apical structures and containing transversely septate, ellipsoidal to fusiform ascospores (Döbbeler, 1997, Gueidan et al., 2014).

Of the families of Chaetothyriales, the Herpotrichiellaceae species exhibit highly diversified life styles and show recurrent infection of a variety of vertebrate hosts (de Hoog 2014). Often opportunistic behaviour in human patients is partly explained by a saprobic behaviour combined with thermotolerance, as in Mucorales where resistance to high temperatures – often associated with other types of extremotolerance – is classically viewed as a prime virulence factor (Scholer ). Opportunistic species often possess dynamic and versatile pathways to sequester carbon from a wide range of substrates in the environment. By chance, when an opportunistic pathogen colonises its host, the abundance and diversity of genes associated with acquiring energy from particular carbon sources might be an advantage. Thus, metabolic plasticity combined with tolerance of adverse conditions could be considered as virulence factors in opportunistic fungi.

In Herpotrichiellaceae, warm- as well as cold-blooded vertebrates with intact immunity are commonly affected (Seyedmousavi ), suggesting the presence of intrinsic virulence factors that are independent from temperature. This led us to perform a comparative genome approach in order to comprehend the general background of ecology-driven traits, adaptation to harsh and toxic environments, and association with vertebrate hosts. The phylogeny of the family Herpotrichiellaceae has been intensively investigated for several years by multi-locus sequence analyses based on ITS, TEF1, BT2, and ACT1, and occasionally with other genes. de Hoog recognised six approximate clades, which showed somewhat different ecological trends (Fig. 1). The europaea-clade located in the basal position has recently been upgraded to family level as Cyphellophoraceae (Réblová ). The jeanselmei-clade is basal to the Herpotrichiellaceae s.s. and contains several clinically relevant species, next to species which were often derived from environments rich in toxic monoaromatic hydrocarbons (Zeng ). The dermatitidis-clade contains thermophilic Exophiala species from hot, low-nutrient water systems, sometimes causing disseminated infections in humans (de Hoog ). The salmonis-clade harbours mainly mesophilic water-borne Exophiala species, often infecting aquatic animals such as fish and amphibians, but rarely humans (de Hoog ). The two remaining clades comprise the major agents of phaeohyphomycosis and chromoblastomycosis, but species can also be found in the environment on plant debris (Salgado et al., 2004, Vicente et al., 2008). The carrionii-clade harbours some species that consistently cause chromoblastomycosis and which may perhaps be regarded as primary human pathogens (de Hoog ). The same pattern is observed in the bantiana-clade, which harbours Fonsecaea and Cladophialophora, with an abundance of species causing serious human diseases (de Hoog et al., 2011, Najafzadeh et al., 2011a, Najafzadeh et al., 2011b, Sun et al., 2012) as well in Rhinocladiella mackenziei. The trends in all clades are approximate since pathogenic species are often flanked by free-living species. Also herpotrichiellacean asexual and sexual morph genera are polyphyletic, but as yet molecular phylogeny is too unstable to replace morphology-based taxonomy (Untereiner, 1995, Untereiner and Naveau, 1998, Haase et al., 1999, de Hoog et al., 2011).

The origin of Chaetothyriales is estimated at approximately 229 MYA during the Middle Triassic (Gueidan ). It has been suggested that the Permian–Triassic (P–T) mass extinction, which deeply affected terrestrial and marine ecosystems, led to the development of a thermotolerant life-style on rock, possibly in association with toxin-producing lichens. After this, a rapid diversification of Chaetothyriales took place. In this vision, extremo- and thermotolerance, and an efficient metabolism of carbon sources are atavisms from this period (Gueidan ). The five families proposed in Chaetothyriales all contain a number of basal rock-inhabiting species with epiphytic or epilithic growth, suggesting a common origin of these life styles (Gueidan ).

The dark colouration of chaetothyrialean mycelium is determined by the high production of melanin pigments, which was shown to contribute to the above discussed ecological niches as well contributing to resistance against host immune responses (Schnitzler et al., 1999, Zhang et al., 2013). The presence of melanin alone is not sufficient to explain pathogenicity as these polymers are known to be present in many Pezizomycotina, and additional factors discussed above may be involved to explain the pathogenic status of these fungi. The virulence of opportunistic black yeasts has been suggested to have evolved from adaptations to extreme environments, e.g. melanisation (Schnitzler et al., 1999, Feng et al., 2001), meristematic growth (Mendoza et al., 1993, Karuppayil and Szaniszlo, 1997), and general extremotolerance (Liu ). Application of concepts of “focused” virulence and “dual ecology” may be considered for chaetothyrialean fungi to explain their ability to infect vertebrate hosts (Casadevall ). Although the source of many black fungal infections are plant-debris and occasionally living plants, their association as common degraders of plant biomass could be a misconception. In order to understand the basic biology of Herpotrichiellaceae, their phenomenal adaptation to extreme environments, and mechanisms associated with infection of vertebrate hosts we sequenced the genomes of 23 BY type species and compared them to related pathogens in Eurotiales and Onygenales (Fig. 1). The general genomic characteristics (i.e., genome size, synteny, gene content, repetitive elements), phylogenomic tree, transposable elements, sex-related genes, gene family expansions and contractions, evolution of protein- and carbohydrate-degrading genes, and secondary metabolism were deeply investigated in order to understand processes of adaptation of Chaetothyriales to multiple environments.

Materials and methods

rDNA LSU phylogeny

Phylogenetic assessment was carried out for all 172 black yeast fungal strains deposited at the CBS-KNAW Fungal Biodiversity Centre (CBS), Utrecht, The Netherlands (Table S1). LSU rDNA sequences were retrieved from GenBank and aligned by means of MAFFT v. 7.273 (Katoh & Standley 2013). Isolates and GenBank accession numbers are listed in Table S1. Phylogenetic analyses using Maximum Likelihood (ML) and a Neighbour-Joining (NJ) were performed by MEGA v. 6 (Tamura ) with Kimura 2-parameter model and statistical bootstrapping procedure involving 500 replicates.

Strains, DNA and RNA extraction

A set of 23 black fungal ex-type strains was obtained from CBS-KNAW Fungal Biodiversity Centre and cultivated in Malt Extract Broth (MEB) for 7 d with shaking at 150 rpm at 25 °C (Table S2). DNA extraction was performed via a cetyltrimethylammonium bromide (CTAB)-based method and phenol-chloroform/isoamyl alcohol purification (Möller ). Total DNA was purified with Qiagen Genomic Buffer Set and the Qiagen Genomic-tip 100/G. Total RNA was isolated with RNEASY Mini kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. The additional strain Cl. carrionii KSF (dH 23894) DNA was obtained from 7-d-old mycelia cultured on Sabouraud Glucose Agar (BBL™) at 25 °C. DNA was extracted using the DNeasy Plant Mini Kit (Qiagen) according to manufacture protocols.

Genome assembly and gene prediction and annotation

The genome of Cl. carrionii KSF (dH 23894) was pyrosequenced using the platform 454 GS FLX (Roche). Shotgun and 3Kb paired-end libraries were sequenced using the GS FLX Titanium XLR70 chemistry (∼450 bp reads). This genome was assembled using the NEWBLER software combining both paired-end and shotgun libraries (Margulies ). The P. attae genome was sequenced and annotated as previously described (Moreno ). For the other 21 species, the genomes were sequenced using Illumina technology. The genome of E. dermatitidis was previously described (Chen ). For the 20 remaining species, genomic DNA was used to construct two libraries with approximate insert size of 180 bp and 3 kb; for F. multimorphosa only a 180 bases-insert library was constructed. Each library was sequenced on an Illumina HiSeq 2000 to generate 101 base paired-end reads. All sequence was assembled using Allpaths (version R48559 for most assemblies); assemblies were inspected for regions of aberrant coverage, %GC, or sequence similarity using GAEMR (www.broadinstitute.org/software/gaemr) and contaminating sequence including was removed.

Genes were predicted and annotated by combining calls from multiple methods. A training set was generated using Genewise (Birney ) and Genemark (Lomsadze ), and then GlimmerHmm (Majoros ), Snap (Korf 2004) and Augustus (Stanke & Waack 2003) was used to generate ab initio gene models. For seven species, strand-specific libraries were constructed from total RNA using the Illumina TruSeq RNA Library prep. For each species, paired 76 base reads were generated on an Illumina HiSeq 2000. RNA-Seq was assembled using Trinity (Grabherr ) (version r20140413p1) in genome-guided mode (with parameters genome_guided_max_intron 10000 – SS_lib_type RF – trimmomatic – min_kmer_cov 2). All assembled transcripts were aligned to the genome using PASA (Haas ) and used to update gene models, predict alternatively spliced transcripts, and add UTR predictions. In addition, any ORF present in the PASA transcripts that did not overlap a gene prediction was used to recover missed genes. The best gene model at a given locus was selected from these data sets using EVidenceModeler (EVM) (Haas ); conserved genes missing in gene sets were identified using OrthoMCL (Li ) and combined with the EVM set (Haas ). All raw sequence data, assemblies, and annotations were submitted to NCBI (Finn ) (Table S2).

Annotation of transposons

In order to ensure a robust detection of repeat element, we used inverted repeat finder (IRF) (Warburton ) and Repeat Modeler (http://www.repeatmasker.org/RepeatModeler.html). IRF was set to identify pairs of repeats within a given of 20 kb. False positives candidates were filtered using the reference Pfam profile (using pfam_scan.pl with E-value threshold 0.00001) and RPS-BLAST against CDD profiles (with E-value threshold 0.001) (Finn et al., 2010, Marchler-Bauer et al., 2011). Multiple overlapping hits, were removed by cd-hit (Fu ) clustering with sequence similarity threshold set to 100 and query coverage set to 99 % of the shorter sequence. The resulting customized reference was merged with RepBase and used as input for Repeat Masker searches (Jurka ). All resulting sequences were translated in six frames and searched against a fixed list of reference Pfam HMM (Hidden Markov Model) profiles (using pfam_scan.pl with E-value threshold 0.01) and RPS-BLAST against CDD profiles (with E-value threshold 0.001). Transposon classification was curated manually based on the encoded protein domains.

Annotation of CYP genes

Identification of Cytochrome p450 monooxygenases (CYPs) were carried out by HMMR v. 3.1 (Finn ) which was used to perform sequence-profile HMM searches with the PFAM (Finn ) profile PF00067 (downloaded from the PFAM protein families database, http://pfam.xfam.org/, last accessed September 16, 2014) against all 23 black yeast proteomes. Proteins that achieved the cut-off 1e−03 were submitted to BLASTP searches against the fungal p450 CYPs database (Nelson 2009) (http://blast.uthsc.edu). The predicted CYPs p450 were assigned to family and subfamily types based on their BLASTP sequence identity. As recommended by the International P450 Nomenclature Committee, the cut-off of sequence identity was set at 40 % for family and 55 % for subfamily levels. Partial CYP p450 sequences (BLASTP identity >40 % and coverage <40 %) were classified as potential pseudogenes.

Annotation of transporter genes

Transporter gene classification was achieved with best match BLASTP (E-value threshold 1e−05, and at least 50 % alignment-length coverage) to transporter sequences available at Transporter Classification Database (TCDB) (Saier ).

Single-copy orthologue extraction and species tree inference

Clustering of single-copy orthologues across multiple fungal species was performed using ORTHOMCL (Li ) version 1.4 with a Markov inflation index of 1.5 and a maximum e-value of 1 × 10−5. Individual amino-acid sequences were aligned with MUSCLE (Edgar 2004) and poorly aligned regions were automatically removed using TRIMAL (Capella-Gutierrez ) under the “-automated1” setting. The sequences were concatenated with FASCONCAT (Kuck & Meusemann 2010) v. 1.0 and species trees were inferred by maximum likelihood RAxML (Stamatakis 2006) using PROTGAMMABLOSUM62 and 1000 bootstraps was used to infer branch support. Beyond the 23 herein analysed black yeast-like fungi, the following outgroups from the orders Eurotiales and Onygenales were applied: Trichophyton rubrum, Coccidioides immitis, Paracoccidioides brasiliensis, Aspergillus nidulans and A. fumigatus (Fig. 1).

Genome-scale chaetothyrialean phylogeny and divergence times

The phylogenomic position of Chaetothyriales was inferred based on 264 single-copy orthologous protein clusters identified among 53 fungal species as mentioned above. Concatenated amino-acid sequences were aligned using MUSCLE (Edgar 2004). In order to select the most-reliable positions in the alignment, TRIMAL (Capella-Gutierrez ) was used to eliminate poorly aligned regions (-automated1 option) resulting in 124 693 amino acid positions in the final alignment. Phylogenetic tree and branch lengths were inferred by Maximum Likelihood via a stochastic algorithm implemented in IQ-TREE software (Nguyen ). Best-fit amino acid model selection was assessed using an automatic model selection (MODELFINDER) and also considering the FREERATE model (-m TESTNEW option), which assesses the fit of multiple of multiple mixture GTR within the same model, in many cases having a better fit when compared to models that use a single parametric distribution (Soubrier ). Phylogenetic branch support was inferred by the ultrafast bootstrap approximation approach (UFBOOT), a measure that is better correlated to the actual probability of existence of a branch than the usual bootstrap (Minh ). Divergence times were inferred using the RELTIME method (Tamura ) implemented in the MEGA 7 software (Kumar ) using the LG model (Le & Gascuel 2008). Batrachochytrium dendrobatidis (Chytridiomycota) was used as outgroup and three calibration constraints were considered for divergence time estimations: (1) Basidiomycota/Ascomycota split: 390–1 490 MYA; (2) Pezizomycetes crown: 230–970 MYA; and (3) Sordariomycetes stem: 210–890 MYA. Those calibrations were based on conservative intervals considering both primary (fossil) and secondary calibrations discussed in Lucking . A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories; +G, parameter = 0.6655).

Functional domain prediction gains and losses

To identify functional domain gains and losses, INTERPRO (Mitchell ) domains were predicted using INTERPROSCAN (Jones ) in all 23 black yeast-like species and in the outgroups with previously released genomes: Trichophyton rubrum, Coccidioides immitis, Paracoccidioides brasiliensis (order Onygenales), Aspergillus nidulans and Aspergillus fumigatus (order Eurotiales). Gene family evolution was estimated with CAFÉ (De Bie ) v. 3.0 using a significance family-wide p-value threshold of <0.05 and VITERBI p-values of <0.001. To search for birth (λ) values we ran the program with the “-s” option. Two files were used as input in CAFÉ analyses: a table containing the organism's number of copies of each INTERPRO domain and an ultrametric tree constructed based on the species tree using a custom R script.

Results

Phylogeny

The aligned LSU dataset of 172 black fungal strains was used to determine a phylogenetic tree of the entire order Chaetothyriales; the LSU gene was sufficiently conserved to allow confident comparison over the entire dataset. Both Maximum Likelihood and Neighbour-joining analyses produced corresponding trees in which the same clades were supported (Fig. 2). Moreover, the tree topology was congruent with previously reported phylogenies of Chaetothyriales (Réblová et al., 2013, Gueidan et al., 2014), supporting the presence of five distinctive families: Chaetothyriaceae, Cyphellophoraceae, Epibryaceae, Herpotrichiellaceae, and Trichomeriaceae. In the most represented family Herpotrichiellaceae, the species were resolved in six clades with different ecological preferences as reported by de Hoog . Overall, this family included several clinically relevant fungi as well as species isolated from a variety of environmental sources, especially sites contaminated with toxic monoaromatic hydrocarbons. Two sub-clades at family level resolution were identifies within Herpotrichiellaceae: The upper clade harbour most of the Exophiala and Rhinocladiella asexual morphs while the lower clade is overrepresented by the genus Fonsecaea and Cladophialophora. Similarly, the family Cyphellophoraceae, which forms a supported monophyletic group, harbours both saprobic and medically important species responsible for mild opportunistic infections in human and animals. In contrast, the majority of the isolates belonging to the family Trichomeriaceae have an inert surface-inhabiting life style, while several are epiphytic. Arthrocladium fulminans seems to be the unique isolate causing a fatal disseminated human disease clustering in this family. The family Epibryaceae is located in the basal position, forming a distinct clade from other Chaetothyriales, but at relatively long branches (Fig. 2). Most of the isolates of this family are living plant associated.

Fig. 2

Phylogenetic analysis of members of Chaetothyriales (Class Eurotiomycetes). The Maximum likelihood tree, based on 172 LSU sequences, was determined using MEGA v. 6 with Kimura 2-parameter model with default settings and statistical bootstrapping procedure involving 500 replicates. Bootstrap values above 70 % are shown at the nodes. Family boundaries are indicated with coloured blocks. The tree was rooted to Verrucula granulosaria AFTOL-ID 2304.

Genome assembly and annotation

The assemblies were highly contiguous, with 12 consisting of 19 or fewer scaffolds, suggesting that many correspond to complete chromosomes. Genome assembly size varied from 25.8 Mbp for Capronia coronata CBS 617.96 to 43.3 Mbp for Cladophialophora immunda (CBS 834.96; Table S2, Fig. S1). Repetitive element identification was considered particularly low (ranging from 0.03 % to 5.2 %; Table S3) compared to other fungal species (Galagan et al., 2003, Martinez et al., 2012, Teixeira et al., 2014). This suggests that repeat content might not play an important role in determining genome size in black yeast-like fungi.

Genes were predicted combining de novo reconstruction of transcriptomes from RNA-seq data for some species and with ab initio and sequence homology based gene models. Corresponding with genome assembly sizes, high gene counts were found in Capronia coronata (9 231 predicted genes) and Cladophialophora immunda (14 033 predicted genes) (Table S2, Fig. S1). However, we did not observe a phylogenetic correlation between genome size and total gene number in the species examined (Fig. S1). Species of the jeanselmei- and bantiana-clades mostly experienced an increase in genome size as well as in predicted Open Reading Frames (ORFs) compared to ancestral populations (Fig. S1). Exceptions were E. aquamarina CBS 119918 with 41.7 Mb, while E. sideris CBS 121828 had a size of 29.5 Mb, as small as members of the carrionii-clade and similar to those of the dermatitidis-clade. In contrast, members of the dermatitidis-clade experienced a notable decrease in genome size and gene content (Fig. S1). Within Ascomycota, BY genomes had the highest percentages of G+C content reported to date, i.e., varying from 49 % in E. aquamarina to 54.3 % in Cl. carrionii, which could contribute to their thermotolerance (Nishio ). This corresponds with low single dinucleotide repetitions found in BY genomes.

Using ORTHOMCL clustered proteins, we determined the protein core families that were conserved in all black yeasts under investigation and other related fungi. This resulted in 4 031 genes per genome in the core set conserved in all species (Fig. 3). The KOG annotation for these amino-acid sequences revealed that proteins responsible for housekeeping functions, particularly for translation and RNA processing, were more represented in the core set (Fig. S2). We also assessed the proteins specific to each clade. We considered as clade-specific proteins those proteins that were present in orthologous groups found in a unique clade but were absent from all others. Non-core proteins may provide insight into specific processes and may be indicative of certain ecological preferences. For example, enzymes related to metabolism of carbohydrate (G) were found to be over-represented in the jeanselmei-clade (p-value = 1e−04; Fisher's exact test). Similarly, enzymes associated to secondary metabolites (Q) were found to be enriched in the bantiana- (p-value = 3e−13; Fisher's exact test), salmonis- (p-value < 1e−08; Fisher's exact test) and jeanselmei-clades (p-value = 7e−08; Fisher's exact test).

Fig. 3

Distribution of orthology classes in black yeasts and closely related fungi: core genes found in all genomes are shown in green, shared genes present in more than one but not all genomes in blue and genes that were unique to only one of the 28 analysed genomes in yellow.

On the other hand, the dermatitidis-clade proteins were under-represented for these functions (G: p-value 2–e01; Q: p-value 9–e02; Fisher's exact test) suggesting a reduced secondary metabolite producing capacity (Fig. S2).

Transposable elements

The members of the families Herpotrichiellaceae and Cyphellophoraceae have low content of transposable elements (Fig. S3, Table S3). Prevention of accumulation of transposable elements in BY genomes might be driven by the hyper-mutation process of repeat-induced point mutation (RIP). The scarcity of transposable elements results in decreased abundance of transposon encoded proteins such as reverse transcriptase (RT domain -IPR00477). Despite the low incidence of repetitive elements in BY genomes, we detect several TEs in the bantiana- and jeanselmei-clades, especially of the DNA Transposons LINE and the LTRs retrotransposons when compared to remaining clades (Fig S3). Rhinocladiella mackenziei, not assigned to any clade, also contained a higher number of elements with some specific expansions, such as the Helitron class (Table S3). The E. aquamarina genome presented the highest number of TEs (5.2 %), possibly reflecting its relatively significant genome expansion compared to other BY genomes (Fig. S3). It has been described previously that eukaryotic genomes of moderate sizes tend to have a linear correlation between complexity and genome size (Metcalfe & Casane 2013). Black yeast moderate genome sizes correlate well with the scarcity of repeats.

Within the Onygenales, which are generally related to animal hosts either as saprophytes or pathogens, there are organisms with small compact genomes and others with expanded complex genomes. The transposon-rich Ajellomycetaceae (Blastomyces, Histoplasma, and Paracoccidioides) and Onygenaceae (Coccidioides) compared to dermatophyte Arthrodermataceae (Trichosporon, Arthroderma, Microsporum), which have streamlined genomes with single repeats. The opportunistic Onygenales seem to have a more diverse TE landscape whereas specialised dermatophytes reduced their genomes. Blastomyces transposons have expanded up to 63 % of the genomes in a low GC genomic environment (Muñoz ). Lower GC is expected to favour mobile element integration (Wicker ). Additionally during genome duplication often mobile elements proliferate (Ma ).

Vma1-a inteins reveal a new evolutionary trend

We detected the presence of self-cleaving parasite proteins of the MEROPS N09 family, nested within Asparagine Peptide Lyases among some of the BY genomes. The N09 domain is commonly found within intein-containing V-type proton ATPase catalytic subunit A in several species of yeasts and genera of Archaea, i.e. Thermoplasma and Pyrococcus (Perler 2002) (vma-1a and vma-1b inteins, respectively) (Perler et al., 1994, Liu, 2000). This mobile element is spliced out from host protein sequences (or exteins) after its translation through an autocatalytic process. This parasite domain, which was suggested to have been acquired by a process of horizontal gene transfer (HGT), has a high sequence similarity with Archaea / Bacteria because of its convergent evolution along the fungal tree of life (Goddard and Burt, 1999, Poulter et al., 2007, Swithers et al., 2013). We detected the presence of vma-1a class intein in the R. mackenziei (Z518_00231) and F. pedrosoi (Z517_06303) genomes, suggesting a broader distribution within Ascomycetes (Fig. S4). We extended our analysis throughout other Pezizomycotina using INBASE (Perler 2002), and the vma-1a intein was also found in members of Sordariomycotina, i.e. Sporothrix schenckii, S. brasiliensis, and Stachybotrys chartarum genomes, bringing the total number of non-yeast species with inteins herein to five. In contrast, this self-cleaving protein is widely distributed among Saccharomycotina. Those five Pezizomycotina species are nested in a monophyletic clade apart from remaining species of Saccharomycotina, which may represent a new class of this element. The Host V-type proton ATPase protein, splicing and DOD homing endonuclease motifs were all identified and conserved with Candida glabrata vma 1-a (Fig. S4). However, the DOD homing endonuclease motif blocks D and E do not seem to be conserved with vma 1-a in Saccharomycetales. On the other hand, motif blocks D and E appear to be highly conserved with those presented in the PRP8 intein among Pezizomycotina (data not shown). Here we report on additional vma-1a intein, showing that they are found in more diverse fungal species within Pezizomycotina. V-ATPases are in general responsible for acidification of a variety of intracellular compartments, especially the vacuolar membrane vesicles of Eukaryotes. These mobiles genetic elements are self-spliced due stress adaptation (Senejani et al., 2001, Topilina et al., 2015, Novikova et al., 2016) and may play an important role in the regulation of extremotolerance of many BYs.

Origin of black yeast species and divergence times

A multigene phylogenetic species tree for a broad panel of 53 fungal species was generated using 264 single-copy orthologues. Representatives of chaetothyrialean families other than Herpotrichiellaceae and Cyphellophoraceae are still missing. With this limitation, two groups are sufficiently remote to conclude that as yet the order Chaetothyriales harbours two monophyletic families, Herpotrichiellaceae and Cyphellophoraceae which ancestry is found around 130 MYA during the cretaceous period (Fig. 4). Based on morphological and molecular methods with conserved genes, Coniosporium apollinis has previously been placed early in the Chaetothyriales. However, using a genome-scale phylogenetic tree, we demonstrated that this fungus is more closely related to the order Botryosphaeriales in the Dothideomycetes (Fig. 4). The Chaetothyriales are close to Verrucariales and Phaeomoniellales. These orders are displayed as paraphyletic branches and compose, along with Onygenales and Eurotiales, the subphylum Eurotiomycetes. Judging from the calibrated phylogenetic tree, the early and major BY lineages of Herpotrichiellaceae and Cyphellophoraceae are contemporaneous and emerged around 75–50 MYA during/after the Cretaceous-Paleogene (K-Pg) extinction event (Fig. 4). It is worth noting that the radiation of Chaetothyriales took place more recently than that of Onygenales.

Fig. 4

Genome-scale of chaetothyrialean phylogeny and divergence times. Calibration points are highlighted in blue and were used to infer the divergence times for Chaetothyriales (upper panel). The red node displays the divergence dates of Chaetothyriales and the red asterisk bolded area highlights a common era for both Cyphellophoraceae and Herpotrichiellaceae. The bottom scale presents the main geological and periods and eras.

Gene family evolution in black yeast

INTERPROSCAN was used to identify protein domains in all 23 black yeasts and in related fungi in Eurotiomycetes, including Trichophyton rubrum, Coccidioides immitis, Paracoccidioides brasiliensis (order Onygenales), and Aspergillus nidulans and Aspergillus fumigatus (order Eurotiales). The species tree was inferred by maximum likelihood RAxML (Stamatakis 2006) based on the concatenated amino acid sequences of 4 031 single-copy orthologous genes shared by all 23 species. Domain gain events (expansions) and domain loss events (contractions) were estimated with CAFÉ (De Bie ) in each black yeast and ancestral node of the species tree. The dynamic evolution of protein domain families in black yeast is shown in Fig. 5A.

Fig. 5

Dynamic evolution of protein families. (A) Phylogenetic tree showing the relationship between species and altered protein families. Pie diagrams and numbers at the nodes represent the abundance of contractions (red) expansion (blue) and No change (black) of 1771 protein families during evolution of black yeasts. (B) Heatmap showing expansion and contractions of protein families found in species belonging to the bantiana- and carrionii-clades, respectively. Domains are grouped by category similarity. All domains shown are significantly changed, and were identified using CAFE with cut-off of family p-values <0.05 and Viterbi p-values <0.01.

Overall, we observed 46 genomic novelties associated with protein domain expansions and contractions, which arose early in the evolution of these fungi and that were present in the common ancestor of Chaetothyriales examined (Table S4). We speculate these expanded domain families have provided selective advantage and extremotolerance for adaptation to ecological niches that are subjected to environmental stress. Black yeasts are known for their extremotolerance and are able to grow and thrive in hostile habitats, such as those containing toxic compounds, high and low temperature, scarcity of nutrients, or conditions of dryness (Gostincar ). We assessed a correlation between the seven domain family expansions and the ecological preferences in herpotrichiellaceous black yeast. These functional domains are likely to be involved in metabolic processes with oxidoreductase activity. Among them, four domains are related to Aldehyde dehydrogenases (ALDHs), which catalyse the oxidation of different aldehydes to their corresponding carboxylic acids (Perozich ). Since several aldehydes are toxic at low levels, this vast repertory of ALDHs present in BYs is likely to play a role in diverse reactions supporting extremotolerance. Three domains are related to zinc-containing alcohol dehydrogenase (Adh), which catalyse the oxidation of alcohols to their corresponding acetaldehyde or ketone. The IPR013154 and IPR013149 correspond, respectively, to the N-terminal and C-terminal portions of this enzyme and IPR011032 represents an oligomeric molecular chaperone associated with the N-terminal region involved in the folding protein process (Walter 2002). Adh are thought to be proteins prone to evolutionary changes following gene duplication due to their ability to assume new functions as consequence of their broad spectrum of substrates (Piskur et al., 2006, Conant and Wolfe, 2008). Furthermore, a physiological role of Adh has been reported in many biochemical pathways including stress tolerance, pathogenicity, detoxification, and substrate specificity (Piskur et al., 2006, Grahl et al., 2011). As expected for a fungal family with many members tolerating extreme conditions, the expansion of alcohol dehydrogenase domains in the black yeast from a common ancestor may have determined the diversification of these organisms in a range of ecological niches. Another important domain expansion verified in the common ancestor of all black yeasts analysed was the trichothecene efflux pump (IPR010573), which might have been important in black yeast to colonize sites contaminated with this class of compounds.

Analysis focussing on individual clades revealed that the bantiana- and carrionii-clades, which have pronounced trends in vertebrate infection, evolved in opposite directions (Fig. 5B). Several domains expanded in the bantiana-clade appeared to be reduced in the carrionii-clade. This would suggest that specific expansions in the bantiana-clade are attributed to ecological preferences in these organisms. However, the clades contain Fonsecaea pedrosoi and Cladophialophora carrionii, respectively, which cause the same disease, chromoblastomycosis, with the same invasive form, the muriform cell, and which thus do not share specific domains.

We did not observe expansions exclusive to the dermatitidis-clade. Previously domain expansions attributed to Exophiala dermatitidis (Chen ), such as IPR002656 and IPR020843, were also found expanded in members of the jeanselmei- and salmonis-clades. Unlike truly pathogenic fungi possessing a specialized thermosensitive tissue phase (Sharpton ), we did not observe massive functional domain loss compared to the ancestral black yeasts.

Cytochrome p450 expansion and diversification

Cytochrome p450 genes (CYPs) play a fundamental role in primary, secondary, and xenobiotic metabolism (van den Brink ). Due to their participation in a large number of detoxification reactions as well as in the metabolism of specific xenobiotics which may be co-assimilated as carbon source, CYPs are thought to be critical for the colonization of new ecological niches (Moktali ). In fungi, point mutation and overexpression of CYP family-specific genes have been found to be responsible for drug resistance (Lamb et al., 1997, Lupetti et al., 2002, Ma et al., 2006). The evolution of fungal pathogenesis is thought to be associated with CYP family expansion and diversification through gene duplication. Our CYP prediction analysis revealed an extraordinary p450 repertoire in black yeast-like fungi ranging from 231 predicted CYPs in Cladophialophora psammophila to 60 predicted CYPs in Capronia coronata (Table 1). Notably, Cl. psammophila was found in a hydrocarbon-polluted environment (Badali ), while Ca. coronata is a coloniser of decorticated wood in nature (Müller ). A comparison of the predicted number of CYPs to those of other species in the Fungal Cytochrome P450 Database (FCPD) (Park ) showed that some black yeasts are among the Ascomycota species with the highest number of CYPs (Fig. 6).

Table 1

Overview of Cytochrome p450 in black yeasts.

Clade	Species	Strain	# CYP	# Family	# Subfamily	# CYP not assigned
jeanselmei-clade	Exophiala xenobiotica	CBS 118157	164	62	39	41
	E.spinifera	CBS 89968	122	56	28	30
	E.oligosperma	CBS 725.88	131	52	30	30
	E.sideris	CBS 121828	97	40	23	25
dermatitidis-clade	E.dermatitidis	CBS 525.76	62	24	27	9
	Capronia epimyces	CBS 606.96	99	40	32	15
	C.coronata	CBS 617.96	60	25	19	16
Rhinocladiella mackenziei-clade-clade	Rhinocladiella mackenziei	CBS 650.93	161	56	46	44
carrionii-clade	Cladophialophora carrionii	CBS 160.54	101	37	31	29
	C.yegresii	CBS 114405	88	34	26	26
	C.semiimmersa	CBS 27337	109	44	28	36
bantiana-clade	Fonsecaea pedrosoi	CBS 271.37	164	70	38	38
	F.multimorphosa	CBS 102226	165	67	44	40
	C.immunda	CBS 834.96	144	51	38	37
	C.bantiana	CBS 173.52	175	68	42	48
	C.psammophila	CBS 110553	231	85	52	57
salmonis-clade	E.aquamarina	CBS 119918	179	68	36	51
	E.mesophila	CBS 402.95	75	39	19	19
Cyphellophoraceae	Phialophora europaea	CBS 101466	117	49	28	33
	P.attae	CBS 131958	135	59	32	37
	Coniosporium apollinis	CBS 100218	77	37	25	8
Outgroup	Verruconis gallopava	CBS 437.64	84	39	23	14

Fig. 6

Distribution of CYP p450 genes in Ascomycota. TOP 10 fungal species with highest CYP p450 numbers in ascomycetous genomes, based on search against the Fungal Cytochrome P450 Database (FCPD).

A total of 2740 CYP sequences were clustered in 131 families (Table S5) and 175 subfamilies according to their amino acid sequence identity against the Fungal p450 CYPs database (Nelson 2009). One hundred and nine partial CYP p450 sequences (BLASTP coverage >40 %) were classified as potential pseudogenes due to the occurrence of premature stop codons or presence of frameshifts (Table S6). These sequences are shorter than their homologues in other fungi. Potential pseudogenes were not included in downstream analysis. Comparative analyses revealed striking differences and expansions across the black yeast-like fungi in a range of CYP p450 families. We observed notorious CYP family expansions, mainly, but not exclusively, in species belonging to the bantiana-clade (CYP530, CYP682, CYP504, and CYP52) (Table S5). These CYP families potentially affect the metabolism of phenolic compounds and aromatic hydrocarbons (Olivera et al., 1994, Cox et al., 1996, Lin et al., 2011, Moktali et al., 2012, Zhang et al., 2012). Our findings are consistent, to some extent, with previous studies showing that some black yeasts appear to have adapted to grow in environments polluted with aromatic hydrocarbons (Woertz et al., 2001, Prenafeta-Boldú et al., 2002, Prenafeta-Boldú et al., 2006, Zhao et al., 2010a). Particularly important due to its abundance in some black yeast species, CYP530 is thought to participate in the degradation of several fatty acids and hydrocarbons (Moktali ). This CYP was found ranging from 13 copies in Cladophialophora psammophila and Fonsecaea pedrosoi to complete loss in Cl. yegresii (Table S7). The phylogenetic tree of CYP530 revealed multiple recent duplications and expansions. In addition, we observed two monophyletic clades likely correspond to distinctive subfamilies of CYP530 (Fig. S5). This gene redundancy observed might have been used to guard the above described critical functions as was shown in other fungi (Skamnioti ). To the best of our knowledge, the 13 copies is the highest rate of CYP530 reported in the fungal Kingdom (Table S7). Since Cladophialophora yegresii was only isolated from thorns of living cactus and was able to grow as an endophyte in cactus tissue (Zeppenfeldt , de Hoog ), it might be speculated that the absence of genes involved in secondary metabolism, such as CYP530, may implicate a biotrophic lifestyle where the organism obtains essential nutrients from its host.

At the subfamily level, we verified that the housekeeping genes CYP51F (encoding lanosterol 14α-demethylase) and CYP61A (encoding sterol delta22-desaturase), which are implicated in sterol biosynthesis (Yoshida and Aoyama, 1984, Podust et al., 2001, Lepesheva et al., 2008, Park et al., 2011) comprise one of the most conserved subfamilies across black yeast-like fungi. Azole antifungal agents interacting with CYP51, lead to growth inhibition and the death of fungal cells due to an ineffective conversion of lanosterol to ergosterol (Yoshida, 1988, Kelly et al., 1997). It has been demonstrated that additional copies of, as well as point mutations in, the CYP51 gene may lead to acquisition of resistance in fungi (Sanglard et al., 1998, Jones et al., 2014a). Our analyses revealed that most species have two CYP51F copies, whereas members of the dermatitidis-clade, Rhinocladiella mackenziei and the outgroups Coniosporium apollinis and Verruconis gallopava, have a unique CYP51F gene (Fig. S6). The important Y136F mutation (Mullins et al., 2011, Jones et al., 2014a) associated with CYP51 copy number variation and involved in azole resistance was not identified in these genes. This suggests that the tolerant allele, responsible for the azole resistance, is acquired only in the presence of azole fungicides. CYP61A was found in a single copy in all genomes studied.

Aromatic compound metabolism

Comparative analyses revealed that the genes PHA and HGD are organized in a syntenic cluster with at least six additional conserved genes (Fig. S7). We verified that this gene cluster organisation was retained by natural selection in most Herpotrichiellaceae. Besides the PHA and HGD genes, this cluster includes a variable number of genes coding for hypothetical proteins, an MFS transporter, a trehalose-6-phosphate hydrolase (T6P-hydrolase), and a fumarylacetoacetase (Fig. S7). T6P has been linked to diverse roles, such as energy source, protectant against stress of heat, freezing, starvation, dehydration, and desiccation (Wiemken, 1990, Iturriaga et al., 2009), and is important in fungal pathogenicity (van Dijck et al., 2002, Petzold et al., 2006, Ngamskulrungroj et al., 2009). The presence of PHA, HGD, and fumarylacetoacetase in this cluster overlaps the styrene degradation pathway, which might support the involvement of these genes in the degradation of aromatic compounds (Fig. S7). The MFS transporter may be involved in energy production transporting simple sugars across the mitochondrial membrane. As the synteny of these genes is highly conserved in several black yeast-like fungi, we hypothesize that the cluster configuration was probably acquired by their common ancestor, and subsequent gene rearrangement resulted in the current gene order and orientation in the extant species.

Secondary metabolism

Fungal secondary metabolites (SMs) are natural products important for the colonization of specific ecological niches. Despite their wide variation, all secondary metabolites are produced by a few common biosynthetic pathways and classified according to the enzyme involved in their biosynthesis: polyketides (PKS), non-ribosomal peptides (NRPS), terpenes and indole alkaloids (Keller ). We identified a large number of potential gene clusters for secondary metabolite present in black yeast (Table 2). The majority of these biosynthetic clusters correspond to PKS I/III (101 clusters), terpene (91 clusters) and NRPS (61 clusters), although it was verified that some species possess hybrid clusters (Table 2). In addition, the PKS III cluster was found only in Chaetothyriales since Coniosporium apollinis and Verruconis gallopava lack such gene cluster.

Table 2

Summary of secondary-metabolite gene classes in black yeast.

Species	Terpene	III PKS	I PKS	NRPS	Terpene/Indole/PKs I	NRPS + terpene	Phosphonate	Lantipeptide	I pks/terpene	I PKS/NRPS	NRPS + indole	Indole
Capronia coronata	4	1	2	4	0	0	0	0	0	0	0	0
Exophiala dermatitidis	4	1	2	5	0	0	0	0	0	0	0	0
C.epimyces	6	1	6	2	0	0	1	0	0	0	0	0
Cladophialophora psammophila	3	1	2	3	1	0	0	0	0	0	0	0
E.aquamarina	6	1	6	4	0	0	1	0	0	0	0	0
C.carrionii	5	1	4	3	0	0	0	0	0	0	0	0
Phialophora attae	4	1	3	1	0	0	0	0	0	0	0	0
P.europaea	4	0	2	2	0	0	0	0	0	0	0	0
C.semiimmersa	4	0	4	4	0	0	0	0	0	0	0	0
E.xenobiotica	4	1	6	2	0	0	1	0	0	0	0	0
E.oligosperma	4	1	3	4	0	0	0	0	0	0	0	0
E.spinifera	3	1	4	3	0	0	1	0	0	0	0	0
Verruconis gallopava	3	0	5	2	0	0	0	0	0	0	0	0
E.mesophila	4	1	1	2	0	0	0	0	0	0	0	0
E.sideris	4	2	2	2	0	0	2	1	0	0	0	0
Coniosporium apollinis	5	0	3	1	0	1	0	0	0	0	0	0
Fonsecaea pedrosoi	4	1	3	3	0	0	0	0	0	0	0	0
Rhinocladiella mackenziei	4	1	9	5	0	0	0	0	0	0	0	0
C.bantiana	3	1	2	2	1	0	0	0	0	0	0	0
F.multimorphosa	4	1	5	2	0	0	0	0	0	0	0	0
C.immunda	4	2	5	2	0	0	0	0	0	0	0	1
C.yegresii	5	1	2	3	0	0	0	0	0	0	0	0

Melanin synthesis

Fungi may produce melanin via distinct pathways: the eumelanin via the DHN and DOPA pathways, and the pyomelanins via l-tyrosine degradation pathway (Langfelder ). Recently, homologues of these three pathways have been identified in the pathogenic black yeast Exophiala dermatitidis and in other filamentous fungi (Youngchim et al., 2004, Chen et al., 2014). Similarly, we found that members of Herpotrichiellaceae possess several melanin-associated genes, suggesting they would be able to produce melanins using all different pathways, as was also suggested for Fonsecaea monophora (Li ). Unlike other filamentous fungi, where the melanin genes are frequently encoded in biosynthetic gene clusters (Kimura and Tsuge, 1993, Woo et al., 2010), we did not verify this organisation in black yeast-like fungi.

The dark polymer 1,8-dihydroxynaphthalene (DHN) melanin is produced via the DHN-melanin pathway and believed to be the best characterised fungal melanin biosynthetic pathway. Comparative analyses between previously released melanin-associated genes (Chen ) and our dataset revealed that many, but not all black yeasts possess homologues for the production of melanin by the DHN pathway (Table S8). The equally dark-pigmented outgroups Verruconis gallopava and Coniosporium apollinis outside of or basal to the Chaetothyriales showed the highest number of missing genes, including the known multicopper oxidases (MCOs) required for melanin biosynthesis. This suggests that the DHN-melanin pathway has been better conserved among the Herpotrichiellaceae and Cyphellophoraceae. However, MCOs with low similarity to known and well-characterised enzymes have been reported in fungi (Tamayo-Ramos ) and additional knowledge about their enzymatic properties is required to elucidate the DHN-melanin pathway in these species.

Similar to the DHN pathway, DOPA-melanin pathway homologues were identified across black yeast-like species. Of particular interest was the high number of tyrosinases and laccases found in Herpotrichiellaceae, but not in the outgroups Verruconis gallopava and Coniosporium apollinis (Table S8). The presence of multiple laccases only in Herpotrichiellaceae supports a diversification of this enzyme that occurred late in the evolution of black yeasts. A possible explanation for the presence of multiple laccase genes would be the various functions that have been attributed to this enzyme other than pigmentation, such as degradation of organic pollutants and lignin, and stress tolerance (Baldrian, 2006, Rodriguez Couto and Toca Herrera, 2006).

Protein degradation

The overall counts of the main MEROPS (Rawlings, 2010, Rawlings et al., 2014) peptidase families revealed the abundance of serine- (S) and metallo- (M) peptidase families in Chaetothyriales (Table S9). With the exception of the dermatitidis-clade, members of both Herpotrichiellaceae and Cyphellophoraceae presented specific and significant number of S09 (S09X sub-family), S33, and M38 families according to CAFÉ analysis (Fig. S8). S09X and S33 families appear to be significantly depleted in Eurotiales and Onygenales, while these proteins might play an important role in the ecology of Chaetothyriales (Muszewska ). Cluster analysis of sequences from the S09X family revealed that most BY protein expansions were found in two different clusters (Fig. S8). Protein sequence classification showed that S09X corresponds to alpha/beta hydrolase fold-3 (IPR013094/PF07859) and proteins containing a carboxylesterase type B (IPR002018/PF00135) domain (Fig. S9). Gene tree reconstruction showed that main gene duplication events were at the basis of the M38 IPR002018/PF00135 domain expansion in BY, while several losses in Eurotiales and Onygenales explain the relative accumulation of IPR013094/PF07859 proteins in BY (Fig. S9). According to several authors, the expansion of metalloproteases M35 and M36 could be associated with mammal-host association (Sharpton et al., 2009, Martinez et al., 2012, Whiston and Taylor, 2015). According to MEROPS classification, the M35 and M36 protein families are depleted in Cyphellophoraceae and absent among Herpotrichiellaceae (Fig. S9). On the other hand, we detected an expansion of M38 proteins in BY, which may be associated with β-aspartyl dipeptidase acting in the release of iso-aspartate residues from peptides as characterised for bacteria (Borek ). Cluster analysis of the M38 family revealed that most of the BY protein expansions were found to be enriched in the highlighted three clusters (Fig. S10). Protein sequence classification of those three clusters revealed that M38 BY enrichment corresponds to an amidohydrolase/metal-dependent hydrolase (IPR011059, IPR006680/PF01979) domain containing proteins. In cluster III, beyond these domains, we detected the presence of a tryptophan synthase (IPR001926) domain expanded in Herpotrichiellaceae and Cyphellophoraceae (Fig. S10).

Carbohydrate-active enzymes

Carbohydrate-active enzymes (CAZymes) are responsible for the degradation, modification, and biosynthesis of carbohydrates and glycoconjugates (Cantarel ). The family classification system based on amino-acid sequence and structure similarities has been used to group the CAZymes into five classes of enzyme activities and one associated module: glycoside hydrolases (GHs), glycosyltransferases (GTs), polysaccharide lyases (PLs), carbohydrate esterases (CEs), auxiliary activities (AAs), and the associated module carbohydrate-binding modules (CBMs) (Cantarel ). In this study we determined the CAZymes composition distributed across the black yeasts and compared this with other ascomycetes. In total, 154 CAZymes families were identified in the predicted protein sets. Generally, the highest and the lowest number of CAZymes were found in members of the jeanselmei- and dermatitidis-clades, respectively, although the variation observed between species within clades was considered low. Some CAZymes appeared to be clade-specific. For example, the GH62 family was only found in the europaea-clade (Cyphellophoraceae). On the other hand, several CAZymes were identified in all species examined (Table S10).

Some striking depletions were verified in CAZyme families involved in the degradation of plant material. Plant cell wall polysaccharides are subdivided into three categories: cellulose, hemicellulose (including xylan, xyloglucan, glucogalactomannan, galactan, and respective side chains), and pectin (composed of galacturonan, rhamnogalacturonan, and respective side chains) (Amselem ). Most black yeast-like fungi lack the pectinases PL1, PL3, PL4, PL7, PL9, and PL10 (Table S10). Comparable depletions have been reported in species of Onygenales (Desjardins ), and Sporothrix (Teixeira ) while they are present in Eurotiales. The β-1,4-glycosyl hydrolase family 28 (GH28) is another family linked to the breakdown of pectins. This enzyme is absent in the dermatitidis-clade, jeanselmei-clade, and salmonis-clade, but present in the bantiana- and europaea-clades. Similarly, pectin methylesterase family CE8 and pectin acetylesterase family CE12 are absent in Herpotrichiellaceae. Comparable patterns are found in the xylan-associated enzyme family GH11 (endo-β-1,4-xylanase), present only in the carrionii- and europaea clades, as well as in Eurotiales, and CE 1 (acetyl xylan esterase) missing in Onygenales and in all black yeasts examined.

Depletions were also verified in chitin-related enzymes, a critical component of the fungal cell wall (Latgé 2007). Black yeasts on average have 5 members of the chitinase family GH18 per species. In contrast, Onygenales and Eurotiales have 10 and 21 members per species on average, respectively (Table S10). Moreover, the Carbohydrate-Binding Module Family 18, which is often found attached to a number of chitinase catalytic domains, is depleted in black yeasts. These comparisons suggest that the breakdown of chitin is likely reduced compared to other filamentous fungi.

The family AA4 contains vanillyl-alcohol oxidase (VAO), which is missing in several other ascomycete fungi and is well represented in BY. VAOs catalyse the oxidation of a wide range of phenolic compounds and are abundant in black yeast genomes ranging from 10 copies in Cladophialophora psammophila to two copies in Capronia coronata. This finding is consistent with the ability of many black yeasts to degrade aromatic compounds (Isola ).

Cell-wall biosynthesis

The cell wall is an essential structure involved in protective functions against osmotic pressure and environmental stress (Bowman & Free 2006). The three major fungal cell wall constituents are chitin, mannan, and β-glucan. These components have been implicated in fungal virulence and represent targets for immune surveillance mechanisms (Bulawa ). In agreement with previously published data (Chen ), BY genomes encode an arsenal of genes involved in chitin synthesis. All 7 proposed classes of chitin synthase genes (CHS) previously described in fungi (Roncero 2002) were present in BY, except Class VI, which is missing in Rhinocladiella mackenziei. This species is recognised as an important causative agent of cerebral phaeohyphomycosis (Li & de Hoog 2009); mutants in CHS-VI are viable and less virulent (Bulawa ). Proteins linked to the regulation and exportation of chitin synthase are conserved in BYs (Table S11). In contrast, comparative analysis of chitin degradation genes showed that black yeasts lack chitosanase, which is conserved in Saccharomyces cerevisiae and Schizosaccharomyces pombe (Table S11). Additionally, BYs have fewer chitinase proteins belonging to the family GH18 compared to other filamentous fungi, as described above (Table S11). Chitin deacetylases, which are believed to be secreted exclusively during modification of chitin in the cell wall (Zhao ), are missing in the carrionii-clade and in R. mackenziei.

Investigation of the genes related to synthesis and processing of 1,3-α-glucan revealed they are altered significantly in the species analysed. The Ags family of 1,3-α-glucan synthase is absent only in Exophiala dermatitidis, but the Agn family of 1,3-α-glucanases is absent from the dermatitidis-, jeanselmei-, and salmonis-clades, and in Rhinocladiella mackenziei, even though these families are present with multiple copies in Aspergillus. Furthermore, BYs possess a putative α-amylase believed to be involved in the formation and/or modification of α-glucans (Table S11).

Transporters

Black yeasts like other filamentous Ascomycota possess a large proportion of genes associated with transporter activity. Our InterProScan analysis revealed that the most abundant protein domain verified in several BYs contains several families of transporters, particularly the Major Facilitator Superfamily (MFS). To better understand the transporter mechanisms in BYs, we annotated transporter subfamilies across all 21 species based on their best match to the curated transporter database TCDB (Saier ). Overall, black yeasts possess more MFS transporters than species of Onygenales and Eurotiales.

The most abundant transporter subfamily found in BYs is a potential nicotinate permease. It has 27 candidate genes in the outgroups Verruconis gallopava and Coniosporium apollinis, but up to 93 candidate genes in Exophiala aquamarina. This transporter belongs to the family of the Anion:Cation Symporter (ACS) (TC 2.A.1.14.11) of the major facilitator superfamily. Another MFS subfamily with a remarkably high number of predicted members is the trichothecene efflux pump (TC 2.A.1.3.47) of the Sugar Porter (SP) family. Since toxin efflux pumps are responsible for mediating both intrinsic and acquired resistance to toxic compounds, this result provides genomic insight into the known extremotolerance of black yeast-like fungi. Moreover, this finding is consistent with the expansion of the trichothecene efflux pump protein domain (IPR010573), as described above.

At family level, the Sugar Porter (SP) Family, the Anion:Cation Symporter (ACS) Family, and the Drug:H+ Antiporter-1 (12 Spanner) (DHA1) Family are the most abundant in BYs. Interestingly, among the family DHA1 we verified that the subfamily 2.A.1.2.77, which confers phenylacetate resistance, is well-represented in the majority of the species examined. Other families verified in all BYs analysed comprise the Ferroportin (Fpn) Family (TC 2.A.100), the Proton-dependent Oligopeptide Transporter (POT/PTR) Family, and the Equilibrative Nucleoside Transporter (ENT) Family.

Sexual and parasexual reproduction

Fungi exhibit a wide diversity of reproductive modes, including sexual, asexual, and parasexual cycles. Recombination is an important and needed process in any fungal life-cycle, and may alter virulence traits, increase fitness in new ecological niches, and eliminate deleterious mutations (Heitman, 2006, Lee et al., 2010, Ni et al., 2011). We used models of sexual and parasexual cycles of Aspergillus nidulans and Neurospora for BY comparisons (Glass et al., 1990, Paoletti et al., 2007, Debuchy et al., 2010, Zhao et al., 2015). We first identified the mating-type idiomorph within each assembled genome. Homothallism of Capronia coronata and Ca. epimyces was confirmed by identifying both MAT1-1 and MAT1-2 Aspergillus homologues closely clustered in a single assembled scaffold (Fig. 1). With the exception of outgroup species Verruconis gallopava (Venturiales), the remaining analysed 20 genomes of asexual species harboured a single mating type idiomorph (either MAT1-1 or MAT1-2) within each assembly, confirming that these fungi are heterothallic (Fig. 1).

We analysed the MAT locus organisation within the main groups of Herpotrichiellaceae and Cyphellophoraceae using genomic information from the MAT flanking genes. Among Eurotiomycetes, the flanking genomic regions of the MAT locus harbours APN2, SLA2, APC5, and COX13 genes, which are conserved and organized in synteny (Coppin et al., 1997, Fraser et al., 2007, Paoletti et al., 2007). We first aligned and compared the gene models from both dothideaceous species V. gallopava and Coniosporium apollinis. The APN2, COX13, APC5, and CIA30 genes appear to be conserved in synteny and preserved in the right MAT flanking region (Fig. 7). However, the SLA2 gene was not found in the MAT locus in these species, but is not genomic linked. Coniosporium apollinis is inferred to be a heterothallic species since it harbours a single copy of the MAT1-1 gene in its genome. In the left flanking site of the MAT locus we found a homologous protein that was syntenically conserved between the two dothideomycete species (PV09_01802/W97_06799); in the genomic alignments it is adjacent to the MAT1-1 genes of Coniosporium apollinis W97_06800 and W97_06801 (Fig. 7). Thus, the MAT locus of V. gallopava harbours two different ORFs: PV09_01800 and PV09_01801, which are not present in the MAT locus of the Coniosporium apollinis genome. According to the protein classifications and annotation, the ORF PV09_01800 is unique to V. gallopava and no putative domains were found. On the other hand, the ORF PV09_01801 encodes a homeodomain-like (HD) protein that carries a DNA-binding homeodomain motif (Fig. 7). This domain is found within the mating types 1 and 2 genes (MAT/MTLα2, Pi and MAT/MTLa1) in yeasts of Saccharomycotina and Taphrinomycotina in Ascomycota, as well in Basidiomycota (Martin ). According to Lee , this domain was lost during speciation of Pezizomycotina, but our analysis of additional species revealed that the HD domain was recognized as a potential mating regulator in Venturiales (Fig. 7). On the other hand, we confirmed the lack of the HD in Eurotiomycetes (including BYs) once an α-box and HMG were found in the MAT locus.

Fig. 7

Mating type locus structure of ancestral lineages C. apollinis and V. gallopava (top-right panel), the Cyphellophoraceae (bottom-left panel) and the salmonis-clade of Herpotrichiellaceae (top-left panel). Sexual loci for each fungal species are displayed in each respectively scaffold and the corresponding genes and accession numbers are displayed to each gene.

We detected the MAT1-1 (α-box) and MAT1-1-5 genes within the mating type 1 locus and/or the MAT1-2 (HMG) gene in Chaetothyriales (Fig. 1). The function of MAT1-1-5 in mating is not well established, and appears poorly conserved with MAT1-1-4 gene among Onygenales (Mandel et al., 2007, Burmester et al., 2011). As reported previously for some ascomycete species (Yun et al., 1999, Tsui et al., 2013), we obtained indirect evidence of a truncated version of the MAT1-1 gene within the MAT1-2 idiomorph, potentially driven by unequal recombination at the MAT locus in an ancestor of Chaetothyriales (Fig. 8, Fig. 9, Fig. 10, Fig. 11, Fig. 12). The loss of a functional α-box domain suggests that the truncated MAT gene might have diverged under low selective pressure after unequal recombination, or was silenced due to interference if both HMG and α-box domains were present. The COX13 gene appears not to be conserved among Chaetothyriales in the flanking regions of the MAT locus as usually observed in Eurotiomycetes (Coppin et al., 1997, Debuchy et al., 2010, Lee et al., 2010).

Fig. 8

Mating type locus structure of heterothallic species R. mackenziei (lower panel), E. dermatitidis and closely related homothallic Ca. coronata and Ca. epimyces (dermatitidis-clade upper panel). Mating type genes are represented in each corresponding assembled scaffold. Accession numbers are displayed to each gene.

Fig. 9

Mating type locus organization of heterothallic species from jeanselmei-clade, Herpotrichiellaceae. Mating type genes are represented in each corresponding assembled scaffold. Accession numbers are displayed to each gene.

Fig. 10

Mating type locus structure of heterothallic species from carrionii-clade, Herpotrichiellaceae. Mating type genes are represented in each corresponding assembled scaffold. Accession numbers are displayed to each gene.

Fig. 11

Mating type locus organization of heterothallic species from bantiana-clade, Herpotrichiellaceae. Mating type genes are represented in each corresponding assembled scaffold. Accession numbers are displayed to each gene.

Fig. 12

Distribution of heterokaryon (het) containing genes in 23 black yeast-like fungi and related Ascomycota. (A) Total counts of het containing genes (IPR010730 domain) for each species. (B) Pairwise similarity graphs generated by clustering analysis of 1439 het-containing proteins from Chaetothyriales and related species.

The Cyphellophoraceae species Cyphellophora europaea and Phialophora attae presented a rather conserved MAT locus structure compared to other Eurotiomycetes. The right MAT flanking domain harbouring the genes SLA2, APC5, and SAICAR5 appeared to be conserved, and at the opposite side in the left flanking area of both species the APN2 and other hypothetical proteins were organized in synteny (Fig. 7). Gene content within the MAT locus diverges in the Cyphellophoraceae: Cyphellophora europaea has MAT1-1 and MAT1-5 configuration, while Phialophora attae harbours the MAT1-2 gene.

The structure of the MAT locus of Herpotrichiellaceae deviates from that of most other members of Eurotiomycetes. We observed an expansion or collapse of the canonical MAT structure compared to model species in, for example, Aspergillus. The flanking site of the MAT genes of some BY was inflated with the accumulation of novel genes or was even unrelated to the MAT locus in other ascomycetes, which suggests a low selective pressure in this important genomic region within the family (Fig. 7). Exophiala aquamarina in the salmonis-clade had a heterothallic MAT locus structure with the MAT1-1 gene, as well as flanking genes SLA2, VPS13, and APN2 conserved in synteny with other Eurotiomycetes. The heterothallic MAT locus structure of E. mesophila lacked this structure. The right flanking area of E. mesophila showed homology and structural conservation with SLA2 and VPS13 genes of E. aquamarina, but lacked synteny in the left flanking region of the MAT locus (Fig. 7). No homology at the left flank of the MAT locus was detected between the two species. In addition, the APN2 gene was located in another scaffold of E. mesophila, unrelated to the MAT locus. Within the dermatitidis-clade we detected an expansion of the MAT locus, which followed a speciation process of the three members of this clade. Exophiala dermatitidis is placed as the basal taxon of the dermatitidis-clade and presents a well-conserved, heterothallic MAT structure with other Eurotiomycetes (Fig. 8). On the other hand, we detected a chromosomal expansion at the right flanking site of both Capronia homothallic MAT loci, which was followed by gene inflation at this locus. This locus has some peculiar features. First, we identified a novel MAT gene that is found within the MAT locus only in those three species (HMPREF1120_08861/A1O3_06090/A1O1_07968). Second, within the canonical SLA2-APN2 MAT locus structure, unique genes were detected within each of the three species, with the highest frequency in Ca. epimyces, since this has a larger SLA2-APN2 genomic range. We also detected an expansion of the MAT locus that is followed by a speciation process with acquisition and inflation of genes in the Exophiala species of the jeanselmei-clade, E. xenobiotica, E. spinifera, and E. oligosperma (Fig. 9). Frequent appearance of new and family-specific genes is observed throughout the Herpotrichiellaceae along the mating type genes, which might be a source of adaptive novelties of mating regulators.

Exophiala sideris is the most basal species in the jeanselmei-clade and its MAT locus structure followed the classical SLA2-APN2 configuration. However, we detected a fused MAT1-1/MAT1-2 gene configuration in this specie (Fig. S11). Protein classification analysis revealed that both α-box and HMG domains are present within a single mating regulator gene, leading us to hypothesize this as an unusual gene fusion event potentially giving a homothallic status to this species. The protein was blasted against the Conserved Domain Database (CDD) (Marchler-Bauer ) in order to achieve MAT gene configuration common to fungi, as found e.g. in the homothallic ascomycetes Curvularia homomorpha, Bipolaris luttrellii, and Penicillium rubens (Fig. S11). The disposition of both α-box and HMG domains varies across the species panel analysed, either being separated along the gene or fused, where mostly HMG binding sites are found within the α-box domain (Fig. S11). The latter case led us to speculate that HMG insertions could be an atavism from an ancient homothallic state, since the majority of the gene sequence is related to MAT1-1, or it could be a product of gene fusion and unequal crossing over of two opposite mating type strains. This last scenario confirms earlier reports where cryptic homothallism was proven to occur in Curvularia homomorpha and Bipolaris luttrellii (Yun ). Possibly the fused MAT1-1/MAT1-2 also plays a role in cryptic homothallism of our species under study. Homothallism as an ancestral state was demonstrated in the carrionii-clade, which differed from what was observed in the jeanselmei- and dermatitidis-clades (Fig. 10): Capronia semiimmersa harbours both MAT1-1 and MAT1-2 genes in a single haploid genome. In addition, these species exhibits an apparent expansion of the MAT locus in that the classical SLA2-APN2 configuration was not found, while new genes were detected down/upstream the MAT genes (Fig. 10). On the other hand, the heterothallic species Cl. carrionii and Cl. yegresii displayed the SLA2-APN2 structure, and new carrionii-clade-specific genes were detected along with the MAT genes. In addition, we detected a Cl. carrionii-specific gene acquisition within the MAT locus as represented by the yellow boxes in Fig. 10.

The most degenerated MAT locus structure within the Herpotrichiellaceae was found within the bantiana-clade. Cladophialophora bantiana (MAT1-1) and Cl. psammophila (MAT1-2) shared a large degree of synteny at the left side of the MAT genes extending to the SLA2 gene (Fig. 11). We also detected unique and shared genes between this specific Herpotrichiellaceae clade, represented by yellow boxes in Fig. 11. At the right flank of the MAT genes, APN2 are poorly conserved between these two species. The APN2 gene is assembled in a scaffold different from that of the MAT genes. The remaining species of the bantiana-clade, Fonsecaea pedrosoi, F. multimorphosa, and Cladophialophora immunda did not share any synteny at the flanking regions of the MAT genes.

Overall, we detected a low selective pressure within the MAT locus structure of Chaetothyriales, compared to other Eurotiomycetes (Coppin et al., 1997, Debuchy et al., 2010, Lee et al., 2010, Burmester et al., 2011). The species-specific, non-characterised genes and gene duplications near the MAT genes are unique to Chaetothyriales. Despite the presence of components of the mating/pheromone-response pathway and their respective domains among Chaetothyriales, we hypothesize that due to the low selection of the MAT locus, those domains could represent an alternative system for generating diversity via parasexuality. The parasexual cycle has not been explored or characterised in any species of Chaetothyriales. Parasexuality is a process triggered by cell-cell fusion and ploidy reduction through random chromosome loss and has been reported in various filamentous fungi including Aspergillus nidulans, Neurospora crassa, and Podospora anserina (Pontecorvo, 1956, De Serres, 1962, Labarere and Bernet, 1977). Undifferentiated vegetative cells undergo haploid cell fusion and produce heterokaryons, as with aberrant ploidy due to mitotic crossing-over (Tolmsoff 1983). Cell compatibility is dependent on a combination of loci known as heterokaryon (het) incompatibility loci (allorecognition loci) (Saupe and Glass, 1997, Zhao et al., 2015). When incompatible, the vegetative cells undergo genetically programmed cell death (Glass ). There is strong evidence that recombination can be triggered by a parasexual cycle in some of the consistently asexual BYs: (1) Based on gene family evolution analysis of BY, we identified a dramatic increase of het-containing proteins in most BY genomes. (2) The expansion of het proteins in BYs exceeds the number of allorecognition loci in N. crassa or Aspergillus species (Table S12). (3) Random gene duplication across BY genomes, which could be linked to impaired mitotic chromosomal reduction via aneuploidy.

With few exceptions, het-containing proteins are expanded in most BY genomes (Table S12). While species of Eurotiales harbour 2–12 het-containing proteins, in the Chaetothyriales the numbers range from a single copy in E. dermatitidis and Cyphellophora europaea to 134 copies in Cl. psammophila. Clustering analysis via CLANS of het-containing proteins (PF06985/IPR010730 – Pfam/InterPro) of 23 black yeasts, Onygenales, Eurotiales, and N. crassa did not show any significant sub-cluster within PF06985/IPR010730 containing genes. Alignment of the 1439 het-containing proteins in order to access the phylogenetic distribution did not return well-conserved alignment blocks due to high sequence dissimilarity. In the attraction graphs of the PF06985/IPR010730 family (Fig. 12), four clusters within het-containing proteins were visually defined in order to narrow down alignment discrepancies. We extracted the sequences related to each of the four identified groups and protein alignment was performed for phylogenetic analysis. All predicted het containing proteins of N. crassa (Zhao ) were included, among which were three het loci, viz. tol, het-6, and pin-c. Phylogenetic analysis shows that most het-expanded BY harbour orthologues of well-characterised Neurospora tol, het-6, and pin-c (Fig. S12). We also identified a cluster of proteins that were related to the HET-E/D/R loci in Podospora anserina (Fig. S12). Beyond the het domain, this particular group of proteins display a GTP binding site followed by WD40 repeats, which play an important role in specificity of vegetative incompatibility in the P. anserina parasexual cycle (Espagne ). BY genomes display a vast repertoire of heterokaryon incompatibility proteins and mechanisms of sexual or parasexual recombination, which needs further investigation.

Discussion

Black yeasts and similar fungi in the order Chaetothyriales are known for their morphological plasticity, asexual diversity, and divergent habitat choice. To date (01-01-2016), 23 species of the order have had their genomes sequenced (Fig. 1). The order contains four or five families according to Réblová who upgraded the “europaea-clade” to family level (Fig. 1). Fig. 2 shows a phylogenetic tree of all available species based on LSU-sequences, which is as yet the only parameter alignable over the entire order. Two ecological, highly speciose groups affiliated to the Chaetothyriales have recently been detected but as yet have not been formally named: asexual species with ant-associated ecology and those parasitizing on lichens (Muggia et al., 2013, Muggia et al., 2016). In addition, many described Capronia species (Barr, 1972, Aptroot et al., 1997, Etayo et al., 2013) have not been genotyped. Some of the lichen-pathogens produce well-recognisable sexual morphs with single, small setose cleistothecia with pale to brown, septate ascospores resembling Capronia or Trichomerium, underlining that the sexual morphs in the order, quite in contrast to the asexual morphs, show little variation. Judging from the above, the Chaetothyriales in the phylogenetic tree of Fig. 2 show severe defects in taxon sampling. Therefore it is difficult to reconstruct a reliable phylogeny with ancestral and derived families in the correct position. Phylogenomic data are thus far available for Herpotrichiellaceae and Cyphellophoraceae only.

Species of Chaetothyriales can be found on phanerogam leaf litter and plant debris, but in contrast to common litter fungi such as Alternaria and Cladosporium, decomposed, tannin-rich material is mostly preferred. Species selectively grow with the presence of aromates and etheric oils, e.g., on babassu coconut shells, which contain a remarkable diversity of black yeasts (Nascimento ). Many species even seem to prefer artificial, human-made habitats, such as gasoline tanks (Isola ) and toxic mine waste rich in heavy metals (Seyedmousavi ). Several studies have indicated that herpotrichiellaceous black yeasts and their filamentous counterparts are potent degraders of toxic monoaromates and are frequently found inhabiting industrial bio-filters (Middelhoven et al., 1989, Middelhoven, 1993, Cox et al., 1997, Prenafeta-Boldú et al., 2006). This property is observed in members of nearly all clades (Woertz et al., 2001, de Hoog et al., 2006, Badali et al., 2011, Seyedmousavi et al., 2011). Phenolic and indolic compounds are substrate units for melanin formation, and the production of this pigment ultimately results from their oxidative polymerization (Jacobson, 2000, Plonka and Grabacka, 2006, Vavricka et al., 2010). Numerous additional studies have reported the presence of clinical and non-clinical species in oil-related environments (Phillips et al., 1998, Prenafeta-Boldú et al., 2001, Prenafeta-Boldú et al., 2006, Sterflinger and Prillinger, 2001), and noted preference of creosoted wood over untreated wood (Seyedmousavi et al., 2011, Döğen et al., 2013a, Döğen et al., 2013b, Gümral et al., 2014). The fungi become nearly the sole colonizers when the material contains toxic hydrocarbons, e.g., on creosoted telephone poles and railway sleepers (Gümral ).

Several possible routes of aromatic hydrocarbon metabolism have become known in black yeast-like fungi (Cox ). Among these the degradation of benzene derivatives via phenylacetate and homogentisate seems to be one of the most important pathways in the family Herpotrichiellaceae (Gunsch ). Overexpression of the genes 2-hydroxy phenylacetate (PHA) and homogentisate 1,2-dioxygenase (HGD) when the fungus is grown in the presence of ethylbenzene, supports the involvement of these enzymes in organic compound degradation (Gunsch ). It has also been demonstrated that genes of catabolic pathways, which may enhance survival under different environmental conditions, are physically clustered preserving gene order and orientation (Keller & Hohn 1997). This organization might favour the co-inheritance and the co-expression of multiple enzymes, which handle toxic intermediate compounds (Takos and Rook, 2012, McGary et al., 2013) and suggests an essential role for this ability in Herpotrichiellaceae. Regarding the phenolic compound catabolism, homologues of the styrene pathway were found present in the species studied, except in Capronia coronata: this species lacks the genes coding to 4-hydroxyphenylpyruvate dioxygenase and maleylacetoacetate isomerase. This suggests that across the analysed black yeasts, only Capronia coronata would not be able to synthetize pyomelanin via the accumulation of homogentisate.

Remarkably, the same species described above can also be found in clear water environments that are very poor in nutrients. For example, Exophiala dermatitidis on the one hand occurs on creosoted wood (Döğen ) and in gasoline (Isola ), but it is also common in steambath facilities (Matos ), hot springs (Sudhadham ) and dishwashers (Zalar ). This strongly suggests that the natural competitive abilities of these fungi are very low, so that they evade confrontation with other mircoorganisms and escape in hostile environments (de Hoog 1993). Gostincar classified E. dermatitidis as a polyextremotolerant fungus.

The above description holds true for the derived family of Chaetothyriales, the Herpotrichiellaceae. The origin of this ecology might be found in the life style of members of ancestral families, Chaetothyriaceae, Epibryaceae, and Trichomeriaceae. In the general phylogeny of Fig. 1, numerous members are epiphytic on plants, occurring as sooty moulds (Chomnunti ) or as rock-inhabiting fungi (Isola ). These habitats require similar abilities as described above to cope with conditions that suppress most competitors, e.g., lack of nutrients, dryness, or extreme and changing temperatures (Middelhoven, 1993, Vicente et al., 2008, Zhao et al., 2010a). A number of species have been isolated from green plants and have been regarded as host-specific plant pathogens, e.g., Cladophialophora hostae (Crous ), Exophiala eucalyptorum (Crous ), or Metulocladosporiella musae (Crous ). CAZyme families, involved in the degradation of plant material, do not seem to play a major role in the black yeast-like fungi, as most of them lack the pectinases PL1, PL3, PL4, PL7, PL9, and PL10 (Table S10). Comparable depletions have been reported in species of Onygenales, which contain numerous obligate pathogens of humans and other vertebrates, while they are present in Eurotiales, with e.g. Aspergillus fumigatus essentially being a degrader of plant debris. Only the Epibryaceae, which have recently been associated with Chaetothyriales (Gueidan ) after they had been regarded as members of Dothideomycetidae for decades (Stenroos ), are known biotrophs. However, ancestral clades in Chaetothyriales including Epibryaceae show long branches (Fig. 2) and their affiliation with the order may be due to incomplete taxon sampling.

Voglmayr isolated a large number of undescribed black fungi from ant domatia and cartons which appeared to be closely related to the Chaetothyriales. Ants and fungi have close symbiont relationships. The ecological roles of chaetothyrialean symbionts appear to be different between domatia and carton associations (Voglmayr ). In the domatia they play important roles such as in recycling ant waste, and there is also evidence that the ants feed on the fungi, contributing to rapid recycling of nitrogen in the tripartite symbiosis of ant, plant and fungus (Little and Currie, 2008, Defossez et al., 2011, Blatrix et al., 2012). In the carton association, ant tunnels or nests are largely made of black fungal material. It has been hypothesized that the fungus serves as building material (Lauth ) and that the fungal layer on the carton walls could act as mechanical protection against radiation, humidity and microbial decomposition (Mayer and Voglmayr, 2009, Zakharova et al., 2013), enhancing the durability. Commonly several black species co-occur on the same carton, forming complex associations, which rely on constant maintenance and care by the ants. Most carton species lack conidiation. In domatia a rather specific association with the host may be observed, the fungi producing a dense mat inside the domatium. These species are hyaline to light brown and show conidiation. Ants possess a great arsenal of exocrine glands (Moglich et al., 1974, Attygalle et al., 1989, Poulsen et al., 2002, Fernandez-Marin et al., 2006) secreting organic compounds that are effective against fungi and bacteria (Schlüns & Crozier 2009). Some species of leaf-cutting ants exude antimicrobial flavonoid and tannin-like compounds and communicate by aliphatic hydrocarbons (Brandt ). Chaetothyriales have been shown to use toxic compounds such as aromatic hydrocarbons as unique nutritional carbon sources (Prenafeta-Boldú et al., 2006, Zhao et al., 2010a), which thus act as key selective agents promoting the dominance of Chaetothyriales in ant nests. The evolutionary origin of Formicinae dates back to Cretaceous times, around 100 MYA. Given the large extant diversity of Formicinae, as well as of associated Chaetothyriales, it seems possible that ants and black yeast-like fungi have diversified in concert. We estimated that the common ancestor of Herpotrichiellaceae–Cyphellophoraceae emerged around 75–50 MYA, shortly after the Cretaceous-Paleogene (K-Pg) extinction event. In a similar study, Gueidan calculated the ancestral groups of Chaetothyriales with a rock-inhabiting life-style and lichenised Verrucariales around 250 MYA. Verrucariales might have been the first to colonise barren rock after the meteor impact that marked the transition from Perm to Trias. Early Chaetothyriales were thought to be hyperparasites on lichens (Gostincar ) lacking an algal component. Lichens produce large amounts of toxic metabolites, and thus the pathogens must have been able to cope with harsh climatic conditions as well as with life in toxic habitats. This may have been the period where early Chaetothyriales developed stress tolerance, nutritional oligotrophism, and physiological versatility to survive the wide array of toxic secondary metabolites produced by the lichens (Gostincar ).

With these ecological and evolutionary speculations, two factors are of particular interest: melanin and action of the cytochrome p450 enzymes. The presence of eumelanin produced via the DHN and DOPA pathways, and pyomelanins via l-tyrosine degradation pathway (Alviano et al., 1991, Sun et al., 2011, Eisenman and Casadevall, 2012, Li et al., 2016) are fundamental to the obligatory melanisation of the cell wall to enhance stress tolerance. Melanised fungi are resistant to environmental challenges found particularly in extreme habitats, including irradiation, nutrient depletions, and high temperature (Rosas & Casadevall 1997); this matches well with the conditions of a rock-inhabiting lifestyle. The melanised fungi are even able to survive in radioactive environments. Fungi growing on surfaces with direct sunlight exposure are highly adapted to cope with ionizing radiation via the constitutive presence of melanin, which acts as energy transduction molecules (Dadachova & Casadevall 2008). Melanised Exophiala dermatitidis cells exposed to ionizing radiation grow faster than non-exposed cells, suggesting a critical role of irradiated melanin and its conversion as a source of energy in herpotrichiellaceous fungi (Dadachova ). Their energy transduction is multifunctional, also playing an essential role in resistance to oxidative burst of vertebrate phagocytes (Henson et al., 1999, Jacobson, 2000, Eisenman and Casadevall, 2012).

Melanins can also act as scavengers of free and oxidative radicals, are cross-linked to fungal cell wall carbohydrates, and also interact with surrounding molecules. These black pigments may therefore have a protective function against natural predators, such as amoebae (Nosanchuk and Casadevall, 2003, de Almeida-Paes et al., 2012). In the most derived family, the Herpotrichiellaceae, the amoebic counterpart may have become the mammalian host's innate immune phagocyte, conferring protection against free radicals of host immune cell oxidative burst (Little and Currie, 2008, Defossez et al., 2011, Blatrix et al., 2012). Melanins thus have become a potential virulence factor. Indeed, we witness the largest number of systemic opportunistic species on humans and cold-blooded animals in the Herpotrichiellaceae, while this behaviour is nearly absent from other families, with Arthrocladium fulminans in the Trichomeriaceae as the only exception (Nascimento ).

Very little is known about the cell wall organization in BY and its dynamic composition. Loss of cell wall alpha-glucan synthase appears to be specific to Exophiala dermatitidis and the two outgroups Verruconis gallopava and Coniosporium apollinis. Other enzymes involved in 1,3-alpha-glucan processing are missing from remaining black yeast species studied. In contrast, the presence of the α-amylases, believed to be involved in the formation and/or modification of α-glucans (van der Kaaij ), suggests that BY carry a cell wall that deviates from that of other filamentous fungi. Importantly, some recent studies have shown that the chronicity observed in chromoblastomycosis is due to a failure of pattern recognition receptor (PRR) costimulation (van der Kaaij ). For example, innate recognition of F. pedrosoi is mediated by C-type lectin receptors (CLRs), but not by Toll-like receptors (TLRs), triggering an inadequate protective inflammatory response and leading to chronic infection (van der Kaaij ). TLRs can recognize pathogen-associated molecular patterns (PAMPs), such as α-glucans in fungi (Takeda ). We therefore speculate that BY genetic variation associated with cell wall composition, the main source of PAMPS, might affect the pattern of recognition of invading pathogens by the host. Indeed, enzymatic treatment to modify α-glucans from the conidial surface has been reported to decrease the phagocytic index in other ascomycetes (Bittencourt ). Further studies will provide a better understanding of the role of α-glucans in the induction of innate immune response as well as its structural organization in black yeasts. No specific virulence differences were found between closely related, pathogenic versus environmental siblings, e.g. between Cladophialophora bantiana/Cl. psammophila, and between Cladophialophora carrionii/Cl. yegresii. It may be concluded that pathogenicity of black fungi is primarily of opportunistic nature, enhanced by the combination of extremotolerance and assimilative abilities of aromatic compounds, similar to pathogenicity in Cryptococcus.

Genetic diversity plays an important role in the adaptation of fungal populations to changing environments (Milgroom, 1996, Heitman, 2006, Giraud et al., 2010). Increased genotypic variation allows some individuals in a given population to inherit variations of alleles that are more suitable for a specific condition such as virulence (Engering ). Genetic diversity is favoured by recombination that positively eliminates deleterious alleles that have accumulated during clonal development (Heitman, 2006, Barton, 2010). The two main sources of recombination in ascomycetes are sexual and parasexual reproduction (Calo ). Sexual reproduction in those fungi is orchestrated by the mating type locus (MAT1-1 and/or MAT1-2), which codify for transcription factors (α-box and HMG, respectively) that regulate genes required for mating and meiosis (Coppin ). Two mating recognition systems are well described in Ascomycota: homothallic (self-fertile) species carry both mating-type alleles (MAT1-1/MAT1-2 genes) in a single haploid genome, and sexual reproduction takes place in a single individual. Heterothallic (self-sterile) fungi are species that carry a single allele (MAT1-1 or MAT1-2) per haploid genome and need an opposite mating-type for sexual reproduction (Debuchy et al., 2010, Lee et al., 2010, Ni et al., 2011). The other main source of recombination in fungi is the parasexual cycle, where meiosis and development of sexual structures remain absent (Pontecorvo 1956). This non-sexual mechanism is governed and controlled by the het locus (heterokaryon incompatibility) that confers vegetative recognition of some filamentous ascomycetes (Glass ). This process is initiated by cytoplasmic fusion in which different nuclei and cytoplasmic organelles co-occupy the same cellular space. Nuclear fusion takes place producing a diploid nucleus (karyogamy), which is, however, unstable and produces segregants by recombination, evoking the mitotic crossing-over followed by haploidisation.

There are two direct indications that recombination takes place in BY. First, the Capronia sexual morph is a homothallic (self-fertile) genus and is polyphyletic, distributed over Herpotrichiellaceae and Cyphellophoraceae (Untereiner, 1995, Untereiner and Naveau, 1998, Untereiner, 2000), and the Trichomeriaceae have, morphologically, a very similar sexual morph in Trichomerium (Chomnunti ). Second, by in silico methods, recombination has been reported to occur in Cladophialophora carrionii (Deng ). However, as most species of Chaetothyriales lack a known sexual morph, only asexual morphs have been recognised. This leads us to hypothesize the following options: (1) species may lack a genomic apparatus for sexual reproduction, (2) species may be purely heterothallic but with cryptic sex taking place under specific conditions, or (3) species may have genomic signatures for a parasexual cycle under het locus control. In addition, we have identified a potential fused MAT1-1/MAT1-2 gene configuration in Exophiala sideris disposing both α-box and HMG domains (Fig. 9, Fig. S11). Cryptic homothallism was proven to occur in Curvularia homomorpha and Bipolaris luttrellii, which are ascomycete species presenting similar configurations as herein reported for Exophiala sideris (Fig. S11). Another remarkable finding was the occurrence of a mating regulator that codifies for a transcription factor carrying a DNA binding homeodomain motif in V. gallopava (Fig. 7). This domain is largely found in yeasts in Ascomycota and Basidiomycota (Martin ) and, according to (Lee ) this domain was lost during speciation of Pezizomycotina. We herein hypothesize that the HD domain is a potential mating regulator in Venturiales (Fig. 7) and suggest this domain was maintained more recently in Eurotiomycetes speciation.