Population study of C. tropicalis

The original reference genome sequence of C. tropicalis MYA-3404 was sequenced in 2009, resulting in a genome assembly consisting of 23 supercontigs totaling 14.6 Mb with 6,258 annotated genes [16]. We used Illumina data from resequencing of the reference strain to assemble the 23 supercontigs into 16 scaffolds, called Assembly B (see Materials & Methods). The assembly was subsequently further improved as described by Guin et al [39].

77 unique C. tropicalis isolates from different geographical locations were collected and sequenced using Illumina technology. For convenience, we named these strains ct01 to ct78, including only one of two isolates with very similar sequences (S1 Table). Most isolates came from clinical sources from the USA, Spain and Ireland. Twelve environmental isolates were included, eleven collected from soil or compost in the USA and Ireland, and one from coconut water in India. The reference strain C. tropicalis MYA-3404 (ct11), which was previously sequenced by Sanger sequencing [16], was also resequenced, as were three engineered auxotrophic derivatives in two genetic backgrounds [40,41].

Variants were identified by mapping reads to C. tropicalis MYA-3404 Assembly B and calling variants with the Genome Analysis Toolkit (GATK) [42]. Analysis of the distribution of allele frequencies in heterozygous biallelic SNPs showed that the majority of isolates are diploid, i.e. the ratio of reference to non-reference allele frequency is 50:50. However one isolate, C. tropicalis ct66 is triploid (peaks of allele frequency at 0.33 and 0.66), and another isolate, C. tropicalis ct26, appears to be octaploid (peaks of allele frequency at approximately 0.5, 0.12 and 0.87) (S1 Fig). In addition, we observed single-chromosome aneuploidies in four isolates (S1 Fig). C. tropicalis ct06 and C. tropicalis ct18 each have three copies of scaffold 8, and C. tropicalis ct14 and C. tropicalis ct15 both have three copies of scaffold 4 (trisomy). C. tropicalis ct14 (CAY3764) and C. tropicalis ct15 (CAY3763) were both derived from C. tropicalis AM2005/0093 and were used as the background to generate gene deletions [41]. Generating gene deletions has been found to induce aneuploidies in C. albicans [43].

Smaller copy number variants (CNVs) were identified in several scaffolds. The largest of these is a duplication of approximately 253 kb on scaffold 7 of isolates C. tropicalis ct04 and ct33, from approximately 350 kb to 603 kb (S2 Fig). Several shorter CNVs were also identified. These include a reduction in copy number of a 35 kb region from ~974 kb to 1 Mb on scaffold 4 in seven isolates: ct12, ct14, ct15, ct26, ct33, ct36 and ct69. Three large copy number variants (ranging from 23 to 235 kb in length) described by Guin et al. [39] in chromosomes 4, 5, and R (i.e. scaffolds 3, 8 and 4, respectively) were not observed in the C. tropicalis isolates described here.

Most isolates have approximately 2–6 heterozygous variants (including SNPs and indels) per kilobase similar to the type strain [16]. This is comparable to the level of heterozygosity seen in C. albicans (2.5–8.6 SNPs per kilobase) [25,26]. One isolate (C. tropicalis ct20) is extremely homozygous, with 0.84 heterozygous variants per kilobase. This isolate also has a higher proportion of homozygous variants compared to the reference (83% of total variants are homozygous, compared to an average of 41% in other isolates). However, six isolates have exceptionally high levels of heterozygosity (Fig 1A). These include one clinical isolate from Spain (C. tropicalis ct25), and five environmental isolates from soil, one from the USA (C. tropicalis ct42) and four from Ireland (C. tropicalis ct75, ct76, ct77 and ct78. Ct77 and ct78 were isolated from the same soil sample). These isolates have 36–49 heterozygous variants per kilobase. Phylogenetic analysis shows that most isolates cluster together (Cluster A in Fig 1B). However, the six heterozygous isolates are extremely divergent (Cluster B, Fig 1B). These six isolates separate into two groups, one containing C. tropicalis ct25, ct42, ct75 and ct76, and a second containing C. tropicalis ct77 and ct78.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Identification of novel isolates of C. tropicalis.

(A) Genome variation among C. tropicalis isolates. Variants were identified using the Genome Analysis Toolkit HaplotypeCaller and filtered based on genotype quality (GQ) scores and read depth (DP). Variants for all 77 isolates are shown according to variant type. Isolates are labelled on the X-axis by strain ID. One isolate (C. tropicalis ct20) has mostly homozygous variants, and six isolates have very high levels of heterozygous variants. (B) Six isolates of C. tropicalis are highly divergent. Variants were called as in (A). For heterozygous SNPs, a single allele was randomly chosen using RRHS [93] and for homozygous SNPs, the alternate allele to the reference was chosen by default. This process was repeated 100 times and 100 SNP trees were drawn with RAxML using the GTRGAMMA model [94]. The best-scoring maximum likelihood tree was chosen as a reference tree and the remaining 99 trees were used as pseudo-bootstrap trees to generate a supertree. Pseudo-bootstrap values are shown as branch labels. The six divergent isolates (Cluster B) are labelled according to their country of origin (see 1C). (C) SNP phylogeny of isolates from Cluster A indicates that clade structure is not associated with geography. The phylogeny of cluster A is shown in detail. Pseudo-bootstrap values are shown as branch labels. Isolates are labelled according to their country of origin, and environmental isolates are indicated with an asterisk. The reference strain, C. tropicalis MYA-3404, is labelled. Five putative clades are highlighted with colored bubbles. These clades are supported by principal component analysis (PCA) (S3 Fig). A sixth group was also identified by PCA, encompassing the remainder of isolates in the tree (S3 Fig).

