Whole genomes from 118 Cryptococcus gattii isolates, representing both the Pacific Northwest (PNW) subtypes and the known global VGII molecular genotype diversity, along with isolates of each of the other major C. gattii molecular types (VGI, VGIII, and VGIV), were sequenced and assembled de novo, using data from multiple next-generation sequencers (all Illumina technology). We obtained at least 10× coverage across an average of 98.27% of the genome (mode, 98.64%; range, 90.63% to 99.43%) compared to the assembled R265 reference (VGIIa), with an average depth of 111× (range, 24× to 532×). The average assembly size was ~17.55 Mb (range, 17.36 to 17.83 Mb).

We determined the presence of single nucleotide polymorphisms (SNPs) across and between subtypes. A total of 1,282,876 SNPs were identified from all genomes, including the VGI, VGII, VGIII, and VGIV molecular types; 544,881 SNPs were parsimony informative, that is they distinguished between more than one strain and were not autapomorphic (Fig. 1). Within the 115 VGII genomes, 310,969 SNPs were found, with 221,248 being parsimony informative (Fig. 2). Consistent with recently emerged and low-diversity clones, only 717, 246, and 528 SNPs were identified within the VGIIa, VGIIc, and VGIIb populations, respectively (Fig. 3, Fig. 4, and Fig. 5, respectively). A combined analysis of only the three PNW clones identified little to no homoplasy (data not shown), as previously described (17).

Population-level genome analyses were initiated with an SNP-based assessment of phylogenetics and population structure, followed by gene presence/absence comparisons to understand the genetic relationships between and among the VGII subtypes in the PNW and other geographically diverse subtypes. A detailed description of the data analysis flowchart is shown in Fig. S1 in the supplemental material.

Isolates chosen for phylogenetic analysis were representative of the global diversity of C. gattii VGII based on their multilocus sequence type (MLST), including a number of isolates from each of the PNW lineages (VGIIa, VGIIb, and VGIIc). Additionally, isolates of the other three C. gattii populations were analyzed. Based on whole-genome SNP typing (WGST), the VGII strains were clearly separated from other molecular types (i.e., VGI, VGIII, and VGIV), with greater than 500,000 SNPs between VGII and any other major molecular type (Fig. 1).

Using maximum-parsimony analysis for WGST, the VGII isolates separated mostly along the lines of their MLST (Table 1; see Fig. S2 and S3 in the supplemental material). A number of isolates with distinct sequence types (STs) were clearly more related to each other and likely represent members of the same clone with single or limited nucleotide variation within the MLST markers (e.g., VGIIc includes ST6 and ST49). Multiple clonal lineages are seen in the whole VGII phylogeny, including the PNW lineages (Fig. 2 and Fig. S2), and the long terminal branches indicate substantial genetic distance between most lineages. The maximum-parsimony tree (Fig. S2) had a relatively low consistency index of 0.4, indicating either a high degree of homoplasy from multiple mutations at the same sites or that VGII is not evolving in a strictly tree-like manner. To better visualize the relationships between lineages, we developed a neighbor-joining phylogenetic split tree network, or “neighbor net” (Fig. 2). The network contained many phylogenetic splits, which could not be condensed into a single tree phylogeny. This finding may indicate that sexual recombination events between distinct subtypes are contributing to the genetic diversity of VGII, as these events will cause additional splits due to different regions of the genome being best described by different trees.

When the phi test for recombination was performed using the SNP data, the test indicated that recombination was present with a very high level of significance (P = 0.0). We then examined the phi value across the genome to see whether recombination could be localized to a particular subset of the genome. Using nonoverlapping windows of 500 concatenated SNP positions, we tested for recombination separately in 673 “windows” across the genome (encompassing variable genome window sizes due to the nonuniform density of SNPs in the genomes). Recombination was detected at a significant level in 631 of these windows when the P value cutoff was adjusted for making 673 comparisons (P < 7.2⁻⁵), indicating that the recombination signal was common across most of the genome and not localized to a single region.

Due to the evidence for recombination among isolates, a model-based Bayesian approach (24) for identifying population structure and individual assignment to populations was added to our phylogenetic analyses. fineStructure analysis placed the 32 distinct VGII subtypes, as identified by WGS and MLST (see Fig. S2 and S3 in the supplemental material), into 11 “populations” or groups related by shared genomic regions, with some of these containing only a single representative (Fig. 6). For several of these groups, members shared many genome regions with members from other groups, indicating incomplete or recent separation between groups. This finding is also supported by the presence of isolates, such as B9764 and WM 05.525, which apparently share genetic material with several populations and may represent admixed individuals. Importantly, all three of the PNW subtypes were placed into separate fineStructure “populations,” all of which include members from South America. The neighbor net network (Fig. 2), maximum parsimony-based WGST (Fig. S2), and MLST (Fig. S3) analyses all demonstrate broad dispersal of subtypes from South America throughout the VGII population.

