Trichosporonales species have fused HD and P/R mating type loci

Since 2015, genome sequences have been published for Trichosporonales species from the genera Apiotrichum, Cutaneotrichosporon, Takashimella, Trichosporon, and Vanrija [28, 48–54] (Table 1). We analyzed the genomes of 29 isolates that belong to 24 Trichosporonales species for the organization of the MAT loci, and found that all of them contain fused MAT loci, i.e. the MAT loci are physically linked on the same contig, with the mating-type determining genes for the P/R locus (STE3 and the pheromone precursor genes) and the HD locus (SXI homeodomain genes) located ~55 kb apart from each other. The gene content between the P/R and HD loci is nearly identical in all species except for the early-branching genus Takashimella (Fig 2). Except for four strains described below, all analyzed MAT loci comprise either the combination of the HD gene SXI1 with the pheromone allele STE3a, or the combination of the SXI2 gene and STE3α allele (Table 1, Figs 1 and 2). Each MAT allele carries a single pheromone precursor gene in the vicinity of the STE3 gene (Fig 2, S1 Text). The SXI/STE3 combinations found in the majority of Trichosporonales are different from the MAT alleles in C. neoformans, where the SXI1 gene is combined with STE3α in the MATα allele, and SXI2 is combined with STE3a in the MATa allele [24, 38]. To avoid confusion with the C. neoformans nomenclature and to be compatible with allele designations in other basidiomycetes, we named the STE3α-containing MAT allele A1, and the STE3a-containing MAT allele A2 (Figs 1 and 2), with each allele containing physically linked HD and P/R loci.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. MAT loci of Trichosporonales.

Classic mating type-defining genes (SXI1 or SXI2, STE3, and the pheromone precursor genes) are shown in red, other genes that are part of the mating-type locus or flanking the mating-type locus (gene g6341) of C. neoformans are shown in green, and genes not present in the C. neoformans mating type locus are shown in white. The MAT allele is indicated on the right, with A1 and A2 MAT alleles shown above and below the dashed line, respectively. For simplicity, the A1* and A2* alleles are not included in this Figure, and are instead depicted in S5 and S6 Figs. The proposed MAT locus region is enclosed in a yellow box in each strain. Orthologs are connected by grey or orange bars when in the same or opposite orientations, respectively.

https://doi.org/10.1371/journal.pgen.1008365.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Genome assemblies that were analyzed in this study.

https://doi.org/10.1371/journal.pgen.1008365.t001