https://doi.org/10.1371/journal.ppat.1009138.g001

The remaining isolates (Cluster A) are shown in more detail in Fig 1C. There is evidence of some population structure, with at least five well-supported clades identified by principal components analysis (colored ovals in Figs 1C and S3 and S4 Table) and many lineages outside these clades. However, there is little obvious correlation between phylogeny and geography. Two clades contain only isolates from the USA, but this likely reflects the overrepresentation of isolates from the USA in our collection. In addition, although some of the environmental isolates cluster together, others are closely related to clinical isolates (Fig 1C). There is therefore no clear distinction between clinical and environmental isolates.

Origins of the heterozygous C. tropicalis isolates

The levels of heterozygosity in the six divergent C. tropicalis isolates are similar to those observed in the hybrid species C. metapsilosis and in hybrid isolates of C. orthopsilosis [26–28]. This suggests that these C. tropicalis isolates may also be hybrids, that is, they may have at least one different parent to most C. tropicalis isolates. Hybrid genomes are characterized by regions of heterozygosity, resulting from differences between the homeologous chromosomes, alternating with regions of homozygosity. This results in distinct bimodal patterns of subsequences (k-mers) in sequencing reads, which represent the heterozygous and homozygous regions of the genome. Such bimodal k-mer patterns are characteristic of heterozygous genomes and have been observed in hybrid isolates of C. orthopsilosis, C. metapsilosis, C. inconspicua and C. albicans [25,38].

We find that the k-mer frequency distribution of four of the six divergent C. tropicalis is also bimodal, with one peak at approximately 100X (the average genome-wide coverage) and one at approximately 50X (half the average genome-wide coverage) (Fig 2A). The full and half coverage peaks represent homozygous regions and heterozygous regions respectively. Approximately half of the heterozygous k-mers (i.e. k-mers that map to heterozygous regions of the genome) are not represented in the reference genome sequence, because it is a collapsed haploid reference sequence from a non-hybrid isolate (C. tropicalis MYA-3404). For the remaining two divergent isolates (C. tropicalis ct25 and ct42), a bimodal k-mer distribution was not observed, possibly because the sequence coverage was too low. In these two isolates, the peak of k-mer multiplicity was less than 20, lower than any other isolates, which may obscure the signal (S4 Fig). This analysis suggests that at least four of the divergent isolates are hybrids, resulting from mating between two related, but distinct, parents. For all four isolates, the heterozygous peak is considerably higher than the homozygous peak, indicating that the hybridization event(s) are recent, and very little loss of heterozygosity (LOH) has occurred.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Novel C. tropicalis isolates result from hybridization.

(A) Analysis of k-mer distribution profiles reveals hybrid genomes. K-mer analysis of sequencing readsets was performed with the k-mer Analysis Toolkit (KAT [82]). For each of four divergent isolates, the number of distinct k-mers of length 27 bases (27-mers) is displayed on the Y-axis and k-mer multiplicity (depth of coverage) is displayed on the X-axis. K-mers that are present in the reference genome are shown in red, and k-mers that are absent from the reference genome are shown in black. There are two distinct peaks of k-mer coverage at approximately 50X and 100X. This pattern implies that most of the genomes are heterozygous (k-mers at 50X coverage) with few homozygous regions (k-mers at 100X coverage). Approximately half of the heterozygous k-mers in the readsets are not represented in the reference sequence. This pattern has been observed in hybrid isolates from other yeast species [25]. (B) Analysis of phased variants identifies two distinct haplotypes in divergent isolates of C. tropicalis. Variants were phased using HapCUT2 [44] into blocks covering 10–13 Mb of the genome. For each phased block, percentage difference from the reference strain in each haplotype was calculated as the number of variants divided by the length of the block. For 84–87% of the blocks, one haplotype is <0.3% different to the reference sequence and one haplotype is >1% different to the reference sequence. All phased blocks for each of the six hybrid isolates are shown as haplotype pairs, with the member of the pair more similar to the reference (haplotype A) shown in blue and the member of the pair less similar to the reference shown in orange (haplotype B) or purple (haplotype C). Percentage difference to the reference sequence is displayed on the Y-axis and position in the genome (chromosome, position (bp)) is displayed on the X-axis.

https://doi.org/10.1371/journal.ppat.1009138.g002

To further investigate the origins of the six divergent isolates, we attempted to separate the haplotypes of the two parental chromosomes. Approximately 500,000–700,000 heterozygous sites were identified per isolate. The heterozygous sites were placed in phased blocks, using HapCUT2 [44]. On average, 86% of the variants in each isolate were successfully phased, with a total phased span in base pairs of approximately 10–13 Mb (Table 1). The proportion of phased variants was slightly higher in isolates that were sequenced using longer read lengths. For example, 88 and 86% of variants were phased in C. tropicalis ct25 and ct42, which were sequenced using 250 bp fragments, compared to 85% in the other isolates, which had read lengths of 150 bp. Longer read lengths likely facilitate more accurate mapping across short homozygous regions.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Results of haplotype phasing.