The members of the PNW subtypes (i.e., VGIIa, VGIIc, and VGIIb) showed strict clonality within their respective subtype, with consistency indices (CI) of 1.0 for each subtype tree (Fig. 3 to 5), and significant separation from each other, as described previously (17). Bioinformatic attempts to determine a date for the most recent common ancestor of each subtype did not yield useful estimates (analysis not shown), which may be due to the small amount of SNP data available or potentially the nonclocklike evolution of C. gattii.

The VGIIa subtype contained two distinct clades, a North American-dominated group and a South American group (Fig. 3). The North American VGIIa clade showed little topological structure with a 15-branch polytomy consistent with a radiation emanating from a recent common ancestor. One subclade, joined by 13 shared SNPs, contained a set of three isolates previously described as having lower virulence (13, 20), one of which was isolated in Brazil in 1981 and was genomically separated from the others by 261 additional SNPs. The other two members of this lower-virulence subclade were isolated in Washington State in 1972 and California in 1990 and had SNP branch lengths of 0 and 52, respectively. The VGIIc population was also strictly clonal (CI = 1.0) with a total of 55 parsimony informative SNPs (Fig. 4). This subtype phylogeny contained a 5-branched polytomy but with several distinct subclades. The branch lengths were short (<16 SNPs); notably, all of these isolates were collected between 2008 and 2012.

Unlike VGIIa and VGIIc, the VGIIb subtype is dispersed beyond the Western Hemisphere (6, 13, 15, 25). This subtype contained 528 SNPs (147 parsimony informative) and multiple distinct clades (Fig. 5). North American isolates were found in multiple locations across the tree, including in one large derived clade where they shared membership with isolates that were collected in Asia and Australia. The Asian-Australian/North American subclade shows polytomy, which is consistent with the presence of a genetic bottleneck and subsequent radiation, or may be a reflection of limited sampling. A single isolate from Florida shared membership in a separate basal clade with an isolate from Australia.

Ab initio gene prediction using the pretrained Augustus software identified a median of 6,276 genes per genome (range, 6,213 to 6,359). The gene content analysis initially focused on differences and similarities between the PNW subtypes. The use of BLAT score ratio (BSR) for this analysis allowed us to make gene predictions from the genome sequences, where exon-coding DNA sequences were concatenated into single contiguous sequences to create the in silico coding DNA (is-cDNA), the presence or absence of which could then be scored for the PNW genomes. The gene sequences from the three PNW subtypes clustered into 6,302 is-cDNAs. The gene sequences from the remaining 29 non-PNW subtypes clustered into 6,541 is-cDNAs. The individual and shared gene differences in the PNW subtypes were scored against representative genomes from the other known global VGII subtypes. Additionally, we surveyed for the presence or absence of specific genes from any identified PNW clone (e.g., PNW-VGIIa) within the representative genomes from the other C. gattii major molecular types (VGI, VGIII, and VGIV), as well as within the C. neoformans var. grubii and C. neoformans var. neoformans reference genomes (Table 2). A heat map was developed (Fig. 7) to visualize the gene presence, absence, partial presence, and/or presence of genes with low sequence similarity. We identified several gene content similarities and differences in all of the comparisons (Table 2). Many of the identified putative genes, or is-cDNAs, were located in contiguous blocks of assembled genome sequence, labeled clusters 1 to 7 in Fig. 7. When the PNW clones were compared with each other, VGIIa had the largest number of PNW clone-specific gene content differences. Additionally, the PNW VGIIa-specific genes were typically absent in the other VGII subtypes and non-VGII genomes. Conversely, the PNW VGIIb-specific and VGIIc-specific genes (i.e., those genes found only in their respective subtype in the PNW) were commonly found in the other VGII subtypes and non-VGII C. gattii genomes. A large number of genes that were present in both VGIIa and VGIIc were absent from VGIIb. These genes, most of which were located adjacent to each other in clusters 1 and 2, were found in nearly all other C. gattii genomes, indicating a likely gene loss event in the PNW VGIIb population. An analysis of genes common to all three PNW populations but absent in other VGII subtypes identified a three-gene cluster (“cluster 6”) that was found to be absent in most C. gattii genomes but was present, with a high degree of similarity, in the C. neoformans var. grubii genome (Fig. 7).

A comparison of phylogenies based on linked SNPs and gene content showed dissimilar relationships among VGII subtypes (data not shown). This suggests that differential gene content has not become fixed within the VGII lineages. As such, placement of an isolate in a phylogeny or within a subtype using SNP data is not necessarily a good indicator of its gene content, and conversely, gene content is not a good indicator of phylogenetic placement.