The fact that the genomes of 24 Trichosporonales species carry fused MAT loci with mostly identical gene content in two allelic combinations distributed throughout the phylogenetic tree (Fig 1, see also Fig 2 and the next results section) suggests that the event physically linking the two MAT loci predates the diversification of the Trichosporonales clade. However, at the beginning of this study, only one MAT allele was known for each of the Trichosporonales species. To verify that both MAT alleles, A1 and A2, can be found in different strains within a single species, we analyzed additional Vanrija humicola strains (Table 2). V. humicola is a member of a phylogenetically early-branching lineage within the Trichosporonales, and the previously sequenced strain (JCM1457) carries the A2 allele (Figs 1 and 2) [28]. Several other strains obtained from culture collections were analyzed by RAPD (Rapid Amplification of Polymorphic DNA) genotyping, which revealed banding patterns similar to those from the type strain V. humicola JCM1457 (= CBS571, MAT A2) (S1 Fig). These strains were analyzed further by PCR for the mating type genes, and based on the analysis of PCR amplicons of the HD gene-containing regions, strains CBS4282 and CBS4283 were preliminarily identified as MAT A1 strains. To assess this further, we sequenced the genome of CBS4282 using Illumina sequencing. The genome was assembled into 21 scaffolds with 5089 predicted genes. At 22.63 Mb, the assembly size is similar to that of strain JCM1457 (A2) (22.65 Mb, [28]). A k-mer analysis of the Illumina reads for CBS4282 showed a single peak as expected for a haploid genome (S1 Fig). Analysis of the MAT region revealed fused MAT loci within a ~50 kb region of scaffold 02 carrying the SXI2 gene and STE3α allele, confirming that this strain is indeed an A1 strain, and thus showing that alternative alleles can be found in different strains of a single Trichosporonales species (Figs 2 and 3, S2 and S3 Figs, Table 2). This was further confirmed by an analysis of the recently published genome of a second Apiotrichum porosum strain [53], with the two A. porosum genome sequences available now representing the A1 and A2 alleles (Figs 2 and 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Phylogenetic analysis of proteins encoded by genes within and around the MAT loci in Trichosporonales.

A. The A1 and A2 MAT alleles of V. humicola and A. porosum are shown for reference of the positions of the analyzed genes within the MAT region (depicted by black arrows). B. Individual genealogies of the selected genes across species. In the core MAT region (enclosed in yellow in panel A), a deep trans-specific polymorphism is only seen for the STE3 gene. In V. humicola, most of the analyzed genes within this region branch by mating type. The MAT alleles of the analyzed strains are given next to the phylogeny. Genus name abbreviations and other color features are given as in Figs 1 and 2.

https://doi.org/10.1371/journal.pgen.1008365.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Strains used in this study.

https://doi.org/10.1371/journal.pgen.1008365.t002

Alternative MAT alleles in Trichosporonales show significant chromosomal rearrangements, but each displays highly conserved structure across species

The gene order within the core MAT region between the SXI2 and STE3α/MFα genes is very well conserved between V. humicola CBS4282 (A1) and the A1 alleles of the other Trichosporonales species analyzed, whereas it differs from the V. humicola A2 strain JCM1457 by two inversions, one of which encompasses most of the core MAT region (Fig 2). The same is true for the A2 alleles, where the core MAT region is also rather conserved between V. humicola JCM1457 (A2) and other Trichosporonales (Fig 2).

In addition to the key mating-type determinants, the MAT loci of the Trichosporonales contain other genes previously shown to be required for mating and filamentation in other species (e.g. STE11, STE12 and STE20), and their common presence in the MAT loci of Tremellales suggests these genes were anciently recruited to the MAT locus (Fig 2). In the Trichosporonales, STE11 and STE20 are found within the core MAT region (between the HD and P/R genes), and STE12 is in the vicinity of the core MAT region in most species (Fig 2). One exception turned out to be STE12 in the Cutaneotrichosporon lineage. In C. oleaginosum strains IBC0246 (A2) and ATCC20509 (A2), STE12 is located on a different scaffold than the MAT locus. While STE12 in IBC0246 (A2) is located at the end of its respective scaffold, its location in ATCC20509 (A2) is around position 1.1 Mb of a 2.5 Mb scaffold, indicating that even if this scaffold was linked to the MAT scaffold, the linkage would not be very tight. To confirm that STE12 is indeed not linked to the MAT locus in C. oleaginosum, we sequenced the genome of the third of the three known strains of this species, ATCC20508. Analyses of the ATCC20508 (A2) and ATCC20509 (A2) genome assemblies showed that both the MAT locus and STE12 are located within large co-linear scaffolds in both strains, providing further evidence that STE12 is not linked to the MAT locus in this species (S4 Fig). Consistent with this, in the other Cutaneotrichosporon species analyzed, STE12 is either on the same scaffold as MAT, but at a distance of 400 to 750 kb (S5 Fig), or is present on different scaffolds than MAT. Thus, while this suggests that the location of STE12 in the vicinity of MAT likely represents the ancestral state in the Trichosporonales, its relocation to a different genomic location in the Cutaneotrichosporon species indicates that STE12 might not be fully linked to MAT in the Trichosporonales.

Another difference in MAT organization that can be observed in the genus Cutaneotrichosporon is the combination of HD genes and STE3 alleles. While all other strains of each species carry either the A1 (STE3α+SXI2) or A2 (STE3a+SXI1) MAT loci, Cutaneotrichosporon dermatis carries an allele that has the overall genomic organization of A1, but has instead a SXI1 gene (S5 Fig). Because the A1 vs. A2 designation is based on the STE3 variant, this allele was called A1*. Interestingly, A1* alleles can also be found in Apiotrichum veenhuisii, whereas a corresponding A2* allele, in which STE3a is linked to the SXI2 gene, can be found in Apiotrichum gracile (Fig 1 and S5 Fig). Furthermore, both A1* and A2* alleles can be found in Cutaneotrichosporon mucoides, one of several hybrid species within the Trichosporonales (S6 Fig, S2 Text) [28, 52]. These MAT alleles seem to represent a derived state and might have been independently generated by recombination involving crossing-over within the MAT region adjacent to the SXI1/SXI2 genes that is co-linear between the A1 and A2 alleles (S5 Fig).

Potential ancient origin of the fused MAT loci in Trichosporonales

The nearly identical gene content and overall high degree of conservation in gene order and allele combinations of HD and P/R loci in the Trichosporonales suggests that the fusion of the MAT loci most likely occurred in the common ancestor of the Trichosporonales lineage, although other hypotheses are possible, which are discussed below. To determine at what time this fusion might have occurred, we estimated divergence times of the basidiomycete lineages included in Fig 1 using three calibration points (see S7 Fig and Materials and Methods). Based on this analysis, the Trichosporonales and Tremellales have a common ancestor dating back approximately 179 million years ago (MYA), whereas the earliest split within the Trichosporonales occurred approximately 147 MYA (S7 Fig). Thus, under the assumption of a MAT fusion in the common ancestor of the Trichosporonales, this fusion would have occurred between 179 and 147 MYA (S7 Fig). While our estimates of divergence times for the Microbotryum and Trichosporonales lineages agree well with recent analyses [8, 56], it should be noted that previous estimates for the last common ancestor of the pathogenic Cryptococci resulted in earlier divergence times ranging from 40 to 100 MYA [57–60]. Nevertheless, even if the divergence time of the ancestor of the pathogenic Cryptococci was underestimated in our analysis, this would still place the last common ancestor of the Trichosporonales, and thus the likely MAT fusion in this group, at a much earlier time point.

Analysis of recombination suppression in the fused MAT loci in V. humicola

Suppression of recombination is a hallmark of sex chromosomes of animals and plants and the mating-type chromosomes of algae and fungi [6, 10, 29, 35, 38, 61, 62]. In basidiomycetes, recombination cessation between P/R and HD loci is often established following fusion (i.e. the physical linkage) of the two loci, and sometimes expanding far beyond the physically linked MAT loci [8, 29]. One consequence of recombination cessation can be the accumulation of transposable elements as well as increased genetic differentiation between allelic sequences [30]. This was observed in C. neoformans, where the MATa and MATα alleles differ significantly in gene organization, and the MAT locus contains more remnants of transposons and other repeat sequences than other genomic regions except for the centromeres and rDNA repeats [24, 63]. Furthermore, both C. neoformans MAT alleles are highly rearranged in comparison with their Cryptococcus gattii counterparts [38].

In contrast, the MAT A1 and A2 alleles of V. humicola are overall co-linear apart from two inversions, and this is generally the case throughout the Trichosporonales (Figs 2 and 3A). An analysis of repetitive sequences shows that there is no accumulation of repeat regions within the MAT alleles of V. humicola (S1 Table). Thus, while the two inversions might impair meiotic pairing and therefore recombination in this region, the conserved gene order and lack of repeats made us consider if the V. humicola MAT loci, and more generally the Trichosporonales MAT loci, are regions of suppressed recombination.

To test this, we first analyzed phylogenetic trees for several genes present within or adjacent to the MAT loci of Trichosporonales (Fig 3). For genes in regions undergoing meiotic recombination, alleles associated with alternative mating types are expected to display a species-specific topology in a phylogenetic tree, whereas genes in MAT regions should cluster by mating type if recombination suppression predates speciation (i.e. with the A1 alleles of the different species branching together rather than each of the alleles clustering with the A2 allele from the same species) [34, 38]. In the Trichosporonales, only STE3 clearly shows a mating type-specific pattern with an ancient trans-species polymorphism in Trichosporonales and Tremellales (Fig 3B). None of the other genes tested showed a mating type-specific phylogenetic pattern for the Trichosporonales at such a deep phylogenetic level, indicating that recombination was not suppressed at the base of the clade. However, for the genes STE20, IKS1, STE11, and MYO2 within the core MAT locus, sequences from two V. humicola A2 strains (JCM1457 and UJ1) group separately from the A1 sequence (CBS4282) (Fig 3B). Therefore, these genes might have evolved mating type-specific alleles at the species level, although more strains from more species would have to be investigated to further confirm this finding (Fig 3).

The possibility of mating type-specific alleles is further supported by BLASTN analyses comparing alleles in the three available V. humicola genomes, which showed that within the core MAT region, alleles from the A2 strains (JCM1457 and UJ1) are more similar to each other than to the A1 strain CBS4282, whereas outside of this region this is not the case (Fig 4A and 4B). This suggests that the MAT locus of V. humicola might be a region of recent recombination suppression.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. MAT alleles of V. humicola strains.

A. Two strains (JCM1457 and UJ1) are MAT A2, while strain CBS4282 is a MAT A1 strain. Genes shown in the same color (orange within the core MAT region, and purple outside of the MAT region) were the most similar pair in a pairwise BLASTN analyses among the three alleles of the different V. humicola strains. B. Schematic view of BLASTN results. For each analyzed gene, it is shown if alleles of the most similar pair came from the strains with the same (A2/A2) or different (A1/A2) MAT alleles. Within the core MAT region, higher sequence similarity was observed between the two A2 alleles than between A1 and A2 alleles. C. Synonymous divergence between MAT alleles in different V. humicola strains. Synonymous substitutions per synonymous site and standard errors (d_S +/- SE) are shown for pairwise comparisons between strains CBS4282 (A1), JCM1457 (A2), and UJ1 (A2) for genes in the region containing the mating type genes. For STE3 genes, for which no d_S value could be calculated in comparisons of strains with different mating types, values are shown as undefined (Un) at the top of the diagram. The inset below shows d_S values in the region between the core MAT genes, gene order as in strain CBS4282 (A1), STE3 is left out for clarity. D. Boxplots of distribution of d_S values shown in C (only genes for which values could be calculated were used in this comparison). The plots show the distribution of d_S values, with median value as a horizontal line in the box between the first and third quartiles. Outliers were left out for better visibility. Significantly different values are indicated (p-values from Student‘s t-test).

https://doi.org/10.1371/journal.pgen.1008365.g004

To further test whether recombination is suppressed between the MAT alleles in V. humicola, we analyzed levels of synonymous divergence (d_S) between alleles on the MAT-containing scaffold as well as on three other scaffolds of CBS4282 (A1) and the two V. humicola A2 strains (S8 Fig). In the four longest scaffolds, genes within large regions of the scaffolds are more similar between CBS4282 (A1) and the JCM1457 (A2) than between the two A2 strains, where stretches of low d_S values occur much less frequently (S8 Fig). Interestingly, we also noted regions of high similarity between CBS4282 (A1) and JCM1457 (A2) that are interrupted by stretches containing more divergent alleles. This pattern is consistent with ongoing recombination in natural V. humicola populations, even though sexual reproduction has not yet been observed.

In contrast to other genomic regions, including the co-linear regions outside of the core MAT region where overall divergence is lowest between CBS4282 (A1) and JCM1457 (A2) (Fig 4C and S8 Fig), most of the genes in the core MAT region show slightly higher levels of divergence between the A1 strain and the two A2 strains than between the two A2 strains (Fig 4C). This is consistent with an absence of genetic exchange between the A1 and A2 alleles, as one would expect if recombination in this region is suppressed. These findings might be explained by a reduced recombination rate in the core MAT region carrying the inversions in alternate MAT alleles, which could lead to accumulation of mutations in the two MAT alleles and thus to elevated synonymous divergence between A1 and A2 alleles. However, divergence between A1 and A2 strains within the core MAT region is only moderately elevated for most genes, and the difference between the average d_S values is statistically significant (P< 0.047, t-test statistics 1.8101) only when comparing the analysis of CBS4282 (A1) vs. UJ1 (A2) with the analysis of JCM1457 (A2) vs. UJ1 (A2) (Fig 4D). One possible explanation could be that genetic exchange within the inverted regions may not be completely inhibited, because exchange via non-crossover gene conversion or double crossover can still occur within the inverted regions [64]. This should result in a so-called suspension bridge pattern with divergence in the middle of an inverted region lower than towards the inversion breakpoints [65]. To test this, we performed BLASTN analyses using sliding windows of genomic sequences of the MAT locus and adjacent regions (S9 Fig). The results show less sequence similarity in the regions of putative inversion breakpoints between A1 and A2 strains compared to regions within the inversions and outside of the MAT region, consistent with the hypothesis that a certain amount of recombination is occurring within the inverted regions.

The incomplete suppression of recombination (and thus the possibility of ongoing genetic exchange) might be a reason why the gene order within each A1 and A2 allele in different species is surprisingly well conserved (Fig 2). To test if this degree of gene order conservation extends beyond the MAT locus, we compared the MAT-containing scaffolds of Trichosporonales. As expected for a range of species separated by millions of years of evolution, synteny is conserved between closely related species, but not between species that diverged long ago (S10 Fig), except for the MAT region. This suggests that the MAT region of Trichosporonales is an ancient cluster of tightly linked loci (also known as 'supergenes') that seem to segregate as a stable polymorphism within the populations of each species. Thus, the Trichosporonales appear to contain very stable MAT alleles with respect to gene order, combined with (slightly) suppressed recombination between different alleles.

The MAT loci in the Trichosporonales have a significantly lower GC content compared to the overall genomic GC content

One curious observation was found during the analysis of the Trichosporonales MAT region, namely a lower GC content in the MAT regions compared to the surrounding regions for the analyzed Vanrija, Cutaneotrichosporon, and Trichosporon strains (S10 and S11 Figs). One possible explanation for a lower GC content could be an accumulation of AT-rich transposable elements, but an analysis of repeats in the three V. humicola strains showed that there are only a few (strain CBS4282, A1) or no (strains JCM1457 and UJ1, both A2) repeats present within the MAT region of these strains (S1 Table). Another explanation might be a lower density of coding regions, which tend to have a higher GC content than non-coding regions. However, an analysis of the GC content only in the coding sequences within and around the MAT region of strain CBS4282 showed the same pattern of lower GC content within the MAT region (S11B Fig). Another possible explanation for the lower GC content might be the accumulation of mutations due to reduced recombination. Under the (simplistic) assumption that mutation frequencies for all nucleotide exchanges are similar, this would drive the GC content towards 50%. In genomes with an average GC content of more than 50%, this would appear as a region with lower GC content, and this would apply to the analyzed strains with an average genomic GC content of 58 to 63% (S11 Fig). It has been shown that spontaneous mutations tend to be AT-biased in many species including several fungi [66–70], and thus an increased accumulation of mutations would generally lead to a lower GC content in the corresponding regions. Both explanations can be supported by an analysis of the codon usage within the MAT region of the three V. humicola strains (S2 Table). Among the amino acids that can be encoded by more than one codon, there is a trend for GC-rich codons to be used less frequently in the MAT region compared to the rest of the genome, consistent with AT-biased (or GC content-equalizing) mutations combined with selection for conserved protein sequences due to functional constraints. Additional studies of species with a GC content of less than 50% might be useful to test these hypotheses.

Fused MAT loci evolved several times in basidiomycetes

In this study, we analyzed MAT loci from the Trichosporonales, and found that all of the analyzed species harbor fused MAT loci with a single HD gene, an arrangement that among the Tremellomycetes has so far only been found in the pathogenic Cryptococci [18]. Physically linked MAT loci have been identified previously in other basidiomycetes (Malassezia, Microbotryum, Sporisorium, Ustilago), but these MAT loci carry the ancestral arrangement of two HD genes except in Microbotryum violaceum caroliniana, in which the HD1 gene was specifically lost from the A2 mating type [8, 18, 20–23, 71, 72]. Fused MAT loci in fungi can be considered supergenes, which are linked genetic loci that facilitate the co-segregation of beneficial allele combinations [61]. Supergenes have evolved in a number of species, and can regulate different traits, e.g. mimicry in butterflies, social organization in ants, or mating in fungi [8, 73, 74]. Our findings add to the growing number of instances of supergene formation for fungal MAT loci, even though the molecular mechanisms maintaining supergenes, e.g. suppressed recombination, might not be identical in all cases as discussed below.

An analysis of several strains from V. humicola showed that there are two main MAT alleles (i.e. two main allelic combinations at the fused MAT loci), and both alleles are distributed throughout the Trichosporonales phylogenetic tree, although additional strains would have to be analyzed to establish if these two alleles are the most prevalent in this group. Nevertheless, if the mating system is predominantly inbreeding, rare alleles have no advantage and can be gradually lost by genetic drift. Thus, fusion of MAT loci, which can evolve by selection because it is beneficial in selfing mating systems, is predicted to lead to a reduction in the number of MAT alleles [4, 75]. In addition, theoretical modelling predicts that a combination of facultative and rare sexual reproduction, low mutation rates, and a small effective population size should also lead to a reduction in the number of MAT alleles [76]. Population sizes and mutation rates are, however, not known for Trichosporonales species, and no sexual development has been observed so far making it possible that asexual reproduction is the predominant form of propagation for many Trichosporonales. Another possibility could be that only two alleles per MAT locus existed in a tetrapolar ancestor due to other factors, such as linkage of the MAT loci to their respective centromeres and concomitant loss of allele diversity, as recently proposed for tetrapolar Microbotryum species [77]. Analyses of species branching at the base of the Trichosporonales lineage could help to clarify this hypothesis. It should be noted that the HD and P/R loci are located at the two ends of the MAT region with an inversion between opposite alleles. Therefore, an odd number of crossovers within the region between HD and P/R would result in genetically imbalanced progeny, which would be inviable, whereas an even number of crossovers would not result in allele exchange at the HD or P/R loci. Thus, recombination could occur, but would most likely not change the HD and P/R allele combinations. In the case of the A1* and A2* alleles, recombination might have occurred in the rather small co-linear genomic region between the SXI1/SXI2 genes and the adjacent inversion, and thus might be a rare event.

Convergent evolution of fused MAT loci in Trichosporonales and Tremellales

Overall, our finding of fused MAT loci with a single HD gene in all of the Trichosporonales species investigated is most consistent with the hypothesis that a fusion event of the HD and P/R loci followed by loss of one HD gene per MAT allele occurred independently in the Trichosporonales and the pathogenic Cryptococci lineages. The alternative hypothesis that the fusion of the HD and P/R loci occurred in the common ancestor of the Trichosporonales and pathogenic Cryptococci is less parsimonious as it would imply multiple independent reversions to tetrapolarity in the non-pathogenic Cryptococci. In addition, while several mating-associated genes can be found in the MAT loci of both groups, the majority of genes within the MAT locus of the Trichosporonales is not found in the MAT locus in the Tremellales, consistent with independent fusion events in the two lineages. This hypothesis is also supported by the combination of the SXI1/SXI2 genes and the STE3a/α alleles in the A1 and A2 alleles of the Trichosporonales, which is different from those of the MATa and MATα alleles of the pathogenic Cryptococci. A model is proposed in Fig 5 to explain the current situation in the Tremellomycetes, where recruitment of several mating-associated genes into the P/R loci of a tetrapolar ancestor occurred first as these genes can also be found within the P/R loci of extant Tremellales with unlinked MAT loci [33–35]. After the split of the Tremellales and Trichosporonales lineages, fusion of the HD and P/R loci as well as the loss of one HD gene per MAT allele occurred in the ancestor of the Trichosporonales, whereas within the Tremellales this happened only in the ancestor of the pathogenic Cryptococci.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Model for the evolution of the MAT loci in Tremellomycetes.

Genes at the core MAT loci include the homeodomain transcription factor genes (HD locus, shown in red/pink) and the pheromone and receptor genes (P/R locus, shown in blue). Genes involved in sexual development, but not originally part of a MAT locus are shown in green, other genes are shown in grey. Only genes from the C. neoformans MAT locus that are also linked to the core MAT genes (STE3, HD genes) in the Trichosporonales are shown (STE11, STE12, STE20, IKS1, MYO2, RPL22). Other genes present at the MAT loci are left out for clarity. A trend towards integrating other developmental genes into the MAT loci is reflected in the recruitment of the STE12 gene and an ancestral STE11/STE20 cluster into the P/R locus. This resulted in a tetrapolar arrangement with (at least) two alleles for each MAT locus. In the Tremellales, tetrapolarity is retained in the extant Tremella, Kwoniella, and non-pathogenic Cryptococcus lineages, whereas in the pathogenic Cryptococci, the MAT loci fused and one HD gene was lost from each MAT allele. Fusion of the MAT loci occurred independently and earlier in the ancestor of the Trichosporonales, also with loss of one HD gene, resulting in two known MAT alleles with combinations STE3a-SXI1 and STE3α-SXI2 that differ from those in the known Cryptococcus alleles. The STE12 gene is shown in outline only to indicate that it is not part of the MAT locus in all species; STE12 may have been recruited to the P/R locus in the Tremellales-Trichosporonales common ancestor and eventually evicted from the MAT locus in the Trichosporonales possibly associated with the fusion event of P/R and HD loci, or instead STE12 was recruited to MAT only in the Tremellales. The phylogenetic relationships and inferred dates of speciation (numbers in black on tree nodes) are depicted according to S7 Fig. The blue bars correspond to 95% confidence intervals.

https://doi.org/10.1371/journal.pgen.1008365.g005

An independent fusion might suggest that similar selective pressures were acting to result in similar evolutionary trajectories in the Trichosporonales and pathogenic Cryptococci. In the latter group, it has been hypothesized that the pathogenic lifestyle might make it difficult to find a mating partner while associated with a host. In this case, linkage between the two MAT loci is expected to be favored by selection, as this increases the odds of compatibility between the gametes derived from a single diploid zygote [4, 8, 9, 18, 29]. Extant members of the Trichosporonales can be associated with hosts as commensals or pathogens, but this lineage also includes many soil- or water-associated saprobes [40, 41]. It is possible that the ancestor in which MAT fusion occurred was associated with a host, and that the physically linked MAT loci remained stable during evolution of species with saprobic lifestyles. However, a scarcity of compatible mating partners might also occur under other conditions, e.g. if a population is derived from few progenitor cells that propagated through mitotic cell divisions during favorable conditions and switched to sexual development when nutrients were depleted.

An alternative hypothesis to the fusion of the MAT loci occurring in the common ancestor of the Trichosporonales would be multiple independent fusions in different lineages within the Trichosporonales. However, in this case, one would expect different genomic arrangements, similar to what has been observed in several cases of convergent evolution of fused MAT loci in the genus Microbotryum [8]. It is, in principle, possible that there is strong selective pressure that favors the specific genomic configurations observed in the Trichosporonales MAT loci for the A1 and A2 alleles, but the nature of this selective pressure, which would have to act on multiple species in different environments, is not clear. Another hypothesis could be that the fusion of MAT loci is ancient, but the inversion, which may have contributed to limiting or repressing recombination, occurred recently and independently in several Trichosporonales lineages. This would explain the lack of trans-species polymorphisms observed for most of the genes at the fused MAT loci. Such a scenario of multiple independent inversions would also require strong selective pressure in multiple species to yield not only functionally similar, but essentially the same genomic configuration. Under the hypotheses of more recent, independent fusions or inversions, one might expect to find species within the Trichosporonales lineage that have not undergone MAT fusion (i.e. physical linkage) or inversions, similar to the genus Microbotryum, where tetrapolar and bipolar species are intermingled in the phylogenetic tree [8]. So far, such species have not been identified in the Trichosporonales. Therefore, the hypothesis of a fusion event in the last common ancestor of the Trichosporonales appears to be the most parsimonious at present, despite several unusual genomic features discussed below.

In addition to the fused MAT loci, the presence of only one HD gene per MAT allele in all examined Trichosporonales is similar to pathogenic Cryptococci. This was noted in the initial analysis of genome sequences for two Trichosporonales strains now designated as MAT A2 [46], and the current study shows that this is a conserved feature throughout the Trichosporonales. Outside of the Tremellomycetes, the majority of analyzed basidiomycete species harbor at least two HD genes per MAT allele. In the case of two HD genes, one belongs to the HD1 class, and the other to the HD2 class, whereas in the Agaricomycetes, species with more than one HD1 or HD2 gene per MAT allele are known [18]. The HD1 and HD2 genes from the same MAT allele are not compatible with each other, irrespective of whether the mating system is tetrapolar or bipolar [16, 18]. An exception is the genus Wallemia, where fused MAT loci exist with only one HD gene, but until now, only one MAT allele harboring a SXI1 homolog is known in this genus, whereas another mating type appears to not contain any SXI gene [18, 78, 79]. It might be that mating in Wallemia species does not rely on HD compatibility. Functionally, two compatible HD1 and HD2 genes from different MAT alleles are necessary and sufficient for sexual development not only in C. neoformans, where only one HD1/HD2 gene pair is present after mating [80, 81], but also in U. maydis and Coprinopsis cinerea, where two (U. maydis) or more (C. cinerea) compatible HD1/HD2 combinations can be present after mating [82–84]. The presence of multiple HD gene paralogs or increased HD allele diversity within a species is advantageous under outcrossing, because it allows more frequent dikaryon viability after mating when the P/R locus also exists in multiple alleles [16]. However, if a species is predominantly selfing, a single MAT locus with two alleles will provide the highest percentage of compatibility between gametes from a single tetrad [4, 16, 18]. In such cases, the loss of one HD gene per MAT locus should not be a problem from either a functional or evolutionary point of view, and thus this genomic configuration might be observed in other species with fused MAT loci, unless the HD genes have other functions unrelated to mating. Another possible reason for the scarcity of HD loci with only a single HD gene might be that having two different HD genes occupying allelic positions results in a genomic region that might not be able to properly pair during meiosis. In the ascomycete Neurospora crassa, unpaired regions are known to trigger MSUD (meiotic silencing of unpaired DNA) through an RNAi-based mechanism, although interestingly the mating type genes of N. crassa are immune to silencing in their normal genomic location [85–87]. Thus, it might be possible that MAT loci with a single HD gene have properties at the DNA level that are independent of the functions of the encoded proteins.

One common feature of fungal MAT loci and sex chromosomes in other eukaryotes might be the recruitment of sex-associated genes. The presence of the STE11 and STE20 within the core MAT loci of Trichosporonales and Tremellales suggests an ancient linkage of these genes to the MAT locus. In C. neoformans, STE20 is required for the formation of proper heterokaryotic filaments and basidia after mating, and STE11α is required for mating and filament formation [88, 89]. STE20 displays MAT-specific alleles [33, 34] and is also located within the P/R locus of other basidiomycete species, e.g. Leucosporidium scottii and red yeasts in the Pucciniomycotina [17]. However, the functions of many genes within the MAT loci are not known, and further experiments will be needed to determine if these genes have functions in mating or sexual reproduction.

Genomic signatures that distinguish the MAT loci from surrounding genomic regions in Trichosporonales

A common feature of fungal MAT loci and the sex-determining regions in other eukaryotic groups is that recombination is usually suppressed in these regions. In several fungi, suppressed recombination is observed not only within the MAT locus itself, but spreading out from the MAT locus along the MAT-carrying chromosome resulting in so-called evolutionary strata of stepwise recombination suppression, similar to findings in sex-chromosomes of mammals and plants. Examples are the ascomycete Neurospora tetrasperma and the basidiomycete genus Microbotryum [8, 10, 29, 32, 90–94]. In the latter, it was recently shown that the linkage of HD and P/R loci evolved independently at least five times within the genus, followed in each case by further stepwise expansions of the regions of suppressed recombination beyond the fused MAT loci, possibly through neutral rearrangements [8, 29]. In contrast, the analysis of V. humicola did not yield indications for suppressed recombination beyond the core MAT region that lies between the SXI1/SXI2 and pheromone genes. Furthermore, the gene content and gene order within the MAT loci and the HD/STE3 allele combinations in the analyzed Trichosporonales species are largely conserved, suggesting that the fusion of HD and P/R loci occurred once in a common ancestor. Regions outside of the MAT locus have generally undergone substantial rearrangements (e.g. translocations) in different species as seen by the varying genomic locations of STE12 in the genus Cutaneotrichosporon and the lack of synteny outside of the MAT region in more distantly related Trichosporonales.

In Microbotryum, it was shown that evolutionary strata ranging from mating-type- to species-specific patterns can be distinguished along chromosomes harboring recently linked MAT loci [8, 29]. These analyses were performed based on alleles from several strains per species and of both mating types. One challenge for the phylogenetic analysis in the Trichosporonales is that most species are currently represented by a single strain. Analyzing additional alleles from more strains of each species will be necessary to fully evaluate the phylogeny of the genes located within and surrounding the MAT locus.

In the Trichosporonales, the two MAT alleles, A1 and A2, differ in general by two inversions. If the MAT fusion is relatively ancient, i.e. occurred in the last common ancestor of the Trichosporonales, a remaining question is why the genetic differences in the two known alleles are not more pronounced. One explanation could be that recombination within the inverted regions might still be possible via non-crossover gene conversion or double crossover [65]. This would be consistent with the pattern observed in V. humicola of lower divergence between A1 and A2 alleles within the inverted regions compared to divergence at the inversion breakpoints. In C. neoformans, gene conversion occurred in a GC-rich intergenic region within the MAT locus and was proposed as a mechanism for maintaining functionality of those genes within the MAT locus that are essential [39]. Gene conversion has also been observed in regions with suppressed recombination in the mating-type chromosomes of N. tetrasperma, in the mating type locus of the green algae Chlamydomonas reinhardtii, and in sex chromosomes of animals [95–101].

Gene conversion might be associated with an increase in GC content as it has been observed in a number of eukaryotes that gene conversion can be biased towards GC [102–104], which would not be compatible with the observed lower GC content in the MAT region of Trichosporonales. However, a recent analysis of two fungal species found no evidence of a GC-bias in gene conversion [105], and therefore at present gene conversion as an explanation for low divergence between MAT alleles in Trichosporonales remains possible. It is tempting to speculate that the lower GC content observed in the MAT region of Trichosporonales is connected to the incomplete suppression of recombination and/or the (yet unknown) mechanisms that lead to the low levels of divergence between the MAT alleles. In C. neoformans, two regions of higher GC content outside of the MAT locus are associated with increased recombination [106]. A correlation between higher GC content and increased recombination has been demonstrated in humans, although it is not clear in this case what is cause and what is effect, whereas in the yeast Saccharomyces cerevisiae, a correlation was described at the kilobase range, but no long-range correlation was found [107–109].

Another hypothesis to explain the low levels of divergence between the MAT alleles could be same mating-type mating ("unisexual reproduction", i.e. sexual reproduction without the need of a partner with a different mating type), which occurs in several pathogenic Cryptococci and the human pathogenic ascomycete yeast Candida albicans [110, 111]. Similar to reproduction after fusion of MATa and MATα cells, "unisexual reproduction" in Cryptococci also entails diploidization and meiosis including genetic exchange [111–113]. Diploidization can occur without cell fusion through endoreplication, or through fusion of two cells carrying the same MAT allele [114]. In the latter case, this would allow recombination within the MAT region, which could prevent degeneration of the MAT locus despite suppressed recombination when paired with a different MAT allele. Recombination within the MAT locus during "unisexual reproduction" was shown for C. neoformans [115]. So far, no sexual cycle has been observed in the Trichosporonales, and therefore it is not known what forms of sexual reproduction occur in this group. This is an important open question that requires detailed investigation in future studies.

Strains and growth conditions

C. oleaginosum and V. humicola strains used in this study are given in Table 2. Strains were grown on YPD medium at 25°C on solid medium or in liquid culture with shaking (250 rpm) for preparation of nucleic acids.

CHEF (contour-clamped homogeneous electric field) electrophoresis analysis of the C. oleaginosum isolates

CHEF plugs of the C. oleaginosum isolates were prepared and the electrophoresis was carried out as described in previous studies with slight modification [33, 35]. Specifically, the plugs were run using two different sets of switching time to get better resolution of the larger (S4A Fig top, 20–30 minutes linear ramp switching time) and smaller (S4A Fig bottom, 120–360 seconds linear ramp switching time) chromosomes, respectively.

Genotyping and analysis of mating type genes in V. humicola strains

A collection of V. humicola isolates were genetically screened by modified RAPD using primer 5’-CGTGCAAGGGAGCACC-3’ with 48°C as annealing temperature (S1A Fig). The V. humicola strains showing identical genotypic profile as the type strain CBS571 (Table 2) were further analyzed for the presence of the SXI1, SXI2, STE3, and MYO2 genes using PCR with oligonucleotides given in S3 Table. PCR fragments were either sequenced directly or cloned into pDrive (Qiagen, Hilden, Germany) and sequenced.

Genome sequencing, assembly, and annotation of V. humicola CBS4282 and C. oleaginosum ATCC20508

Genomic DNA samples were extracted using a modified CTAB protocol as previously reported [35, 116]. Specifically, to enrich samples with high molecular weight, after precipitation, the genomic DNA was picked out from the solution instead of spun down, and the samples were checked by CHEF for their sizes and integrity, following manufacturer’s protocol (BioRad, Hercules, CA, USA). Sequencing of the ATCC20508 and CBS4282 genomes was carried out at the Sequencing and Genomic Technologies Core Facility of the Duke Center for Genomic and Computational Biology, using large insert library (15–20 kb) and PacBio Sequel (2.0 chemistry) for ATCC20508, and generating 151 nt paired-end Illumina reads for CBS4282. The PacBio sequence reads for C. oleaginosum ATCC20508 were assembled using the HGAP4 assembly pipeline based on the Falcon assembler [117] included in the SMRT Link v5.0.1 software by Pacific Biosciences. The Illumina sequence reads of V. humicola CBS4282 were trimmed using Trimmomatic (v0.36) [118] with the following parameters to remove adapter contaminations: ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10:1:TRUE LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:40. Trimmed reads were assembled with SPAdes (v3.11.1) [119], and contig sequences were improved using pilon [120] based on the Illumina reads mapped to the assembly using Bowtie2 (v2.2.6) [121]. k-mer frequencies were analyzed based on the CBS4282 Illumina reads as described previously [122, 123]. Gene models for the newly sequenced strains as well as for strains for which genome assemblies, but not annotation was available from GenBank were predicted ab initio using MAKER (v2.31.18) [124] with predicted proteins from C. oleaginosum as input [46].

Analysis of mating type regions, synonymous divergence, and repeat content

MAT regions in Trichosporonales genomes were identified by BLAST searches [125] against the well-annotated MAT-derived proteins from C. neoformans [38], and manually reannotated if necessary. The short and not well conserved pheromone precursor genes were not among the predicted genes, and were identified within the genome assemblies using custom-made Perl scripts searching for the consensus sequence M-X(15–60)-C-[ILMVST]-[ILMVST]-X-Stop.

Synteny analysis of the genomes of V. humicola strains and C. oleaginosum strains was done with nucmer from the MUMmer package (v3.23) [126]. Synteny plots were drawn with Circos [127]. Synteny between MAT regions of different Trichosporonales species was based on bidirectional BLAST analyses of the corresponding predicted proteins. For the d_S plots, alleles in V. humicola strains were identified by two-directional BLAST analysis, and MUSCLE (v3.8.31) [128] was used to align the two alleles per gene per strain pair. Synonymous divergence and standard errors were estimated with the yn00 program of the PAML package (v4.9) [129]. Analysis of transposable elements and other repeats in V. humicola genomes was performed with RepeatMasker (Smit AFA, Hubley R, Green P. RepeatMasker Open-4.0. 2013–2015, http://www.repeatmasker.org) based on the RepbaseUpdate library [130] and a library of de novo-identified repeat consensus sequences for each strain that was generated by RepeatModeler (Smit AFA, Hubley R. RepeatModeler Open-1.0. 2008–2015, http://www.repeatmasker.org/RepeatModeler/) as described [123].

Linear synteny comparison along the MAT-containing scaffolds (S10 Fig) was generated with Easyfig [131] using a minimum length of 500 bp for BLASTN hits to be drawn.

Species tree and gene genealogies

Of the 24 Trichosporonales species for which genome assemblies were available we selected only 21 for phylogenetic analysis as the remaining three species (Trichosporon coremiiforme, Trichosporon ovoides and Cutaneotrichosporon mucoides) were shown to be hybrids in previous studies [28, 52]. Additionally, we selected other well-studied bipolar and tetrapolar representatives belonging to the three major Basidiomycota lineages, plus two ascomycetes as outgroup. To reconstruct the phylogenetic relationships among the selected members, the translated gene models of each species were clustered by a combination of the bidirectional best-hit (BDBH), COGtriangles (v2.1), and OrthoMCL (v1.4) algorithms implemented in the GET_HOMOLOGUES software package [132] to construct homologous gene families. The Cryptococcus neoformans H99 protein set was used as reference and clusters containing inparalogs (i.e. recent paralogs defined as sequences with best hits in its own genome) were excluded. A consensus set of 32 protein sequences was computed out of the intersection of the orthologous gene families obtained by the three clustering algorithms. Protein sequences were individually aligned with MAFFT v7.310 [133] using the L-INS-i strategy and poorly aligned regions were trimmed with TrimAl (-gappyout). The resulting alignments were concatenated to obtain a final supermatrix consisting of a total of 21,690 amino acid sites. We inferred a maximum-likelihood phylogeny using the LG+F+R5 model of amino acid substitution in IQ-TREE v1.6.5 [134]. Branch support values were obtained from 10,000 replicates of both ultrafast bootstrap approximation (UFBoot) [135] and the nonparametric variant of the approximate likelihood ratio test (SH-aLRT) [136].

For phylogenetic analysis of selected genes within and outside the MAT region in Trichosporonales, protein alignments were generated and trimmed as above and subsequently used to infer maximum likelihood phylogenies in IQ-TREE. Consensus trees were graphically visualized with iTOL v4.3.3 [137].

Divergence time analysis

Divergence times were estimated in MEGA-X [138] using the RelTime approach [139]. Contrary to other methods, RelTime does not require assuming a specific model for lineage rate variation and was shown to be as accurate as other approaches using relaxed and strict molecular clock models [140]. The reconstructed species tree, with branch lengths in the units of number of substitutions per site, was used as input and transformed into an ultrametric tree with relative times. The final timetree was obtained by converting the relative node ages into absolute dates by using three calibration constraints: 0.42 million year (MY) corresponding to the divergence between Microbotryum lychnidis-dioicae and Microbotryum silenes-dioicae [141]; 41 MY for the Ustilago—Sporisorium split; and 413 MY representing the minimum age of Basidiomycota. The latter two calibration points were obtained from the Timetree website (http://www.timetree.org/), which should be referred to for additional information and references.

Data availability statement

The V. humicola CBS4282 genome sequence (BioProject PRJNA475686) has been deposited at DDBJ/ENA/GenBank under the accession QKWK00000000. The version described here is version QKWK01000000. Illumina reads have been deposited in the NCBI SRA database under accession SRP150316. The C. oleaginosum ATCC20508 genome sequence (BioProject PRJNA475739) has been deposited at DDBJ/ENA/GenBank under the accession QKWL00000000. The version described here is version QKWL01000000. Pacific Biosciences reads have been deposited in the NCBI SRA database under accession SRP150334.