https://doi.org/10.1371/journal.ppat.1009138.t001

Variants in the phased genomic regions or “blocks” were assigned to one of two haplotypes. For each phased block greater than 1 kb, the difference of each haplotype in that block to the reference sequence was calculated as the number of variants assigned to the haplotype divided by the length of the block. For the majority of blocks (84–87%), one haplotype (which we refer to as haplotype A) has >99.7% identity to the reference and the second haplotype is more than 1% different to the reference (Fig 2B). The alternative haplotypes were constructed by substituting all variant sites in the reference sequence with alleles that had been assigned to the alternative haplotype. The alternative haplotypes of all six isolates are 4.0–4.6% different from the reference strain. The alternative haplotypes of four of these isolates, C. tropicalis ct25, ct42, ct75 and ct76, which we refer to as haplotype B, are approximately 1% different from each other. The alternative haplotypes of the other two, C. tropicalis ct77 and ct78, called haplotype C, are approximately 3% different in sequence to the B haplotypes in the other four isolates (and less than 1% different in sequence from each other).

These analyses strongly suggest that the six novel isolates originated from mating or hybridization between related parents, one of which is very similar to the C. tropicalis reference, and others that are > 4% different. The second parent is not the same for the six divergent isolates. We therefore refer to most C. tropicalis isolates as AA diploids, to four isolates as AB diploids, and to two isolates as AC diploids. All AB and AC isolates contain only one rDNA locus (D1/D2 region), which is 99% identical to the reference haplotype A. The rDNA sequences were confirmed by PCR amplification and Sanger sequencing (GenBank accession numbers MW584905—MW584910).

Loss of heterozygosity (LOH) analysis in C. tropicalis hybrid isolates suggests three parental strains

Loss of heterozygosity (LOH) describes tracts of the genome that are essentially homozygous, most likely due to gene conversion or mitotic recombination. We observe a pattern of heterozygous regions alternating with homozygous (LOH) regions in all C. tropicalis isolates (Fig 3A). We defined heterozygous regions of the genome as regions of at least 100 bp in length containing at least two heterozygous variants; all remaining regions of the genome were classified as homozygous, or LOH, regions, as long as they were at least 100 bp in length.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Loss of heterozygosity in C. tropicalis isolates.

(A) Hybrid and non-hybrid isolates differ in the extent of LOH across the genome. The eight largest scaffolds in the reference genome are displayed horizontally from left to right and labelled from 1 to 8. LOH blocks are shown in pink and heterozygous (“HET”) blocks are shown in green. Centromere positions are indicated with “C”, telomere positions are indicated with “T” and the rDNA locus is indicated with “R”. Isolates are labelled on the left-hand side. The re-sequenced reference strain C. tropicalis MYA-3404 (labelled as “Ref”) is shown as a representative of the non-hybrid (AA) isolates. The genomes of the AA isolates consist mostly of LOH blocks. The AA isolate C. tropicalis ct20 has undergone extensive LOH, covering >99% of the genome. In contrast, in the AB/AC isolates, the majority of the genome consists of heterozygous blocks. (B) LOH is limited to short tracts of the genome in hybrid isolates. The histograms show the frequency of LOH blocks of different lengths in the six hybrid isolates and two AA (non-hybrid) isolates the re-sequenced reference strain C. tropicalis MYA-3404 (labelled as “Ref”) and C. tropicalis ct20. Frequency is shown on a log scale on the Y-axis while length in kilobases (kb) is shown on the X-axis, with a bin width of 1000 bp. The average length of LOH blocks in the hybrid isolates ranges from 286–416 bp. A similar pattern is observed in all six hybrid isolates, i.e. a predominance of short LOH blocks, with very few long tracts of LOH. In the non-hybrid isolates (e.g. C. tropicalis MYA-3404), LOH blocks are generally longer. C. tropicalis ct20 has the longest average LOH block length (~10 kb).

https://doi.org/10.1371/journal.ppat.1009138.g003

Only 4% on average of the non-hybrid (AA) genomes are heterozygous. Heterozygous regions in the AA genomes have a mean length of 208 bp and a maximum length of approximately 7.6 kb (S5 Table). In C. tropicalis ct20 only 0.37% of the genome is heterozygous, with a mean block length of 213 bp. In contrast, on average, 69% of the six hybrid genomes consists of heterozygous regions, with a mean length of approximately 900 bp, and a maximum length of approximately 13.8 kb.

