Genome-scale phylogeny and contrasting modes of genome evolution in the fungal phylum Ascomycota

To infer a genome-scale phylogeny of the Ascomycota, we used 1107 publicly available genomes from species belonging to Ascomycota (Saccharomycotina, 332; Pezizomycotina, 761; and Taphrinomycotina, 14) and six out-groups from the sister fungal phylum Basidiomycota. The numbers of species examined represent ~332 of 1000 (~30%) of all described species of Saccharomycotina, ~761 of 82,000 (~1%) of Pezizomycotina, and ~14 of 140 (~10%) of Taphrinomycotina. All genomes were retrieved from the National Center for Biotechnology Information (NCBI) GenBank database, ensuring that only one genome per species was included (tables S1 and S2). Analysis of genome assembly completeness reveals that 1021 of 1113 (~92%) genomes have more than 90% of the 1315 full-length Benchmarking Universal Single-Copy Orthologs (BUSCO) genes (fig. S1) (31).

The 1315 BUSCO genes from 1107 Ascomycota species and six out-groups were used to construct a phylogenomic data matrix (see Materials and Methods). After constructing the multiple amino acid sequence alignment and trimming ambiguous regions for each of these 1315 BUSCO genes, we kept only the 815 BUSCO genes that had taxon occupancy of ≥50% for each subphylum (i.e., ≥7 Taphrinomycotina, ≥166 Saccharomycotina, and ≥381 Pezizomycotina) and whose amino acid sequence alignments were ≥300 sites in length. In the final set of 815 BUSCO genes, alignment lengths range from 300 to 4585 amino acid sites (average, 690), and numbers of taxa range from 851 to 1098 (average, 1051) (table S3). The final data matrix contains 1107 taxa, 815 genes, and 562,376 amino acid sites.

Inference using concatenation- and coalescent-based approaches yielded a robust, comprehensive phylogeny of the Ascomycota phylum (Fig. 1). The vast majority of internodes in both the concatenation-based (1103 of 1110, 99%) and the coalescent-based phylogeny (1076 of 1110, 97%) received strong (≥95%) support and were congruent between the phylogenies inferred using the two approaches; only 46 of 1110 (4%) internodes were incongruent between the two phylogenies (figs. S2 and S3).

Our higher-level phylogeny of Ascomycota is generally more congruent with previous genome-scale phylogenies (2, 17, 18) than with few-genes-from-many-taxa phylogenies (12–15), particularly with respect to relationships among the nine classes in the subphylum Pezizomycotina. For example, genome-scale studies, including ours, consistently favor a clade consisting of Pezizomycetes and Orbiliomycetes as the sister group to the rest of the Pezizomycotina (17, 18), while studies based on a few genes recovered either Orbiliomycetes (14–16) or Pezizomycetes (13) as the sister class to the rest of the Pezizomycotina (Fig. 2A). Our phylogeny also strongly supported the placement of the class Schizosaccharomycetes, which includes the model organism S. pombe, as the sister group to the class Pneumocystidomycetes, which contains the human pathogen Pneumocystis jirovecii (Fig. 2B). A recent genome-scale study of 84 fungal genomes showed that our result is consistent with the phylogeny inferred using an alignment-free composition vector approach but not with the phylogeny inferred using maximum likelihood (ML), which instead recovered Schizosaccharomycetes as the sister group to Taphrinomycetes (17). Last, both concatenation- and coalescent-based approaches supported the placement of the subphylum Saccharomycotina as the sister group to the subphylum Pezizomycotina (Fig. 2C). This result is consistent with most previous studies that analyze multiple sequence alignment data (15, 17, 19, 32, 33), but not with a recent study that analyzed genomic data with an alignment-free method and placed the subphylum Saccharomycotina as the sister group to the subphylum Taphrinomycotina (34).

