Selection of Orthologous Genes for Construction of a Highly Resolved Phylogenetic Tree and Clarification of the Phylogeny of Trichosporonales Species

The order Trichosporonales (Tremellomycotina, Basidiomycota) includes various species that have clinical, agricultural and biotechnological value. Thus, understanding why and how evolutionary diversification occurred within this order is extremely important. This study clarified the phylogenetic relationships among Tricosporonales species. To select genes suitable for phylogenetic analysis, we determined the draft genomes of 17 Trichosporonales species and extracted 30 protein-coding DNA sequences (CDSs) from genomic data. The CDS regions of Trichosporon asahii and T. faecale were identified by referring to mRNA sequence data since the intron positions of the respective genes differed from those of Cryptococcus neoformans (outgroup) and are not conserved within this order. A multiple alignment of the respective gene was first constructed using the CDSs of T. asahii, T. faecale and C. neoformans, and those of other species were added and aligned based on codons. The phylogenetic trees were constructed based on each gene and a concatenated alignment. Resolution of the maximum-likelihood trees estimated from the concatenated dataset based on both nucleotide (72,531) and amino acid (24,173) sequences were greater than in previous reports. In addition, we found that several genes, such as phosphatidylinositol 3-kinase TOR1 and glutamate synthase (NADH), had good resolution in this group (even when used alone). Our study proposes a set of genes suitable for constructing a phylogenetic tree with high resolution to examine evolutionary diversification in Trichosporonales. These can also be used for epidemiological and biogeographical studies, and may also serve as the basis for a comprehensive reclassification of pleomorphic fungi.


Introduction
The order Trichosporonales was proposed by Fell et al. [1] and, in "The Yeasts, A Taxonomic Study" 5 th ed., it includes the genera Trichosporon (37 species) [2] and Cryptotrichosporon (1 species) [3], as well as some species of highly polyphyletic genera, Bullera (3 of 41 species listed) [4] and Cryptococcus (10 of 70 species listed) [5]. Since the type species of the genera Bullera and Cryptococcus were assigned to the Tremellales, species within these two genera in the Trichosporonales are misleading taxonomically.
The genus Trichosporon, as it is defined at this time, is polyphyletic and divided into five clades; namely, Brassicae, Cutaneum, Gracile, Ovoides and Porosum. Each is traditionally characterized based on a combination of serological properties (serogroup) and the major ubiquinone (as the phenotype), as shown in Fig 1a. In the Cryptococcus species mentioned above, Weiss et al. [6] recently reported the amended concept of the genus Vanrija [7] by following several discussions on this group [2,3,5], and proposed the use of Vanrija as the genus name for Cryptococcus humicola and phylogenetically related species.
Clinically, biotechnologically and agriculturally important species are included in this order: approximately one-third of species in the genus Trichosporon are associated with infections or allergies in humans [2]; Trichosporon asahii (Ovoides clade) is the causative agent of trichosporonosis and summer-type hypersensitivity pneumonitis (SHP) [8]. Trichosporon mycotoxinivorans (Gracile clade) detoxifies mycotoxins, such as ochratoxin A and zearalenone [9], and Trichosporon porosum (Porosum clade) is thought to function in the global carbon cycle because it possesses CbhI genes for cellobiohydrolase [10] and oleaginous activity [11]. Additional oleaginous species have been reported among the Trichosporonales, including Cryptococcus curvatus [12], Vanrija albida [13], Vanrija musci [13], Trichosporon cacaoliposimilis (Gracile clade) [14] and Trichosporon oleaginosus (Cutaneum clade) [14]. Some species show cellulase activity, and most species can assimilate various compounds including aromatic and aliphatic compounds [2]. Since the expression (as shown above) and implied potential activities have been characterized in this order, understanding why and how this evolutionary diversification occurred is important. We have performed various studies on the genus Trichosporon and related species, including intra-and inter-species diversity; however, the phylogenetic relationships and evolutionary processes in this order remain unclear [2].
This study was performed to clarify the phylogenetic relationships among Trichosporonales species, especially basal species and clusters within this lineage. The draft genomes of 17 Trichosporonales species were determined, 30 orthologous protein-coding DNA sequences (CDSs) were selected, and phylogenetic trees were constructed based on the respective genes and a concatenated alignment. When selecting CDSs suitable for this analysis, the draft transcriptomes of T. asahii and T. faecale were used to identify coding regions for these species. Each CDS of T. asahii and T. faecale was aligned with that of Cryptococcus neoformans (outgroup), and those of other species were added and aligned based on codons, because the intron position of the respective genes is not conserved in this order.