Analysis of the population genome structure of the PNW VGII clones provides additional significant evidence of multiple and continued emergence of apparently regionally fit lineages from within VGII originating in South America, as has been proposed previously (9, 26). As the apparent level of recombination across VGII negates the use of a single SNP-based tree to understand the true phylogenetic relationships among all VGII members, we conducted a neighbor net analysis, providing a better understanding of relationships within the phylogenetic network. For example, the Australian ST33 subtype, which appears most genetically distant on the WGST SNP tree (see Fig. S2 in the supplemental material), still appears most genetically distant in the network; however, there are several distinct splits uniting this subtype with other subtypes (see splits A and B on Fig. 2). The ST33 lineage, therefore, may have continued to exchange genetic material with these taxa after its derivation from the bulk of currently sampled VGII taxa.

From within the genomic diversity of VGII, we can see the emergence of the three low-diversity and highly fit clones (i.e., VGIIa, VGIIb, and VGIIc) that have become prevalent in North America. The emergence of these three clonal subtypes clearly represents at least three distinct events. Our phylogenetic analyses suggest that the phenotypically differentiated VGIIa and VGIIc lineages evolved in South America and spread to the PNW. Additionally, the widely disseminated VGIIb lineage and the VGIIa subtype were likely exported, possibly out of South America, on multiple occasions, and the VGIIb subtype has emerged in multiple geographic locations.

The likely most recent translocation of the VGIIa lineage, as inferred from the very short branch lengths within its respective clade (Fig. 3), resulted in its distinct emergence on Vancouver Island, spreading to the Canadian mainland and into Northwest United States. There is a distinct separation of the PNW dominant VGIIa clade and a clade of South American origin (i.e., Brazil, Argentina, and Japan [in a Brazilian emigrant]) strains. The two VGIIa clades likely have a recent common ancestor, and the lack of its presence elsewhere suggests that this subtype originated in South America. Additionally, WGST allowed us to identify genetic separation between the strains identified during the recent emergence on Vancouver Island and the subclade composed of the older, but less virulent, North American VGIIa strains: NIH-444 (1972; Washington) and CBS7750 (1990; California) (16). This lower-virulence subclade also contained a 1981 Brazilian isolate (ICB-107). It is possible that the latter isolate represents the originating lineage in Brazil and that the two North American isolates represent two distinct translocation events of that same lineage to North America prior to the eventual PNW emergence. Strain ICB-107 was isolated from a patient in Recife, a major port city in Brazil, representing a possible point of export. A lack of epidemiological information associated with that particular isolate prevents further understanding at this time.

The VGIIc subtype has been found only in the PNW. Its most recent splits (i.e., closest genetic relatives) in the phylogenetic network were found exclusively in South America (see split A in Fig. 2). The very short branch lengths and lack of distinct separate clades within this subtype (Fig. 4) are indicative of an individual emergence of this subtype, likely from South America, where all of the most closely related subtypes were identified.

VGIIb is the most diverse of the VGII lineages found in the PNW and likely has been translocated to North America and other parts of the world on multiple occasions. The origin of PNW VGIIb strains is not clear. A VGIIb subclade defined by the putative loss of a 12-gene 10-kb contiguous genomic region (see Fig. S4 in the supplemental material; cluster 1 in Fig. 7) contains most, but not all, PNW VGIIb strains along with strains collected in the previous decade in Australia and Asia (6, 25). These strains are the most commonly found VGIIb strains in the PNW and they reportedly are less virulent than the other VGIIb strains found in the PNW (16). This 10-kb contig is absent from the derived clade yet occurs in all other, more-basal VGIIb strains, including two additional isolates found in PNW, which strongly suggests a single gene loss event in a parental strain. This event, again, provides evidence for multiple independent introductions of VGIIb to North America, which are discernible by both SNP states and genetic content. Further studies are necessary to assess the phenotypic impact of the loss of this genomic component.

The maximum possible consistency indices (CI = 1.0) for the VGIIa, VGIIb, and VGIIc phylogenetic trees demonstrate a lack of observable recombination occurring within or among the North America VGII isolates. These data, taken with the fact that all PNW strains contain only the mating type α locus, indicate that there is limited to no sexual reproduction occurring in these strains. Same-sex mating was previously hypothesized as an evolutionary driver of these subtypes (13). As no phylogenetic inconsistencies were identified among strains within each PNW subtype, there is no obvious evidence of same-sex mating in the region as well, although the progeny of such events may not be readily identified due to the genomic similarities of the parental strains. As sexual reproduction is typically necessary for the development of infectious propagules, it is possible that individual isolates undergo unisexual reproduction (selfing) or haploid fruiting, which have been previously described in C. neoformans (27, 28), to induce a yeast-to-hyphal transition.

