An additional set of 26 genomes were downloaded from the Saccharomyces Genome Database (SGD) and prepared similarly to what is described above for the learning set. These genomes were investigated for the number of CDSs contained and for coverage in terms of the candidate gene markers obtained from the learning set. The total number of CDSs in each genome and the presence of the gene markers obtained from the learning set were investigated. Unfortunately, we found that 12 of these new genomes lacked an important amount of CDSs (less that 4000 versus the expected more than 6000 in the S288C reference strain) possibly indicating some technical issues with the sequencing/annotation processes. Consequently, the total number of CDSs shared with the learning set of strains was highly reduced, seriously compromising the creation of a genome-wide phylogeny. In addition, six genomes sequenced in Borneman et al. (14) were heterozygous and the presence of a high number of ambiguous bases led to the impossibility of including them in the validation set. After this quality control step, we were forced to discard 18 out of 26 available genomes (Validation Table V1 in Supplementary Data). The remaining set of eight new genomes shared 4766 CDSs with the previous strains in the learning set and therefore was suitable for comparison.

The 13 genes we described as classifier (see ‘Results’ section) were actually reduced to eight given that 5 out of 13 genes could not be found in the sequences of the new genomes. For the validation we proceeded as follows: (i) we reproduced the full-genome-based tree as for the learning genome set, containing the already evaluated 39 strains plus the eight additional ones (to map their positioning compared to the original five clusters); (ii) we considered as part of a ‘validation cluster’ only those genomes whose distances from the other members of the clusters were comparable to the cluster definition according to Liti et al.; only four out of the eight genomes fell into one of the previously defined Liti clusters; (iii) we created trees using each of the eight candidate genes independently; (iv) we compared the trees obtained with each candidate gene with the genomic tree to evaluate whether the topology was the same for all the new strains. We therefore considered as validated all the candidate genes that were able to assign at least all but one of the newly analyzed strains to the Liti's genome-wide clustering topology.

The list of genes satisfying the preliminary quality control described above was used to assemble a series of gene-specific fasta formatted files containing all the sequence of that gene in each strain. Such files were subject to multiple sequence alignment by ClustalW version 2.0.10 (15,16) using the fast heuristic pairwise algorithm (quicktree option). After alignment, sequences were optionally clipped at 3′ and/or 5′ to cope with alternative start codons or premature stop codons. The aligned sequences were then processed to extract SNPs/indels (i.e. all zero-entropy alignment columns were removed, meaning that all the positions occupied by gaps in the alignment and representing insertion or deletion events of any extension were retained in the SNPs/indels sequence). Both full-length alignments and SNP/indel-based alignments were subject to a systematic phylogenetic analysis.

We developed an automated pipeline to construct and test complete phylogenetic analyses of genes or gene combinations using both full-length sequences or distilling only the polymorphic/indel sites (SNPs/indels). Specifically, we tested the phylogenetic relationships of 39 available strains of S. cerevisiae. Given a set of complete known genomic sequences, the script was designed to (i) generate multi-fasta files for each of all different genes, (ii) calculate the different combinations of genes (one, two or three), (iii) perform the nucleotide alignments, (iv) identify the SNPs/indels, (v) use distance and neighbor-joining methods to compute the phylogenetic tree using the SNPs/indels or the full-length gene or gene combination, (vi) perform the bootstrap analysis, and eventually (vii) identify the phylogenetic trees that satisfy the evolutionary relationships proposed by the sequence of the entire genomes (see a schematic representation of the pipeline in Figure 1).

We introduced a computational procedure that automated the process of screening the huge number of trees generated from this approach, searching for those actually representing the clusters defined of interest a priori. Basically, we implemented a script that, given a tree, evaluated the leaves descendant from each of its nodes, counting the nodes that encompassed exactly the members of a given cluster. In this work, we investigated the five main clusters identified by Liti and coworkers (12) (named ‘Reference’, ‘West African’, ‘Malaysian’, ‘North American’ and ‘Wine/European’) as those representing the phylogeny of the 39 S. cerevisiae strains. When the evaluation of all the nodes was completed, a total count of five in the tree flagged it as ‘matching’, since it contained all the strains in the same clusters defined by the genome-wide SNPs/indels-based-analysis. In addition, terminal branch lengths were taken into account to evaluate the resolution of the trees at the strain level.

We used all the genes to reproduce the phylogenies proposed by Liti and coworkers and to verify if our procedure using coding sequences could generate the same results. To this aim, the 5850 genes were used to generate a neighbor-joining tree based on pairwise SNPs/indels distances. In this way a total of 226961 SNPs/indels were identified and used. As expected, the phylogenetic tree obtained was superimposable on that of Liti (Figure 2A), confirming that the procedure we developed was consistent and reliable.

