Whole Genome Analysis of 132 Clinical Saccharomyces cerevisiae Strains Reveals Extensive Ploidy Variation.

Budding yeast has undergone several independent transitions from commercial to clinical lifestyles. The frequency of such transitions suggests that clinical yeast strains are derived from environmentally available yeast populations, including commercial sources. However, despite their important role in adaptive evolution, the prevalence of polyploidy and aneuploidy has not been extensively analyzed in clinical strains. In this study, we have looked for patterns governing the transition to clinical invasion in the largest screen of clinical yeast isolates to date. In particular, we have focused on the hypothesis that ploidy changes have influenced adaptive processes. We sequenced 144 yeast strains, 132 of which are clinical isolates. We found pervasive large-scale genomic variation in both overall ploidy (34% of strains identified as 3n/4n) and individual chromosomal copy numbers (36% of strains identified as aneuploid). We also found evidence for the highly dynamic nature of yeast genomes, with 35 strains showing partial chromosomal copy number changes and eight strains showing multiple independent chromosomal events. Intriguingly, a lineage identified to be baker's/commercial derived with a unique damaging mutation in NDC80 was particularly prone to polyploidy, with 83% of its members being triploid or tetraploid. Polyploidy was in turn associated with a >2× increase in aneuploidy rates as compared to other lineages. This dataset provides a rich source of information on the genomics of clinical yeast strains and highlights the potential importance of large-scale genomic copy variation in yeast adaptation.

Yeast species classified under the Saccharomyces sensu stricto species complex, in particular S. cerevisiae, have played vital cultural and economic roles in human history (see Sicard and Legras 2011 for review). Indeed, yeast is probably one of the first species domesticated by man. Ancient evidence for the role of Saccharomyces in fermentation dates as far back as the Neolithic period, including wine jars dating to 7000 B.C. China (McGovern ), 5400–5000 B.C. Iran (McGovern ), and ∼3150 B.C. Egypt (Cavalieri ). Another Egyptian study revealed some of the oldest evidence for the active use of yeasts in bread and beer making around 1300–1500 B.C. (Samuel 1996). Yeasts are widely used in industrial applications, with yeast fermentation being harnessed in the development of biofuel technologies that produce ethanol as a renewable source of energy (e.g., Farrell ).

With extensive human contact over nearly 7000 yrs, it is not surprising that S. cerevisiae is often found as an asymptomatic human gut commensal (Zerva ; McCullough ; Posteraro ; Salonen ; Erdem ). However, it has become clear that not all human–S. cerevisiae interactions are harmless. Indeed, S. cerevisiae has been identified with increasing frequency as the causal invasive species in fungal infections, often but not necessarily associated with the use of S. cerevisiae var. boulardii (Perapoch ; Lherm ; Wheeler ; Munoz ; Enache-Angoulvant and Hennequin 2005; de Llanos ,b). While most fungal infections are caused by C. albicans and C. glabrata, S. cerevisiae has been estimated to be responsible for up to 3.6% of all known cases of fungemia, the most severe clinical manifestation of fungal infections (Piarroux ; Smith ; Munoz ). It is also known to cause or be associated with a multitude of other symptoms that contribute to poor medical outcomes, including pneumonia, peritonitis, and liver abscesses (Dougherty and Simmons 1982; Aucott ; Doyle ; Mydlik ; Munoz ). In extreme cases, yeast appears to be the direct cause of mortality, usually by inducing sepsis (Hennequin ; Piarroux ). Because of such reports, S. cerevisiae is now considered a human pathogen, though one of relatively low virulence (Malgoire ).

Unlike species that are fully adapted to a pathogenic lifestyle, it is unclear what factors play a role in the opportunistic virulence of a usually benign commensal. There is clear epidemiological evidence supporting impaired host immunity as an important factor. For example, a large proportion of identified cases were of low virulence and occurred in patients who were immunocompromised, had underlying medical conditions, underwent recent antibiotic treatments, used intravenous catheters, or a combination thereof (Enache-Angoulvant and Hennequin 2005; Munoz ). Infections appear to occur frequently through inadvertent exposure to environmental sources such as contaminated catheters, probiotic treatments, or close contact with other infected individuals (Perapoch ; Lherm ; Cassone ; de Llanos ). This is supported by evidence that S. cerevisiae cannot cross the epithelial barrier in vitro in the absence of epithelial compromise (Perez-Torrado ).

However, even given fortuitous exposure and host condition, there is some evidence that at least two inherent phenotypic traits of yeast strains, the ability to survive at high temperatures and the propensity for pseudohyphal formation, are highly associated with virulence (McCusker ). Virulent isolates are more likely to produce multiple distinct colony phenotypes as compared to avirulent isolates, which may contribute to their ability to survive in complex host environments (Clemons ). A study using mice infected with clinical and nonclinical strains showed that some clinical strains appeared more virulent than nonclinical strains, but virulence phenotypes generally overlap between the two groups, suggesting that preexisting virulence phenotypes play a complementary role to opportunistic exposure (Clemons ). It should be noted, however, that another study found no appreciable differences between clinical and nonclinical strains (Klingberg ), although this result does not contradict the possibility that a subset of ‘safe’ commercial yeast strains are phenotypically more virulent (de Llanos ), and may be more likely to cause infections given the opportunity.