In order to identify genetic content associated with niche adaptation to the PNW, we explored the population genome content for genes associated with each of the North American VGII subtypes compared to each other and to other C. gattii molecular types. We were able to identify a subtype-specific gene present (compared to other PNW subtypes) in each of the three subtypes, as well as shared genes in VGIIa and VGIIc compared to the VGIIb subtype. These genomic differences could be the basis for specific niche adaptations and possible differences in virulence and tissue tropism (18). None of the PNW VGIIa, VGIIb, and/or VGIIc clone-defining genes or gene regions was completely unique to North American populations; all were identified in at least once in other global VGII genotypes analyzed.

The gene regions specific to VGIIa in the PNW were rarely seen in other global subtypes and may represent greater likelihood of novel genes. The majority of these VGIIa-specific regions were in a contiguous gene cluster (gene cluster 3, which contains is-cDNA 21 to 29, in Fig. 7) that was seen only in three other global subtypes, originating from Brazil, Colombia, and Australia. This gene cluster contains several genes associated with sugar transport and alcohol dehydrogenase and is possibly important for regional plant sugar metabolism (Table 2). This cluster also includes a collagen-binding domain, likely from a collagenase enzyme, which has been identified in other pathogenic fungi (29). Collagen is a component of all major lung cells (30), and collagen binding in the lung has been associated with colonization of Aspergillus in the asthmatic lung (31) and increased virulence in other pulmonary infections (32). The presence of a collagen-binding domain, therefore, may be relevant to binding of cryptococci to host lung cells and contribute to increased virulence or lung tissue tropism or both in the PNW VGIIa infections (18).

Interestingly, the VGIIa cluster 3 gene cluster is also absent from the other major C. gattii types (VGI, VGIII, and VGIV), except for variants of two of the sugar transport/metabolic domains that were present in the VGI genome, possibly an indication of a historical gene transfer event. Additionally, this gene cluster was absent from the C. neoformans var. grubii and C. neoformans var. neoformans reference genomes (Fig. 7).

An additional gene complex found in the VGIIa group, a DEAD box RNA helicase (is-cDNA 33), was rarely seen in VGII subtypes outside the PNW. Semihomologous genes or gene fragments were identified in VGIIc and four additional VGII subtypes and in the VGIII genome, as well as in the C. neoformans var. grubii genome. The DEAD box RNA helicase has been found previously in multiple eukaryotic and prokaryotic genomes and was recently established as an important virulence factor in C. neoformans, with reported impacts on environmental sensing, pathogenesis, and RNA metabolism (33). Additionally, RNA helicases have been associated with fungal cold response (34). As such, its presence/partial presence in VGIIa and VGIIc and relatively few other VGII lineages may suggest a linkage to differential virulence and/or temperate climate phenotype. Its absence from most other C. gattii genomes suggests one or more intertype or interspecies gene transfer events.

An analysis of the genome content of the three North American VGII subtypes against all other VGII genomes in the current study revealed an ~10-kb region, containing three genes, that was consistently found in all VGIIa, VGIIb, and VGIIc genomes. This cluster (cluster 6, which contains is-cDNA 44 to 46, in Fig. 7) was present in only a limited number of other VGII subtypes and was absent from the representative VGI, VGIII, and VGIV strains. BSR analysis revealed that the region had high similarity with a homologous region of the same size in the C. neoformans var. grubii reference genome but was missing from C. neoformans var. neoformans. The three genes in this cluster were identified as encoding rhamnosidase, avenacinase, and an unspecified sugar transporter. Rhamnosidase has been previously identified in Cryptococcus albidus (35), Aspergillus spp. (36), and other basidiomycetes and ascomycetes (37). The presence of rhamnosidase may enable more efficient use of rhamnose as a carbohydrate source. Of particular note, rhamnose is a free sugar produced by Rhamnus (buckthorn), a plant with a large distribution in the PNW. Rhamnose is also found in members of the Uncaria family, a genus found throughout the tropical region. It is a component of many plant and bacterial polysaccharides, as well as plant glycosides. Additionally, rhamnosidase in soft-rot fungi has been shown to have lignocellulolytic activity, which is involved in the degradation of hardwoods (38). Its presence in the PNW population may also provide an adaptive advantage to old-growth forests of the region, as C. gattii been recovered from rotting and decaying wood in other geographic areas (39, 40).

