Evolution of pathogenicity and sexual reproduction in eight Candida genomes.

Candida species are the most common cause of opportunistic fungal infection worldwide. Here we report the genome sequences of six Candida species and compare these and related pathogens and non-pathogens. There are significant expansions of cell wall, secreted and transporter gene families in pathogenic species, suggesting adaptations associated with virulence. Large genomic tracts are homozygous in three diploid species, possibly resulting from recent recombination events. Surprisingly, key components of the mating and meiosis pathways are missing from several species. These include major differences at the mating-type loci (MTL); Lodderomyces elongisporus lacks MTL, and components of the a1/2 cell identity determinant were lost in other species, raising questions about how mating and cell types are controlled. Analysis of the CUG leucine-to-serine genetic-code change reveals that 99% of ancestral CUG codons were erased and new ones arose elsewhere. Lastly, we revise the Candida albicans gene catalogue, identifying many new genes.

Four species, C. albicans, C. glabrata, C. tropicalis and C. parapsilosis, together account for ∼95% of identifiable Candida infections1. While C. albicans is still the most common causative agent, its incidence is declining and the frequency of other species is increasing. Of these, C. parapsilosis is a particular problem in neonates, transplant recipients, and patients receiving parenteral nutrition; C. tropicalis is more commonly associated with neutropenia and malignancy. Other Candida species, including C. krusei, C. lusitaniae and C. guilliermondii, account for <5% of invasive candidiasis. Almost all Candida species, with the exception of C. glabrata and C. krusei, belong in a single Candida clade (Fig. 1) characterized by the unique translation of CUG codons as serine rather than leucine2. Within this, haploid and diploid species occupy two separate subclades (Fig. 1).

Figure 1

Phylogeny of sequenced Candida and Saccharomyces clade species. Tree topology and branch lengths were inferred with MrBayes (see supplementary methods S5). Posterior probabilities are indicated for each branch. The * marks a branch that was constrained based on syntenic conservation40.

To determine the genetic features underlying their diversity of biology and pathogenesis, we sequenced six genomes from the Candida clade (Fig. 1). These include a second sequenced isolate of C. albicans (WO-1) characterized for white-opaque switching, a phenotypic change that correlates with host specificity and mating3, 4. We also sequenced the major pathogens C. tropicalis and C. parapsilosis, Lodderomyces elongisporus, a close relative of C. parapsilosis recently identified as a cause of bloodstream infection5, and two haploid emerging pathogens, C. guilliermondii and C. lusitaniae. We compared these to the previously sequenced C. albicans strain (SC5314)6-8, Debaryomyces hansenii9, a marine yeast rarely associated with disease, and nine species from the related Saccharomyces clade (Fig. 1). These species span a wide evolutionary range and show large phenotypic differences in pathogenicity and mating, allowing us to study the genomic basis for these traits.

Genome sequence and comparative annotation

We found enormous variation in genome size and composition between the Candida genomes sequenced (Table 1). Each genome assembly displayed high continuity, ranging from nine to 27 scaffolds (Supplementary table 1). Scaffold number and size largely match pulsed-field gel electrophoresis estimates for all genomes, and telomeric repeat arrays are linked to the ends of nearly all large scaffolds (Supplementary text S2). Genome size ranges from 10.6 to 15.5 Mb, a striking difference of nearly 50%, with haploid species having smaller genomes. GC content ranges from 33% to 45% (Table 1). Transposable elements and other repetitive sequences vary in number and type between assemblies (Supplementary text S6). Regions similar to the Major Repeat Sequence (MRS) of C. albicans were found only in C. tropicalis, suggesting that MRS-associated recombination could contribute to the observed karyotypic variation between C. tropicalis strains10.

Table 1

Candida genome features