Unlike phenotypic traits, the genetic factors that may influence or drive opportunistic yeast infections, if any, have been more difficult to pin down. Some gene candidates have been reported to be associated with virulence phenotypes (Muller ; Strope ), but most such studies were conducted on haploid derivatives from clinical isolates, whereas a microsatellite study showed that clinical yeasts are often highly aberrant in their ploidy (Muller and McCusker 2009a). This variability in ploidy and chromosomal copy number is unsurprising yet important to our understanding of yeast adaptation. Historically, yeast has had to adapt to new, stressful, and drastically different environments repeatedly in the process of domestication and commercialization. Some of the best examples of how yeast strains accomplished this come from our understanding of how S. cerevisiae adapted to commercial fermentation (see Dequin and Casaregola 2011 for review). Commercial yeasts are extremely tolerant of large-scale copy number changes and both bread and brewing strains are often found to have >2n genomes (Albertin ). Certain strains of sherry yeasts have stable aneuploid chromosomes and a few ale strains have 4n alleles at several loci (Legras ; Albertin ). Large structural rearrangements and copy number variants are widespread and may be associated with selective forces in stressful environments (Bidenne ; Carro ; Wei ; James ; Schacherer ; Dunn ). Aneuploidy, whole chromosomal copy number changes that typically have large fitness effects (Sunshine ), are also common in wild as well as clinical strains (Hose ; Gasch ). Ploidy changes and aneuploidy are thus two large-scale genomic changes known to result in substantial fitness effects in yeast (Cromie and Dudley 2015), and both have been associated with adaptation (Selmecki , 2015; Pavelka ; Rancati and Pavelka 2013; Thorburn ; Zörgö ). However, a large-scale, whole genome survey of ploidy in clinical yeast strains has not been conducted. In order to investigate the true extent of ploidy variation and aneuploidy abundance in clinical yeast strains, we sequenced 144 strains of yeast in their natural ploidy states, the majority of which are clinical yeasts isolates, to an average genomic coverage of 80× per strain. We found that 32% of strains assayed contained >2 copies of each chromosome (3n–5n), and 35% contained aneuploidies at one or more chromosomes. This highlights the importance of whole genome and chromosomal copy number variation in yeast population genomics and adaptation.

Materials and Methods

Library prep and sequencing

DNA extraction and sequencing was carried out in conjunction with strains prepared as in Lang . Briefly, DNA was extracted using a modified glass bead lysis protocol and libraries were constructed with TruSeq DNA library kits. One hundred base pair paired-end Illumina sequencing libraries were prepared at The Genome Institute, Washington University School of Medicine, and run on Illumina GAII with 12 samples pooled per lane (TruSeq adapters AD001-AD012), to an average depth of 80× coverage per sample.

Mapping and SNP calling

The mapping pipeline was similar to that described in Zhu . Briefly, fastq files were mapped to a reference genome panel (consisting of S. cerevisiae, S. paradoxus, S. mikatae, S. kudriavzevii, and S. uvarum) with bwa v0.5.9 (Li and Durbin 2009), sorted and indexed with samtools v0.1.18 (Li ), and assigned strain ID with Picard tools v1.55. Duplicated read pairs were removed and remaining reads locally realigned with GATK v2.1-8 (McKenna ). The UnifiedGenotyper was used to call candidate variants across all samples simultaneously. The resulting VCF file was filtered for variants with coverage depth >8X, >2 reads, and >15% of reads supporting alternative variant. Around 600 kb of the genome – these regions were annotated in the SGD database as simple repeats, centromeric regions, telomeric regions, or LTRs (SGD project – http://downloads.yeastgenome.org/curation/chromosomal_feature/SGD_features.tab, downloaded 4th August 2012) were excluded from analysis due to their susceptibility to mismapping and associated miscalls. For a full list of commands, refer to Supplemental Material and Zhu ).

Aneuploidy calling, gene copy number, and GO term analysis

The aneuploidy calling process was similar to the method used in Zhu and carried out by computing coverage depth of nonoverlapping 1 kb windows. Whole chromosome aneuploidies were called if average coverage of a chromosome differs >35% from other chromosomes in the same line (likelihood P < 0.001, χ2-test), and visually confirmed through coverage traces (see Supplemental Material figures for examples). Gene copy number was estimated using normalized gene coverage against genome-wide distribution, assuming binomial distribution of read depth. GO term analysis was conducted with g:Profiler (Reimand et al. 2011) and AmiGO 2 (Carbon ).

Neighbor-joining tree and Structure

Variants for the Liti strains (Liti ) were identified according to an older version 2 of the yeast genome. Only sequenced single nucleotide polymorphics (SNPs) from the Liti strains (not those inferred by imputation) were converted to version 3 coordinates using LiftOver scripts from the UCSC genome browser website. SNPs called in both datasets with coverage in at least half of the strains within each dataset were retained for neighbor-joining (NJ) tree construction. The NJ tree was constructed using PhyML-3.1 (http://www.atgc-montpellier.fr/phyml), with missing data included as ‘N’. All tree figures were constructed with the online program Interactive Tree of Life (Letunic and Bork 2006). Structure (version 2.3.2) analysis was run with independent allele frequencies under an admixture model, for a 15,000 burn in period and 5000 iterations of sampling.

Data availability

Sequencing data are available through the National Center for Biotechnology Information Sequence Read Archive, number SRP072079.

Results

High coverage whole genome sequencing

One hundred and forty-five S. cerevisiae strains were cultured and sequenced in their natural ploidy states. Thirteen of these strains were from nonclinical sources; the remaining 132 came from clinical sources. Eight strains were recorded as recently derived from strains already present in the collection, and were expected to be genetically very similar. The clinical strains were collected from at least 12 different geographical locations and isolated from at least 14 human tissue types or procedures. Some of the strains have been previously studied in one or more earlier publications: YJM798 was one of the best-studied clinical yeast strains with a fully sequenced genome (Wei ), 14 strains were part of a clinical yeast pathogenesis study (Clemons ), 4 strains were part of a clinical yeast phenotypic study (de Llanos ), 85 strains were part of a genetic diversity by microsatellite study (Muller and McCusker 2009a), 13 strains were part of a GWAS study (Muller ), and 45 strains were sporulated and a single segregant fully sequenced in Strope ) (Supplemental Material, Table S1). Unlike earlier studies, such as Liti and Strope ), we sequenced the clinical strains as is, without sporulation. While extensive heterozygosity, coupled with polyploidy and aneuploidy, makes de novo genome assembly a challenging task, it does provide easily accessible information on chromosomal copy number and heterozygous sites.