The genes found in VGIIc, but not VGIIa or VGIIb (is-cDNA 34 to 41), relate primarily to sugar transport, metabolism, and gene expression and are commonly seen in other VGII subtypes and non-VGII cryptococcal genomes, suggesting that the lack of homologous genes in VGIIa and VGIIb represents independent gene loss events.

Among the three gene regions found only in the North American VGIIb strains is a haloacid dehalogenase-like hydrolase (is-cDNA 18). This gene was rather infrequent in other VGII subtypes, although a highly homologous region was seen in VGI, VGIII, and VGIV. A haloacid dehalogenase was found to be a virulence factor in Paracoccidioides brasiliensis, necessary for lung epithelial cell adhesion (41). No function for this gene region has been established in Cryptococcus.

The majority of all other between-PNW subtype genetic differences include genes or domains that are seen relatively frequently in other VGII subtypes. The presence or absence of these genes will need to be further explored to understand any functional role associated with the PNW populations.

The 12-gene cluster that was absent in the derived clade of VGIIb (cluster 1, which contains is-cDNA 1 to 12, in Fig. 7; see Fig. S4 in the supplemental material) was present in all other VGIIb isolates, including the separate clade found in North America, and was found, wholly or partially, in all other Cryptococcus isolates. This gene complement is composed of structural, regulatory, transport, and repair genes. The significance of this deletion remains unknown, but it is a marker of the derived lower-virulence clade of VGIIb. Phenotypic studies are necessary to determine whether there are other differences between the derived clade and the rest of the VGIIb subtype.

Beyond the clear differences of gene content between and among VGII subtypes, multiple lines of evidence demonstrate the likelihood of “genome sharing” throughout the VGII population. In addition to the “all VGII” tree consistency index of ~0.4 (meaning a high degree of homoplasy) (see Fig. S2 in the supplemental material), and a significant phi test for recombination (P = 0.0), it is evident that gene content does not dictate phylogeny. We assume that recombination throughout the VGII population accounts for the high degree of homoplasy, which may be accomplished through larger genome content sharing events (i.e., bisexual or same-sex mating [19]). While the neighbor net phylogenetic network recovers a mostly star-like shape, there are apparent relationships between lineages (Fig. 2). The fineStructure analysis conducted here assigned lineages to specific “populations” and assessed levels of relatedness between genomes using adjacent linked SNPs to identify the nearest neighbor relationships for portions of the genome bounded by recombination events. Evidence of recombination was identified in South American isolates. For example, isolates WM 05.275 (originating in Colombia) and WM 05.529 (originating in Brazil) share many of their linked loci, although they are distinct subtypes and are separated by WGST, and they have distinct gene content by BSR. The same pattern is seen with isolates WM 05.419 and WM 05.547, both Brazilian clinical isolates. The VGIIa (R265 [British Columbia, Canada]) and VGIIb (B9563 [Washington]) representatives showed clear evidence of limited genome sharing with most isolates in the analysis. Conversely, the VGIIc genome showed almost no sharing with the vast majority of other VGII subtypes, except its nearest phylogenetic relatives from Brazil (samples WM 05.462 and IAL-3234). Taken together, these results indicate a relatively weak population structure within VGII or a recent subdivision of populations that have not yet allowed gene absence and presence to become fixed within recombining groups. Again, the high consistency indices for whole-genome SNP phylogenies of each PNW subtype (Fig. 3 to 5) indicate a lack of observable recombination occurring within these clones.

The presence of the three-gene cluster with a high degree of homology to C. neoformans var. grubii (cluster 6) provides evidence of a possible previous hybridization event between C. gattii and C. neoformans var. grubii, resulting in the possible introgression of this cluster in either a limited number of C. gattii subtypes or in the C. neoformans var. grubii genome. More whole-genome sequencing of the latter will be needed to ascertain the direction of introgression. Natural C. gattii and C. neoformans var. grubii hybrids have been recently identified from clinical samples in Colombia and Brazil (42). Additionally, genetic remnants of previous hybridization events, also referred to as “identity islands,” have been reported in the case of a nonreciprocal exchange of a multigene cluster from C. neoformans var. grubii to C. neoformans var. neoformans (43). The existence of hybrids and subsequent introgression are an additional mechanism of genome evolution (beyond recombination, SNP mutation, and gene gain, loss, and duplication) and may allow continued successful emergence of novel phenotypes.

No clear progenitor strains have been identified for VGIIa, VGIIb, or VGIIc, although the first two have been found in South America. The phylogenetic evidence shown here is consistent with multiple translocations of VGIIb to the PNW and elsewhere, as well as possible multiple introductions of VGIIa into North America. The VGIIc subtype lacks representation outside the PNW; however, its consistent SNP and gene content relationship with subtypes from Brazil provides compelling evidence of a likely South American origin for VGIIc as well. The results of the analysis also confirm that South America is a major source of VGII diversity and, with numerous subtypes and extensive recombination, it appears to represent an evolutionary accelerator for the VGII subtype and clones emerging from it.