Gene selection
A summary of the draft genome sequences is shown in Table 1. The obtained genomic size of T. asahii was similar to other recently published data (ALBS01000000 and AMBO01000000). The range of genomic sizes of Trichosporonales species, estimated based on the total length of contigs in this study, ranged from 18.2 Mb (C. curvatus) to 42.2 Mb (Trichosporon coremiiforme). We did not observe a trend between genome size and phylogenetic position or morphology in this study. by a grant of the Innovative Cell Biology by Innovative Technology (the Cell Innovation Program) from MEXT, Japan to Yoshihide Hayashizaki. In silico biology, inc. provided support in the form of salaries for author AO but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The specific roles of this author is articulated in the 'author contributions' section.' Competing Interests: Akira Ohyama is the CEO of In silico biology, inc. There are no patents, products in development or marketed products to declare. This does not alter the authors' adherence to all the PLOS ONE policies on sharing data and materials, as detailed online in the guide for authors.
The candidate genes for phylogenetic analysis selected in this study are shown in Table 2. These 30 genes were selected as described in the Materials and Methods without consideration of function. They are orthologous genes shown in the Microbial Genome Database for Comparative Analysis (MBGD) (http://mbgd.nibb.ac.jp). Notably, we found that the genomes of T. coremiiforme and T. ovoides occasionally include two paralogs for some of the studied genes, whereas those of other species contained only one ( Table 2). Even including the two paralogs from T. ovoides and T. coremiiforme, all phylogenetic trees based on a single gene showed that the former constituted a subclade with T. inkin and the latter with T. asahii and T. faecale (data not shown). This topology was the same as reported previously [2], indicating that the paralogous genes were not transferred from phylogenetically distant species. Therefore, these two species were excluded from the robustness analyses of the tree topology.

Phylogeny of Trichosporonales
First, we obtained sequence data for three regions (RNA polymerase II subunits I (RPB1) and II (RPB2) and translation elongation factor 1-α (TEF1-α)) that have been used for the phylogenetic analyses in Ascomycotina and Basidiomycotina [15][16][17][18][19]. The sequences of these regions  Phylogeny of Trichosporonales of 12 ascomycetous yeasts (S1 Table) were obtained from the database and aligned to extract the corresponding regions from our data. The conserved domains of genes predicted by Microbial Genome Annotation Pipeline (MiGAP) (http://www.migap.org/) showed good alignment based on the codons, and we found that Trichosporonales species possessed spliceosomal introns and that their position was strain-specific. Intron positions also differed from those in C. neoformans, which was used as a model species for gene prediction using HMM-based gene finder AUGUSTUS [20] (Fig 2 and S2 Table). While the intron positions of Agaricomycetes species are sufficiently conserved to distinguish and characterize the higher taxa [16], those of Trichosporonales species are not conserved. Therefore, we identified and aligned CDS regions of T. asahii and T. faecale by referring to their mRNA data and then added genes from other species, as indicated in the Materials and Methods. The sequences used in this study covered almost the complete CDS region; 72,531 nucleotides (24,173 residues) among 30 concatenated genes. Gene alignments are available upon request.
Phylogenetic trees for the 16 species based on nucleotide and amino acid sequences using the 30 concatenated genes are shown in Fig 1b and 1c, respectively. The trees based on a single gene are presented in (S3 Table). Trees inferred from the concatenated alignments based on nucleotide (Fig 1b) and amino acid (Fig 1c) sequences had the same topology and increased reliability compared with the tree based on the D1/D2 region of the LSU rRNA gene (Fig 1a) and other single-gene trees. A phylogenetic tree inferred from the amino acid sequence alignment including consensus sequences of T. coremiiforme and T. voides (see Materials and Methods) showed an identical topology (S1 Fig). Thus, they currently represent the best estimate of the phylogenetic relationship in Trichosporonales.
To evaluate the robustness of the topology of the tree inferred from the concatenated alignments, we identified all nodes recovered during the single-gene analyses and recorded bootstrap values when the nodes were present in a tree inferred from the concatenated alignment (Table 3). Although T. ovoides should be included in the phylogenetic analyses of the Trichosporon genus since this is the type species of the genus, we excluded this species from the analysis, as mentioned above. Nonetheless, we continue to use the name of the clade, Ovoides, in this report.
In  Table 3, and hereafter]. Regarding the Cutaneum clade, the two species examined (T. cutaneum and T. dermatis) formed a cluster in trees inferred from the concatenated alignment, which was also observed in single-gene trees [nodes c (27/30) and p (26/30)]. The Ovoides and Cutaneum clades formed a monophyletic group, which was supported with 100% and 72% bootstrap values for nucleotide and amino acid sequences, respectively. The tree based on the D1/D2 region of the LSU rRNA gene did not support this cluster (Fig 1a), and single-gene trees showing more than 70% bootstrap values were 4 of 30 [nodes e] and 1 of 30 [node r] based on nucleotide and amino acid sequences, respectively.
The Porosum and Gracile clades formed a cluster, and the connection of the Brassicae clade was shown in the tree based on the D1/D2 region of the LSU rRNA gene (Fig 1a). In trees inferred from the concatenated alignment, the Brassicae and Gracile clades were clustered and connected to the Porosum clade with 100% bootstrap values for both nucleotide and amino acid sequences. The bootstrap values for single-gene trees were similar to those for trees inferred from the concatenated alignment: two-thirds supported the connection of the Brassicae and Gracile clades [nodes i (19/30) and v (20/30)], and one-third supported their connection to the Porosum clade [nodes j (13/30) and w (9/30)].
Cryptococcus curvatus is in an adjunct position in the Cutaneum clade with 100% bootstrap values for both nucleotide and amino acid sequences (Fig 1b and 1c). We previously evaluated the phylogenetic position of the species based on rRNA genes and spacer regions, but were    Table 3.     Fig 1b and 1c. Numerals show the bootstrap value of each node which corresponded to that of the concatenated tree [boldface, more than 70%; (), less than 70%].-indicates that the node on the concatenated tree is not shown in respective CDS tree. * 2 , Nodes showing more than 70% bootstrap value (boldface) were treated as "positive". * 3 , The number of base substitutions per site from averaging over all sequence pairs are shown. Analyses were conducted using the Kimura 2-parameter model [47]. All ambiguous positions were removed for each sequence pair.
* 4 , The number of amino acid substitutions per site from averaging over all sequence pairs are shown. Analyses were conducted using the JTT matrix-based model [45]. All ambiguous positions were removed for each sequence pair. doi:10.1371/journal.pone.0131217.t003 unable to obtain good resolution [21][22]. In this study, the single-gene trees did not frequently support the node [more than 70% bootstrap support nodes for d (6/30) and q (7/30)]. Vanrija humicola was located outside of the Trichosporon species, which has been reported previously [2,5,6,21,23]. As described above, Weiss et al. [6] proposed the use of genus Vanrijia as the genus name, and we follow their suggestion in this report. The genus Vanrija sensu Moore [7] is polyphyletic; thus the emendation and subsequent reclassification including new combinations is required. We believe that our genes will contribute to the reclassification; namely, delineation for the species at the basal position between Vanrija and other genera.
The Bullera koratensis and Cryptococcus tepidarius clades occurred at remote positions in the Trichosporonales. In addition, the intron positions were closer to those of C. neoformans than to those of other Trichosporon species (S2 Table). Okoli et al. [3] proposed the genus Cryptotrichosporon based on the D1/D2 region of the LSU rRNA plus ITS regions and the 18S rRNA gene, and noted that B. koratensis and C. tepidarius were phylogenetically closely related to their species. These species might be pivotal in clarifying the phylogenetic relationships between Tremellales and Trichosporonales, since they obtain intermediate positions.

Phenotypic properties and phylogeny
Trichosporon has been divided into five clades based on a combination of serological properties and major ubiquinone as described above [2,24]. Based on morphological aspects, species in the clades Brassicae, Cutaneum, Gracile, and Ovoides produce abundant arthroconidia, whereas those in the Porosum clade rarely produce arthroconidia or not at all [2]. Fig 1b and  1c showed clustering of the Brassicae and Gracile clades and a connection to the Porosum clade with 100% bootstrap values for both nucleotide and amino acid sequences, indicating that this morphological characteristic supports the linkage between Brassicae and Gracile clades. Recently, Nagy et al. [25] reported that Zn-finger transcription factor families are related to switches between yeast and filamentous forms. Since the production of arthroconidium has been reported in the genera Guehomyces and Tausonia (Cystofilobasidiales) and species in the family Dipodascaceae (Saccharomycetales), genes associated with arthroconidium production may be found when the genomic data are accumulated.
Regarding serological properties, five clades; namely, Cutaneum, Ovoides, Brassicae, Gracile and Porosum, are characterized by serogroup I, II, III, I-III and I-III (reacting with both factors I and III), respectively [26]. The clustering topology of clades Brassicae (serogroup I) and Gracile (I-III), and the connection to clade Porosum (I-III), indicated that these three clades can be characterized by positive reactions with III as common antisera (S2 Fig). This grouping, serogroups I, II and III, was proposed previously [27] and confirmed in this study. The epitopes of serological group II and III have been reported [28][29]; thus, this characteristic would be a useful property for differentiation.
For the major ubiquinone in ascomycetous yeasts, ubiquinone isoprene chain length ranged from 6 to 10 (namely Q-6 to Q-10) and showed a strong relationship with phylogeny of higher taxa with some exceptions, and was used as a phenotypic characteristic [15,30]. In basidiomycetous yeasts, most species contain Q-10 and few genera are characterized by different chain lengths: Mrakia and Cystofilobasidium (both in Cystofilobasidiales) are characterized by Q-8, and Q-9 was observed in Trichosporon (excluding the Cutaneum clade), Vanrija, Guehomyces, and Itersonilia in Agaricomycotina. In this study, Q-10 species C. curvatus clustered with the Cutaneum clade, which has been characterized by Q-10 (S2 Fig), and this characteristic supports this linkage. Although the taxonomic significance of ubiquinone isoprene chain length remains unclear due to the occurrence of both Q-9 and Q-10 in the same genus; e.g., in Leucosporidium [31], this could be used as a distinctive characteristic with serological properties when the Cutaneum clade together with C. curvatus is proposed as a separate genus. We believe that further studies will clarify how the isoprene chain length is controlled by prenyltransferase and related mechanisms.
Clarification of the topology of basal species in a lineage by obtaining a highly resolved tree Recently, the number of fungal phylogenetic trees based on genomic data have increased due to the greater amount of genomic data available e.g., [32][33][34][35][36]. However, studies on fungi are still fewer compared with bacteria and archaea, mainly because of the amount of available genomic data. Based on the NCBI BioProject (http://www.ncbi.nlm.nih.gov/bioproject/), only 84 and 230 species were available for yeasts and filamentous fungi, respectively, as of May 15th, 2014. With the available fungal genomic data, Rokas et al. [32] demonstrated that a minimum of 20 concatenated genes in an analysis were sufficient to represent the phylogeny based on all genomic data for Saccharomyces species. Kuramae et al. [33] showed that 40-45 concatenated proteins would be required to resolve the tree of life of the fungal kingdom by cophenetic correlation analysis. Based on these reports, we believe that concatenated analyses of genes could be used for robust studies on the fungal tree of life, although the number of genes required would depend on the range of the target groups.
Rokas et al. [32] used 106 genes and compared a total length of 127,026 nucleotides (42,342 residues), which corresponded to roughly 1% of the genomic sequences and 2% of the predicted genes. According to this estimation, the use of "a concatenated analysis of 20 or more genes", corresponding to roughly 0.2% of the genomic sequences, is sufficient to construct an accurate phylogeny. Our data based on 72,531 nucleotides and corresponding to 0.24% (for Trichosporon laibachii, genome size 30.6 Mb) to 0.40% (for C. curvatus, genome size 18.2 Mb) of the total genomic sequences of Trichosporonales support the results of Rokas et al. [32].
While the number of genes required to represent a tree with a resolution comparable to the resolution of a genome tree has been discussed [32,34], an example of such a set of genes has not been reported. In this study, we propose a set of genes suitable for constructing a phylogenetic tree with high resolution to examine the evolutionary diversification in Trichosporonales. Table 3 indicates that the use of more genes does not improve the resolution of the tree of this taxonomic group. In addition, among the nucleotide sequences used in this study, those for four genes (phosphatidylinositol 3-kinase TOR1, glutamate synthase (NADH), malate synthase and phosphoenolpyruvate carboxykinase) showed better resolution than those typically used for multigene analyses, such as TEF1-α, RPB1 and RPB2. For the tree based on amino acid sequences of phosphatidylinositol-3-kinase TOR1 (even if used alone), the same topology was obtained as for the tree inferred from the concatenated alignment. The tree inferred from the concatenated alignment based on RPB1, RPB2 and TEF1-α (regions corresponding to those used by Kurtzman and Robnett [15]) revealed the existence of five clades in Trichosporon (nodes b, c, f, g and h) with 100% bootstrap values and a clear relationship among the Brassicae, Gracile and Porosum clades (nodes i and j; S3 Fig). This tree had better resolution than the respective single-gene trees, as shown in (S4 Table), but did not clarify the relationship between the Ovoides clade and the other four Trichosporon clades (nodes e and r), or the positions of V. humicola (nodes k and x) and C. curvatus (node d). Obviously, the use of these three genes is not enough to clarify the phylogenetic relationship of the Trichosporonales. This was our original motivation for the present study.
In yeasts, resolution of the phylogenetic relationships has particularly been increased for the ascomycetous teleomorphic genera using partial sequences of TEF1-α and cytochrome oxidase II, in addition to rRNA genes. The later addition of RPB1 and RPB2 improved the resolution for the Saccharomycetales yeast species [15]. Recently, Koufopanou et al. [37] reported primer sets for 14 genes to clarify the topology of basal species among the Saccharomycetales. Of our 30 genes, glutamate synthase and adenosyl homocysteinase were used for multiple gene analysis by Koufopanou et al. [37], and CDC19 was selected as a candidate during the first screening stage but was not used in the published analysis. Gene duplication in the Tor pathway was also reported in some fungi by Shertz et al. [38]. Careful attention must be paid when conducting multigene analyses in other taxa using the Tor gene; however, this gene can be used for Trichosporonales, since we did not detect a Tor paralog in our experiments.

Conclusions
In this study, phylogenetic trees inferred from concatenated alignments (nucleotide and amino acid sequences) of 30 genes extracted from draft genome data showed greater resolution than in previous reports, indicating that they currently represent the best estimate of the phylogenetic relationship in Trichosporonales. In addition, we suggest that the genes identified in this study are good candidates for constructing a multi-locus sequence analysis system for both inter-and intra-species analyses. Specific primers are required to adopt these newly discovered promising genes. These genes can also be used for epidemiological and biogeographical studies. A set of genes suitable for constructing a phylogenetic tree to examine evolutionary diversification. The relationship between the closely related Trichosporonales and Tremellales in the Agaricomycotina might be clarified using these genes [6,17,39]. Moreover, our study provides a platform for future studies on a comprehensive reclassification of pleomorphic fungi.

Genome sequencing
Strains and DNA and RNA extraction. Strains used in this study were received from the Japan Collection of Microorganisms at the RIKEN Bioresource Center (Table 1). Cultures grown on YM agar (Difco, Detroit, MI; BD, Franklin Lakes, NJ, USA) were collected and freeze-dried, and DNA was prepared according to a previous method [40]. The integrity and quality of the extracted genomic DNA were examined by agarose gel electrophoresis using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The concentration of DNA was determined using a Quant-iT PicoGreen dsDNA Assay Kit (Life Technologies, Gaithersburg, MD, USA). Total RNA was prepared from T. asahii and T. faecale cells cultivated at 37°C or 28°C using an RNeasy Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions after disruption by glass beads. The integrity of the RNA was verified by agarose gel electrophoresis. The RNA concentration was determined using a Nano-Drop 2000 spectrophotometer.
Genome sequencing. Draft genomes of 17 selected species in the Trichosporonales were determined using a 454 GS FLX (Roche Diagnostics, Mannheim, Germany) and/or Illumina (Illumina Inc., San Diego, CA, USA) next-generation sequencer. Since karyotyping data is not available for the Trichosporonales, we were unable to estimate the genome sizes of these species before the beginning of this study. Thus, the genome size of C. neoformans (approximately 20 Mb) was used as an initial estimate.
DNA library preparation and sequencing were performed according to the manufacturer's instructions. Briefly, a paired-end DNA library was prepared from 1 μg genomic DNA, following DNA fragmentation into approximately 350-bp average sized fragments using an S2 ultrasonicator (Covaris, Woburn, MA, USA) and then sequenced using a HiSeq 2000 sequencer (Illumina Inc.; S5 Table). The T. asahii genome was also sequenced using the GS FLX Titanium platform (Roche Diagnostics). A DNA library was prepared from 0.6 μg T. asahii genomic DNA, followed by fragmentation into approximately 900-bp average sized fragments using a Covaris S2 ultrasonicator and then sequencing using GS FLX pyrosequencing technologies. De novo assembly was performed using the CLC Genomics Workbench (version 4.5 or later) and CLC Genomics Server (version 3.6 or later) (CLC Bio Co., Aarhus, Denmark) with the default settings, including the "Update contigs" option, which edited contig sequences based on the reads mapped to the contig. The de novo assembly produced contigs with an average N50 value of 37,614. Key attributes of the draft genome sequences are summarized in Table 1. Library generation, sequencing and de novo assembly were performed at the Genome Network Analysis Support Facility, RIKEN CLST (Yokohama, Japan).
Transcriptomic sequencing. An RNA-Seq library was prepared as described previously [41]. Briefly, polyA(+) RNA was extracted from 10 μg of total RNA using Dynabeads Oligo (dT)25 (Life Technologies), according to the manufacturer's protocols, and was fragmented by heating at 70°C for 3.5 min in 0.5× fragmentation buffer (Life Technologies). Fragmented RNA was purified using the RNeasy MinElute Kit (Qiagen) and then dephosphorylated using Antarctic phosphatase (New England Biolabs, Beverly, MA, USA). The dephosphorylated RNA was phosphorylated using T4 polynucleotide kinase (New England Biolabs) with 10 mM ATP and was purified again using the RNeasy MinElute Kit. Purified RNA was concentrated and ligated at both ends with 3' and 5' adapters from the TruSeq Small RNA Sample Prep Kit (Illumina Inc.) using T4 RNA ligase 2 and T4 RNA ligase 1 (New England Biolabs), respectively. Adapted RNA was reverse transcribed using PrimeScript Reverse Transcriptase (Takara Bio, Kyoto, Japan) and amplified using Phusion High-Fidelity DNA Polymerase (New England Biolabs), as well as forward and reverse primers from the TruSeq Small RNA Sample Prep Kit. cDNA libraries with a size range of 160-320 bp were selected by removing PCR primers using 1.2 volumes of AMPure XP beads (Beckman Coulter, Miami, FL, USA) and subsequent gel extraction using Pippin Prep (Sage Science, Beverly, MA, USA). The four resulting cDNA libraries were pooled and then subjected to 101-bp single-read sequencing on a HiSeq 2000 sequencer. The sequencing produced 97 million and 88 million reads for T. asahii cultivated at 37°C and 28°C, respectively, and 80 million and 117 million reads for T. faecale cultivated at 37°C and 28°C, respectively. The reads were trimmed of adapter, low-quality and ambiguous sequences using the CLC Genomics Workbench with default parameters and were then assembled de novo with the CLC Genomics Workbench and CLC Genomics Server. The draft transcriptome assemblies of T. asahii at 37°C and 28°C and of T. faecale at 37°C and 28°C consisted of 26,569 (N50 = 1,307), 26,769 (N50 = 1,382), 24,105 (N50 = 1,362) and 25,349 (N50 = 1,389) transcripts, respectively. Library generation, sequencing, and de novo assembly were performed at the Genome Network Analysis Support Facility, RIKEN CLST.
Auto-annotation. After the de novo assembly, auto-annotation was performed using MiGAP (http://www.migap.org/), a publicly available annotation pipeline system of microbial and fungal genome draft sequences operated by the National Institute of Genetics. Briefly, MiGAP provides consecutive gene searching and annotation functions and consists of prediction programs. Exons and introns were predicted using the HMM software AUGUSTUS 2.5.5 [20]. The HMM model species used for this prediction is cryptococcus_neoformans_neofor-mans_JEC21. After gene prediction, each ORF was subjected to a homology search against the following protein sequence databases using NCBI BLAST: the fungi division of RefSeq, whole entries of TrEMBL and whole entries of NR.

Alignment and phylogenetic analyses
Gene selection and alignment. For selection of genes, we used genomic data from C. neoformans (AE017341) and C. gattii (CP000286) as references. Since our target genes were those that can be used for genus-level discrimination, they were required to be conserved in the same genus. A total of 1,234 CDSs were selected from C. neoformans, C. gattii and T. asahii genomic data under the conditions of identity = 60 and coverage = 60, using in-silico molecular cloning (IMC; In Silico Biology, Inc., Yokohama, Japan). CDS homologs showing more than 97% similarity between C. neoformans and C. gattii were listed. Candidate T. asahii genes were selected using IMC under the criteria that they include approximately 300 or more amino acids and do not include many introns, since our purpose was to develop signature/marker sequences. CDS homologs of other species were then selected using IMC. Genes used in this study are listed in Table 2 with DDBJ/GenBank/EMBL accession numbers.
The coding regions of genes from T. asahii and T. faecale were determined using a BLAST search against the respective mRNA data. Multiple alignment of respective genes was first performed using the CDSs of T. asahii and T. faecale using those of C. neoformans as an outgroup. Subsequently those of other species were added and aligned based on codons. In some cases, part of a gene could be found in one contig and the other parts in different contigs. In this case, manual alignment was performed.
When making an alignment file, the following treatment was performed for two CDSs beforehand: An obvious insertion (1,557 nucleotides) in CND03540 (DNA-dependent RNA polymerase II (RPB2)) homolog was found in T. asahii but not in other species, and the region was eliminated from the analyses. Since the same was observed in the published T. asahii genome data (for both ALBS01000000 and AMBO01000000), this insertion was assumed to be distinct to this species or these strains in Trichosporonales. For CNF03740 (phosphatidylinositol-3-kinase TOR1) homolog, the length of CDSs of 9 species studied, viz. Trichosporon brassicae, T. cutaneum, T. dermatis, Trichosporon scarabaeorum, B. koratensis, C. curvatus, C. tepidarius, V. humicola and C. neofomans (outgroup) were 7,167 nucleotides (2,388 amino acids), whereas those of other species were found to be longer and varied. Since they aligned well from the 5' end to position 7,167 by codon, 7,167 nucleotides (2,388 amino acids) were used for the analyses; the downstream sequences were not used. Since CDSs were orthologous genes shown in the Microbial Genome Database for Comparative Analysis (MBGD) (http:// mbgd.nibb.ac.jp), CDSs were concatenated using SeaView [42].
Alignment including the consensus amino acid sequence of T. coremiiforme and T. ovoides. As shown in Table 2, we found paralogous genes in T. ovoides and T. coremiiforme. In this case, we created a consensus sequence from paralogs (from a multiple alignment) with X indicating a non-identical amino acid. When a partial sequence of gene was found only in one of the two sequences, this partial sequence was included in the consensus sequence.
Alignment of frequently used regions of RPB1, RPB2 and TEF1-α in fungi. Sequence alignments of RPB1, RPB2, and TEF1-α that are often used for the phylogenetic analyses of fungi were generated by referring to sequences of 12 ascomycetous yeasts [15]. The reference sequences used for this extraction are shown in (S1 Table).
Alignment of the D1/D2 region of LSU rRNA gene. All sequences of the D1/D2 region of LSU rRNA gene were obtained from GenBank/EMBL/DDBJ database (S1 Table). A multiple alignment was performed using MEGA6 [43] and edited manually.
Phylogenetic analyses. Trees (a phylogenetic tree based on each gene and a tree inferred from the concatenated alignment) were constructed using the maximum likelihood method in MEGA6 [43]. When a tree was constructed based on nucleotide sequences, the Tamura-Nei model [44] was used for the analyses with options Rates and Patterns = uniform rates. Initial trees for the heuristic search were obtained automatically by applying the neighbor-joining and BIONJ algorithms to a matrix of pairwise distances estimated using the maximum composite likelihood (MCL) approach and then selecting the topology with the highest log likelihood value. When a tree was inferred based on amino acid sequences, the JTT matrix-based model [45] was used for analyses with option Rates and Patters = uniform dates. Initial trees for the heuristic search were obtained automatically by applying the neighbor-joining and BIONJ algorithms to a matrix of pairwise distances estimated using a JTT model and then selecting the topology with the highest log likelihood value. All sites were used in all analysis with 100 times bootstrap analyses [46].
Average base substitutions per site and average amino acid substitutions per site for each gene were conducted using the Kimura 2-parameter model [47] and the JTT matrix-based model [45], respectively.
Supporting Information S1 Table. Sequence accession numbers retrieved from the GenBank/EMBL/DDBJ databases in this study. (XLSX) S2 Table. Intron distribution in translation elongation factor 1-α (TEF1-α), DNA-directed RNA polymerase ii largest subunit (RPB1) and DNA-dependent RNA polymerase II (RPB2). (XLSX) S3 Table. Phylogenetic trees based on amino acid and nucleotide sequences. Table. Comparison of topology between the 30 gene concatenated tree and single and concatenated tree of the RPB1, RPB2 and TEF1-α genes. (XLSX) S5 Table. Basic statistics of whole-genome shotgun sequencing. (XLSX) S1 Fig. Phylogenetic tree inferred from the amino acids sequences including consensus sequence of T. coremiiforme and T. ovoides. The amino acid consensus sequences for T. coremiiforme and T. ovoides were created (see Materials and Methods) and labeled "consensus" in parenthesis following the species name. The tree was constructed using the maximum likelihood method based on a JTT matrix-based model [45]. Numerals on each node represent the percentages from 100 replicating bootstrap samplings [46]. There were a total of 24,173 amino acids in the dataset, and the tree with the highest log likelihood (-219231.8349) is shown. B., Bullera; C., Cryptococcus; T., Trichosporon; V., Vanrija. The trees in (a) and (b) were constructed using the maximum likelihood method based on the Tamura-Nei model [44] for (a) and a JTT matrix-based model [45] for (b). Numerals on each node represent the percentages from 100 replicating bootstrap samplings (frequencies of less than 60% are not shown) [46]. Letters in brackets after bootstrap values indicate the node name (see Table 3