Library sequencing, mapping, and analysis protocols are detailed in full under Materials and Methods. Briefly, libraries were constructed for each strain, and 101 bp, paired-end Illumina reads were generated, with an intended average genome-wide read depth of ∼80X. As it is also known that both hybridization and introgression are not uncommon in clinical budding yeasts (Muller and McCusker 2009b; Strope ), the fastq reads were mapped to a combined Saccharomyces sensu stricto genome that included sequences from S. cerevisiae, S. paradoxus, S. mikatae, S. kudriavzevii, and S. uvarum. The S288C reference sequence was downloaded from the Saccharomyces Genome Database (Sheppard ); all other species were downloaded from the Saccharomyces sensu stricto Resources site (www.saccharomycessensustricto.org). One strain was subsequently removed from analysis due to a likely contamination during library construction, with the vast majority of reads failing to map entirely (Materials and Methods). The coverage depth for the other 144 individual strains ranged from 57X to 215X within the S. cerevisiae sequences. Reads were then deduplicated, locally realigned, and recalibrated following best recommended practices on the Genome Analysis Toolkit (GATK) website. SNPs were filtered by coverage depth, mapping quality, and minimum number of reads supporting variant allele. SNPs falling within known repeat sequences and telomeric regions were subsequently excluded from analysis.

No clinical specific point mutations found

Within S. cerevisiae genomic sequences, a total of 477,512 SNPs were called across all non-S288C strains. Individual strains each carry between 36,786 and 97,758 variants (0.3–0.8%) as compared to S288C. Although the number of nonclinical strains was small, allele frequencies at diallelic positions could be compared between clinical and nonclinical strains. Correlation between the clinical and nonclinical populations was high, with an r2 value of 0.9346 (Figure S1). Fst values were calculated for every position, but values did not exceed 0.15, signifying that no SNPs were differentiated between clinical and nonclinical strains (Figure S2). Thus, no further analysis to identify clinical specific alleles was performed.

Population phylogeny recapitulates earlier studies

An NJ tree was constructed with the newly sequenced strains as well as 35 sequences previously published in Liti ) (Figure 1). The two NA soil strains, the three Malaysian strains, the three sake strains, the two baker’s strains, and the lab strains W303 and S288C cluster into independent clades, as would be expected (Strope ). The clade consisting of YJM789, a strain isolated from lung, also included other strains isolated from lung, paracentesis fluid, and blood. Six clinical strains appeared to be closely related to Fleischmann’s baking yeast. Twenty-three clinical strains fell with four commercial yeast strains from Italy on the NJ tree (although it was unclear what precisely the commercial niche was). Four of these 23 strains were documented as clinical isolates from nearby Greece. Three clinical strains further down the tree appeared to be nearly identical to a probiotic strain labeled as S. boulardii. Ten strains were closest to two known vaginal isolates, although information regarding their sites of isolation was not available. A handful of strains fell with RM11_1A, a derivative of a natural vineyard isolate from California. Finally, a couple of clinical strains appeared most similar to baker’s yeast strains.

Figure 1

Neighbor-joining tree of 144 sequenced strains and 35 Liti ) strains, based on SNP differences within the S. cerevisiae genome. Groups of strains that fall within a cluster containing noncommercial strains from a single known source are highlighted in color blocks. Clades containing only clinical strains without a known source are collapsed, with the full tree available in Figure 4.

Figure 4

Heatmap of chromosomal copy number (0n–5n) for each of 144 strains aligned next to each strain on the phylogenetic tree. The 35 Liti 2009 strains do not have chromosomal copy number information panels (in white). The 16 chromosomes of S. cerevisiae are sorted in increasing order in an outwards direction. Panel color from lightest to darkest green represent chromosomal copy number in increasing order.

Population structure identifies multiple origins of clinical strains

Overall, the clinical strains were scattered throughout the clades of nonclinical strains on the phylogenetic tree, suggesting multiple independent origins. To obtain a more explicit grouping of ancestral source population, we used Structure (Pritchard ; Falush , 2007; Hubisz ) on the same set of SNPs used to create the NJ tree (Materials and Methods). The clinical strains were distinguished into five main populations in our analysis. The nonclinical strains that corresponded to these populations included baking strains, sake/soil strains, wine strains, and unspecified commercial strains from Italy. A fifth population consisted of only clinical strains, including YJM789, similar to previous studies that found an unknown source population for clinical strains (Liti ; Strope ). The remainder of the clinical strains showed varying amounts of admixture from various population sources (Figure 2).

Figure 2

Structure assignment of source populations for all 180 strains. Strains are aligned from left to right (x-axis) and sorted by their percentage composition of each population (y-axis) as computed by Structure. Source populations are annotated with text above the figure. Clinical or nonclinical strains are annotated with arrows below the figure. Structure (version 2.3.2) analysis was run with correlated frequencies, admixture model based on base pair distance from preceding SNP, 15,000 burn in period, and 5000 iterations of sampling.

Introgressions from S. paradoxus and S. kudriavzevii