Analysis of heterozygous regions in the six hybrid isolates reveals further support for the hypothesis that they originated from different hybridization events involving different parent strains (B and C). If we assume that the hybrid isolates were derived from a mating event between two parental isolates, we can expect that the heterozygous regions of the genome in the hybrid isolates should be derived equally from the two parent strains. Therefore, if two hybrids originated from hybridization between the same parental strains, the heterozygous regions of their genomes should carry the same variants. However, if two hybrids originated from hybridization between different parental strains, the variants in common heterozygous regions will be different. Shared heterozygous regions were defined as regions of heterozygosity that are common to all isolates. For partially shared heterozygous regions, the portion that was common to all isolates was extracted and analyzed. Shared heterozygous regions in all six hybrid isolates cover 5.8 Mb, with approximately 45% of all heterozygous positions (totaling 485,018 variants) in these regions present in all six. The relatively low proportion of shared variants observed indicates that the six hybrid isolates did not all originate from the same parental strains. However, there is a much higher degree of conservation of variants in shared heterozygous regions among the four AB isolates; 94% of 419,440 heterozygous variants in 6.7 Mb are present in all four. Similarly, the two AC hybrids share 98% of 620,569 variants across 9.6 Mb. This further indicates (in line with our previous analyses) that the four AB isolates share a common origin, and that the two AC isolates share a common origin that is separate from the origin of the AB isolates.

There is extensive LOH in the non-hybrid isolates, covering on average 95% of the genome (S5 Table). In C. tropicalis ct20, >99% of the genome is in LOH blocks. The average length of LOH blocks across all non-hybrid isolates (excluding C. tropicalis ct20) is approximately 1.8 kb with a maximum length of 238 kb. In contrast, limited LOH is observed in the six hybrid (AB/AC) isolates, with an average of 13,139 LOH blocks of at least 100 bp, covering between 25 and 42% of the genome. The average length of LOH blocks in the AB/AC isolates is 330 bp, but can be as long as 112 kb (Fig 3B). Only 1.6% of LOH blocks (equating to 731 LOH blocks or 0.5% of all LOH length in bases) is conserved among all six isolates. There are more shared LOH regions in the four AB isolates; 17% of LOH blocks (equating to 5,131 LOH blocks) in these isolates are identical. In the AC isolates, 55% of LOH blocks are identical (equating to 8,807 LOH blocks). There is a large LOH block at the start of scaffold 4 (equivalent to Chromosome R [39]) covering approximately 400 kb, that is shared between four of the hybrid isolates (C. tropicalis ct25, ct42, ct77 and ct78). The LOH block extends from the telomere to the rDNA locus, although the exact end point differs, and it is interrupted by some small heterozygous regions. A larger LOH block, encompassing this region and extending to the centromere, was identified in a complete, chromosome-scale assembly of C. tropicalis and in the related species Candida sojae [39]. Two of the AB hybrids (C. tropicalis ct75 and ct76) are different, in that only the rDNA locus itself has undergone LOH. The same results were obtained when comparing to the updated C. tropicalis reference genome from Guin et al. [39] (S1 Text and S5 Fig and S6 Table).

We considered the possibility that the homozygous isolate C. tropicalis ct20 might represent one parent of the hybrid isolates. We therefore compared it with both haplotype A and haplotypes B and C of the six hybrid isolates by computationally reconstructing both subgenomes of each hybrid strain. We constructed a putative A haplotype from C. tropicalis ct20 by substituting bases in the reference with homozygous variants identified in this isolate. For the hybrid isolates, the A haplotype was constructed by substituting variants that were originally assigned to haplotype A during haplotype phasing (see Materials & Methods, subsection Haplotype splitting). Similarly, B and C haplotypes were constructed by substituting variants that were assigned to either B or C. The A haplotypes from the hybrids share, on average, approximately 8% of variants with C. tropicalis ct20 (i.e. approximately 8% of variants identified in C. tropicalis ct20 and a given hybrid isolate are identical). There is even less similarity between the B and C haplotypes and C. tropicalis ct20; only 1% of variant sites in C. tropicalis ct20 and the hybrid haplotypes B or C are identical. C. tropicalis ct20 therefore has an A haplotype, but it is unlikely that it is a parent, or closely related to a parent, of the hybrid isolates.

Mating type-like loci (MTL) in C. tropicalis isolates

Most AA isolates (45) are heterozygous at the MTL, similar to previous reports [20,21] (S1 Table). Fifteen are homozygous for MTLa/a and seven are homozygous for MTLα/α. Three isolates have three copies of the MTL. The triploid isolate C. tropicalis ct66 is MTL a/a/α. C. tropicalis ct18 is trisomic for scaffold 8, which carries the MTL, and is MTL a/α/α. C. tropicalis ct06 is also trisomic for scaffold 8, and has three copies of MTLα. The MTL idiomorphs of the octaploid isolate, C. tropicalis ct26, could not be definitively determined by assembling the Illumina data or by PCR, but it appears to have 7 copies of MTLα and one copy of MTLa (S6 Fig).

All six AB and AC isolates contain both MTLa and MTLα idiomorphs. In the AB isolates, the MTLa idiomorphs are >99% identical to that of the reference strain (haplotype A) with only three nucleotide changes across the entire locus (8,180 bp). These include synonymous and nonsynonymous substitutions in PAPa and PIKa. In addition, one isolate (C. tropicalis ct42) has a nonsynonymous substitution in MTLa1. Apart from this, the MTLa idiomorphs in the AB isolates are identical. The MTLa idiomorph therefore likely originated from the A parent. The MTLα loci are >99% identical in all four AB isolates, and ~7% different to the reference strain, indicating that it was donated by the B parent. All AB isolates therefore most likely resulted from mating between the same parents, an MTLa parent similar to the reference strain (parent A), and an MTLα parent which is approximately 4% different (parent B).