To evaluate whether our genome-scale data matrix robustly resolved the three historically contentious branches discussed in the previous paragraph, we quantified the distribution of phylogenetic signal for alternative topologies of these three phylogenetic hypotheses at the level of genes and sites using an ML framework presented by Shen et al. (35). First, we found that the phylogenetic support for each of the three branches stemmed from many genes, i.e., it was not dominated by a small number of genes with highly disproportionate influence (table S4). Second, we found that the topology recovered by both concatenation- and coalescent-based approaches in our study had significantly higher frequencies of supporting genes and supporting sites (G-test) than two alternative competing hypotheses, ranging from 0.45 to 0.65, in all three branches examined (Fig. 2, A to C, and tables S4 and S5). Neither of the two alternative conflicting phylogenetic hypotheses for each of the three branches received frequencies of supporting genes and supporting sites that were equal to or greater than one-third (0.33), the value expected if the relationships among the taxa were represented by a polytomy. In summary, we found that very few branches (<5%) were incongruent between concatenation- and coalescent-based phylogenies (figs. S2 and S3). We also found that individual genes and sites exhibited robust phylogenetic signal for specific, historically contentious branches (Fig. 2, A to C). These results suggest that genome-scale amounts of data from a comprehensive sampling of taxa will robustly resolve major lineages of the tree of life.

We next used the robust phylogeny, a relaxed molecular clock approach, and six widely accepted time calibration nodes (see Materials and Methods) to infer the time scale of evolution of Ascomycota. We inferred the origin of the phylum to have taken place 563 Ma ago (95% CI, 631 to 495 Ma ago), the origin of the subphylum Saccharomycotina 438.4 Ma ago (CI, 590 to 304 Ma ago), the origin of the subphylum Pezizomycotina 407.7 Ma ago (CI, 631 to 405 Ma ago), and the origin of Taphrinomycotina crown group 530.5 Ma ago (CI, 620 to 417 Ma ago). Notably, the taxonomic placement of all budding yeasts into a single class, Saccharomycetes, whose origin coincides with the origin of the subphylum Saccharomycotina, means that the last common ancestor of this sole class of budding yeasts is much more ancient than those of any of the nine classes (based on current taxon sampling) in the subphylum Pezizomycotina (fig. S4 and table S6). For example, the most ancient class in Pezizomycotina is Pezizomycetes, whose origin is dated 247.7 Ma ago (CI, 475 to 193 Ma ago) (fig. S4 and table S6). The other outlier, albeit with much larger confidence intervals, is class Neolectomycetes in Taphrinomycotina, which we estimate to have originated 480.4 Ma ago (CI, 607 to 191 Ma ago) (fig. S4 and table S6).

Comparison of our inferred dates of divergence to recent studies using a few genes (11, 13, 14) shows that our estimates are generally younger. For example, our estimate for the phylum Ascomycota is around 563 Ma ago, whereas the origin of the Ascomycota was dated at least 533 Ma ago by Hyde et al. (9). However, recent studies that used 18S and 28S ribosomal RNA markers from 111 fungal taxa (11) and four gene markers from 145 Ascomycota taxa (13) inferred the origin of the Ascomycota to have taken place around 606 and 588 Ma ago, respectively. This result is consistent with findings of previous studies [e.g., (36)], where inclusion of large numbers of genes was found to also result in younger estimates of divergence times, perhaps because of the influence of larger amounts of data in decreasing the stochastic error involved in date estimation. In summary, generation of a genome-scale timetree for more than 1000 ascomycete species spanning the diversity of the phylum provides a robust temporal framework for understanding and exploring the origin and diversity of Ascomycota lifestyles.