In addition to reads mapping to the S. cerevisiae genome, three strains also had some sequence reads that mapped to other Saccharomyces species’ genomes.

YJM264 is a known hybrid (Muller 2009). Reads covered complete chromosomes from both S. cerevisiae and S. kudriavzevii species. Coverage plots across the genomes showed the presence of at least one copy of all S. kudriavzevii chromosomes except V, XI, XIV, and XVI, and likely two copies of S. kudriavzevii chromosomes IX and XV. S. kudriavzevii chromosomes XII and XIII were partially present, with reads mapping to the first 435 kb of XII and the first 364 kb of XIII. Coverage across the S. cerevisiae chromosomes showed complementary patterns, including lower copy numbers for the partial segments on chromosomes XII (first 450 kb) and XIII (first 390 kb), such that combined across both genomes, there are three copies of all chromosomes. The only exceptions were chromosome V, which appeared to have four copies, and chromosome II, which appeared to have a segmental duplication spanning the first 200 kb (Figure S3).

CBS2909 and CBS2910, which were isolated from human feces, are located immediately adjacent on the phylogeny. Both strains showed highly concordant, though not perfectly identical, regions of high coverage along the S. paradoxus genome, ranging between 2 kb and 40 kb in size. Eight hundred and eighteen such 1 kb windows were observed in CBS2909 and 778 such 1 kb windows were observed in CBS2910. This mapping pattern was not observed in any other strain. Mapping results suggest that they contain mainly S. cerevisiae genomes with S. paradoxus or S. paradoxus-like fragments interspaced across the genome (Figure S4). Removing regions that may be prone to mismapping artifacts (artificial coverage in >50% of all strains), 47 remaining regions consisting of 292 S. paradoxus genes were suspected of introgression into CBS2909 (279 genes) and CBS2910 (274 genes) (Table S2). Two hundred and sixty-one genes were shared by both strains. Other yeast strains have also been previously found to contain introgressions from S. paradoxus (Strope ). S. paradoxus is often found in the same habitats as S. cerevisiae, and such results may indicate ancient hybridization (Sampaio and Gonçalves 2008; Sniegowski et al. 2002).

No other clinical strains showed significant coverage across non-S. cerevisiae sequences (Table S3).

Ploidy variation is common in clinical strains

In S. cerevisiae strains, ploidy variation is common. For 96 of the strains that have been included in previously published work, sporulation patterns indicated that 11 strains were haploids, 56 were diploids, 15 were triploids, and 14 were tetraploids (Muller and McCusker 2009a). However, because a number of generations had elapsed since the strains were first made or studied, we assigned ploidy state to each of our sequenced strains de novo, using observed SNP frequency distributions within the sequence reads. Haploid strains should contain only SNP frequencies of 1; diploids should contain frequencies of both 0.5 and 1, as some SNPs will be heterozygous, while others will be homozygous, and so on (Figure 3). We note, however, that homozygosing events will eliminate heterozygous SNPs, while not changing ploidy. Thus, we were unable to unequivocally rule out euploidy in homozygous strains in the absence of aneuploidy or previous published work. All 11 haploid strains matched expectations. Three of the diploid strains, YJM1104/YJM1105/YJM436, appear to have been recently homozygozed, possibly during isolation and propagation of the strain—their SNP distributions had no visible clusters of SNPs at P = 0.5 as compared to other diploid strains. YJM525 and YJM954 were listed as diploid strains yet had profiles similar to that expected of triploids, while YJM464, YJM1094, YJM1099, and YJM1114 were listed as triploid strains, yet had SNP frequency distributions that indicated diploidy. YJM1146 was also listed as triploid but appeared to be tetraploid based on our SNP read frequency data. Finally, YJM950 and YJM952 were listed as tetraploids, but resembled triploids in our data. Many of these strains possessed aneuploid chromosomes that may have affected sporulation efficiency or viability, leading to the discrepancy between sporulation patterns and our SNP frequency profiles. Across all 144 strains, there were 16 (11%) haploids, 82 (57%) diploids, 23 (16%) triploids, and 23 (16%) tetraploids. Most strains that were 3n or 4n fell within two clades on the phylogeny, the Italian commercial clade and the baker strain clade, suggesting a common recent ancestry (Figure 4).

Figure 3

Top left: allele frequency plot of a typical haploid with SNPs at allele frequencies at 1. Top right: allele frequency plot of a typical diploid with SNPs at allele frequencies at 0.5 and 1. Bottom left: allele frequency plot of a typical triploid with SNPs at allele frequencies at 0.33, 0.67, and 1. Bottom right: allele frequency plot of a typical tetraploid with SNPs at allele frequencies at 0.25, 0.5, 0.75, and 1. Chromosomes are lined up in increasing order (x-axis), with raw SNP allele frequencies along the genome plotted (y-axis).

Aneuploidy is found in more than one third of clinical strains surveyed

Inconsistencies in allele frequency distributions across chromosomes within a single strain were also observed. This could have been due to recent homozygozing of individual chromosomes, or more likely due to individual chromosomes being aneuploid. We assigned the presence of aneuploidies in each strain by examining coverage for 1 kb nonoverlapping windows across all 16 S. cerevisiae chromosomes. Fifty-three of 144 strains (36%) were observed to contain aneuploidies, a number nearly five times higher than the 7.5% previously observed in newly isolated segregants (Strope ). Some chromosomes were more often observed in an aneuploid state than others, possibly correlated with chromosomal size (Figure 5). The most commonly observed whole chromosomal aneuploidy was chromosome IX, with 22 observed events across all strains, as compared to the 1–13 events observed in all other chromosomes. When the rate of observed aneuploidies was compared between strains with ploidy states that are thought to be more stable (haploids and diploids) and strains with ploidy states that are thought to be unstable (≥3n), there was a significantly higher aneuploidy rate in the latter (1n/2n – 27.3%, ≥3n – 56.5%, Fisher’s exact test P = 0.0009). There were also 35 instances of inconsistencies within a chromosome due to a partial chromosomal copy number change (Figure 6). Partial chromosomal copy number changes were most often seen in chromosomes III and VII. There appeared to be no correlation with chromosomal size (Figure 5).