In the two AC isolates, the MTLα idiomorphs are also identical to each other, and they are >99% identical to the reference strain. The MTLa idiomorphs are identical to each other, and approximately 96% identical to the reference strain. The MTLa idiomorph in the AC isolates therefore originated from the C parent, and the MTLα idiomorph originated from the A parent.

Analysis of phenotypic variation in C. tropicalis

To measure the phenotypic diversity within C. tropicalis, the growth of 68 AA isolates was tested in 61 different conditions, including alternative carbon sources, stressors (e.g. calcofluor white, congo red), heavy metals (e.g. zinc, cobalt, cadmium) and antifungal drugs (e.g. fluconazole, ketoconazole, caspofungin) (S7 Fig). Because nitrogen and carbon metabolism are important virulence attributes in fungi [45], the ability of C. tropicalis isolates to use different sole nitrogen sources (e.g. amino acids, gamma-aminobutyric acid (GABA)) was also tested (S7 Fig). The AB and AC isolates and the engineered lab isolates C. tropicalis ct13, ct14 and ct15 were excluded from the analysis.

The C. tropicalis isolates show wide variation in their growth characteristics (S7 Fig). We attempted to identify genome variants that are associated with specific growth defects. For this analysis, only conditions that resulted in a growth defect of at least 70% compared to the control condition in at least one strain were included (i.e. 25 conditions using YPD as a base media, and 10 conditions using different nitrogen sources). Reduced growth was scored as 1, and growth similar to the control was scored as 0. Predicted genomic variants were annotated with SnpEff [46] to identify those that were likely to have a major impact on protein function. 390,321 variant sites were identified in total across 68 isolates. The majority of variants (~75%) were SNPs, with the remainder consisting of small insertions and deletions (indels) (S8 Fig). Most variants are found in intergenic regions, or are silent or missense mutations. Only variants that were predicted to have a high impact, including frameshifts, gene fusion events, loss or gain of a stop codon, or variation at splice donor or acceptor sites (9,261 variants, S8 Fig), were included in the genotype-phenotype correlation analysis.

One clinical isolate, C. tropicalis ct04, identified by cosine similarity analysis [47], has impaired growth when valine or isoleucine (branched chain amino acids) are provided as the sole nitrogen source (Fig 4A). Compared to other isolates, C. tropicalis ct04 also grows poorly on 2% sodium acetate, 2% starch and in the absence of a carbon source. There are 40 variants unique to this isolate that are predicted to have a high impact on protein function (S7 Table). One of these is a heterozygous deletion of two bases in CTRG_06204 (BAT22), an ortholog of the S. cerevisiae BAT1/2 genes that encode a branched-chain amino acid aminotransferase (BCAT). Amino acid metabolism and acquisition can directly affect pathogenicity, and disruption of branched-chain amino acid metabolism is associated with loss of virulence (reviewed in [48]). BCATs catalyze the final step of biosynthesis and the first step in the degradation of the branched chain amino acids valine, isoleucine and leucine [49]. The deletion results in a frameshift which introduces a premature stop codon at amino acid Gly30 of the Bat22 protein (Fig 4B). We determined if introducing an equivalent change into other genetic backgrounds using CRISPR/Cas9 [50] would result in the same phenotype. A repair template was designed to delete two bases and also to destroy the target of the guide RNA to prevent recutting. The gene was edited in three different C. tropicalis isolates ct09, ct44 and ct53. All edited strains can no longer use valine or isoleucine as sole nitrogen sources (Fig 4C). However, unlike C. tropicalis ct04 they have no growth defect on sodium acetate, starch or in the absence of carbon sources, indicating that another variant, or combination of variants, is responsible for these phenotypes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Disrupting BAT22 prevents growth of C. tropicalis on branched chain amino acids as a sole nitrogen source. (A) Phenotype analysis of C. tropicalis isolates.

Growth of C. tropicalis ct04 is shown on solid media. Strains were grown in 2x2 arrays; two biological replicates (top and bottom rows), with two technical replicates each (left and right columns), of each strain were tested. C. tropicalis ct04 replicates are outlined with red boxes. C. tropicalis ct04 cannot utilize valine or isoleucine as a sole nitrogen source and also exhibits a growth defect on solid media with 2% starch or 2% sodium acetate as the sole carbon source, or on solid media without a carbon source provided. (B) Editing of BAT22. Plasmid pCT-tRNA-BAT22 was generated to edit the wild type sequence of BAT22 (CTRG_06204) using CRISPR-Cas9. The sequences of the reference C. tropicalis BAT22 (CtBAT22 (wt)), BAT22 from C. tropicalis ct04 (CtBAT22 (ct04)) and edited BAT22 (CtBAT22*) are shown. The guide sequence is highlighted with a black box, the PAM sequence is shown in bold, and the Cas9 cut site is indicated with a red scissors. C. tropicalis isolates ct44, ct09 and ct53 were transformed with pCT-tRNA-BAT22 and a repair template (RT_BAT22_2bpDel_SNP) generated by overlapping PCR using RT_BAT22_2bpDel_SNP-TOP/BOT oligonucleotides. The repair template contains two 60 bp homology arms and deletes two bases in BAT22 resulting in the same frameshift observed in C. tropicalis ct04. (C) Edited strains have defects in branched-chain amino acid metabolism. 5-fold serial dilutions of C. tropicalis ct04, ct09(wt; bat22*), ct44 (wt; bat22*) and ct53 (wt; bat22*) in the same conditions tested in (A). The edited strains cannot use valine or isoleucine as sole nitrogen sources.