To begin to understand the similarities and differences in the modes of genome evolution between subphyla, we assessed seven different genomic properties for each species. Because the number of genomes from Taphrinomycotina taxa is substantially smaller (14) than Saccharomycotina (332 genomes) and Pezizomycotina (761) (Fig. 3), we focused on examining the evolution of different genomic properties between Saccharomycotina and Pezizomycotina. Specifically, we found that Saccharomycotina exhibited a 1.6-fold higher evolutionary rate (on average, 1.80 substitutions per site in Saccharomycotina versus 1.12 substitutions per site in Pezizomycotina), 1.24-fold lower guanine-cytosine (GC) content (40% versus 50%), 3-fold smaller genome size (13 Mb versus 39 Mb), 1.9-fold lower number of protein-coding genes (5734 versus 10,847), 1.3-fold lower number of DNA repair genes (41 versus 54), 1.2-fold higher number of tRNA genes (179 versus 146), and 1.3-fold smaller estimates of nonsynonymous to synonymous substitution rate ratio (d_N/d_S) (0.053 versus 0.063), compared with Pezizomycotina (Fig. 3A, Table 1, and table S7). In addition, we found that gene density (that is, the number of genes divided by genome size) in Saccharomycotina is much higher than that in Pezizomycotina (on average, 451 genes per megabase in Saccharomycotina versus 290 genes per megabase in Pezizomycotina), which is in agreement with previous work showing that gene density is inversely correlated with genome size (37).

Analysis of standard Pearson’s correlations among the seven genomic properties revealed that two pairs exhibited statistically significant contrasting patterns between Saccharomycotina and Pezizomycotina. Specifically, evolutionary rate shows negative correlation with GC content in Saccharomycotina but positive correlation in Pezizomycotina, and GC content shows negative correlation with number of DNA repair genes in Saccharomycotina but positive correlation in Pezizomycotina (Fig. 3B). These correlations are largely consistent before (i.e., standard Pearson’s correlations) and after (i.e., phylogenetically independent contrasts) accounting for correlations due to phylogeny (table S8).

For each of the seven properties, we used our genome-scale phylogeny (Fig. 1) to infer the ancestral character states and reconstruct their evolution in the Saccharomycotina ancestor and the Pezizomycotina ancestor. Comparison of ancestral states along branches on the Saccharomycotina part of the phylogeny to those on the Pezizomycotina part of the phylogeny showed that all genomic properties, except the number of tRNA genes, exhibited different modes of evolution (Fig. 4 and Table 1). For example, most Saccharomycotina branches exhibit evolutionary rates of at least 1.0 amino acid substitutions per site, whereas those of Pezizomycotina exhibit evolutionary rates between 0.7 and 1.4 substitutions per site (Fig. 4A). However, the inferred values for these properties of the last common ancestors of Saccharomycotina and Pezizomycotina are quite similar. For example, the inferred state values for the Saccharomycotina last common ancestor and the Pezizomycotina last common ancestor are 1.1 and 0.9 substitutions per site for evolutionary rate and 43 and 47% for GC content (Table 1), respectively. The same trends are also observed across lineages, such as Lipomycetaceae, which is the sister group to the rest of the Saccharomycotina, and the clade consisting of Pezizomycetes and Orbiliomycetes, which is the sister group to the rest of the Pezizomycotina (Fig. 4, A and B).

Comparison of the trait values for the seven genome properties between extant Saccharomycotina and Pezizomycotina branches to those of the Saccharomycotina and Pezizomycotina last common ancestors showed that evolutionary rate, GC content, genome size, and number of protein-coding genes were the properties with the highest amounts of evolutionary change (Figs. 3 and 4 and Table 1). Ancestral state reconstruction also enabled inference of the direction of evolutionary change for each of the evolutionary properties. For example, the Saccharomycotina and Pezizomycotina last common ancestors, as well as branches in Lipomycetaceae and branches across Pezizomycotina, exhibit similar evolutionary rates, whereas the rest of the nodes and branches in the Saccharomycotina part of the phylogeny exhibit much higher evolutionary rates. This pattern suggests that the higher levels of genomic diversity in Saccharomycotina stem from an acceleration of evolutionary rate that occurred within the subphylum, after the divergence of Lipomycetaceae from the rest of the Saccharomycotina (Fig. 4, A and B).