Figure 5

Top: correlation between whole chromosome aneuploid events observed (y-axis) and chromosomal size (x-axis). Bottom: correlation between partial chromosomal copy number changes observed (y-axis) and chromosomal size (x-axis).

Figure 6

Histogram of the total number of complete and partial aneuploid events observed within all 144 sequenced strains, by chromosome.

Combined information from SNP frequency distributions and aneuploidy calls suggested that multiple chromosomal scale copy number events might be present in some strains. One diploid strain YJM1098 appeared to have a 3n aneuploid chromosome XII, but SNP frequencies also suggested that chromosome VIII has been recently homozygozed (Figure 7). Another triploid strain, YJM466, had four chromosomal copies of both chromosome VI and IX, but chromosome VI showed SNP frequency distribution of a tetraploid, while chromosome IX showed that of a diploid (Figure 8). YJM525, YJM674, YJM810, YJM811, YJM813, and YJM815 also showed such indications of multiple events.

Figure 7

Coverage and allele frequency plots for YJM1098. Chromosomes are lined in increasing order (x-axis), with total coverage in 1 kb windows (top panel) or raw SNP allele frequencies (bottom panel) plotted (y-axis). YJM1098 showed a copy number change in chromosome XII with corresponding shift in SNP allele frequency distribution, as well as no copy number change in chromosome VIII, but an allele frequency distribution suggesting recent chromosomal homozygosing event.

Figure 8

Coverage and allele frequency plots for YJM466. Chromosomes are lined in increasing order (x-axis), with total coverage in 1k b windows (top panel) or raw SNP allele frequencies (bottom panel) plotted (y-axis). YJM466 is a triploid strain that showed a copy number increase at chromosomes VI and IX. Chromosome IX showed a corresponding shift in SNP allele frequency distribution to that of 4n, but chromosome VI showed a shift in SNP allele frequency distribution that resembles 2n.

Damaging mutation in NDC80: possible genetic basis of high ploidy/aneuploidy rates in the baking lineage

Strains that were polyploid or carried aneuploid chromosomes were not randomly distributed across the population. One particular clade, containing 36 members from CBS1464 to YJM811 on the phylogeny tree, contained an unusually high number of polyploid and aneuploid strains (Figure 4). Twenty-nine out of 36 (81%) strains from the clade made up of baker’s-like and Italian commercial strains were 3n/4n, whereas only 17/98 (17%) of all other nonsegregant strains were (Fisher’s exact test P = 2.05e-11). Similarly, 24/36 (67%) of strains from the baker’s/Italian clade contained aneuploid chromosomes, whereas 32/109 (29%) of all other strains were (Fisher’s exact test P = 1.3e-04).

We thus investigated if there were any mutations that might explain the high levels of aneuploidy (and ploidy) observed in this particular lineage. Hundreds of genes have been associated with genomic stability and faithful replication and segregation during the cell cycle (Ouspenski ; Stirling ; Cheng ). In particular, a number of genes involved in sister chromatid separation and segregation in yeast have been well studied and characterized. We focused on the genes encoding the kinetochore proteins (), , (YIL144W), , and the chromosome segregation genes , , (Goh and Kilmartin 1993; Biggins ; Cho and Harrison 2011). We used SIFT to look for damaging SNP mutations that are likely to affect the function of these genes (Materials and Methods).

Most of these genes are essential and had very few damaging mutations in the strains sampled (: 0 strains, : 3 strains, : 1 strain, : 4 strains). However, a much larger number of strains carried copies of damaging mutations in (41 strains), (57 strains), and (27 strains). Of these, one particular proline to glutamine mutation in at ChrIX:78597, mostly in the heterozygous form, was only found in the baker/commercial lineage, affected 32/36 strains, and may be part of the explanation for the high rate of aberrant chromosomal copies observed in this clade. Additionally, all 36 baker/commercial lineage strains contained damaging mutations in either (89%) or (89%), with half also carrying affected copies of (54%). These numbers are significantly different from strains outside the lineage, where only 29/98 strains carried any damaging mutations in these genes, and the relative frequencies are just 7% in (Fisher’s exact test P = 2.2e-16), 26% in (P = 4.5e-11), and 10% in (P = 0.035). Strikingly, almost all other damaging mutations in found outside of the baker/commercial lineage were in polyploid strains, with just two exceptions in two closely related segregant strains that shared a unique mutation not found in any other strain.

Gene loss and copy number changes in clinical strains

To look for gene copy number changes that may be adaptive in all clinical strains, we investigated recent gene losses and gains for verified genes in the 144 sequenced clinical strains. Genes lost or gained in a few strains may be recent sporadic events, but a common gene loss or gain found in most clinical strains could indicate an adaptive advantage.

Conservatively, genes were classified as lost only if no reads mapped to any position within the coding sequence of a gene, and the Bonferroni corrected likelihood that a gene of that length is not sequenced given each strains’ average genomic coverage was <0.05. Ninety-five genes passed the filter (Table S4). We focused on recent gene losses – those found only in <5% (1–7 strains) of the sampled population. The 10 most significant GO terms associated with the 52 deleted genes suggested a depletion of deletion in genes involved in cellular component organization and biogenesis (Table 1). There was enrichment in vitamin B6 related biosynthetic and metabolic pathways, largely due to a four-gene deletion (, , , ) in one strain from chromosomes VI and XIV. These genes are found in the subtelomeric region, and are known to vary in copy number (Stambuk ). Conversely, common gene deletions (7–86 strains) showed enrichment in deletion of genes involved in transmembrane glucoside and saccharide transporter activities (Table 1).