https://doi.org/10.1371/journal.ppat.1009138.g004

Analysis of genotype-phenotype correlations becomes more challenging when analyzing variants in intergenic regions or variants with a predicted low or intermediate impact. By restricting our analysis to only variants that were predicted to have a high impact, we were able to narrow our search to a small number of variants, which could be manually validated. We also identified ~69,000 variants in upstream regions (occurring within 300 bp upstream of an annotated gene), which is too many to verify experimentally.

Strain collection and growth

C. tropicalis isolates were collected from a variety of clinical (anonymized) and environmental sources (S1 Table). Yeast were isolated from soil samples as described in Sylvester et al. [79], or following three passages in YPD (1% yeast extract, 2% peptone, 2% glucose) with chloramphenicol (3% [wt/vol]) and ampicillin (10% [wt/vol]). For phenotype analysis, isolates were inoculated as 2x2 arrays (two independent cultures with one technical replicate of each) into 200 μl of YPD broth in 96-well plates and incubated at 30°C for 24 h. Stocks were diluted in 96-well plates containing 200 μl of water by dipping a 12x8 pin bolt replicator (V&P Scientific) three times in the culture and then transferring it to the water. Once diluted, the cultures were pinned onto 85 unique media on solid agar plates and incubated at 30°C for 48 h (S2 Table). For 60 conditions, the base media was YPD, with 2% agar including 2% glucose as a carbon source. Glucose was substituted with different carbon sources where indicated, or compounds were added at the indicated concentrations (S2 Table). To test the ability to use specific nitrogen sources (24 conditions), the base media was 0.19% of YNB (Yeast Nitrogen Base) without ammonium sulfate or amino acids, 2% glucose and 2% agar. Nitrogen sources were added as indicated (S2 Table). Spider media was tested as the 85th condition (S2 Table). Plates were photographed and growth was measured using SGAtools [80]. SGAtools was designed to analyze synthetic genetic interactions and assumes that average growth on a plate does not vary. This was not true for several media, where many strains grew poorly. We therefore compared the growth of each strain on the test media to the growth of the same strain on YPD, or on YNB with ammonium sulfate, as a control, using the raw data extracted from SGAtools. For each strain in each analyzed growth condition, the SGAtools scores (ranging from 0 to 1.8) were converted to a binary score where a growth ratio above 0.3 (no growth defect) was assigned 0, and a ratio below or equal to 0.3 (major growth defect) was assigned 1. These scores were chosen to be very stringent—only conditions which resulted in reducing growth to approximately 30% of that under the control conditions were judged as a defect. We found that SGAtools could not reproducibly identify enhanced growth in these conditions. The raw data for the image analysis is available at https://doi.org/10.6084/m9.figshare.13128839.v1.

Genome sequencing

For most C. tropicalis isolates, genomic DNA was isolated by phenol-chloroform extraction followed by purification using the Genomic DNA Cleanup and Concentration kit from Zymo Research (catalogue number D4065). For three isolates (C. tropicalis ct76, ct77 and ct78), genomic DNA was extracted and purified using the QIAamp DNA Mini Kit from QIAGEN (catalogue number 51304). For most isolates, library preparation and sequencing was performed at the Earlham Institute, Norwich, UK using the LITE method (Low Input Transposase-Enabled), a custom Nextera-based system. These isolates were sequenced on two lanes of an Illumina HiSeq 2500 generating 2x250 bp paired-end reads. For five isolates (C. tropicalis ct51, ct75, ct76, ct77 and ct78), library preparation and sequencing was performed by BGI, Hong Kong, generating 2x150 bp paired-end reads, on an Illumina HiSeq 4000. Our genome sequences of two isolates (C. tropicalis ct20 and ct21) were almost identical. These may represent independent isolates of the same strain, or one isolate may have been accidentally sequenced twice. We therefore included only one of these (C. tropicalis ct20) in subsequent analysis.

For the 72 unique isolates sequenced using the LITE method, Nextera adapters were removed using TrimGalore v0.4.3 with the parameters “—paired” “—length 35” “—nextera” and “—stringency 3”. Custom adapters and low-quality bases were trimmed using Skewer v0.2.2 with the parameters “-m pe” “-l 35” “-q 30” “-Q 30” [81]. For 5 isolates sequenced by BGI, adapters were removed by the sequencing provider and reads were quality trimmed using Skewer. K-mer distribution profiles were analyzed using the k-mer Analysis Toolkit v2.4.2 using the default k-mer length of 27 bases [82]. All genomes were assembled using SPAdes v3.9.1 with parameters “—careful” “-t 12” “-m 60” [83]. Assembly statistics were assessed using QUAST v4.4 [84]. To confirm the species identity of hybrid isolates, the D1/D2 domain of the large subunit of the ribosomal DNA was amplified using standard universal primers NL-1 and NL-4 (S3 Table).

Mating type-like locus analysis