Our extensive whole-population genome analysis enabled an in silico search for genetic factors associated with differing phenotypes in the PNW populations; similar in silico genome-wide association studies (GWAS) have been useful for studying adaptations in other widely dispersed fungi (34). The initial SNP and gene analyses based on this population-level study provide an important perspective on C. gattii diversity and will serve as a baseline for future studies. The identification of variant sugar transport and metabolism genes, as well as previously unexplored virulence genes (e.g., collagen adhesion gene in VGIIa), provides excellent candidates to pursue in functional studies that will help explain the PNW outbreaks and other C. gattii emergence events. The comprehensive SNP analyses will also be useful for future mutation-based functional analyses. Given the obvious gene content differences seen within VGIIb, even demonstrably clonal subtypes are not uniform and may require intrasubtype comparative analyses.

We recognize that limitations of this study include sampling bias and hence possible overrepresentation of isolates from some geographic regions and underrepresentation from others. There could be further diversity of VGII in other regions that were not sampled or that were sparsely sampled. Additionally, given the finding of C. gattii VGII strains in undisturbed Amazon rainforest and the amount of diversity already seen in South America (7), there is doubtless much more unsampled diversity on that continent. Additionally, this study was designed to identify important genomic features from whole-genome sequence data and was not designed to identify important mutations impacting gene regulation or expression. Nevertheless, our findings provide a significant advance in our understanding of the emergence of C. gattii in the PNW though the importance of the gene differences is unknown and requires further functional genomic analyses. As has been seen in recent years, C. gattii continues to spread to new locales, and a comprehensive understanding of its population genomics may help researchers, public health officials, and clinicians to better respond to the impacts of its continued evolution and emergence.