Table 1

GO enrichment terms for genes that showed complete deletion, copy number gain, and copy number loss

	#	Expected	Fold Enrichment	+/−	P Value
Deleted genes (1–7 strains)
Biological processes
GO biological process complete
Vitamin B6 biosynthetic process (GO:0042819)	4	0.06	>5	+	1.00E-03
Pyridoxine metabolic process (GO:0008614)	4	0.06	>5	+	1.00E-03
Vitamin B6 metabolic process (GO:0042816)	4	0.07	>5	+	1.60E-03
Cellular macromolecule metabolic process (GO:0044260)	4	19.42	0.21	−	1.95E-03
Pyridoxal phosphate biosynthetic process (GO:0042823)	4	0.09	>5	+	3.53E-03
Pyridoxal phosphate metabolic process (GO:0042822)	4	0.09	>5	+	3.53E-03
Aldehyde biosynthetic process (GO:0046184)	4	0.09	>5	+	3.53E-03
Cellular component organization or biogenesis (GO:0071840)	3	17.35	<0.2	−	3.83E-03
Macromolecule metabolic process (GO:0043170)	5	20.23	0.25	−	4.30E-03
Molecular function
GO molecular function complete
Amine-lyase activity (GO:0016843)	4	0.05	>5	+	1.71E-04
Pyridoxal 5′-phosphate synthase (GO:0036381)	4	0.05	>5	+	1.71E-04
Lyase activity (GO:0016829)	8	0.88	>5	+	2.55E-03
Carbon-nitrogen lyase activity (GO:0016840)	4	0.11	>5	+	4.84E-03
Cellular component
GO cellular component complete
Intracellular (GO:0005622)	24	38.94	0.62	−	5.74E-03
Intracellular part (GO:0044424)	24	38.81	0.62	−	6.98E-03
Macromolecular complex (GO:0032991)	4	16.95	0.24	−	1.40E-02
Intrinsic component of membrane (GO:0031224)	24	11.14	2.15	+	3.96E-02
Deleted genes (>8 strains)
Biological processes
GO biological process complete
Maltose metabolic process (GO:0000023)	3	0.06	>5	+	4.34E-02
Molecular function
GO molecular function complete
Oligosaccharide transmembrane transporter activity (GO:0015157)	3	0.01	>5	+	2.92E-04
Disaccharide transmembrane transporter activity (GO:0015154)	3	0.01	>5	+	2.92E-04
Maltose transmembrane transporter activity (GO:0005363)	3	0.01	>5	+	2.92E-04
α-Glucoside transmembrane transporter activity (GO:0015151)	3	0.02	>5	+	6.91E-04
Glucoside transmembrane transporter activity (GO:0042947)	3	0.02	>5	+	6.91E-04
Sugar transmembrane transporter activity (GO:0051119)	4	0.08	>5	+	1.56E-03
Carbohydrate transmembrane transporter activity (GO:0015144)	4	0.1	>5	+	2.71E-03
Carbohydrate transporter activity (GO:1901476)	4	0.1	>5	+	2.71E-03
Cadmium ion binding (GO:0046870)	2	0.01	>5	+	3.31E-02
Cellular component
GO cellular component complete
Plasma membrane (GO:0005886)	11	2	>5	+	1.08E-03
Ccell periphery (GO:0071944)	11	2.86	3.85	+	3.32E-02
Duplicated genes (1–2 strains)
Molecular function
GO molecular function complete
Thiosulfate sulfurtransferase activity (GO:0004792)	3	0.05	>5	+	1.45E-02
Duplicated genes (>2 strains)
Biological processes
GO biological process complete
Transmembrane transport (GO:0055085)	13	2.8	4.65	+	4.08E-03
Cellular component organization or biogenesis (GO:0071840)	3	15.01	<0.2	−	4.30E-02
Molecular function
GO molecular function complete
Transporter activity (GO:0005215)	14	3.04	4.61	+	8.68E-04
Sucrose α-glucosidase activity (GO:0004575)	3	0.04	>5	+	9.93E-03
β -Fructofuranosidase activity (GO:0004564)	3	0.04	>5	+	9.93E-03
Transmembrane transporter activity (GO:0022857)	11	2.39	4.61	+	1.69E-02
Glucosidase activity (GO:0015926)	4	0.19	>5	+	4.03E-02
Cellular component
GO cellular component complete
Macromolecular complex (GO:0032991)	3	14.67	0.2	−	2.34E-02
Organelle (GO:0043226)	16	29.44	0.54	−	2.71E-02
Intracellular organelle (GO:0043229)	16	29.37	0.54	−	2.97E-02
Organelle part (GO:0044422)	7	19.76	0.35	−	3.36E-02
Plasma membrane (GO:0005886)	12	3.21	3.74	+	3.47E-02
Intracellular organelle part (GO:0044446)	7	19.68	0.36	−	3.70E-02
Loss of copy genes (1–2 strains)
Biological processes
GO biological process complete
Cellular component organization or biogenesis (GO:0071840)	1	11.68	<0.2	−	2.25E-02
Cellular macromolecule metabolic process (GO:0044260)	2	13.07	<0.2	−	3.27E-02
Molecular function
GO molecular function complete
Oligo-1,6-glucosidase activity (GO:0004574)	3	0.03	>5	+	2.68E-03
Acid phosphatase activity (GO:0003993)	3	0.05	>5	+	1.54E-02
α-Glucosidase activity (GO:0090599)	3	0.06	>5	+	2.79E-02
Cellular component
GO cellular component complete
Organelle (GO:0043226)	6	22.9	0.26	−	3.84E-06
Intracellular organelle (GO:0043229)	6	22.85	0.26	−	4.30E-06
Membrane-bounded organelle (GO:0043227)	5	21.18	0.24	−	1.46E-05
Intracellular membrane-bounded organelle (GO:0043231)	5	21.13	0.24	−	1.62E-05
Intracellular (GO:0005622)	11	26.21	0.42	−	5.17E-05
Intracellular part (GO:0044424)	11	26.12	0.42	−	6.29E-05
Organelle part (GO:0044422)	4	15.37	0.26	−	2.42E-02
Intracellular organelle part (GO:0044446)	4	15.31	0.26	−	2.63E-02
Loss of copy genes (>2 strains)
Biological processes
GO biological process complete
Transmembrane transport (GO:0055085)	20	5.1	3.92	+	1.91E-04
Cellular component organization or biogenesis (GO:0071840)	8	27.36	0.29	−	1.11E-03
Maltose metabolic process (GO:0000023)	5	0.16	>5	+	1.28E-03
Macromolecule metabolic process (GO:0043170)	12	31.91	0.38	−	2.51E-03
Regulation of biological process (GO:0050789)	4	19.5	0.21	−	7.50E-03
Cellular macromolecule metabolic process (GO:0044260)	12	30.63	0.39	−	9.19E-03
Cellular component organization (GO:0016043)	7	23.83	0.29	−	1.00E-02
Ion transport (GO:0006811)	15	3.92	3.83	+	1.30E-02
Regulation of cellular process (GO:0050794)	4	18.43	0.22	−	2.28E-02
Thiamine-containing compound biosynthetic process (GO:0042724)	4	0.15	>5	+	3.07E-02
Thiamine biosynthetic process (GO:0009228)	4	0.15	>5	+	3.07E-02
Thiamine metabolic process (GO:0006772)	4	0.16	>5	+	4.19E-02
Molecular function
GO molecular function complete
Transmembrane transporter activity (GO:0022857)	20	4.35	4.6	+	7.55E-06
Transporter activity (GO:0005215)	21	5.53	3.8	+	7.92E-05
Substrate-specific transmembrane transporter activity (GO:0022891)	17	3.92	4.34	+	2.74E-04
Substrate-specific transporter activity (GO:0022892)	17	4.75	3.58	+	3.80E-03
Ion transmembrane transporter activity (GO:0015075)	13	2.86	4.55	+	5.03E-03
Sugar transmembrane transporter activity (GO:0051119)	5	0.25	>5	+	5.51E-03
Oligosaccharide transmembrane transporter activity (GO:0015157)	3	0.04	>5	+	7.76E-03
Disaccharide transmembrane transporter activity (GO:0015154)	3	0.04	>5	+	7.76E-03
Maltose transmembrane transporter activity (GO:0005363)	3	0.04	>5	+	7.76E-03
Sodium-exporting ATPase activity, phosphorylative mechanism (GO:0008554)	3	0.04	>5	+	7.76E-03
Cellular component
GO cellular component complete
Plasma membrane (GO:0005886)	25	5.85	4.27	+	1.88E-07
Cell periphery (GO:0071944)	27	8.38	3.22	+	1.35E-05
Intrinsic component of membrane (GO:0031224)	40	17.56	2.28	+	2.51E-05
Protein complex (GO:0043234)	2	19.2	<0.2	−	6.80E-05
Integral component of membrane (GO:0016021)	38	16.82	2.26	+	9.03E-05
Membrane part (GO:0044425)	40	20.68	1.93	+	2.37E-03

