All strains sequenced in this study were obtained from the Japan Collection of Microorganisms (JCM), the RIKEN BioResource Center (BRC), Tsukuba, Japan, and the National Biological Resource Center (NBRC), National Institute of Technology and Evolution (NITE), Kisarazu, Japan (Table 1). Strains were incubated for one to five days in YM liquid medium (10 g glucose, 3 g yeast extract, 3 g malt extract, and 5 g peptone per litre) at 25 °C, with shaking at 100 rpm.

From R. toruloides strains JCM 10020, JCM 10021, JCM 10049, JCM 24501 and L. creatinivora (JCM 10699), cells were freeze-dried and crushed with a mortar and pestle, genomic DNA was extracted according to Raeder and Broda (1985) and purified with Genomic-tip 100/G columns (Qiagen, Netherlands); Illumina libraries with insert sizes of approximately 240 bp were constructed with a TruSeq DNA PCR-free (Illumina, USA) preparation kit, and the libraries were paired-end sequenced using the HiSeq 2500 (Illumina, USA). In the case of the other 31 strains (Table 1), genomic DNA was extracted after crushing the cell wall mechanically with a mortar and pestle and, following a protocol optimized for Westase enzyme (Takara, Japan), and genomic libraries with insert sizes of 350 bp were constructed with an Illumina DNA Prep, (M) Tagmentation kit (Illumina, USA), and the libraries were paired-end sequenced on a NextSeq 500 (Illumina, USA). For all libraries, sequences of 151 bp were obtained.

The genome sizes were estimated with GenomeScope1.0 (Vurture et al. 2017) following k-mer (k = 21) counting with Jellyfish 2.3.0 (Marcais and Kingsford 2011). Draft genomes were assembled with Allpath-LG v52155 (Gnerre et al. 2011) for JCM 10020, JCM 10021, JCM 10049, JCM 10699, and JCM 24501, and with SPAdes genome assembler v3.15.3 (Prjibelski et al. 2020) for the other strains. Completeness of assembled genomes was assessed with BUSCO 5.4.2 (Manni et al. 2021) with fungi_odb10 (n = 758) selected as the reference database.

Genes were predicted with MAKER 2.31.8 (Cantarel et al. 2008) using the combination of GeneMark-ES 4.21 (Lomsadze et al. 2005) and Augustus 3.0.3 (Stanke et al. 2006) for JCM 10020, JCM 10021, JCM 10049, JCM 10699, and JCM 24501, and with BRAKER 2.1.6 (Bruna et al. 2021) using the combination of GeneMark-ES 4.69 (Lomsadze et al. 2005) and Augustus 3.4.0 (Stanke et al. 2008). The protein dataset of RefSeq R. toruloides NP11 (GCF_000320785.1, Zhu et al. 2012) was used as an Augustus hint.

For the construction of a concatenated BUSCO single-copy tree, conserved single-copy orthologs were searched from protein sequences with BUSCO 5.4.2 (Manni et al. 2021) with fungi_odb10 (n = 758) selected as the reference database. The sequences of 331 BUSCO orthologs found as complete single copies in all strains were extracted with SeqKit 2.2.0 (Shen et al. 2016). Concatenated sequences were aligned with MAFFT 7.508 (Katoh et al. 2013). The phylogenetic tree was constructed with IQ-tree 2.2.0.3 (Nguyen et al. 2015) with “-m MFP -B 1000 --alrt 1000 --abayes -T 8” options. The gene concordance factors (gCF) and site concordance factors (sCF) were calculated using the WAG amino acid model and executed with the default commands, following the instructions on the IQ-tree website. For the construction of a tree integrating trees for all orthologous genes, OrthoFinder 2.5.4 (Emms and Kelly 2019) was used. For this analysis a species tree was obtained with STAG (Emms ane Kelly 2018) while the integrated gene tree was rooted with STRIDE (Emms and Kelly 2017), both of which were included in the OrthoFinder execution. All applications were run with default options. The distance-based phylogenetic tree with Balanced minimum evolution (BME) support values was constructed with JolyTree 1.1b (Criscuolo 2019; Criscuolo 2020) using default options.