Why do Saccharomycotina exhibit higher evolutionary rates compared with Pezizomycotina? Studies in other lineages, such as vertebrate animals (38), have previously shown that evolutionary rate is positively associated with generation time. Assuming that mutation rates are equal, species with shorter generation times will replicate their genomes more frequently, accruing more mutations per unit time. While the generation times of most fungi in our phylogeny are unknown, the generation times of model organisms in Saccharomycotina are thought to be shorter than those in Pezizomycotina. For example, the doubling time of the budding yeasts S. cerevisiae and C. albicans under optimal conditions is 90 min (39), while that of the filamentous fungi Aspergillus nidulans and Neurospora crassa is between 2 and 3 hours (40, 41). It is also well established that absence or loss of DNA repair genes can lead to an increase in genomes’ mutation rates (42). We found that Saccharomycotina have, on average, fewer DNA repair genes (41) than Pezizomycotina (54) (Fig. 3 and Table 1). Thus, an alternative but not mutually exclusive explanation for the higher evolutionary rates in Saccharomycotina may stem at least partially from the presence of a smaller set of genes involved in DNA repair in their genomes than in the genomes of Pezizomycotina. Last, other life history traits (e.g., smaller cell size, faster metabolism, and larger population size) that have been associated with variation in the rate of molecular evolution (43) might also contribute to higher evolutionary rates of the Saccharomycotina.

Variation in genomic GC content has historically been of broad interest in biology. Average GC content values of different genomic regions (e.g., intergenic regions and protein-coding regions) in Saccharomycotina are consistently lower than those in Pezizomycotina (figs. S5 and S6A). Similarly, gene-wise average estimates of GC content showed that all 815 BUSCO genes in Saccharomycotina have lower GC content values than those in Pezizomycotina (fig. S6B). Moreover, we found that the frequencies of amino acids encoded by GC-rich codons in Saccharomycotina are much lower than those of amino acids encoded by GC-rich codons in Pezizomycotina (fig. S7). Ancestral state reconstruction of genomic GC content along branches on the phylogeny shows that the last common ancestors of the Saccharomycotina and Pezizomycotina, as well as branches in Lipomycetaceae and branches in classes Pezizomycetes and Orbiliomycetes, exhibit intermediate GC content of around 45%. In contrast, GC content of most branches within the rest of Saccharomycotina (i.e., all major clades of Saccharomycotina, including extant taxa, except Lipomycetaceae) evolved toward 40%, while GC content within the rest of Pezizomycotina (i.e., all classes, including extant taxa, except Pezizomycetes and Orbiliomycetes) evolved toward 50%. This pattern suggests that the evolution of lower overall GC content in Saccharomycotina occurred after the divergence of Lipomycetaceae from the rest of Saccharomycotina and that the evolution of higher overall GC content in Pezizomycotina occurred after the divergence of the clade consisting of Pezizomycetes and Orbiliomycetes from the rest of Pezizomycotina (Fig. 4, A and B).

Why are Pezizomycotina genomes more GC rich compared with Saccharomycotina genomes? There are two possible explanations. The first one is that mutational biases have skewed the composition of Saccharomycotina genomes toward AT (adenine-thymine) content (44). For example, Steenwyk et al. showed that Hanseniaspora budding yeasts with higher AT content lost a greater number of DNA repair genes than those with lower AT content (42), suggesting that the loss of DNA repair genes is associated with AT richness. Consistent with these results, we found that Pezizomycotina genomes contain a higher number of DNA repair genes than Saccharomycotina (Fig. 3 and Table 1). The second potential, not necessarily mutually exclusive, explanation is that mutational biases have skewed Pezizomycotina genomes toward GC richness. It was recently shown that increasing GC-biased gene conversion (gBGC), a process associated with recombination that favors the transmission of GC alleles over AT alleles, can result in a systematic underestimate of d_N/d_S in birds (45). If this is true for Ascomycota, because of the higher GC content of Pezizomycotina genomes, we would expect that their d_N/d_S would be underestimated due to the higher levels of gBGC compared with Saccharomycotina. Consistent with this expectation, by calculating differences in d_N/d_S before and after accounting for gBGC across 815 codon-based BUSCO genes, we found that the underestimate of d_N/d_S in Pezizomycotina is twofold higher than that in Saccharomycotina (Pezizomycotina: average of differences in d_N/d_S = 0.004; Saccharomycotina: average of differences in d_N/d_S = 0.002) (fig. S8).