To observe patterns in gene copy number changes that are not deletions, we included only the 87 strains with no observed or suspected recent aneuploidies, where effects from those events are likely to dominate. One hundred and nine genes were significant for duplication (Table S5). Rare duplications (63 genes found in 1–2 strains each) were weakly enriched for thiosulfate sulfurtransferase activity, whereas common duplications (45 genes in 3–67 strains each) were enriched for genes in membrane transport (Table 1). One hundred and nineteen genes were significant for copy number loss (excluding haploid strains) (Table S6). Rare losses (35 genes; 1–2 strains each) were depleted for genes involved in organelle components, whereas common losses (82 genes; 3–68 strains each) were enriched for membrane components (Table 1). Overall, 90/215 genes were members in two or more of the gene loss, gene copy number gain, and gene copy number loss lists, and likely represent genes that are highly variable in copy number among yeast strains.

Discussion

There has been a renewal in interest in natural populations of ‘wild’ S. cerevisiae, and the adaptive events that lead to commercial yeast strains (Barrio ; Blondin ; Muller and McCusker 2009b; Sicard and Legras 2011; Wang ; Liti 2015; Peter and Schacherer 2015). A survey extending beyond the classic European strains used in lab or commercial environments showed that there was larger diversity in undomesticated S. cerevisiae strains from New Zealand (Gayevskiy and Goddard 2012). Another study on Chinese yeast strains found twice the amount of diversity as compared to previously surveyed strains, suggesting that Chinese yeasts may be a candidate source population (Wang ). Geography thus plays a large part in wild yeast population structure, in stark contrast with commercial strains that have been spread across far reaching human communities through human activities (Goddard ). To understand what differentiates virulent, commensal, wild, and commercially ‘safe’ strains in S. cerevisiae, it is necessary to understand the genetics of pathogenic yeast within this ever-growing context of nonclinical strains.