The common orthogroup (OG) score was calculated following the description in Takashima et al. (2019), using OrthoFinder 2.5.4 (Emms and Kelly 2019) to identify OGs. The pairwise average amino acid identity (AAI) was calculated using CompareM v. 0.1.2 (https://github.com/dparks1134/CompareM) with default parameters. The shared OG rate was calculated based on the results obtained from OrthoFinder 2.5.4 (Emms and Kelly 2019).

The amino acid sequences of predicted genes were assigned to KEGG orthologs (KOs) using the GhostKOALA web server (Kanehisa et al. 2016, accessed Aug 3, 2023). The enzymes involved in the conversion of each nutrient were identified from KEGG metabolic pathway maps obtained via the KEGG Mapper website (Kanehisa and Sato 2020, accessed Aug 3, 2023). The physiological traits referred to were based on the findings of Kurtzman et al. (2011), Satoh and Makimura (2008), and Huang et al. (2011).

Assembly genomes and gene catalogs used as reference, i.e. R. toruloides NP11 (Zhu et al. 2012), R. toruloides IFO 0559 (Zhang et al. 2016), R. toruloides IFO 1236 (Zhang et al. 2016), R. toruloides IFO 0880 v4.0 (Coragetti et al. 2018), R. toruloides ATCC 204091 (Paul et al. 2014), R. taiwanensis MD1149 (Tkavc et al. 2018), R. graminis WP1 (Firrincieli et al. 2015), L. creatinivorum UCDFST 62-1032 (Mondo et al. 2017), M. intermedium 1389 BM 12 12 (Branco et al. 2017), M. lychnidis-dioicae p1A1 Lamole (Perlin et al. 2015), were downloaded from the JGI MycoCosm web portal (Grigoriev et al. 2014). Details are summarized in Suppl. material 1: table S1.

AAI: average amino acid identity

BRC: BioResource Research Center

ITS: internal transcribed spacer

JCM: Japan Collection of Microorganisms

KEGG: Kyoto Encyclopedia of Genes and Genomes

KO: KEGG ortholog

LSU: large subunit

NBRC: National Biological Resource Center

OG: orthogroup

rRNA: ribosomal RNA

We sequenced 35 yeast strains from 27 species of Sporidiobolales and a Leucosporidium strain, maintained in RIKEN BRC-JCM and NBRC, including type strains of nine Rhodotorula, six Rhodosporidiobolus, and eight Sporobolomyces species (Table 1). Expected genome sizes were determined by k-mer analysis and estimated for Sporidiobolales to be between 17 and 26 Mbp whereas that of L. creatinivorum was approximately 28 Mbp (Table 1). The predicted genome sizes of Rhodosporidiobolus were slightly bigger than those of Rhodotorula and Sporobolomyces. The k-mer histograms showed a single peak for each genome, suggesting that these strains have a haploid karyotype (Suppl. material 1: fig. S1), although for a few strains such as S. pararoseus JCM 5350 and S. roseus JCM 5353 small satellite peaks were observed, which can be caused by sequence errors or aneuploidy.

In line with the predicted genome sizes, the genome assemblies for Sporidiobolales were between 17 and 26 Mbp, and that of L. creatinivorum was approximately 28 Mbp (Table 1). These genome sequences are also similar to those of the reference genome (Suppl. material 1: table S1). Also, genome-assembly sizes of Rhodosporidiobolus were slightly bigger than those of the other Sporidiobolales. The comprehensiveness of the genomes assessed with BUSCO was between 78% and 92% with up to 1.2% duplicated genes, suggesting that these genome assemblies represent haploid cells (Table 1).

Homology-based gene prediction estimated that Sporidiobolales strains have approximately 6,400 to 9,800 genes, while L. creatinivorum has 10,454 genes (Table 2). Consistent with the assembly size, the gene numbers of Rhodosporidiobolus were relatively bigger than those of Rhodotorula and Sporobolomyces. The completeness of the BUSCO assessment was 93% to 97% with up to 1.3% duplicated genes, also suggesting that these gene data are highly comprehensive and indicative for haploid gene sets (Table 2).

We predicted the phylogenetic relationships between these strains from whole-genome data, including Sporidiobolales and other Microbotryomycetes (Leucosporidium and Microbotryum) the data for which were retrieved from the JGI Mycocosm (Grigoriev et al. 2014). We constructed phylogenetic trees using three methods: a tree from concatenated protein sequences of BUSCO genes found as a single copy in all strains (Fig. 1), a tree integrating gene trees of all orthologs (Fig. 2), and an alignment-free tree based on the distance of genomic nucleotide sequences (Fig. 3).

Figure 1. 

A phylogenetic tree based on concatenated protein sequences of 331 BUSCO genes found as single copy in all strains. The numbers beside the nodes represent the reliability based on SH-aLRT, aBayes, ultrafast bootstrap values, gene concordance factors (gCF), and site concordance factors (sCF) from left to right.

Figure 2. 

A phylogenetic tree inferred from gene trees of 3,181 orthologs found in all strains. The root is estimated with STRIDE. The numbers beside the nodes represent STAG support values.

Figure 3. 

A phylogenetic tree based on the similarity of whole-genome nucleotide sequences. The numbers beside the nodes represent balanced minimum evolution (BME) support values.

The BUSCO assessment of the gene catalog for each species found that 331 BUSCO orthologs (fungi_odb10) were present as single copies in all strains (Suppl. material 1: table S2). A phylogenetic tree was based on the concatenated amino acid sequences of these 331 BUSCO genes and can be considered as a phylogenetic tree of marker genes expanded to the whole-genome level. This tree supported the monophyly of each of the genera Rhodotorula, Rhodosporidiobolus, and Sporobolomyces (Fig. 1). Rhodotorula is first divided into two clades: one with R. glutinis, R. graminis, R. babjevae, R. araucariae, and R. paludigena, and the other with R. taiwanensis, R. sphaerocarpa, R. mucilaginosa, R. dairenensis, R. pacifica, and R. toruloides. The second clade is further divided into the R. toruloides clade and others. Sporobolomyces can be divided into two clades: one with S. johnsonii and S. salmonicolor, and the other with S. pararoseus, S. salmoneus, S. carnicolor, S. ruberrimus, S. phaffii, S. koalae, S. roseus, and S. blumeae. The general topology is similar to the phylogenetic tree based on concatenation of seven marker genes in a previous study (Wang et al. 2015a), suggesting the robustness of the topology when relying on the concatenation of conserved genes regardless of whether the used sequences are those for nucleic acids or proteins. The support values of nodes obtained by bootstrapping, aBayes, or the SH-aLRT method were generally high, although they may have been overestimated, a corollary associated with concatenation-based multi-locus phylogenetic trees (Seo 2008). When reliability was assessed using gene concordance factors (gCF) and site concordance factors (sCF), most major nodes showed gCF values over 60 and sCF values over 40. However, not all values were high; even among branches strongly supported by methods such as bootstrapping, aBayes, and SH-aLRT with very high scores, there were often cases where the gCF and sCF values were low.

Orthologs were retrieved with OrthoFinder, which assigned 361,309 genes to 13,277 orthogroups, of which 3,181 orthogroups were found in all strains. A phylogenetic species tree inferred by orthologous gene trees also supported the monophyly of each of the genera Rhodotorula, Rhodosporidiobolus, and Sporobolomyces, as well as the subclades mentioned above (Fig. 2). This topology is highly similar to the concatenation-based tree (Fig. 1).

The alignment-free phylogenetic tree based on the distance of nucleotide similarity, although keeping the monophyly of Rhodotorula and Rhodosporidiobolus, resulted in a somewhat different topology (Fig. 3). The monophyly of Sporobolomyces was not supported according to this outcome, and the clade including S. johnsonii and S. salmonicolor was separated from the clade of other Sporobolomyces.

Based on the predicted gene sets, we investigated the gene diversity of genera in Sporidiobolales. Initially, we utilized the presence-absence patterns of orthologs (PAPO) method (Takashima et al. 2019). Although the small sample size may have had an impact, the common orthogroup (OG) score of Rhodosporidiobolus was the highest among the three genera in Sporidiobolales (Fig. 4A). The common OG score of Rhodotorula and Sporobolomyces was relatively similar and lower than that of Rhodosporidiobolus, suggesting that these two genera may exhibit greater ortholog diversity.

Figure 4. 

A Ratio of conserved orthogroups (OGs) in each genus. The common OG score was calculated as the ratio of the number of OGs maintained in all members of a genus to the total number of OGs, as described by Takashima et al. (2019); B plot showing the pairwise average amino acid identity (AAI) similarity and the rate of shared OGs between strains. Open circles, open triangles, open squares, closed squares, and closed circles represent comparisons within Rhodotorula, within Rhodosporidiobolus, within Sporobolomyces, among different genera in Sporidiobolales, and among different orders in Microbotryomycetes, respectively. The black, green, blue, and yellow markers indicate comparisons within the same species, within the same genus but different species, within the same family but different genera, and among different orders, respectively.

Subsequently, we conducted average amino acid identity (AAI) analysis, following the methodology described by Liu et al. 2024. To minimize the influence of gene expansion in specific OGs, we employed the rate of shared OGs instead of shared genes as the y-axis. Similar to the previous study on Saccharomycetaceae (Liu et al. 2024), the distribution of AAI values more clearly differentiated taxonomic ranks compared to shared OGs (Fig. 4B, Suppl. material 1: table S3). In Sporidiobolales, the AAI values within the same species ranged from approximately 88% to 100%, but the values were calculated based on only 5 species with genome data from multiple strains, including R. toruloides (11 strains), R. graminis (2 strains), S. johnsonii (2 strains), S. pararoseus (2 strains), and S. salmonicolor (2 strains). The AAI values among different species, but within the same genus, ranged from approximately 65% to 88%, and the AAI values among different genera ranged from approximately 60% to 66%. The lower limit of interspecific AAI within Rhodotorula (green open circles in Fig. 4B) was lower than that within Rhodosporidiobolus (green open triangles in Fig. 4B) and Sporobolomyces (green open squares in Fig. 4B), suggesting that Rhodotorula exhibits the most variation in terms of conserved orthologs.

Since yeasts have species-specific physiological traits (Kurzman et al. 2010), we examined metabolic genes that contribute to key steps of energy source utilization. We compared gene catalogs derived from sequenced and reference genomes to orthologs annotated in the Kyoto Encyclopedia of Genes and Genomes (KEGG) and counted the number of genes found for each KEGG ortholog (Suppl. material 1: table S4). A clear correlation was observed between the numbers of metabolic genes related to sugar metabolism, nitrogen metabolism, and vitamin biosynthesis with reported physiological traits (Fig. 5). For example, species without the ability to assimilate sucrose and raffinose (Rhodotorula lusitaniae and Rhodotorula poonsookiae) lack beta-fructofuranosidase (K01193, INV/sacA), which hydrolyzes sucrose and raffinose into monosaccharides. Species lacking L-rhamnose 1-dehydrogenase (K18337, LRA1) are unable to assimilate L-rhamnose, with the exception of L. creatinivorum JCM 10699. R. mucilaginosa, S. pararoseus, and S. carnicolor, which cannot assimilate nitrate, lack both nitrate reductase (K10534, NR) and nitrite reductase (K17877, NIT-6). A clade of Rhodotorula, which includes species unable to grow on vitamin-free media (R. dairenensis, R. mucilaginosa, and R. sphaerocarpa), commonly miss the key enzyme for thiamine biosynthesis (K18278, THI5). Species that are unable to fully utilize glycerol as a carbon source (showing negative, weak, or latent growth phenotypes) have no or only one copy of glycerol 2-dehydrogenase (K18097, GCY1), while many strains that grow on media with only glycerol have multiple copies of this ortholog. S. salmonicolor, a species lacking beta-glucosidase (K05349, bglX), is unable to convert cellobiose.

Figure 5. 

Comparison between assimilation phenotypes and the numbers of metabolic enzymes in the genome for which a correlation is observed. The phenotypes are from Kurtzman et al. eds. (2011), Satoh and Makimura (2008), and Huang et al. (2011), and are indicated as positive (+), negative (-), variable (v), slow (s), latent (l), weak (w), and not available (n), in line with Kurtzman et al. (2011). The phylogenetic tree on the left is the consensus topology of Figs 1–3. Superscripts after strain numbers indicate ‘type strain’ (T), or ‘reference strain’ (*).

For several other carbon sources, however, no correlation was found between the number of genes linked to the associated metabolic pathway and the described phenotypes. This was the case for maltose, starch, galactose, sorbose, xylose, arabinose, ribose, mannitol, and ethanol (Fig. 6). Therefore, not all metabolic traits can be explained by an overlap with KEGG gene sets.

Figure 6. 

Comparison between assimilation phenotypes and the numbers of metabolic enzymes in the genome for which no clear correlation is observed. See legend of Fig. 5 for details.

Given the known ability of yeasts in the order Sporidiobolales to produce carotenoids such as gamma-carotene, beta-carotene, and torulene (Buzzini et al. 2007), we also examined the copy numbers of genes involved in carotenoid biosynthesis pathways. Carotenoid biosynthesis involves two enzymes: AL1/CAR1 (phytoene desaturase, equivalent to al-1, carB, crtI, K15745), which converts phytoene to 3,4-didehydrolycopene via lycopene, and AL2/CAR2 (bifunctional phytoene synthase/lycopene cyclase, equivalent to al-2, carRA, crtYB, K17841), which converts 3,4-didehydrolycopene to torulene and lycopene to beta-carotene via gamma-carotene (Suppl. material 1: fig. S2 upper panel; Kanehisa and Sato 2020). With the exception of the S. johnsonii JCM 5296 strain, which possesses two copies of the AL1/CAR1 gene, all strains of Sporidiobolales contain only one copy of each gene involved in carotenoid biosynthesis (Suppl. material 1: fig. S2, lower panel). Thus, differences in carotenoid quantity and composition may depend on the regulation of gene expression rather than the presence of the gene itself.

Several species, such as R. toruloides and R. glutinis, have been reported to contain high levels of lipids (Chattopadhyay et al. 2021). Therefore, we investigated the lipid synthesis-related genes present in their genomes. By analyzing the number of genes involved in the KEGG pathways for Fatty acid biosynthesis (PATH:ko00061), Fatty acid elongation (PATH:ko00062), and Biosynthesis of unsaturated fatty acids (PATH:ko01040), we found that ACSL/fadD (K01897), which encodes long-chain acyl-CoA synthetase, was slightly more abundant in R. toruloides (Suppl. material 1: table S5). However, no other notable differences in gene counts were observed.

To investigate potential marker genes that characterize taxa, we searched for KEGG Orthologs (KOs) that are unique to specific groups. To assess the relationship between gene presence and taxa, we calculated the correlation coefficient between binarized gene presence and taxa. Only KOs identified with an absolute correlation coefficient of 0.7 or higher were considered as relevant. Initially, we analyzed genes that exhibited a varying presence across three orders of Microbotryomycetes: Sporidiobolales, Leucosporidiales, and Microbotryales (Suppl. material 1: tables S6–S8). When comparing Sporidiobolales with the other two taxa, we identified 43 relevant KOs (Suppl. material 1: table S6), of which 26 were absent in all Sporidiobolales strains but present in all the other strains, suggesting they could be potential markers for distinguishing Sporidiobolales from other taxa. These include several oligosaccharide hydrolases such as alpha-glucosidase malZ (K01187) and beta-glucosidase E3.2.1.21 (K01188). These functions might be carried out by alternative orthologs like beta-fructofuranosidase INV/sacA (K01193) and beta-glucosidase bglX (K05349) in Sporidiobolales (Fig. 5). When comparing two Leucospodiales strains with the others, we found 24 KOs that were relevant (Suppl. material 1: table S7). It is noteworthy that both Leucosporidiales strains lack the AL1 (K15745) and AL2 (K17841) genes, which are important components of the carotenoid biosynthetic pathway. The absence of these genes in Leucosporidium can explain that Leucosporidium colonies have a white appearance while Sporidiobolales and Microbotryales colonies become red when carotenoids accumulate (Wang et al. 2015b). When comparing two Microbotryales strains against all the others, we identified 40 KOs (Suppl. material 1: table S8). Despite the limited genome data available for Microbotryales, these KOs could potentially serve as markers for distinguishing Microbotryales from other taxa.

Next, we searched for genes that are specifically present or absent in one of the three genera of Sporidiobolales. In the comparison between Rhodotorula and the other two genera, 14 KOs were identified (Suppl. material 1: table S9); ten KOs came out of the comparison between Rhodosporidiobolus and the other two genera (Suppl. material 1: table S10), and 25 KOs out of the comparison between Sporobolomyces and the other two genera (Suppl. material 1: table S11). Among KOs that differ by presence between all strains of one genus and all strains of the other two genera was ER membrane protein complex subunit 10 (K23570, EMC10), which was absent in Rhodotorula but present in Rhodosporidiobolus and Sporobolomyces (Suppl. material 1: table S9). Protein TBF1 (K22485) was found in Rhodosporidiobolus but not in Rhodotorula and Sporobolomyces (Suppl. material 1: table S10), and uric acid-xanthine permease (K23887, UAPA_C) was absent in Sporobolomyces but present in Rhodotorula and Rhodosporidiobolus (Suppl. material 1: table S11). These KOs can be markers for distinguishing genera among Sporidiobolales, while the functional importance of these proteins is unknown.

In addition to comparisons among existing taxa, we checked KOs whose presence differs between the two clades of Sporobolomyces. We found 83 relevant KOs, including 4 KOs that are absent in S. salmonicolor and S. johnsonii but present in all other clades of Sporobolomyces, and 21 KOs that are present in S. salmonicolor and S. johnsonii but not in other Sporobolomyces clades (Suppl. material 1: table S12).

We also conducted a search for genes that are indicative of the previous classification, which was based on the presence of ballistospores. Specifically, we examined the presence of certain KOs that exhibited significant differences between the ballistospore-forming group (previously classified as Sporobolomyces and Sporidiobolus) and the ballistospore non-forming group (previously classified as Rhodotorula and Rhodosporidium). Despite identifying multiple KOs with absolute correlation coefficients exceeding 0.7, we were unable to find any KO that distinctly distinguished between the two groups (Suppl. material 1: table S13).

Lastly, we conducted a comparative analysis between taxa by focusing on the number of gene components present in specific metabolic pathways or modules. As a result, it was reconfirmed that Leucosporidium completely lacks the components of the carotenoid biosynthesis pathway, map00906 (Suppl. material 1: table S14). Additionally, several correlations were observed between each taxonomic group and the number of components they possessed in various metabolic systems (Suppl. material 1: tables S14, S15). However, apart from the carotenoid biosynthesis pathway, no other metabolic systems were found to exhibit differences in multiple components significant enough to suggest substantial alterations in their metabolism. While a detailed examination of the differing genes may eventually lead to functional insights, this remains a task for future research.

The authors have declared that no competing interests exist.

No ethical statement was reported.

All the fungal strains used in this study have been legally obtained, respecting the Convention on Biological Diversity (Rio Convention).

This work was supported by the foundation of the Institute for Fermentation, Osaka (IFO, no. K-2020-006 to M.T.) and a Grant-in-Aid for Scientific Research from the MEXT (no. 20K06801 to M.T.). R.M. was supported by a Research Grant from the Ministry of Education, Culture, Sports, Science and Technology of Japan to the RIKEN Center for Integrative Medical Sciences. Genome sequencing was partially funded by the Genome Information Upgrading Program of the National BioResource Project of MEXT, with financial support provided to M.O. and R.M.

Y. Ko., N. T., and M. T. designed and coordinated the study. M. O. and M. T. isolated and cultured fungal strains. Y. Ku. isolated and M. M. sequenced genomic DNA. Y. Ko. and N. T. assembled and analyzed genomic data. Y. Ko., N. T. and M. T. drafted the manuscript. All authors contributed to the assessment of results, manuscript revisions and final approval for publication.

Yuuki Kobayashi https://orcid.org/0000-0003-1061-3639

Minenosuke Matsutani https://orcid.org/0000-0002-6262-3188

Keita Aoki https://orcid.org/0000-0003-2079-4031

Masako Takashima https://orcid.org/0000-0002-7686-8661

All sequence reads and assemblies were submitted to DDBJ/NCBI under BioProject accession numbers PRJDB3686 (JCM 10020), PRJDB3687 (JCM 10021), PRJDB3688 (JCM 10049), PRJDB3691 (JCM 24501), PRJDB3719 (JCM 10699), and PRJDB17412 (all other strains). The detailed accession numbers for sequence reads and assemblies were described in Suppl. material 1: table S16. Data for phylogenetic analyses were uploaded to FigShare (https://doi.org/10.6084/m9.figshare.25592013).