Species*	Genome size (Mb)	%GC	Genes (#)	Gene size avg. (bp)	Intergenic avg. (bp)	Ploidy	Pathogen†
Candida albicans WO-1	14.4	33.5%	6,159	1,444	921	diploid	++
Candida albicans SC5314	14.3	33.5%	6,107	1,468	858	diploid	++
Candida tropicalis	14.5	33.1%	6,258	1,454	902	diploid	++
Candida parapsilosis	13.1	38.7%	5,733	1,533	752	diploid	++
Lodderomyces elongisporus	15.4	37.0%	5,802	1,530	1,174	diploid	−
Candida guilliermondii	10.6	43.8%	5,920	1,402	426	haploid	+
Candida lusitaniae	12.1	44.5%	5,941	1,382	770	haploid	+
Debaryomyces hansenii	12.2	36.3%	6,318	1,382	550	haploid	−

C. albicans SC5314 assembly 21 and gene set dated 28-Jan-2008 downloaded from the Candida Genome Database (www.candidagenome.org); D. hansenii assembly from GenBank9. The remaining assemblies are reported as part of this work, and are available in GenBank and at http://www.broad.mit.edu/annotation/genome/candida_group/MultiHome.html

Relative level of pathogen strength (++, strong pathogen; +, moderate pathogen; −, rare pathogen).

Despite the genome size and phenotypic variation among the species, the predicted numbers of protein-coding genes are very similar ranging from 5,733 to 6,318 genes (Table 1). Even the small differences in gene number are not correlated with genome size; the smallest genome, C. guilliermondii, has more genes than the largest genome, L. elongisporus. Instead, genome size differences are explained by a ∼3-fold variation in intergenic spacing (Table 1). Large syntenic blocks of conserved gene order were detected among the four diploid species and between two of the haploid species, C. guilliermondii and D. hansenii (Supplementary fig. 5). While synteny blocks have been shuffled by local inversions and rearrangements, these have been primarily intrachromosomal as chromosome boundaries have been largely preserved across the diploid genomes (Supplementary fig. 5).

Given the high conservation of protein-coding genes across the Candida clade, we used multiple alignments of the related genomes to revise the annotation of C. albicans. We identified 91 new or updated genes, of which 80% are specific to the Candida clade (Supplementary text S4). We also corrected existing annotations in C. albicans, revealing 222 dubious genes, and also identified 190 likely frame-shifts and 36 nonsense sequencing errors in otherwise well-conserved genes (Supplementary text S4). In each case, manual curation confirmed ∼80% of these predictions.

Polymorphism in diploid genomes

To gain insights into the recent history of C. albicans, we compared the two diploid strains, SC5314 and WO-1, which belong to different population subgroups11. Variation in the karyotype of these strains is primarily due to translocations at MRS sequences (12 and Supplementary fig. 1). The two assemblies are largely co-linear with 12 inversions of 5−94 kb between them, except that in WO-1 some non-homologous chromosomes have recombined at the MRS (Supplementary text S8). We found similar rates of single nucleotide polymorphisms (SNPs) within each strain (1 SNP per 330−390 bases), and twice this level between them, suggesting relatively recent divergence. Polymorphism rates in the other diploids range from 1 SNP per 222 bases in L. elongisporus and 1 SNP per 576 bases in C. tropicalis, to a remarkably low 1 SNP per 15,553 bases in C. parapsilosis, more than 70-fold lower than the closely related L. elongisporus.

Striking regions of extended homozygosity are found in three of the four diploid genomes, which may reflect break-induced replication, or recent passage through a parasexual or sexual cycle. C. albicans, C. tropicalis, and L. elongisporus each show large chromosomal regions devoid of SNPs, extending up to ∼1.2 Mb (Fig. 2, Supplementary fig. 6−8). In contrast, the few SNPs in C. parapsilosis are randomly distributed across the genome (Supplementary fig. 9A). A total of 4.3 Mb in total (30%) of the WO-1 assembly is homozygous for SNPs, approximately twice that found in SC5314 (Supplementary text S7 and 7, 13, 14). There is at least one homozygous region per chromosome, none of which spans the predicted centromeres, and only one of which starts at a MRS (Fig. 2). While nearly all homogeneous regions are present at diploid levels and are therefore homozygous, WO-1 has lost one copy of a >300 kb region on chromosome 3 comprising nearly 200 genes (Supplementary fig. 10). The pressure to maintain this region as homozygous in both strains is apparently high, as it is diploid but homogenous in SC5314.

Figure 2

C. albicans WO-1 is highly homozygous. Red lines show SNPs per kb, normalized by coverage, within WO-1, and blue lines show SNPs per kb between WO-1 and SC5314. While both copies of chromosome 5 have rearranged at the MRS (yellow box) in WO-1, we show this as a single chromosome to allow a haploid reference for polymorphism. Relative to SC5314, chromosomes 1, 4, and 6 are in the opposite orientation (Supplementary figure 6).

Usage and evolution of CUG codons

All Candida clade species translate CUG codons as serine instead of leucine15. This genetic code change altered the decoding rules of CUN codons in the Candida clade: while S. cerevisiae uses two tRNAs which each translate two codons, Candida species use a dedicated tRNACAGSer for CUG codons and a single tRNALAGLeu for CUA, CUC, and CUU, as inosine can base pair with A, C, and U (Fig. 3). This remarkable alteration in decoding rules forced the reduced usage of CUG, and also CUA likely due to the weaker wobble, in Candida genes (Fig. 3A). CUU and CUC codons do not display the same bias for infrequent usage (Fig. 3B). An additional pressure influencing codon usage may be the GC content, as usage of leucine codons in Candida species is correlated with %GC composition (Supplementary table 11).

Figure 3

Evolutionary effects of CUG coding. A. Average percent of genes with 0, 1, or 2 CUG and CUA codons in Candida (grey bars) and Saccharomyces (black bars). Error bars indicate standard deviations. All differences are significant with p ≤ 0.0004 (t-test) except for genes with two CUA codons genes (p=0.02). “Candida” and “Saccharomyces” here refer to the CTG and WGD clades in Figure 1, but including P. stipitis and excluding C. dubliniensis. B. CUN codon usage for all codon counts. C. Decoding rules for CUN codons in Saccharomyces and Candida.

We also examined the evolutionary fate of ancestral CUG codons and the origin of new CUG codons (Supplementary table 12). CUG codons in C. albicans almost never (1%) align opposite CUG codons in S. cerevisiae. Instead, CUG serine codons in C. albicans align primarily to Saccharomyces codons for serine (20%) and other hydrophilic residues (49%). CUG leucine codons in S. cerevisiae align primarily to leucine codons in Candida (50%) and to other hydrophobic residue codons (30%). This suggests a complete functional replacement of CUG codons in Candida.

Gene family evolution

To identify gene families likely to be associated with Candida pathogenicity and virulence, we used a phylogenomic approach across seven Candida and nine Saccharomyces genomes (Supplementary text S10). Among 9,209 gene families, we identified 21 that are significantly enriched in the more common pathogens (Table 2). These include lipases, oligopeptide transporters (Opts), and adhesins, all known to be associated with pathogenicity8, as well as poorly characterized families not previously associated with pathogenesis.

Table 2

Gene families enriched in pathogenic Candida sp.

# Annotation	Pathogengenes	Nonpathogengenes	Pval.	Dup.	Loss	Generate	C.alb	C.tro	C.par	L.elo	C.gui	C.lus	D.han	C.gla	Yeast(avg.)
1 GPI-family 18 (Hyr/Iff-like)	56	10	1.4E-16	52	11	16.2	11	18	17	9	3	7	1	0	0.0
2 Leucine-rich repeat (IFA/FGR38-like)	34	0	4.2E-16	32	5	18.3	33	1	0	0	0	0	0	0	0.0
3 Ferric reductase family	45	10	1.9E-12	30	25	2.5	12	19	7	7	3	4	2	0	0.1
4 Reductase family	43	11	3.2E-11	31	30	2.3	7	12	9	6	13	2	4	0	0.1
5 GPI-family 17 (ALS-like Adhesins)	31	5	4.4E-10	29	4	20.5	8	16	5	4	2	0	1	0	0.0
6 GPI-family 13 (Pga30-like)	34	7	5.0E-10	25	5	14.8	12	14	6	6	1	1	1	0	0.0
7 Unclassified	20	0	9.0E-10	13	3	15.9	9	9	0	0	2	0	0	0	0.0
8 Cell wall mannoprotein biosynthesis	38	18	7.2E-07	19	34	2.1	8	7	8	8	11	4	9	0	0.1
9 Major facilitator transporters	25	7	9.2E-07	14	17	2.0	3	3	7	3	10	2	4	0	0.0
10 Oligopeptide transporters	31	13	2.2E-06	23	11	6.7	6	9	9	4	4	3	1	0	0.9
11 Unclassified	25	9	6.3E-06	15	6	11.1	7	9	3	5	3	1	4	2	0.2
12 Amino acid permeases	27	11	7.7E-06	11	18	1.7	6	6	6	4	6	3	6	0	0.1
13 Sphingomyelin phosphodiesterases	18	5	3.2E-05	11	9	7.4	4	5	4	2	3	2	1	0	0.2
14 FGR6 family (filamentous growth)	12	1	3.3E-05	7	1	14.5	8	1	1	0	1	1	1	0	0.0
15 Secreted lipases	20	7	4.6E-05	17	8	9.6	10	5	4	4	1	0	3	0	0.0
16 Cytochrome p450 family	34	21	5.5E-05	23	22	6.0	6	8	10	7	5	4	6	1	1.0
17 Amino acid permeases	16	4	5.6E-05	14	10	1.5	2	3	6	3	2	3	1	0	0.0
18 Zinc-finger transcription factors	31	18	6.2E-05	17	14	12.3	5	8	7	7	7	4	11	0	0.0
19 Unclassified	13	2	6.3E-05	8	0	8.1	3	1	6	1	2	1	1	0	0.0
20 Predicted transmembrane family	17	5	7.2E-05	9	2	7.5	4	4	5	3	3	1	2	0	0.0
21 Unclassified secreted family	20	8	1.1E-04	7	6	9.3	4	4	6	4	4	2	4	0	0.0

Pathogen genes = total genes in family for C. albicans, C. tropicalis, C. parapsilosis, C. guilliermondii, C. lusitaniae, and C. glabrata. Nonpathogen genes = total genes in family for L. elongisporus, D. hansenii, and all Saccharomyces clade species (Figure 1) expect C. glabrata. Pval.= Pvalue of the hypergeometric test; all families shown above have a false discovery rate less than 0.05 (Supplementary text S10c). Dup., Loss = Duplications and losses (Supplementary Methods 10b). Gene rate = average mutation rate for each family (Supplementary text S10d); the average gene rate across all families is 5.8. Yeast (avg.) = average count for all Saccharomyces species.

Three cell wall families are enriched in the pathogens: Hyr/Iff16, Als adhesins17, and Pga30-like proteins (Table 2, Supplementary text S11). The Als family (family 17 in Supplementary tables 22 and 23) in C. albicans is associated with virulence, and in particular with adhesion to host surfaces18, invasion of host cells19, and iron acquisition20. All these families are absent from the Saccharomyces clade species and are particularly enriched in the more pathogenic species (Supplementary table 18). All three families are highly enriched for gene duplications (Supplementary table 19), including tandem clusters of 2−6 genes, and show high mutation rates (fastest 5% of families) (Supplementary text S10d). This variable repertoire of cell wall proteins is likely of profound importance to the niche adaptations and relative virulence of these organisms.

Als17 and Hyr/Iff proteins frequently contain intragenic tandem repeats (ITRs), which modulate adhesion and biofilm formation in S. cerevisiae21 (Figure 4, Supplementary Fig. 19, 20). The ITR sequence is conserved at the protein level across species (Supplementary fig. 19−20). Interestingly, two proteins contain both an Als domain and repeats characteristic of the Hyr/Iff family.

Figure 4

Conserved domains of Als and Hyr/Iff cell wall families. The N-terminal Als domain (red) and Hyr/Iff domain (green) are shown as ovals. Intragenic tandem repeats (ITRs, see Supplementary text S11) are shown as rectangles, colored to represent similar amino acid sequences.

Candida clade pathogens show expansions of extracellular enzyme and transmembrane transporter families (Table 2 and Supplementary table 22). These families are either not found in Saccharomyces (including amino acid permeases, lipases, and superoxide dismutases), or present in S. cerevisiae but significantly expanded in pathogens (including phospholipase B, ferric reductases, sphingomyelin phosphodiesterases, and GPI-anchored yapsin proteases, which have been linked to virulence in C. glabrata22). Several groups of cell-surface transporters are also enriched (including Opts, amino acid permeases, and the major facilitator superfamily). Overall these family expansions illustrate the importance of extracellular activities in virulence and pathogenicity. Genes involved in stress response are also variable between species (Supplementary text S12).

C. albicans also showed species-specific expansion of some families, including two associated with filamentous growth, a leucine-rich repeat family and the Fgr6−1 family (Table 2). As C. albicans forms hyphae while C. tropicalis and C. parapsilosis produce only pseudohyphae, these families may contribute to differences in hyphal growth.

We identified 64 families showing positive selection in the highly pathogenic Candida species (Supplementary table 32). These are highly enriched for cell wall, hyphal, pseudohyphal, filamentous growth, and biofilm functions (Supplementary text S13). Six of the families have been previously associated with pathogenesis, including ERG3, a C-5 sterol desaturase essential for ergosterol biosynthesis, for which mutations can cause drug resistance23.

Structure of the MTL locus

Pathogenic fungi may have limited their sexual cycles to maximize their virulence24, and the sequenced Candida species show tremendous diversity in their apparent abilities to mate. Among the four diploids, C. albicans has a parasexual cycle (mating of diploid cells followed by mitosis and chromosome loss instead of meiosis25), L. elongisporus has been described as sexual and homothallic (self-mating)26, while C. tropicalis and C. parapsilosis have never been observed to mate. Among the three haploids, C. guilliermondii and C. lusitaniae are heterothallic (cross-mating only) and have a complete sexual cycle, while D. hansenii is haploid and homothallic27 (Supplementary text 14c).

To understand the genomic basis for this diversity, we studied the Candida mating-type locus (MTL), which determines mating type, similar to the MAT locus in S. cerevisiae. In both C. albicans and S. cerevisiae, the mating locus has two idiomorphs, a and α, encoding the regulators a1 and α1/α2 respectively. C. albicans MTLa also encodes a2, and both idiomorphs in this species contain alleles of three additional genes without known roles in mating; PAP, OBP, and PIK28. The MTLα and MTLa genes, alone or in combination, specify one of the three possible cell-type programs (a and α haploid, a/α diploid). In C. albicans, the alpha-domain protein α1 activates α-specific mating genes, the HMG factor a2 activates a-specific mating genes, and the a1/α2 homeodomain heterodimer represses mating genes in a/α cells29, 30.

Despite extended conservation of the genomic context flanking MTL, there is tremendous variability in MTL gene content (Fig. 5). MTLa1 has become a pseudogene in C. parapsilosis 31, likely a recent loss since target genes retain predicted a1/α2 binding sites (Supplementary text S14). MTLα2 is missing in both C. guilliermondii and C. lusitaniae (Reedy, Floyd, and Heitman, submitted). A fused mating type locus containing both a and α genes is found in D. hansenii and Pichia stipitis32, 33.

Figure 5

Organization of MTL loci in the Candida clade. A. MTLa-specific genes are shown in grey, MTLα-specific in black, and other orthologs in color. Two idiomorphs from C. albicans and C. tropicalis are shown. Arrows indicate inversions relative to C. albicans. “X” symbols show gene losses and the C. parapsilosis MTLa1 pseudogene. There is no MTL locus in L. elongisporus; gene order around the PAPa, OBPa and PIKa genes is shown. For C. guilliermondii and C. lusitaniae the genome project sequenced one idiomorph, the second was obtained by J. Reedy. In D. hansenii, PAP, OBP and PIK are separate from the fused MTL locus32. B. Placement of gene losses on the phylogenetic tree.

Most surprisingly, all four mating-type genes are missing in L. elongisporus. It contains a site syntenic to MTLa in other species, but this contains only 508 base pairs of apparently noncoding DNA, a length insufficient to encode a1 or a2 even if they were extensively divergent in sequence. We confirmed this finding in seven other L. elongisporus isolates (not shown). The sexual state of L. elongisporus has been assumed to be homothallic because asci are generated from identical cells5, 26, but the absence of an a-factor pheromone, receptor and transporter (Supplementary table 34) as well as MTL, suggests that it may not have a sexual cycle. Alternatively, mating may require only one pheromone and receptor (α), and L. elongisporus may be the first identified ascomycete that can mate independently of MTL/MAT.

Our discovery that MTLα2 and MTLa1 are frequently absent challenges our understanding of how the mating locus operates. The model derived for mating regulation in C. albicans, including the role of a1/α2 in white-opaque switching3, must differ substantially in the other sexual species. This is particularly interesting because the major regulators of the white-opaque switch in C. albicans (WOR1, WOR2, CZF1, EFG134) are generally conserved in other species, but Wor1 control over white-opaque genes appears to be a recent innovation in C. albicans and C. dubliniensis35. Our data also suggest that the loss of α2 occurred early in the haploid sexual lineage.

Mating and meiosis

To gain further insights in their diversity in sexual behavior, we examined whether 227 genes required for meiosis in S. cerevisiae and other fungi have orthologs in the Candida species (Supplementary text S14). A previous report36 that some components of meiosis, such as the major regulator IME1, are missing from C. albicans, led us to hypothesize that their loss could be correlated with lack of meiosis. Surprisingly however, we find that these genes are missing in all Candida species, suggesting that sexual Candida species undergo meiosis without them. Conversely, even seemingly non-mating species showed highly conserved pheromone response pathways, suggesting pheromone signaling plays an alternate role such as regulation of biofilm formation37. These findings suggest considerable plasticity and innovation of meiotic pathways in Candida.

Moreover, we find that sexual Candida species have undergone a recent dramatic change in the pathways involved in meiotic recombination, with loss of the Dmc1-dependent pathway in the heterothallic species C. lusitaniae and C. guilliermondii (Supplementary text 14c). We also found that mechanisms of chromosome pairing and crossover formation have changed recently in these two species, because they (and to a lesser extent D. hansenii) have lost several components of the synaptonemal and synapsis initiation complexes (Supplementary table 35). They have also lost components of the major crossover formation pathway in S. cerevisiae (MSH4, MSH5), while retaining a minor pathway (MUS81, MMS4)38, 39. Overall, if Candida species undergo meiosis it is with reduced machinery, or different machinery, suggesting that unrecognized meiotic cycles may exist in many species, and that the paradigm of meiosis developed in S. cerevisiae varies significantly, even among yeasts.

The genome sequences reported here provide a resource that will allow current knowledge of C. albicans biology, the product of decades of research, to be applied with maximum effect to the other pathogenic species in the Candida clade. They also allow many of the unusual features of C. albicans – such as cell wall gene family amplifications, and its apparent ability to undergo mating and a parasexual cycle without meiosis – to be understood in an evolutionary context that shows that the genes involved in virulence and mating have exceptionally dynamic rates of turnover and loss.

Methods Summary

The full methods for this paper are described in Supplementary Information. Here we outline the resources generated by this project.

Assemblies, gene sets, and single nucleotide polymorphisms are available in GenBank and at the Broad Candida Database website (http://www.broad.mit.edu/annotation/genome/candida_group/MultiHome.html). The Broad website provides search, visualization, BLAST, and download of assemblies and gene sets. The C. parapsilosis assembly is also available at the Wellcome Trust Sanger Institute website (http://www.sanger.ac.uk/sequencing/Candida/parapsilosis/). The revised annotation of C. albicans (SC5314) is available at the Candida Genome Database (www.candidagenome.org). Gene families can be accessed on the Broad website via individual gene pages or by searching with family identifiers (CF#####).