The MTL idiomorph of a subset of isolates was confirmed by PCR using primer pairs MTLa1F and MTLa1R to amplify the MTLa1 gene and MTLα2F and MTLα2R to amplify the MTLα2 gene, as described in Xie et al. [21]. Colony PCR was performed by boiling single colonies in 5 μl sterile deionized water, then adding 12.5 μl MyTaq Red Mix (2X), 1 μl forward primer (100 μM), 1 μl reverse primer (100 μM) and 5.5 μl deionized water. PCR was run for 1 min at 95°C; then for 30 cycles of 30 sec at 95°C, 30 sec at 57°C, 60 sec at 72°C; and then a final 2 min at 72°C.

C. tropicalis reference genome

The C. tropicalis reference genome annotation was updated using RNAseq data for three C. tropicalis strains downloaded from NCBI under BioProject ID PRJNA290183 [85]. RNAseq data were aligned against the original C. tropicalis reference [16] with HISAT2 v2.0.5 with the parameter “—novel-splicesite-outfile” to predict splice sites in the genome [86]. Predicted splice sites were manually validated by examination of transcripts mapping to predicted splice sites. The reference genome sequence was subsequently scaffolded from 23 supercontigs to 16 supercontigs. Areas of overlap between supercontigs in the original reference assembly were identified using Gepard to generate dot matrix plots [87]. Overlapping supercontigs were merged if this arrangement was supported by synteny with other Candida species, using the Candida Gene Order Browser (CGOB) [78], and by data from Illumina resequencing of the reference strain. The final assembly (also known as Assembly B [39]) contained 16 supercontigs and is available under NCBI accession JAFIQD000000000. The C. tropicalis reference was subsequently further improved as described by Guin et al [39].

Variant calling

For isolates sequenced using the LITE method, trimmed reads were aligned to C. tropicalis MYA-3404 Assembly B with bwa mem v0.7.11 to generate two BAM files per sample (one for each lane used for sequencing) [88]. BAM files were sorted with SAMtools v1.7 [89], and duplicate reads were marked using GenomeAnalysisToolkit (GATK) v3.7 Mark Duplicates [42]. BAM files from separate lanes were combined for each sample and marked for duplicates again using GATK MarkDuplicates. For isolates sequenced at BGI, Hong Kong, trimmed reads were aligned to the updated C. tropicalis MYA-3404 Assembly B with bwa mem v0.7.11 as before, generating only one BAM file per sample (each of these samples was sequenced on only one lane of the sequencer). BAM files were sorted with SAMtools v1.7 [89] and duplicate reads were marked using GenomeAnalysisToolkit (GATK) v3.7 Mark Duplicates [42].

The subsequent steps were applied to all samples. Realignment around indel sites was performed using GATK IndelRealigner and variants were called using GATK HaplotypeCaller in “—genotyping_mode DISCOVERY”. Variants were filtered for quality based on genotype quality (GQ) < 20 and read depth (DP) < 10. For SNP trees, gVCFs were generated using GATK HaplotypeCaller with the parameters “—genotyping_mode DISCOVERY” and “—emitRefConfidence GVCF”. Joint genotyping was performed using GATK GenotypeGVCFs to produce a single multi-sample gVCF. SNPs were extracted from the multi-sample gVCF using GATK SelectVariants with parameter “-selectType SNP”. Variants were filtered based on genotype quality (GQ) < 20 and read depth (DP) < 10. For genotype-phenotype analysis, the presence of a variant at a particular site in each isolate was scored as 1, and absence was scored as 0.

Aneuploidy analysis

To calculate copy number variants based on coverage discrepancies, the C. tropicalis MYA-3404 Assembly B genome was split into 1 kb windows using the “makewindows” command from bedtools v2.26.0, with parameters “-i winnum” (label windows sequentially) “-w 1000” (window size 1 kb) [90]. Mean coverage in each 1 kb window was calculated for each sample using the “coverage” command from bedtools [90]. Average whole genome coverage for each strain was calculated using GATK DepthOfCoverage [42]. Coverage ratios for each 1 kb window were calculated as log₂(window coverage/average whole genome coverage). A value of zero was assigned to windows that had zero coverage. The resultant ratios were visualized using the DNACopy package from Bioconductor in R [91]. Ploidy was also visualized using allele frequencies from heterozygous biallelic SNPs extracted from the VCF files using GATK SelectVariants with parameters “-selectType SNP” and “-restrictAllelesTo BIALLELIC”. Allele frequency was calculated as allele depth (AD) /read depth (DP). Histograms of allele frequency for each scaffold in each sample were visualized in R using ggplot2 [92].

Phylogeny

SNP trees were drawn from filtered variants, using only those SNPs that passed the filters described in “Variant Calling”. To account for heterozygous SNPs, the Repeated Random Haplotype Sampling tool (RRHS) v1.0.0.2 was used to select a random allele at heterozygous SNP sites [93]. This process was performed 100 times to generate 100 SNP profiles for each isolate, thereby encapsulating the full heterozygosity of each isolate. For homozygous variant sites, the alternate allele was chosen by default. 100 maximum likelihood (ML) trees were drawn (one for each SNP profile) using RAxML v8.2.12 [94] with the “GTRGAMMA” model. The best-scoring ML tree was chosen as a reference tree and the remaining 99 ML trees were used as pseudo-bootstrap trees to generate a supertree using RAxML v8.2.12 with options “-f b” (draw bipartition information on a reference tree based on multiple trees (e.g. from a bootstrap)) and the “GTRGAMMA” model. Phylogeny was also examined using principal component analysis (PCA) with the ade4 package in R [95].