To collect the greatest possible set of genome representatives of the phylum Ascomycota as of 14 December 2018, we first retrieved the 332 publicly available Saccharomycotina yeast genomes (https://doi.org/10.6084/m9.figshare.5854692) from a recent comprehensive genomic study of the Saccharomycotina yeasts (2). We then used “Pezizomycotina” and “Taphrinomycotina” as search terms in NCBI’s Genome Browser (www.ncbi.nlm.nih.gov/genome/browse#!/eukaryotes/Ascomycota) to obtain the basic information of strain name, assembly accession number, assembly release date, assembly level (e.g., contig and scaffold), and GenBank File Transfer Protocol (FTP) access number for draft genomes from the subphyla Pezizomycotina and Taphrinomycotina, respectively. For species with multiple isolates sequenced, we only included the genome of the isolate with the highest assembly level and the latest release date. We next downloaded genome assemblies from GenBank data via FTP access number (https://ftp.ncbi.nlm.nih.gov/genomes/). Collectively, we included 332 species representing 77 genera and all 12 major clades of the subphylum Saccharomycotina, 761 species representing 257 genera and nine classes of the subphylum Pezizomycotina, and 14 species representing 6 genera and four classes of the subphylum Taphrinomycotina. Last, we used the genomes of six representatives of the phylum Basidiomycota as out-groups. To check whether our sample of taxa were biased toward particular ecological lifestyles, we obtained the trophic categories for the 1107 Ascomycota species used in this study [symbiotrophs (e.g., ectomycorrhizae and lichens), pathotrophs (e.g., biotrophs, parasites, and pathogens), and saprotrophs (e.g., wood rotters and litter rotters)] from the FUNGuild database (47). We found that Saccharomycotina and Pezizomycotina, the two subphyla with the largest numbers of taxa sampled, contained taxa from all three trophic categories in similar proportions. Specifically, the numbers of Saccharomycotina species classified as symbiotrophs, pathotrophs, saprotrophs, and at least two trophic categories are 2 of 154 (1%), 45 of 154 (29%), 47 of 154 (30%), and 60 of 154 (40%), respectively; the numbers of Pezizomycotina species classified as symbiotrophs, pathotrophs, saprotrophs, and at least two trophic categories are 43 of 702 (6%), 183 of 702 (26%), 176 of 702 (25%), and 300 of 702 (43%), respectively. Detailed information on taxonomy and source of the 1107 Ascomycota and six out-group genomes in our study are provided in tables S1, S2, and S7.

To assess the quality of each of the 1113 genome assemblies, we used the BUSCO version 3.0.2 (31). Each assembly’s completeness was assessed on the basis of the presence/absence of a set of 1315 predefined orthologs (referred to as BUSCO genes) from 75 genomes in the OrthoDB version 9 database from the Ascomycota database, as described previously (27). In brief, for each BUSCO gene, its consensus orthologous protein sequence among the 75 reference genomes was used as query in a tBLASTn search against each genome to identify up to three putative genomic regions, and the gene structure of each putative genomic region was predicted by AUGUSTUS v3.2.2 (48). Next, the sequences of these predicted genes were aligned to the hidden Markov model profile of the BUSCO gene. BUSCO genes in a given genome assembly were considered as single copy, “full-length” if there was only one complete predicted gene present in the genome; duplicated, “full-length” if there were two or more complete predicted genes present in the genome; “fragmented” if the predicted gene was shorter than 95% of the aligned sequence lengths from the 75 reference species; and “missing” if there was no predicted gene present in the genome.

To construct the phylogenomic data matrix, we started with the set of 1315 single-copy, full-length BUSCO genes from 1107 representatives of the phylum Ascomycota and six out-groups. For each BUSCO gene, we first translated nucleotide sequences into amino acid sequences, taking into account the different usage of the CUG (cytosine-uracil-guanine) codon in Saccharomycotina (2). Next, we aligned the amino acid sequences using Multiple Alignment using Fast Fourier Transform (MAFFT) v7.299b (49) with the options “--thread 4 --auto --maxiterate 1000” and trimmed amino acid alignments using the trimAl v1.4.rev15 (50) with the options “-gappyout -colnumbering.” We mapped the nucleotide sequences on the trimmed amino acid alignment based on the column numbers in the original alignment and to generate the trimmed codon-based nucleotide alignment. Last, we removed BUSCO gene alignments whose taxon occupancy (i.e., percentage of taxa whose sequences were present in the trimmed amino acid alignment) was <50% for each subphylum (i.e., <7 Taphrinomycotina, <166 Saccharomycotina, and <381 Pezizomycotina) or whose trimmed alignment length was <300 amino acid sites. These filters resulted in the retention of 815 BUSCO gene alignments, each of which had ≥50% taxon occupancy for each subphylum and alignment length ≥300 amino acid sites.

For each of the 815 BUSCO genes, we first inferred its best-fitting amino acid substitution model using IQ-TREE multithread version 1.6.8 (51) with options “-m TEST -mrate G4” with the Bayesian information criterion. We then inferred best-scoring ML gene tree under 10 independent tree searches using IQ-TREE. The detailed parameters for running each gene were kept in log files (see the Figshare repository). We inferred the concatenation-based ML tree using IQ-TREE on a single node with 32 logical cores under a single “LG + G4” model with the options “-seed 668688 -nt 32 -mem 220G -m LG+G4 -bb 1000”, as 404 of 815 genes favored “LG + G4” as best-fitting model (table S3). We also inferred the coalescent-based species phylogeny with ASTRAL-III version 4.10.2 (52) using the set of 815 individual ML gene trees. The reliability of each internal branch was evaluated using 1000 ultrafast bootstrap replicates and local posterior probability in the concatenation- and coalescence-based species trees, respectively. We visualized phylogenetic trees using the R package ggtree v1.10.5 (53).

We used the RelTime method, as implemented in the command line version of MEGA7 (54), to estimate divergence times. The concatenation-based ML tree with branch lengths was used as the input tree. Six time calibration nodes, which were retrieved from the TimeTree database (55), were used for molecular dating analyses: the S. cerevisiae–Saccharomyces uvarum split (14.3 to 17.94 Ma ago), the S. cerevisiae–Kluyveromyces lactis split (103 to 126 Ma ago), the S. cerevisiae–C. albicans split (161 to 447 Ma ago), the origin of the subphylum Saccharomycotina (304 to 590 Ma ago), the S. cerevisiae–Saitoella complicata split (444 to 631 Ma ago), and the origin of the subphylum Pezizomycotina (at least 400 Ma ago) based on the Paleopyrenomycites devonicus fossil (10, 56).

To explore whether divergence time estimates were robust when using different methods, we also inferred divergence times using a Bayesian MCMCTree analysis (57). As the MCMCTree approach is computationally intractable for very large phylogenomic data matrices like the one we assembled for this study, we pruned the full data matrix to retain 48 representatives of 25 major groups (four classes in Taphrinomycotina, nine classes in Pezizomycotina, and 12 major clades in Saccharomycotina) so that we can estimate the divergence times of 47 nodes on the backbone of the Ascomycota phylogeny. We first estimated branch lengths under a single LG + G4 model with codeml in the PAML 4.9e package (57) and obtained a rough mean of the overall mutation rate. Next, we applied the approximate likelihood method to estimate the gradient vector and Hessian matrix with Taylor expansion (option usedata = 3). Last, we assigned (i) the Gamma-Dirichlet prior for the overall substitution rate (option rgene_gamma) as G(1, 8.36), with a mean of 0.11 (meaning 11 × 10⁻¹⁰ amino acid substitutions per site per year), (ii) the Gamma-Dirichlet prior for the rate-drift parameter (option sigma2 gamma) as G(1, 4.5), and (iii) the parameters for the birth-death sampling process with birth and death rates λ = μ = 1 and sampling fraction ρ = 0. We used the autocorrelated-rate model (option clock = 3) to account for the rate variation across different lineages and used soft bounds (left and right tail probabilities equal 0.025) to set minimum and maximum values for the four calibration splits mentioned above. The Markov chain Monte Carlo (MCMC) run was first run for 100,000 iterations as burn-in, then sampled every 500 iterations until a total of 3000 samples were collected. Our results showed that divergence time estimates for the 47 backbone nodes between the RelTime and MCMCTree methods are very similar, and their correlation is highly significant (Pearson’s correlation coefficient r = 0.9907; P < 0.00001) (fig. S9). Furthermore, the average percentage of deviation (that is the average percentage of the absolute deviation from a mean point) between the two methods is small (∼3%) (table S6).

As the subphylum Taphrinomycotina (no. of species = 14) has a much smaller number of species than the subphylum Saccharomycotina (no. of species = 332) and the subphylum Pezizomycotina (no. of species = 761) in our dataset, we focused our analyses on the comparisons of seven genome properties (evolutionary rate, GC content, genome size, number of genes, number of DNA repair genes, number of tRNA genes, and d_N/d_S) between Saccharomycotina and Pezizomycotina. Specifically, for a given taxon, (i) evolutionary rate is a sum of path distances from the most common ancestor of the subphyla Saccharomycotina and Pezizomycotina to its tip on the concatenation-based ML tree (Fig. 1); (ii) GC content is the percentage of GC nucleotides in genome; (iii) genome size is the total number of base pairs in genome in megabases; (iv) number of genes is the number of protein-coding genes in genome [the gene structure was predicted with AUGUSTUS v3.3.1 (48) on Aspergillus fumigatus and S. cerevisiae S288C trained models for Pezizomycotina and Saccharomycotina, respectively]; (v) number of DNA repair genes was estimated by counting the number of unique protein-coding genes with Gene Ontology terms related to DNA repair using InterProscan version 5 (58); (vi) the number of tRNA genes is the number of tRNA genes inferred to be present using the tRNAscan-SE 2.0 program (59); and (vii) d_N/d_S was estimated by calculating the average of the ratio of the expected numbers of nonsynonymous (d_N) and synonymous substitutions (d_S) along the external branch (tip) of a given taxon across 815 trimmed codon-based BUSCO gene alignments and the concatenation-based ML tree. Note that the concatenation-based ML tree was pruned to match the taxon set in the codon-based alignment if the concatenation-based ML tree contained taxa that were absent from the codon-based alignment. We used the YN98 (F3X4) codon model and the free ratio model using bppml and MapNH in the bio++ libraries, following the study by Bolívar et al. (45).

All statistical analyses were performed in R v. 3.4.2 (R core team 2017). Pearson’s correlation coefficient was used to test for correlations among seven variables. To account for phylogenetic relationships of species in correlation analysis, we used the R package ape v5.1 (60) to compute phylogenetically independent contrasts following the method described by Felsenstein (61).

To reconstruct ancestral character states for each of seven continuous properties, we used the R package phytools v0.6.44 function contMap (62) to infer ancestral character states across internal nodes using the ML method with the function fastAnc and to interpolate the states along each edge using Eq. 2 of Felsenstein (61). The input tree was derived from the concatenation-based ML with branch lengths, which was then pruned to keep the 1093 taxa from the subphyla Pezizomycotina and Saccharomycotina.