Dating the Cryptococcus gattii Dispersal to the North American Pacific Northwest

The Teddy Roosevelt effect.

Given the commonality of just under 100 years for the PNW arrival of the VGIIa, -b, and -c subtypes and given their lack of known presence elsewhere on the North American continent (with three notable exceptions: a single Florida VGIIb isolate [28], the recently discovered Québec cluster of VGIIa [16, 29], which seems to have derived from the PNW VGIIa population, and the single VGIIb isolate from a deer in Nova Scotia [17], likely derived from the PNW VGIIb population), it is logical to hypothesize that these C. gattii subtypes arrived in the PNW region shortly into the 20th century. It is possible that populations arrived elsewhere on the continent and failed to become established; a similar scenario was seen with the concurrent spread of Yersinia pestis by boats from Asia to multiple U.S. ports around the turn of the 20th century that only became enzootic through the ports in California (30). What appears to be clear is that there was an outward spread of the PNW C. gattii populations from the PNW coastal regions to locales further inland (including across the country of Canada). These findings would suggest at a minimum that a translocation event or events occurred between Brazil and the PNW coast approximately 100 years ago. As shipping was minimal between Brazil and Pacific Coastal countries in South America until the opening of the Panama Canal, it follows that this historical event may have allowed for ships to carry cargo, and microbial stowaways, between these distant locations for the first time (Fig. 3).

The Panamanian isthmus acts as a land bridge between the American continents and as a natural barrier between the Pacific and Atlantic Ocean. The isthmus has played a crucial role in the geology, hydrology, biology, and climate in the region since the uprising of this intercontinental land bridge nearly 3 million years ago (31). The isthmus has also played a key role in the dispersal of microbes: notably, the fungal pathogen Coccidioides posadasii likely was translocated from North America to South America across the isthmus via mass movements of mammals during the later Great American Biotic Interchange events that occurred nearly a million years ago (32). After the final rising of the isthmus, a permanent land barrier was formed between the continents, separating the oceans, until just over 100 years ago. In 1914, after years of human toil, the Panama Canal, commissioned by Theodore Roosevelt, was completed. This allowed, for the first time in nearly 3 million years, water movement between the Atlantic and Pacific Oceans, subsequently allowing for cargo ships to move between the east and west coasts of the Western Hemisphere continents without having to go through dangerous waters at the southern point of South America. After the opening of the Panama Canal, shipping between ports of Brazil and western North America began in earnest.

Common goods transported out of Brazil at this time were Brazilian hardwood lumber, minerals, coffee, and rubber. While C. gattii is known to have been transported on (33) and to thrive on (34, 35) live trees and seeds, no contamination of lumber products has been documented previously. While we may be unable to ascertain a specific shipping product that may have carried the fungus from Brazil to the PNW and elsewhere, one hypothesis we propose is contaminated ballast water. Ballast water contamination is a well-understood problem in large vessels, and ballast has been known to contain and transport nonnative animals, algae, and microbes around the world (36). At a vessel’s source port, seawater is taken on by ships and stored in their ballast to provide vessel stability. This ballast water is released en route or at the destination port to correct for fluctuations in cargo weight (37). Cryptococcus has been shown to survive in vitro in saline (8), has been isolated from seawater (38, 39), and is known to cause significant morbidity and mortality in marine mammals (40–42), particularly in the PNW (43–45). Again, while only hypothetical, it is possible that a limited number of ships with contaminated ballast water could have transported these pathogens unknowingly to PNW ports and beyond. It is notable that the dominant PNW clade of VGIIb also shares a very recent common ancestor with multiple lineages in Southeast Asia. Again, a shipping vessel from Brazil carrying goods (and contaminants) to the PNW could also have caused the translocation to Asia at this time, given that the new transcontinental seaway would have allowed shipping to these markets as well.