One hundred eighteen C. gattii isolates were included in the analyses (Table 1). Information for clinical strains was collected with the approval and oversight of the institutional ethic committee. Additionally, the two near neighboring taxa, C. neoformans var. neoformans strain JEC21 and C. neoformans var. grubii strain H99, were obtained as whole-genome assemblies from NCBI (http://www.ncbi.nlm.nih.gov/Genbank/index.html) and the Broad Institute of Harvard and MIT (http://www.broadinstitute.org/annotation/genome/cryptococcus_neoformans_b), respectively. The 118 C. gattii isolates were selected based on previously established MLST sequence types (Table 1) to ensure global genomic diversity. The previously published VGIIa reference genome (R265) (see the Broad Institute link above) was resequenced for this analysis, due to the relatively high number of likely errors in the existing Sanger-based reference genome (17). The genomes of all 118 C. gattii isolates were sequenced using Illumina HiSeq, MiSeq, or GAIIx sequencing technology as previously described (17, 44). High-molecular-weight DNA was extracted using the ZR Fungal/Bacterial DNA MiniPrep kit (catalog no. D6005; Zymo Research, Irvine, CA, USA). DNA samples were prepared for paired-end sequencing using the Kapa Biosystems library preparation kit (catalog no. KK8201; Woburn, MA, USA) protocol with an 8-bp index modification. Briefly, 2 µg of double-stranded DNA (dsDNA) sheared to an average size of 600 bp was used for the Kapa Illumina paired-end library preparation as described by the manufacturer. Modified oligonucleotides (Integrated DNA Technologies, Coralville, IA, USA) with 8-bp indexing capability (45), were substituted at the appropriate step. Prior to sequencing, the libraries were quantified with quantitative PCR (qPCR) on the ABI 7900HT (Life Technologies Corporation, Carlsbad, CA, USA) using the Kapa library quantification kit (catalog no. KK4835; Woburn, MA, USA). Libraries were sequenced to a read length of 100 bp or 150 bp on the Illumina HiSeq system, to 75 bp or 100 bp on the GAIIx system, or to 250 bp on the MiSeq system. All WGS data files have been deposited in the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/bioproject; project PRJNA244927).

All sequenced genomes were de novo assembled using the SPAdes assembler v2.5.0 (46). Raw contigs were filtered by length and coverage, as SPAdes tends to produce many short and low-coverage contigs. For isolates assembled from 100-bp or 150-bp paired-end data, contigs shorter than 200 bp were filtered and removed, while contigs shorter than 500 bp were filtered/removed in genomes assembled from 250-bp paired-end MiSeq data. Additionally, for every genome, contigs with coverage below 20% of the average coverage of the 20 largest contigs were filtered out.

The seven MLST loci, CAP59, GPD1, LAC1, SOD1, URA5, PLB1, and intergene sequence (IGS), were aligned to the contigs according to the consensus MLST scheme for the C. neoformans/C. gattii species complex (10). Allele types (ATs) and sequence types (STs) were assigned according to the MLST database at http://mlst.mycologylab.org, and new ATs and STs were added to this database. A maximum likelihood tree was constructed based on ST locus scores using the program MEGA v. 5.05 (47).

Illumina read data were aligned against the herein newly obtained R265 (high virulent VGIIa Vancouver Island outbreak strain) assembly using Novoalign 3.00.03 (Novocraft Technologies, Selangor, Malaysia), and then SNP variants were identified using the GATK Unified Genotyper v2.4 (48). SNP calls were then filtered (http://tgennorth.github.io/NASP/) to remove positions with less than 10× coverage, with less than 80% variant allele calls, or identified by Nucmer as being within duplicated regions in the reference. In the cases of the VGIIb and VGIIc analyses, the read data were aligned by the same methods using the genome assemblies for B9563 and B8825, respectively. In total, five SNP matrices with different taxa were produced: (i) C. gattii molecular types VGI to VGIV; (ii) C. gattii molecular type VGII only with all previously established MLST sequence types; (iii) C. gattii VGIIa; (iv) C. gattii VGIIb; and (v) C. gattii VGIIc.

To understand relationships between isolates, we conducted whole-genome SNP typing (WGST) as previously described (17, 39, 49). Maximum-parsimony SNP trees, based on each of the above five SNP matrices, were constructed, using PAUP* v.4.0b10 (50) and visualized in FigTree v.1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/). The C. gattii VGII tree was rooted using the large C. gattii tree, while the VGIIa, VGIIb, and VGIIc trees were each rooted using the other VGII subtypes as outgroups. Roots are shown simply as the midpoint on the branch determined by outgroup rooting, because the exact number of SNPs assigned to each branch varies slightly with the inclusion of different taxa, due to the SNP filtering strategies employed. The phi test for recombination was conducted on the whole VGII SNP data set, as well as the same data set broken into 500-bp nonoverlapping windows, using the PhiPack software (51). Due to recombination signal strongly present in the VGII SNP data, a neighbor-joining split tree network or “neighbor net” (52) was drawn to visualize the relationship between isolates in the VGII tree. The neighbor net was drawn using SplitsTree4 (53) with the uncorrected P-distance transformation. The GTR+G model of nucleotide evolution, selected by jModelTest2 (54), was also used to calculate distances for input to the neighbor net algorithm, which produced a network showing the same divisions among samples (data not shown).

In addition to whole-genome SNP phylogenies, fineStructure analysis (25) was performed using SNP matrix 2 (see “SNP variant detection” above) to help identify the population structure within VGII. Using one representative from each of the 32 distinct VGII subtypes (i.e., MLST sequence types) (Table 1), visualized by placement in the SNP phylogeny, whole-genome SNP data were reduced to a pairwise similarity matrix using ChromoPainter, which estimates the number of regions in two genomes that coalesce prior to coalescing with another genome. ChromoPainter uses the linkage disequilibrium present between closely positioned SNP loci to help assign a given region of a genome to its best near neighbor, increasing the sensitivity of population detection over unlinked methods. ChromoPainter was run using the linkage model, assuming uniform rates of recombination per base pair of sequence. Next, the population structure is identified by fineStructure using this similarity matrix, where the underlying model assumes individuals within populations will share more regions of their genome with each other, and have a similar amount of admixture with individuals from different populations. After identifying populations, fineStructure then merges populations one at a time, at each step using the most likely population merge, to relate populations to each other through a tree.

Potential gene content comparison of taxa was performed using BLAT score ratio (BSR) analysis, as previously described for prokaryotic genomes (55), with modifications made to accommodate eukaryotic genes containing introns (described below). Gene sequences were predicted for each assembled genome using the ab initio eukaryotic gene prediction software Augustus v2.6.1 (56). For each gene prediction, the coding DNA region sequences were concatenated into a single contig, removing introns to create the putative in silico coding DNA (is-cDNA). is-cDNA sequences from all isolates were then clustered and dereplicated using USEARCH v5.0.144 (57) with a 90% identity cutoff, creating a nonredundant set of consensus is-cDNA predictions encompassing the entire diversity of samples. Putative is-cDNAs were then aligned back to each genome assembly using BLAT v35x1 (58), which can align the cDNA sequence against a reference allowing for intervening intron sequences, using an 80% identity cutoff. BLAT alignments were converted to scores replicating the calculation used by the web-based BLAT server (https://genome.ucsc.edu/). For each putative is-cDNA, the raw score for an isolate was converted into a ratio by dividing the score of a particular isolate by the top score achieved by any isolate, as previously described (55). These BSRs were then imported into the MultiExperiment Viewer (MEV) to generate a heat map which could be manually inspected to identify genes present and absent across different sets of taxa (Fig. 7). Thus, a value of 1.0 (yellow on the heat map) represents the best hit, and a value of 0 (blue) represents no alignment by BLAT.

Using the BSR method described above, the potential gene content between VGIIa, VGIIb, and VGIIc sequence types was compared using representative isolates from each subtype present in the PNW (Table 2). Representative strains were selected based on assembly quality and being part of the PNW outbreak subtypes (i.e., isolates typed as VGIIa and VGIIb but not from the PNW were not included in these sets; all VGIIc isolates were from the PNW). This set included 9 VGIIa isolates, 6 VGIIb isolates, and 12 VGIIc isolates (Table 1). Putative consensus is-cDNAs showing differential content across the three PNW sequence types were identified using a cutoff BSR score of above 0.9 in the taxa of interest and below 0.9 in the excluded taxa. Multiple assemblies from each subtype were included to help reduce the number of genes that may appear to be absent due to variability of draft genome quality with multiple contigs per chromosome. Genes with a low BSR value in a single member of a subtype almost invariably were due to contig breaks in assembly. Unique putative genes were identified for each of the VGIIa, VGIIb, and VGIIc groups as well as for genes contained in both VGIIa and VGIIc but missing in VGIIb, as this comparison may help explain the previously described phenotype differences seen in the PNW.

Expanding the analysis beyond differences within the PNW sequence types, the putative is-cDNA sequences that define the PNW subtypes were scored in each of the additional 29 global VGII subtypes (Table 1 and Fig. 7), a representative sequence from each of the VGI, VGIII, and VGIV major molecular types, as well as C. neoformans var. grubii (H99) and C. neoformans var. neoformans (JEC21). This allowed for the gene differences found in the PNW subtypes to be examined across the diversity of the Cryptococcus neoformans/C. gattii species complex, as well as for the identification of genes present in all of the PNW isolates but not universally present in other VGII subtypes. In addition, the full complement of is-cDNAs was scored in all 29 global VGII subtypes.

To characterize genes absent in the PNW populations but present in many of the other global VGII subtypes, a second set of putative genes (is-cDNAs) was generated using the gene predictions in non-PNW VGII subtypes. While this would presumably contain many of the same genes found in the PNW gene set, it would also contain any genes that may be absent in all of the PNW subtypes. Using the same cutoff scheme as described above, putative genes were identified that were missing or truncated in all of the PNW subtypes but present in at least half of the other VGII subtypes. is-cDNA sequences identified as having a differential content pattern of interest were then located in the genome assemblies of strains R265 (VGIIa), B9563 (VGIIb), B8825 (VGIIc), depending on which taxa contained the is-cDNAs, to look for blocks of adjacent genes showing similar patterns of presence and absence.

Because a low BSR score may indicate gene loss or truncation as well as problems caused by assembly issues in draft genomes, putative gene sequences identified as differential across taxa of interest were subsequently assessed by raw read data alignment. This provides evidence for confirmation of the differential gene content and, in the case of scores between 0.0 and 0.9, helps to identify the type of gene structure differences present between isolates with different scores. To avoid potential problems in aligning genomic read data against putative is-cDNA sequences, which do not contain introns, full-length gene predictions containing introns were dereplicated with USEARCH to serve as read alignment references. Putative is-cDNAs of interest were aligned against the full-length gene sequences using BLAT to identify the corresponding full-length gene. Read data from each isolate in the PNW analysis was then aligned to the full-length gene. These alignments were inspected for each gene to ensure read coverage across the entire gene sequence for isolates with high BSR scores and to ensure a lack of read coverage for isolates with BSR scores indicating gene absence. Moderate BSR scores occurred in a few cases when limited read data successfully aligned to a putatively missing gene (e.g., is-cDNA 9 in VGIIb), or when read data aligned to truncated genes. In an effort to be conservative in identifying differential gene content using short-read data, gene differences that could not be confirmed with read data are not reported. This most often occurred in the case of putative genes that contained highly repetitive sequences, which cause assembly and definitive alignment issues in short-read data.

is-cDNAs, which were confirmed as variable with read data (see above), were translated into amino acid sequences and then searched against the Pfam database (http://pfam.sanger.ac.uk/) to help identify potential protein functions. Protein families and domains achieving a significant score are reported in Table 2. Additional characterization for selected is-cDNAs was performed using a blastp search of the NCBI nonredundant protein database (http://www.ncbi.nlm.nih.gov/RefSeq/).