This task is hampered by the fact that a strain’s source substrate is a poor indication of whether it is truly a virulent, pathogenic clinical strain. The process of isolating S. cerevisiae colonies is often heavily biased toward strains that survive well in sugar rich conditions, and it is not always clear if S. cerevisiae is the dominant species occupying those substrates, or is even particularly well adapted to them (Goddard and Greig 2015). Earlier surveys of the S. cerevisiae population structure and phenotypic diversity suggested that the shift from nonclinical to clinical yeast occurred multiple times, ruling out a single selective event that would have left clear genetic signatures. Instead, clinical strains isolated from human patients do not cluster into genetically distinct groups (Liti ; Schacherer ; Muller ), and many remain genetically most similar to commercial strains (likely their source populations) rather than other clinical strains (de Llanos ,b; McKenzie ). The fact that multiple origins of clinical yeast strains have been found suggest a scenario where transitions are common, originate from multiple genetic backgrounds, and are unlikely to share a single adaptive strategy, if they are adapted to a pathogenic lifestyle at all. To fully study the transition of nonpathogenic yeast to a pathogenic lifestyle, either models are needed that can select for clinical behavior in nonclinical strains, to determine precise changes that might afford this change in lifestyle, or thousands of clinical and nonclinical yeast strains must be sequenced for a comprehensive analysis. Here, we start the process by resequencing 144 yeast strains, including 132 clinical isolates, from a variety of geographical locations and tissue substrates.

We find that earlier results hold over this larger set of clinical strains. NJ tree and Structure results both confirm that many clinical isolates are more closely related to specific commercial strains instead of each other. We find evidence that there were likely independent pathogenic invasions from each of lab, sake, wine, and baking yeast strains. Within the rest of the clinical strains, there was no clear segregation of substrate tissue or geographical origin, and closely related strains often include isolates from entirely different tissue types and locations. It is possible that some of these strains may have been isolated from different tissues in the same patient, although we do not have access to the clinical data to confirm or refute that. It is also worth noting that due to considerations of comparability across datasets, we only used SNPs found and sufficiently covered in all datasets, which may bias the analysis toward relatively common polymorphisms.

Some clues about how yeast strains have adapted to drastic changes and stresses in new environments come from our understanding of how S. cerevisiae adapted to commercial fermentation (see Dequin and Casaregola 2011 for review). Functionally, whole genome or whole chromosome copy number changes can be beneficial under certain extreme conditions in yeast. Chromosome VIII, which carries the and genes, consistently increases in copy number in a copper-rich environment (Whittaker ), and in general yeast appears to have the ability to maintain individual chromosomal copy numbers (Waghmare and Bruschi 2005). Chromosome IX, the most often observed aneuploid in the clinical strains, is especially common in aneuploid state, even compared to other chromosomes (II, III, and VIII) that are comparably stable as aneuploids (Zhu ). These structural variants have been suggested to provide one-step reversible short cuts when adapting to a harsh new habitat (Hittinger and Carroll 2007; Coyle and Kroll 2008; Rancati ; Querol and Bond 2009; Stambuk , Muller and McCusker 2009b; Berman 2016), and a similar process may be involved in clinical yeast adaptation going by the high frequency of polyploid and aneuploid strains. Association between the two also agrees well with the thinking that polyploid genomes are less stable than diploids/haploids, making them important drivers of rapid adaptation (Hickman ; Selmecki ). More importantly, aneuploidies are highly unstable if not maintained by strong selective forces (Zhu ; Kumaran ; Zhu ), and high aneuploidy rates could be indications that this is a selectable trait that may have been beneficial during the transition to pathogenic lifestyle (Hirakawa ), and would be worth detailed phenotypic studies, especially in clades such as the baker’s/commercial strains.

It must be noted that the majority of the clinical strains represented in our dataset were of North American and European origin. For a significant number of the clinical strains collected, we were unable to confidently identify the source population. Many of these strains are hybrids of known populations, but some may represent the genetic diversity of uncharacterized commercial or environmental S. cerevisiae. In addition, the lack of clinical strains that cluster with strains collected from Malaysia, Hawaii, and West Africa could indicate sampling bias. The role of increased human contact in opportunistic invasions may favor wine and baking strains, or it could indicate genetic predisposition of domesticated strains to transition to pathogenic lifestyles, perhaps precisely due to their genetic predisposition to easily adapt to drastically different carbon sources. More likely, it could be a combination of the two. Clinical strains thus present a selective view into a subset of medically relevant S. cerevisiae. More clinical strains collected from other parts of the world will help understand if clinical strains regularly colonize human hosts from noncommercial sources, or are mainly able to do so as a side effect of domestication prolonging exposure and opportunity.

It is undeniable that there is a large environmental shift when yeasts move into a pathogenic lifestyle. However, the corresponding genetic changes that resulted from these selective pressures are complex. While our knowledge of yeast diversity, complexity, and ecological traits are expanding, there remains a need for the whole genome deep sequencing of more yeast strains (Peter and Schacherer 2015), as well as the better documentation and classification of what exactly defines a strain as pathogenic and virulent. Better means of analyzing genomic data with variable ploidy across and within individuals will also be important in the study of S. cerevisiae isolates. While our dataset is limited in the types of analysis that can be conducted, in contrast to earlier studies on haploid sequencing data such as the 100-genomes from Strope , it provides information on a different scale and is complementary to our understanding of the many levels of genetic variation available to S. cerevisiae. Such information, combined with human host traits and immune information, may provide the comprehensive picture required to tease apart yeast pathogenicity.