An alternate and intriguing hypothesis for dispersal to Malaysia and elsewhere in Southeast Asia is the possible contamination of rubber tree (Hevea brasiliensis) seeds, which were infamously “stolen” and exported to Britain by Henry Wickham in 1876 to establish British-controlled rubber tree plantations in Ceylon (now Sri Lanka) and Malay (modern day Malaysia and Singapore) and, subsequently, to British colonies in Africa and India (46). While not previously shown to be associated with Hevea seeds, C. gattii contamination of other tree seeds, such as Eucalyptus spp., has been well documented (33). Brazil, attempting to maintain global control of the rubber trade, tried to prevent the shipment of such seeds, limiting other possible such exports out of Brazil; however, the Wickham seed transport likely resulted in the crash of Brazilian rubber beginning in 1910 (46). The 1876 seed transport timeline and global movement of Hevea seeds at this time could represent an indirect dispersal mechanism. Along these lines, the Brazil nut (Bertholletia excelsa) has been an additional plant-based export out of Brazil around the world for multiple centuries; however, intense exportations began following the Brazilian rubber crash in 1910 (47), providing an appropriate temporal connection to the C. gattii dispersal event. Extensive analyses of well-known ongoing fungal contamination of Brazil nuts during collection and processing have not identified the presence of Cryptococcus spp. (47), and there is no obvious route to localized environmental contamination in the PNW from these products.

While the prevalence of recombination along with the extent of temporal signal within a sample set plays a large role in estimating genome-wide evolutionary rates (18), these sample sets show no statistical evidence of recombination to hinder these analyses. However, the date randomization permutation testing showed there was little true temporal signal within the VGIIb and VGIIc sample sets. Given the narrow sampling dates and the limited number of samples in each data set, it is possible we are only capturing a glimpse of the temporal signal and, with a wider sampling time and more samples, it would be possible to narrow our MRCA estimates.

Whatever the cause of C. gattii dispersal to the PNW, it is clear that those populations are neither ancient nor very recent (i.e., <25 years ago) arrivals to the region. Dispersal in the last 100 years would strongly suggest anthropogenic causes, directly or indirectly, rather than movement by natural animal migrations, as suggested with Cryptococcus populations (48) and other environmental pathogenic fungi (32). The identification of VGIIa in a pet store bird (sample number B11566) and a separate, unrelated case in pet store worker (sample number B11569), both in Québec, are unexpected occurrences (29), and it is of interest to note that the bird reportedly originated from British Columbia and the infecting C. gattii strain clearly evolved from the PNW VGIIa population, suggesting a linkage to human transport of animals as the cause. Conversely, the appearance of a PNW-derived strain of VGIIb (sample number B11567) in a Nova Scotia yearling white-tailed deer (Odocoileus virginianus) (17) is less explainable. There is no other current evidence of long-distance dispersal of VGIIb across the North American continent, nor is there evidence of transport of animals from the PNW to the region, as the infected deer was part of a local deer population. White-tailed deer are nonmigratory, providing further evidence of the presence of VGIIb in Nova Scotia, and this is suggestive of a unique introduction/dispersal event in the region (17). The province of Nova Scotia is a peninsula that is surrounded by four major bodies of water, and the city of Halifax, throughout its history, has been an active sea port, one of the major international points of entry for shipping into Atlantic, Canada. More investigation of the local environment is necessary to understand the history of this incident and provide estimates of regional endemicity. As the North American populations of C. gattii continue to evolve and disperse, it will be useful to continually apply genomic dating to understand the nature of these events and the expanding impact of these fungi on human and veterinary health.

Samples.

A total of 66 genomes were included in the analysis of the three clonal lineages found in the PNW. The 22 new genomes from this study include 13 C. gattii recent isolates collected either in 2015 or 2016 from Washington, Oregon, and Canada, which were comprised of 7 VGIIa, 1 VGIIb, and 5 VGIIc isolates (Table 1). Previously published isolate genomes (12) were included for each genotype, including 22 VGIIa, 16 VGIIb, and 24 VGIIc genomes (Table 1). The VGII lineage tree was based off of 134 genomes, 112 of which were previously published (12), including the above-mentioned samples (Table 1). The recently described msh2 variant lineage of VGIIa (49) was not included in the genotype-specific genomic dating analyses (described below) due to possible confounding from the hypermutation phenotype, except for the NIH444 genome, which had no identifiable impacts on the total number of SNP mutations (12).