In order to test the minimal number of genes necessary and sufficient to reproduce the most important features of this phylogeny—the five most important clusters (see ‘Materials and Methods’ section)—we performed phylogenetic analyses using single genes and combinations of doublets and triplets of genes, for the 39 strains of S. cerevisiae. Table 1 summarizes the number of trees that were tested and those that could reproduce the phylogenesis of interest.

We evaluated the phylogenetic performances of the SNPs/indels of the 5850 genes taken as a single source of variability (singletons). We considered as valid only those genes generating trees harboring the five clusters obtained with the full-genome analysis, according to the scheme presented above. Interestingly, we observed that 41 genes did not present any variation in their sequences and were therefore excluded from the analysis. The complete description of these genes, defined as interesting because of their extreme conservation among all the 39 strains, is reported in Supplementary Table S3 and discussed later on. Of the remaining 5809 genes, 13 were identified as generating matching trees (Supplementary Figures S2–S14) when SNPs/indels-based analyses were performed. The name, identification and a brief description of these genes are presented in Table 2. Even if these genes produced trees containing the five clusters, we did not have any statistical support with SNPs/indels. Therefore, for each one of these 13 selected genes, we repeated the phylogenetic analysis using full-length alignments and assessing the level of confidence at each node with a 100-bootstrap background (17). This analysis pointed out three genes generating matching trees with bootstrap support >60% (Table 2, marked in bold), namely YPR152C, YJL099W and YJL057C. The trees obtained with these genes are reported in Figure 2B, C and D, respectively.

We then focused our attention on the 13 genes able to recapitulate the genome-wide phylogeny (Table 2). These genes presented a high number of polymorphic sites, ranging from 20 to 628, the latter being a clear evidence of an insertion/deletion event. We initially observed that, though highly informative, none of the 13 genes alone could separate the strains sharing the same cluster. We therefore used the pipeline we created to investigate whether combinations of two or three of these genes could improve the identification at the strain level. Both SNPs/indels and full-length genes were analyzed, the latter with bootstrap support on nodes (see ‘Material and Methods’ section). We found that neither doublets nor triplets could effectively individualize the strains in all five clusters, though many combinations were able to resolve strains in at least three out of five clusters. In particular, we observed that, as expected, the reference cluster (composed of laboratory strains labeled with REF, S288C and W303) was never resolved, since ‘S288C’ was a simple resequencing control of the ‘REF’ strain. The ‘West African’ cluster (composed only by NCYC110 and DBVPG6044 strains) was hardly ever resolved. The remaining three clusters were instead individualized by several combinations of genes. Note that these clusters contain 34 out of 39 strains available. This indicates that the genes we isolated were able to recapitulate the phylogeny observed with the genome-wide analysis, and to classify and fingerprint 87% of the strains, if used in combinations. In SNPs/indels-based analysis, we found that 7 out of the 76 matching trees obtained with doublets had an efficient strain resolution, as well as 50 out of the 286 matching trees obtained with triplets. In contrast, the full-length analysis (with bootstrap support) evidenced that 1 out of 51 doublets and 21 out of 211 triplets had strain resolution. Of particular relevance was the phylogeny obtained with the full YNL161W gene (CBK1, a serine/threonine-protein kinase that regulates pathways of cell morphogenesis), that was sufficiently variable to identify strains when used as the single source of variability (Figure 2E). Such extremely high performance is due to a long insertion/deletion event as evident by observing the alignment of the sequence of this gene. Given that result, it was not a surprise to observe that all the doublets and triplets—that reached the strain resolution—contained the YNL161W gene. All the results obtained with this combinatorial analysis are summarized in Table 3.

During recent years, the variability of some S. cerevisiae genes has been studied as phylogenetic classifiers. We selected 10 genes (Table 4) from published results (5–7) and tested their ability in reproducing classification and strain resolution using the combinatorial approach implemented in our pipeline. Though they presented a number of SNPs/indels ranging from 16 to 71 (Table 4), none of them was able to produce matching trees if used as the sole source of variability using either SNP/indel-based alignment or full-length alignments. The 45 possible doublets tested resulted in six matching trees in SNP/indel-based analysis (involving 7 of the 10 genes) and only two matching trees (involving only four genes) in full-length analysis with bootstrap support and a threshold of 60%. Triplets performances were low: of the 120 tested combinations, only 25 generated matching trees (with eight genes involved) in SNP/indel-based analysis and 11 in full-length analysis (involving the same eight genes). Interestingly, no trees exhibited a separation at the strain level for more than one cluster, indicating a poor resolution at the strain level. All the results obtained with this combinatorial analysis are summarized in Table 5.