Loss of heterozygosity

Loss of heterozygosity (LOH) was calculated in blocks of at least 100 base pairs (bp) across the genome. Heterozygous regions were defined as any region containing at least two heterozygous variants within 100 bp of each other, with a minimum total length of 100 bp. Remaining regions were defined as homozygous, or LOH, regions as long as they were at least 100 bp in length. A similar approach has been used to characterize LOH in other Candida hybrids, e.g. C. metapsilosis [27], C. inconspicua [38] and C. albicans [25]. Commonality of variants in heterozygous regions was examined to determine which isolates originated from the same parental strains. Heterozygous regions shared by all isolates were identified using bedops intersect [96]. In the case of heterozygous regions that were partially shared, the portion that was common to all isolates was extracted and analyzed as a shared heterozygous region. The number of common variants in the shared heterozygous regions was counted as the number of variant sites in these regions with the same genotype in all isolates. Shared LOH regions were defined as LOH blocks with identical start and stop coordinates in the relevant isolates. The analysis was repeated using the updated C. tropicalis genome assembly from Guin et al. [39].

Haplotype splitting

Hybrid haplotypes were phased using HapCUT2 v0.7 [44]. The filtered variants were used as input for the subcommand “extractHAIRS” (extract haplotype-informative reads) to identify “haplotype-informative reads”, i.e. sets of reads that align to the same location in the reference genome but that contain one or more different alleles at variant sites. HapCUT2 was subsequently used to build haplotype blocks from the haplotype-informative reads with parameter “—threshold 30” (Phred-scaled threshold for pruning low-confidence SNPs). Haplotype-informative reads form the basis of these haplotype blocks. Two haplotype-informative reads will be assigned to the same haplotype if they overlap by a certain amount and have matching alleles at variant sites. Thereby, the haplotype block is extended in a continuous manner until no further overlapping reads can be identified. Each block consists of a region of the genome where alleles at variant sites in that region have been assigned to one of two haplotypes. Therefore, each phased block is assigned two sets of variants; one for each haplotype. For each phased block, the percentage difference of each of these two haplotypes was calculated by counting the number of bases in that haplotype that were different from the reference (i.e. the number of variant alleles assigned to the haplotype) and dividing this number by the total number of bases in the block. The difference of each phased block to the reference genome was calculated as the number of SNPs in block/length of block. Haplotypes were subsequently assigned to either the reference haplotype or the alternate haplotype according to their percentage difference; the member of the pair that was more similar to the reference was assigned to haplotype A and the member of the pair that was less similar to the reference was assigned to haplotype B. In the majority of cases, the A haplotype was < 0.3% different to the reference genome and the B haplotype was > 4% different to the reference genome.

Analysis of genotype-phenotype correlation

Variants from non-hybrid isolates were further annotated with SnpEff v4.3t to predict the functional effect of variants [46]. High-impact variants (e.g. variation at splice donor or acceptor sites, variants resulting in a gain or loss of stop or start codon, or frameshifts in genes) were extracted and correlated with phenotypes. Variants were converted to binary scores; 1 for the presence of a variant in a given strain, 0 for the absence. Phenotype scores were coded as 1 for a growth defect (score of 0.3 or less), and as 0 for no growth defect (score above 0.3). For each variant-condition pair, two vectors were generated using the binary scores; the first consists of the scores for every strain with respect to the variant, the second consists of the scores for every strain with respect to the condition. For every variant-condition vector pair, the cosine similarity between the two vectors was calculated as . Any variant-condition pair with a cosine similarity of > 0.85 was selected for further analysis.

Editing BAT22 with CRISPR-Cas9

A 20 bp sequence (guide RNA) targeting C. tropicalis BAT22 (CTRG_06204) was designed using the web tool ChopChop [97]. The guide RNA was generated by annealing of two short oligos (g60BAT22_TOP/BOT, S3 Table), and then cloned into the SapI-digested pCT-tRNA plasmid to generate plasmid pCT-tRNA-BAT22, as previously described in [50]. The repair template carrying the desired modification, including the disruption of the PAM sequence, was generated by primer extension (RT_BAT22_2bpDel_SNP-TOP/BOT) using ExTaq DNA polymerase (Takara Bio, USA). C. tropicalis isolates ct09, ct44 and ct53 were transformed with 5 μg pCT-tRNA-BAT22 and 25 μl of unpurified RT- BAT22_2bpDel_SNP using a previously described method [50]. Transformants were selected on YPD agar plates containing 200 μg/ml nourseothricin (NTC), incubated at 30°C for 48 h. The relevant region was amplified by PCR from two NTC-resistant transformants for each strain using primers bat22_fwd_01/bat22_rev_01 and sequenced using Sanger sequencing. The pCT-tRNA-BAT22 plasmid was cured by growing the cells in the absence of selection on YPD until they failed to grow in the presence of NTC.