Genomic DNA was extracted from the 22 isolates using the ZR fungal/bacterial DNA MiniPrep kit (Zymo Research), following the manufacturer’s instructions. DNA samples were fragmented by sonication and prepared for multiplexed, paired-end sequencing with a 700-bp insert using the library hyper-preparation kit with standard PCR library amplification (KAPA Biosystems) as previously described (32). Libraries were quantified using the KAPA library quantification kit (KAPA Biosystems) and sequenced to a read length of 300 bp on the MiSeq instrument (Illumina, Inc., San Diego, CA).

SNP matrix generation.

Sequenced genomes were assembled de novo using SPAdes version 3.6.0 (50). Reference genomes for phylogenetic analyses were selected based on the assembly quality metric N₅₀ and total assembly length (see Table S1 in the supplemental material). NASP (51) was used to generate an SNP matrix for each of the three C. gattii subtypes, as well as the VGII complex, using the selected reference genomes for each genotype and R265 for the VGII lineage analysis. In brief, sequencing reads were aligned to the reference genome using Novoalign (Novocraft Technologies Sdn Bhd) and SNPs were identified using the GATK Unified Genotyper ToolKit version 2.7-2 (52). NASP filtered SNP loci not present in every sample and with less than 10× coverage and less than a 90% consensus in any sample. Using NUCmer (53), duplicated regions identified within the reference genome were also removed.

Phylogenetic analyses and root-to-tip regression.

IQ-Tree version 1.3.10 (54) was used to infer the best-fitting nucleotide substitution model for each of the three SNP matrices and to produce maximum-likelihood trees with 1,000 nonparametric bootstrap pseudoreplicates for branch confidence. Trees were visualized in FigTree version 1.4.2, and SNP numbers placed on branches using phangorn (55). To test each SNP matrix for evidence of recombination, which can confound divergence-dating analyses (18), the PHI test was conducted using PHIPack (56). In order to assess the temporal signal of the data set, regression analysis implementing root-to-tip genetic distance as a function of the sample collection year was conducted using the software package TempEst version 1.5 (http://tree.bio.ed.ac.uk/software/tempest/). Using the determination coefficient, R², a measure of clocklike behavior was assessed. In an effort to maximize R², the best-fitting root was selected based on TempEst recommendations (18). Additionally, we performed 10,000 random permutations of the sampling dates over the sequences in an effort to evaluate the significance of our regression results (57).

Divergence time analyses.

A Bayesian molecular clock using tip dating was implemented in the BEAST version 1.8.0 software package (58) to infer evolutionary rates and time to most recent common ancestor (TMRCA) for the three VGII genotypes. MEGA7 (14) was used for nucleotide substitution model selection for each genotype, utilizing the Bayesian information criterion results to determine the best-fitting models. As previously described (32), available nucleotide substitution models are limited in BEAST, and therefore, the best-fitting model from MEGA7 that was also available in the BEAST software was implemented. While the SNP matrices only included variable sites, we corrected for the invariant sites by specifying a Constant Patterns model in the BEAST XML file. For each separate analysis, we determined the numbers of constant A’s, C’s, T’s, and G’s and added them to the XML file. Additionally, “path and stepping stone” sampling marginal-likelihood estimators were used in order to determine the best-fitting clock and demographic model combinations (59). The statistical fits of 10 different clock and demographic model combinations were assessed using the log marginal likelihood (Table S2). In BEAST, four independent chains of 4 billion iterations each were run for all molecular clock and demographic model combinations for the VGIIa and VGIIc data sets. Convergence among the four chains for VGIIb analysis completed after 1 billion iterations.

Accession number(s).

All new sequence data files were deposited in the NCBI Sequence Read Archive (BioProject accession number PRJNA388113).