In addition to the five clusters identified by Liti and coworkers on the basis of a genome-wide phylogenetic analysis, several other clusters were evaluated, taking advantage of the pipeline we created. The aim was to assess the capability of some genes or gene combinations to discriminate strains according to different definitions (see Supplementary Table S3 for a complete list of clusters). We initially tested several geography-based clusters (five in total) with different granularity. No gene, gene doublet or triplet generated trees matching any of such cluster definitions. We then tried to investigate whether the source of isolation of the strains could be a possible hub for identifying discriminant genes. Also in this case, no combination resulted in matching trees. We therefore concluded that no gene can be used to assign functional, industrial or ecological roles or provenience to a strain.

In order to investigate whether non-coding regions could be used as a high performance genetic marker, we applied our pipeline to a data set composed of a 400bp region upstream of the first base of the coding sequence of each gene present in all 39 available strains. The choice of this length window for 5′-UTR investigation was based on the analysis of the length of the intergenic regions of the S288C genome resulting in an average of 475±618bp and a median value of 375bp (Supplementary Figure S1). We considered a size of 400bp a good balance, allowing the collection of phylogenetic variations while ensuring not to intersect other coding regions. We investigated a total of 5822 5′-UTR regions using both SNPs/indel-based alignments and full-length alignments. When challenged against the five clusters, none of the trees generated by UTRs proved to be able to match the genome-wide phylogeny, indicating that the variations in gene promoters or other regulatory sequences were not sufficient or at least not well suited to discriminate S. cerevisiae strains at a phylogenetic level. It is interesting to note that a slight variability exists in the 5′-UTR regions of the 41 genes previously described as absolutely conserved among S. cerevisiae strains, reinforcing our finding that the coding sequences of these genes are a feature characterizing the cerevisiae species among the other members of the Saccharomyces genus. As a final test, we repeated the analysis on a 1000bp long 5′-UTR region (without checking for overlaps with preceding coding sequences), evidencing also in this case no satisfactory trees (data not shown).

Careful analysis of the information present in SGD (http://www.yeastgenome.org) indicates that the proposed phylogenetic markers fall into specific gene categories. It is worth noticing that 6 out of the 13 genes enabling the reconstruction of the phylogeny are involved in biological processes related to membrane metabolism. Three genes, namely YNL161W-CBK1 (coding for a serine/threonine-protein kinase located at the top of the growing apical bud in budding cell), YJL099W-CHS6 (coding for a chitin-biosynthesis protein component of the CHS5/6 complex) and YJL051W-IRC8 (coding for a multipass membrane protein located at the cellular bud) play a role in the regulation of budding, cell morphogenesis and proliferation, recombination and cell division. Three other genes are involved in membrane homeostasis such as YAR042W-SWH1 (coding for a lipid-binding protein involved in maintenance of intracellular sterol distribution and homeostasis), YNL125C-ESBP6 (coding for an uncharacterized multipass transporter in the endoplasmic reticulum) and YKL068W-NUP100 (coding for a nucleoporin acting as a component of the nuclear pore complex). Four genes are involved in fundamental nucleic acid related processes: YML056C-IMD4 is involved in the de novo purine biosynthesis, YML080W-DUS1 codes for a tRNA-dihydrouridine synthase that modifies the U16 and U17 positions—and allows the maturation of the tRNA molecule—and, YPR152C-URN1 contributes to the formation of the spliceosome with other splicing factors. Two other genes also take part in processes relevant for yeast metabolism: YJL057C-IKS1, that codes for a putative serine/threonine-protein kinase whose disruption leads to hypersensitivity to copper sulfate (one of the most ancient antifungals, used in agriculture and wine production since Roman times) and could probably be under selection in association with human intervention and YBL052C-SAS3, coding for a histone acetyltransferase that acts as the catalytic component of the NuA3 complex. The acetylation of Lys-14 of histone H3, involved in the silencing the HMR locus, is crucial for mating and therefore is an important process for the evolution of yeast population structure. The last two marker genes are YBR163W-EXO5, the mitochondrial exonuclease V involved in mitochondrial DNA replication and recombination and YOR133W-EFT1, an elongation factor supposed to be highly conserved that form the 80S ribosome and interacts with both RNA and proteins.