Local admixture of amplified and diversified secreted pathogenesis determinants shapes mosaic Toxoplasma gondii genomes

1Most of the diversity of eukaryotic life is contained in early branching, unicellular organisms that differ substantially from model organisms such as yeast, flies, worms and mice. This diversity is illustrated by the protozoan phylum Apicomplexa, estimated to contain more than 5,000 species, most being parasitic on insects and mollusks, while a few cause disease in domestic animals and/or humans. Studies of these few disease-causing agents comprise our limited knowledge of this phylum, which demarcate a deep branching phylogeny that has been estimated to span more than ∼400 my of evolution (Fig. 1a). Over that time frame, it is likely that apicomplexans have adapted to their various vertebrate hosts via multiple independent changes in host range, and yet the molecular mechanisms underlying these adaptations remain largely undefined.

2Although most members of this phylum are adapted to a narrow range of hosts, Toxoplasma gondii stands out as a generalist. The genus is characterized by a single species that enjoys worldwide prevalence in animals including humans. Infections with T. gondii are common, yet they typically only cause disease in immunocompromised hosts, or as a result of transplacental infection. T. gondii is equipped with excellent forward and reverse genetic tools, providing a model for many less-tractable apicomplexan parasites. As a highly successful parasite, T. gondii is positioned to inform us about genomic features that are important for efficient transmission and expansion of host range. Here, we sought to exploit this potential by analysing the composition and diversity of the T. gondii genome in comparison to several closely related apicomplexan parasites.

3T. gondii belongs to the tissue-cyst forming coccidian parasites, which is distinguished from enteric coccidian parasites by having an alternating two-host (heteroxenous) life cycle (Fig. 1a, Table 1). Most tissue-cyst forming coccidian parasites have obligatory heteroxenous life cycles (that is, Sarcocystis spp. and Hammondia spp.), while others share this mode but have evolved additional strategies for transmission (Table 1). Notably, both T. gondii and Neospora caninum can cause congenital infection, while only T. gondii can be transmitted between intermediate hosts by oral ingestion of infected tissues, thus bypassing the sexual phase of the life cycle. These flexible features in the T. gondii life cycle likely aid in transmission through the food chain, thus underlying its broad host range (Table 1). In contrast to our appreciation of differences in life cycle, modes of transmission and host range among these closely related parasites, their molecular bases remain largely unexplored.

4In North America and Europe, the population structure of T. gondii is dominated by three prevalent clonal lineages, which coexist with much more rare, genetically diverse isolates. A fourth clonal lineage is largely confined to North America, where it is more common in wild animals. In contrast, much greater genetic diversity is seen in South America where the population lacks signs of the recent genetic bottleneck and clonal structure seen in the Northern Hemisphere. T. gondii utilizes rodents and birds as natural intermediate hosts, and hence it is particularly well adapted for survival in these niches. Forward genetic mapping studies have identified several families of secretory proteins in T. gondii that are important for thwarting innate immunity and hence facilitating infection in the mouse. Related effectors are conserved in Hammondia hammondi, hence the basis for the dramatic differences in biology of these two parasites remains unclear. Nonetheless, one hypothesis advanced by the study of select laboratory strains in the mouse model is that pathogenicity, and perhaps host range, may depend on the repertoire of such secretory pathogenicity determinants, although this has not been tested on a wider level.

5Here we tested the generality of this hypothesis through genomic analyses of 62 strains of T. gondii in comparison to several closely related parasites. Our findings reveal that expansion and diversification of secretory pathogenesis determinants (SPDs), which are often tandemly clustered, is a prominent feature of the genomes of T. gondii and related tissue-cyst forming coccidians. Furthermore, patterns of block inheritance, due to recent admixture or selective retention, may underlie specific traits that are shared by related lineages of T. gondii containing similar combinations of SPDs. These features define the population structure of T. gondii and have implications for the evolution of transmission, host range and pathogenicity.

6We undertook a comparative genomics approach to understand the population diversity of T. gondii and its relationship to closely related tissue-cyst forming coccidian parasites. First, we generated additional genomic DNA sequence coverage (∼26 × coverage) and RNA-seq data (>1,000 × mean coverage of coding sequence) to improve the assembly and annotation for the reference ME49 strain of T. gondii (Table 2). We also generated a whole-genome sequence for H. hammondi (∼66 × coverage; Table 2) and compared these two closely related parasites to the recently completed genomes of Sarcocystis neurona and N. caninum (Table 2), which cause economically important diseases in horses and cattle, respectively (Table 1). Finally, to provide insight into genetic variation of T. gondii we derived whole-genome sequences for 61 additional isolates that were chosen to span presently known global diversity (Supplementary Data 1). Among the total of 62 T. gondii strains, 16 reference strains representing the major haplogroups were sequenced by both 454 (3 and 8 kb paired-end libraries) and Illumina (300 bp paired-end libraries) technologies and the resulting reads were assembled and annotated separately (∼47 × average sequence coverage) (Supplementary Table 1). The remaining strains were sequenced using Illumina only (∼42 × average sequence coverage) and were aligned to the reference strain ME49 (Supplementary Data 1). Below, we present the comparative analyses of these genomes focusing on three broad themes: (1) comparison of T. gondii to the most closely related tissue-cyst forming coccidian parasites, (2) analysis of the core genome of T. gondii and how it has diversified and (3) examination of how the global population structure of T. gondii has been shaped by local genomic admixture.

7We compared the whole-genome sequences from four related tissue-cyst forming coccidian parasites spanning a range of biological hosts and life-cycle strategies (Tables 1 and 2, Fig. 1). Three of the four organisms have a similar total genome size of 62–65 Mb (N. caninum, T. gondii and H. hammondi), while the S. neurona genome is somewhat larger due to expanded repeats and much larger introns (Table 2, Supplementary Table 1). All four genomes have roughly similar GC compositions and are predicted to encode from 7,000 to slightly more than 8,000 genes located on 14 chromosomes, as verified in T. gondii (Table 2). Similar to other genome sequencing projects, 42–56% of the predicted CDSs (coding DNA sequences) encode genes with a putative functional domain annotation, while 44–58% are hypothetical unknowns. To identify conserved features, we compared the four different genomes to the enteric coccidian Eimeria tenella using OrthoMCL to cluster genes into putative orthogroups. Not surprisingly, more closely related taxa showed a higher degree of shared OrthoMCL clusters (Fig. 1b, Supplementary Data 2). Orthogroups classified by Pfam domains and grouped into the top 20 Gene Ontology (GO) terms (http://geneontology.org/) revealed that all five species share similar orthologous groups for many key biological functions, and that tissue-cyst coccidians are enriched in processes involved in protein modification (Supplementary Fig. 1).

8The most abundant protein domains in tissue-cyst forming coccidian parasites include serine/threonine (S/T) kinases, RNA-binding proteins, PP2C-type S/T phosphatases and calcium-binding motifs (EF-hands) (Fig. 1c). There is a precedent for the importance of S/T kinases, such as the expanded polymorphic family of rhoptry (ROP) kinase virulence determinants in T. gondii and of calcium-binding motifs, including within a family of calcium-dependent protein kinases. However, the abundance of RNA-binding proteins (RMR and RMA motifs) was unexpected, as these have been largely unexplored in T. gondii and closely related parasites. In addition, plant-like AP2 transcription factors are abundant in tissue-cyst forming coccidian parasites (Fig. 1c), consistent with these being major transcription factors in apicomplexans. Also prevalent in N. caninum, T. gondii and H. hammondi are a family of surface antigens (SAG) called the SRS family (Fig. 1c), which are amplified and highly divergent among tissue-cyst forming coccidian parasites. These appear less abundant in S. neurona and absent in E. tenella (Fig. 1c), although this result is likely due to their divergence from canonical SRS domains, as similar families of 6-Cys rich proteins occur in other coccidian parasites and a related family is found in Plasmodium. SRS proteins and related 6-Cys proteins share a common extracellular structural domain, are typically GPI-anchored, and are thought to play diverse roles in cell attachment, invasion and development.

9From the predicted proteomes, we also reconstructed common metabolic pathways, which were highly conserved across T. gondii, H. hammondi and N. caninum, as noted previously. Expanding this analysis to include the 16 reference strains of T. gondii identified paralogues for certain functions, for example, in the pyrimidine and purine metabolic pathways and fatty-acid biosynthesis (Supplementary Fig. 2). Most enzymes involved in energy metabolism were well conserved with few paralogues and/or non-synonymous polymorphisms. Previous studies have established the ability of these pathways to mediate strain-specific growth differences. It is therefore interesting to note that several enzymes with the capacity to modulate flux within these pathways were associated with paralogous expansions and/or significant numbers of non-synonymous polymorphisms (Supplementary Fig. 2, Supplementary Data 3).

10Finally, we compared the position of genes across the chromosomes to establish the extent of synteny (Fig. 1d, Supplementary Fig. 3, Supplementary Table 2). There was a high degree of conservation of chromosomal position of orthologous genes between T. gondii and H. hammondi, and this only slightly decreased when they were compared with N. caninum, as previously reported for comparisons of T. gondii and N. caninum (Fig. 1a). Analysis of the more complete H. hammondi genome provided here revealed that it shares 29 long syntenic blocks with T. gondii harbouring >80% of its genes, with only a few blocks rearranged, most notably a ∼1 Mb reciprocal translocation between chromosomes Ia and IX (Fig. 1d). In contrast, synteny broke down substantially when these three organisms were compared with S. neurona and was completely absent when compared with E. tenella, as has been described previously for a pairwise comparison of T. gondii and E. tenella (Fig. 1d). The loss of synteny since the divergence of enteric from tissue-cyst forming coccidians stands in stark contrast to the conservation of synteny in the kinetoplastidae, fungi and chordates, all groups with greater evolutionary divergence times relative to the coccidians.

11To highlight key features of the T. gondii genome we depicted the coding capacity of the reference ME49 strain as a Circos plot, where the outermost circle indicates the genes encoded by each of the 14 chromosomes (Fig. 2a, Supplementary Fig. 4). By comparing the average sequence read depth across the genome, we identified chromosomal genes with copy-number variation (CNV) (Fig. 2a, second innermost circle). Expanding this analysis to all 62 strains revealed 14 genes that have evidence of CNV in all strains, and 39 genes with CNV in 90% of the strains (Supplementary Fig. 5, Supplementary Data 4). Examination of patterns of CNV also revealed several examples of large segmental duplications or aneuploidy in specific strains (Supplementary Fig. 6), similar to reports from previous genetic crosses. These regions were genetically homogeneous suggesting they arose by duplication events and are not hybrids created by unequal crossing over at meiosis. The significance of these diploid regions is uncertain, although recent studies in yeast indicate that aneuploidy can accelerate evolutionary adaptation.

12Recent comparisons of the draft genomes of T. gondii and H. hammondi, and the published genome of N. caninum, highlighted the expansion of gene families that differ between these otherwise closely related species. Using the newly annotated assembly of the ME49 genome obtained here, and data from 61 additional genomes of T. gondii, we expanded these analyses to examine the distribution of amplified genes, evident either as CNV or tandem arrays in the assemblies. These amplified genes were plotted as concentric coloured circles corresponding to the protein families they belong to and using symbols proportional to their total copy number (Fig. 2a). Many of these amplified genes encode secretory or surface proteins that have been previously implicated in host pathogenesis, referred to here as SPDs. These SPDs include genes encoding secretory proteins found in micronemes (MICs), dense granules (GRA), ROPs, as well as the SRS super family (Fig. 2a, Supplementary Fig. 4). Members of these protein families are known to mediate host cell attachment (MICs), modification of host immunity (GRA and ROP proteins) or adherence and immune evasion (SRS). Within the ME49 reference genome, we detected a total of 57 gene loci with CNV, which contain 176 gene copies, and 95 loci in tandem arrays, which contain 264 gene copies (Supplementary Data 4). Both CNV and tandemly duplicated genes were enriched in SPDs, in particular in genes encoding SAG/SRS, ROP and MICs (Fig. 2b,c, Supplementary Data 4), a pattern also noted previously. For example, SPDs comprised 52 (35%) of the 152 expanded loci and 196 (45%) of the 440 expanded gene copies, despite making up only 375 (4.5%) of the 8,311 total genes in the ME49 genome. Many of the SPDs also show evidence of positive selection, evident in elevated frequencies of non-synonymous (dN) versus synonymous (dS) mutations (Fig. 3a). Among these, the GRA, ROP and SAG genes show some of the highest levels of dN/dS, while metabolic enzymes typically show selection for conservation, as seen by low levels of dN/dS (Fig. 3a, Supplementary Data 5).

13We expanded the analysis of SPDs to examine their diversity among a set of reference genomes representing the 16 major haplogroups (Fig. 3b). OrthoMCL clustering of the SPD families revealed that while most members were represented in all 16 haplogroups, differences in representation and copy number were most evident in the SRS and ROP families (Fig. 3b). Collectively, these analyses reveal that the major difference between T. gondii strains is the diversification of SPD family members.

14Comparison of orthologues for GRA, ROP and SRS genes between T. gondii, H. hammondi and N. caninum revealed substantial differences in clustering by OrthoMCL, suggesting that the divergence among these genes may underlie biological differences between these species (Fig. 3b, Supplementary Fig. 7). In contrast, MIC genes were highly conserved, suggesting these organisms use a similar repertoire of host receptors (Fig. 3b). Comparison of OrthoMCL groupings also identified a number of putative species-specific genes unique to N. caninum, H. hammondi or T. gondii (Supplementary Fig. 8a, Supplementary Data 6). Further analysis indicated that a subset of the putative species-specific genes represent distant orthologues that are classified as separate groups by OrthoMCL (Supplementary Fig. 8b, Supplementary Data 6). Notably, this distant orthologue category is greater when comparing N. caninum to either T. gondii or H. hammondi, versus the pairwise comparison between the later two species. Among these distantly related orthologues, a number encode TgFAM or SRS genes, consistent with the idea that they influence important aspects of the biology (Supplementary Data 6). In contrast, a large number of the genes that differ between T. gondii and H. hammondi show evidence of alternative gene models, including early truncations, premature stop codons and frame shifts (Supplementary Fig. 8b, Supplementary Data 6). In addition, a smaller number of genes were present only in one species and are predicted to be unique, the majority of which were annotated as hypothetical unknowns. Analysis of alternative allele frequencies, RNA-seq data, and sequencing depth coverage failed to find evidence that these predicted differences are due to sequencing or assembly errors and instead suggest that many are genuine (Supplementary Fig. 8 c–f). Consequently, the putative unique gene list provided in Supplementary Data 6 provides a tentative starting point to identify genes that may mediate important biological differences between these closely related species.

15In addition to the previously recognized SPDs, we identified families of genes that are uniquely enriched in the T. gondii genome, referred to here as TgFAM genes (Fig. 3c, Supplementary Data 7), including one previously referred to as Toxoplasma-specific family (TSF), which corresponds to TgFAMC here. Our analysis of multiple T. gondii genomes reveals a much broader set of Toxoplasma-specific families (TgFAMs) (Supplementary Data 7), five of which we have specifically highlighted for their unique domain structures (Fig. 3c). Several TgFAMs are expanded and show evidence of CNV and/or tandem duplication, while others are located at the ends of chromosomes (Fig. 2), as previously noted for the TSF family. This pattern of telomeric clustering has also previously been associated with antigenic variant surface adhesins in Plasmodium. Although T. gondii is not known to undergo antigenic variation, the variable domains of the TgFAMs may represent adaptations to enhance host cell recognition and/or escape immune detection. We have highlighted five of the TgFAMs here based on the fact that they contain conserved signal peptides as well as domain architectures that suggest they may encode surface proteins with extracellular domains that contain conserved protein motifs (Fig. 3c). TgFAM genes are expanded in T. gondii, although they are less common in H. hammondi and N. caninum (Fig. 3b, Supplementary Data 7). In particular, H. hammondi and N. caninum contain far fewer members of TgFAMA and TgFAMB, and TgFAMC appears to be largely absent in N. caninum (Fig. 3b, Supplementary Fig. 7). Notably, many of the TgFAM genes highlighted here are expressed during sexual development in the cat gut or in oocysts that are shed into the environment following the sexual phase (Fig. 3d), suggesting they may play roles during transmission. In addition to the TgFAM genes highlighted here, there are a number of other gene families containing parasite-specific motifs that are expanded in T. gondii, and which may contribute to important biological traits not yet identified (Supplementary Data 7).

16Previous studies have reported the influence of recombination on the global population structure of T. gondii, which shows marked geographic segregation of major haplogroups, although the factors shaping these patterns remain unresolved. To examine the population structure based on genome-wide polymorphism data, we analysed single nucleotide polymorphisms (SNPs) that were defined by comparison of 61 T. gondii strains to the reference ME49 genome and filtered this set to include positions where reliable data were available for all strains (a total of 802,764 positions in each genome (Supplementary Data 8)). Generation of a neighbour network for these data revealed that the 62 strains group closely with haplogroups and major clades that were previously defined by lower resolution genotyping (Fig. 4a). Importantly, similar groupings were defined using admixture (Supplementary Fig. 9) and principal components analysis (Supplementary Fig. 10). Collectively, these findings support a population structure consisting of a small number of clades that show strong geographic segregation, as described previously.

17Although the neighbour network permits visualization of gene flow along several pathways, it does not fully capture the extent or pattern of local genomic admixture among any given pair of strains. To illustrate this more directly, we generated pairwise SNP diversity plots for three of the haplogroups contained in clade D, comparing them to the reference strain ME49 (Fig. 4b). Strains like ARI (haplogroup 12), a sister group of type 2 that is also found in North America, contain large haploblocks that are similar to ME49 (∼60%), interspersed with regions that are divergent (Fig. 4b), consistent with previous findings that these two groups are closely related. In contrast, TgCtPRC2 (haplogroup 13), which is a common clonal genotype in China, shares fewer regions with ME49 (∼40%) and COUG (haplogroup 11), which represents a rare North American lineage found in wild animals, showed almost no conserved regions with ME49 (<1%) (Fig. 4b). Thus although members of a common clade contain distinct genomic patterns that have arisen by different evolutionary paths, it is striking that many share large conserved haploblocks across their genomes.

18To better represent the shared ancestry across strains, we analysed local inheritance patterns using chromosome painting to reveal patterns of local admixture. When strains were aligned by clade, the presence of shared haploblocks across members was evident by common colour patterns (Fig. 4c, Supplementary Fig. 11). These shared regions represent chromosomal haploblocks that show a high degree of shared ancestry, in some cases eroding the boundaries of the clade structure. Noteworthy, this analysis also revealed patterns of local admixture that suggest the occurrence of genetic crosses among strains of different clades, likely favoured by their geographic proximity (Fig. 4c, Supplementary Fig. 11). Similar shared chromosomal haploblocks are also seen in average pairwise plots for SNPs among members of individual clades (Supplementary Fig. 12). Among the most strongly conserved haploblocks is chromosome Ia, which is shared across nearly all clades, with the exception of clades E and F (Fig. 4c,d, Supplementary Fig. 13). The basis for the widespread conservation of chromosome 1a (ref. ) is uncertain, but recent studies suggest that it may be due to the enhanced transmission in domestic cats.

19The analyses presented above suggest that common inheritance of large haploblocks is the major factor in determining the phylogenetic grouping of T. gondii strains. To test this model, we performed two types of analysis to cluster strains using the SNP data. First, we analysed SNPs using the linkage model of ChromoPainter in FineStructure to generate a clustering hierarchy. This model, which combines information across linked markers in a co-ancestry matrix, recreated the clades seen in the previous analysis with several minor exceptions (Supplementary Fig. 14). Separately, we analysed the SNPs using a rolling window method to define an overall similarity index based on how many regions were co-inherited between all pairwise comparisons (1,953 unique comparisons), which produced a highly similar clade structure (Fig. 4e, Supplementary Fig. 15). These analyses reveal that the current population structure is defined by recent genomic admixture, where large chromosomal haploblocks have been inherited in common by members of individual clades. Although recent admixture had previously been suggested by analysing individual regions separately, the present genome-wide analysis of SNPs establishes that this pattern is a defining feature of the population structure of T. gondii.

20To further examine the pattern of long-haploblock inheritance, we compared the ancestry of regions that were conserved with those that were more variable. When SNP diversity was averaged for members of the same clade, it emerged that discrete regions of the genome show very low SNP diversity, while others are highly variable (Fig. 5a). Regions of low average pairwise SNP diversity were observed in all clades, but differed in their frequency and location (Supplementary Fig. 12). To compare the ancestry of different regions of the genomes, we partitioned the genome into two segments based on regions that exhibited low SNP diversity in at least one clade (defined as the union of all regions that were ‘conserved' in at least one clade) versus regions that showed high SNP diversity in all clades (defined as ‘non-conserved') (Supplementary Data 9 and 10). We then compared unrooted phylogenetic trees for the conserved versus non-conserved regions using a Robinson–Foulds distance metric, which measures the degree of difference, or distance, between the two sets of trees. This analysis revealed that the ancestry of the conserved haploblocks was significantly different than that of the non-conserved regions for all 14 chromosomes (Fig. 5b). In addition, we generated neighbour networks based on the conserved versus non-conserved regions (Supplementary Fig. 16). The conserved region network most closely resembled the total SNP network (Fig. 4a) and the network based on the non-conserved regions grouped most strains in similar clades, with several notable exceptions (Supplementary Fig. 16). These findings illustrate the importance of the conserved blocks in influencing the grouping of genotypes into clades.

21To determine the influence of these shared haploblocks on the content of genes found within specific clades, we plotted the distribution of SPDs found within conserved regions shared by members of specific clades (Fig. 5c, Supplementary Data 11). The distribution of SPDs revealed that a number of known pathogenicity determinants were common to conserved regions in specific clades (Fig. 5c). We tested whether these patterns were random, or if they showed specific enrichment of SPDs within conserved regions. When clades A, B, C, D and F were analysed together, the pattern of clustering of SPDs was highly significant (P≤0.005) allowing us to reject the null hypothesis that they were randomly distributed across the genome. SPDs were also significantly clustered in conserved regions when separately analysing clade C (P≤0.005), clade D (P≤0.05), and clade F (P≤0.000001), while clade B was suggestive (P≤0.08) and clade A was not significant (P≤0.2). The failure to observe a significant clustering of SPDs in clade A may be due to its considerable substructure that suggests it may actually be comprised of two or more distinct groups. Nonetheless, it is clear that SPDs are often clustered and are found with increased frequency in conserved, shared regions of the genome. SPDs within these regions also share the recent ancestry of the surrounding conserved regions when analysed using phylogenetic trees (Supplementary Figure 17). These findings are consistent with the hypothesis that recent inheritance of conserved blocks containing specific SPDs is associated with successful expansion of specific lineages and suggest that SPDs impart a selective advantage to members of specific clades.

22The co-inheritance of SPDs within conserved regions provides a tentative list of candidates for further study of genes that may underlie important biological traits shared by specific clades. Among these are a number of SPDs previously implicated in acute virulence in the mouse: for example ROP17 (ref. ) found in conserved regions in clades A and B, ROP5 (ref. ) found in conserved regions in clades A and C, and GRA3 (ref. ) found in conserved regions in A, C, D, and F (Fig. 5c). Low diversity regions also contain a number of SRS genes, encoding immunologically dominant surface proteins, which have previously been implicated in host cell invasion including sporozoite SAG (Sp-SAG also known as SRS28) found in conserved regions in clades A and B, SAG3 (SRS57) found in conserved regions in clades A and C, and SAG2A (SRS34A) (M.E.G., unpublished) found in conserved regions in clades C and D (Fig. 5c, Supplementary Data 11). A second pattern that emerges from this analysis is the presence of clusters of SPDs that are clade-specific, for example clusters of SRS and TgFAM genes on chromosomes IV, V and IX found in clade C, and clusters of various SPDs on chromosomes II, IX, XI and XII in clade D (Fig. 5c, Supplementary Data 11). Although the specific roles of these genes are unknown, they may underlie common traits that distinguish phenotypes characteristic of specific clades.

23Toxoplasma gondii belongs to a diverse and ancient phylum of parasites that antedates the wide range of vertebrate hosts that they currently inhabit. It shares a core set of genes and metabolic processes with closely related tissue-cyst forming coccidian parasites. Despite having similar genomic content, these organisms differ dramatically in their host range, pathogenicity and modes of transmission. We demonstrate here that T. gondii is demarcated from its closest relatives by the expansion of parasite-specific SPDs that are involved in host–pathogen interactions. Diversification of SPDs also highlights key differences among major clades of T. gondii, which are distinguished by common inheritance of large haploblocks in their genomes. Shared inheritance of large haploblocks among related strains reinforces the hypothesis that recombination in the wild, while infrequent, drives important biological adaptations. The distribution of clustered SPDs within conserved regions that show common ancestry identifies a number of candidate genes that may influence both clade specific and more broadly shared traits. Overall, the phenotypic traits of individual strains are likely determined by both their core ancestral genomes, and inheritance of conserved haploblocks, which together comprise their mosaic genomes.

24The mosaic genomic patterns seen in specific clades may underlie differences in population structure that exist in different T. gondii populations between North and South America. Although the common ancestry of conserved blocks among otherwise different genotypes is consistent with recent introgression, the conservation of these regions among members of a given clade may reflect several different mechanisms. Such shared haploblocks may be retained in the face of ongoing recombination in outbreeding populations, suggesting they impart a selective advantage. Alternatively, they may simply reflect recent admixture that has not been eroded due to infrequent recombination, such as in clonal populations. Regardless of their exact histories, strains that inherit conserved haplotype blocks in common will also share clusters of highly related genes, including SPDs that may influence traits such as transmission, host range and pathogenesis.

25Expansion of polymorphic genes that are important in pathogenicity is also a key feature of other pathogen genomes. One feature they share in common is that the amplified genes typically encode surface or secretory proteins that interact directly with the host, either to mediate attachment or immune evasion. Examples include: the expansion of surface antigen variants encoded by VAR genes in the Plasmodium falciparum genome, and the unrelated yet expanded VIR genes in Plasmodium vivax, variant surface glycoprotein encoding genes (VSG) in Trypanosoma brucei, and the expansion of RXLR effectors in oomycetes. It is noteworthy that while gene expansion and diversification are common to each of these examples, the protein families involved are largely distinct and reflect the specialized biology of these diverse pathogens. This pattern suggests that expansion of polymorphic gene families is a common theme that underlies important changes in host range and transmission that characterize the evolution of pathogens in their diverse hosts.

26Sixty-two representative strains of T. gondii were selected from different haplogroups from around the world (Supplementary Data 1). Strains were cultured in human foreskin fibroblast cells, as described previously.

27Sequencing of 16 reference strains of T. gondiii, and of one isolate of H. hammondi, was conducted using a combination of 454, and Illumina PE sequencing technologies (Supplementary Table 1). Sequence reads were screened for contamination and reassembled using Celera Assembler software or Newbler v2.6 (ref. ). Scaffolds were then aligned with MUMmer to T. gondii ME49 chromosome sequences from ToxoDB v8.0 (http://ToxoDB.org) to generate super-scaffolds spanning entire chromosomes. Annotated genomes were deposited into National Center for Biotechnology (NCBI).

28For each of the remaining 46 non-reference strain (Supplementary Data 1), a single Illumina PE barcoded library was prepared from tachyzoite gDNA. Libraries were then pooled into groups of nine samples and sequenced multiplexed in a single lane of an Illumina HiSeq 2000 machine. Sequencing reads were deposited in the Sequence Read Archive repository at NCBI.

29To aid in the curation of ME49 gene models, two tachyzoite-specific Illumina complementary DNA libraries were constructed from mRNA isolated from tachyzoite cultures from ME49 and GT1 strains. Each library was then sequenced in a single lane of an Illumina Genome Analyzer II machine.

30Gene annotations were derived by comparison of the existing ME49 reference genome (http://ToxoDB.org) using a combination of evidence from RNASeq data, cDNA/EST sequences and a variety of software tools to predict potential protein-coding genes using an in-house pipeline at J. Craig Venter Institute (JCVI) (Supplementary Methods). Predicted proteins were run through JCVI's autonaming pipeline, that assign product names based on a number of sequence similarity searches including blastp searches against the previous T. gondii ME49 proteome (ToxoDB v8.0; http://ToxoDB.org) and the GenBank non-redundant protein database, HMM searches against Pfam and TIGRfam databases, and RPS-Blast searches against the NCBI-CDD database. Proteins without any significant hit to other proteins or protein domains were flagged as ‘hypothetical protein'. The final list of product names was then curated by researchers from the Toxoplasma research community before being assigned to working models. ME49 protein-coding genes were assigned similar pub_locus identifiers to the previous genome assembly while newly predicted protein-coding genes were assigned completely new pub_locus identifiers (Supplementary Methods).

31Functional annotation of protein-coding genes in other T. gondii reference strains was performed as above. Genes syntenic to ME49 inherited their product names, GO terms, and Enzyme Commission numbers, while non-syntenic genes acquired their names and other functional annotations from the output of JCVI's autonaming pipeline. Structural and functional annotations of H. hammondi were carried out following a similar approach with slight modifications (Supplementary Methods).

32To identify known protein domains the T. gondii ME49 proteome was searched against Pfam and TIGRfam HMM profiles using HMMER3 ref. ). Proteins matching a particular HMM profile were assigned to that domain and remaining peptides searched against each other using blastp to identify potential novel domains. The top five protein families containing novel para domains (TgFAMs A to E) were analysed using Phobius to identify signal peptides and transmembrane domains. De novo identification of conserved protein domains across members of the same gene family was carried out with MEME. Expression levels for the TgFAMs were obtained from T. gondii Affymetrix Array data available from NCBI GEO records GSE32427 and GSE51780.

33Coding sequences from each cluster of orthologous genes from the 16 T. gondii reference strains were used to estimate dN/dS ratios using a modified version of the Bioperl script bp_pairwise_kaks.pl (http://search.cpan.org/dist/BioPerl/scripts/utilities/bp_pairwise_kaks.pl).

34Illumina reads for each of the 61 other genomes were aligned using Bowtie2—end-to-end against the ME49 reference genome assembly (release date 23 April 2013), identifying a total of 2,342,433 SNPs across all strains. Positions with informative base calls for all 62 strains were identified, generating a final list of 802,764 SNPs that were used for analysis.

35Annotated proteomes were analysed using OrthoMCL v2.0 (ref. ) to define orthologous groups. Clusters of orthologous groups were functionally annotated using GO Slim terms, which are designed to group the many different GO terms into smaller groups of related processes (http://geneontology.org/page/go-slim-and-subset-guide). The proteomes were queried against the Pfam HMM database using HMMER3 to estimate the abundance of Pfam domains.

36SNPs from the 16 T. gondii reference strains that correspond to the 382 proteins in the iCS382 metabolic pathway reconstruction of T. gondii ME49 (ref. ) were downloaded from ToxoDB (http://www.toxodb.org/toxo-release4-0/home.jsp).

37Genome-wide SNPs were saved as FASTA files and directly incorporated into SplitsTree v4.4 (ref. ) to generate unrooted phylogenetic networks using a neighbour-net method and 1,000 bootstrap replicates. Principal components analysis was performed by eigenanalysis of a co-ancestry matrix implemented in fineSTRUCTURE, as described in ref..

38SNP data for ChrIa were plotted as a minimum spanning tree using SplitsTree v4.4 (ref. ) with 2,000 spring-embedded iterations. The 62 strains were clustered into four major groups denoted as monomorphic, divergent, 5′-chimeric and 3′-chimeric chromosome Ia. SNPs present in each cluster were calculated using a custom script over a 10-kb moving window and plotted using Excel.

39The population genetic structure of T. gondii was determined by an unsupervised clustering algorithm, ADMIXTURE with ancestral clusters set from k=1 through 10. The number of ancestral clusters k was determined by estimating the low cross-validation error (CV error) for different k values using five-fold CV.

40We developed a co-ancestry heatmap by using the linkage model of ChromoPainter (http://www.paintmychromosomes.com) and fineSTRUCTURE based on the genome-wide SNP data. The burn-in and Markov Chain Monte Carlo (MCMC) after the burn-in were run for 10,000 iterations with default settings.

41For each strain of T. gondii, the respective.sra files were used to align reads to the 14 ME49 reference chromosomes using Bowtie2 with the end-to-end option. The read depth per base pair, or read bases (RB), across 8,320 chromosomal-mapped genes was determined using samtools mpileup. Plots were generated in R (http://www.r-project.org/). T. gondii gene families organized in tandem arrays were identified with an in-house perl script.

42Genes found to be specific to T. gondii, H. hammondi or N. caninum based on OrthoMCL clustering were further analysed using a combination of sequence alignment tools (Supplementary Methods). Genes were classified based on whether they showed a significant ‘Blastp hit', showed blastn similarity that was either ‘full length' or constituted an ‘alternative gene model', or showed no similarity and were ‘unique'. Further analysis was done to investigate these differences by RNA-Seq, analysis of minimum alternative allele frequency and minimum read depth (Supplementary Methods).

43To determine the extent of recombination and co-inheritance of blocks between strains, low SNP regions (regions of recent co-inheritance or shared blocks) were identified for 10 kb windows for pairwise strain comparisons. A heatmap was generated using the R function heatmap.2 (gplots library (http://www.r-project.org/)) with hierarchical clustering on the % shared blocks value. The number of SNPs per 10 kb window were averaged for all strains within a Clade, and chromosomal regions with low SNP density were identified as above using 10 kb windows that had three or fewer SNPs across a continuous stretch of 10 windows (100 kb), allowing for intermittent outliers.

44We identified genes that belong to the SPD families (that is, MIC, GRA, ROP, SRS and TgFAM) based on the annotation of ME49 accounting for CNV in determining the gene number. We then mapped the position of the SPDs onto the assembled ME49 genome and defined those that fell into conserved or non-conserved regions. To determine if gene type was independent of region type we compared the observed frequency of SPDs and non-SPD genes in conserved versus non-conserved regions of the genome using a χ2-squared analysis. The null hypothesis was that the distribution would be random, and there would be no difference between observed and expected. A P value of ≤0.05 was considered significant cause for rejection of the null hypothesis.

45Phylogenetic trees for the conserved and non-conserved regions were constructed using maximum likelihood as implemented in RAxML version 7.3.0 with the GTR+GAMMA model. Standardized Robinson–Foulds distances were calculated between the conserved and non-conserved trees based on 500 bootstrap replicates. Trees were considered congruent if they had no conflicting branches with bootstrap support of >95%.

46Phylogenetic trees were constructed for the conserved OrthoMCL OG5_0126701 using the Neighbour-Joining algorithm with 1,000 bootstrap replicates as implemented in Geneious ver. 7.1.5 (http://www.geneious.com (ref. )) and visualized with FigTree ver. 1.4.0 (http://tree.bio.ed.ac.uk/software/figtree/).

47The OrthoMCL ortholog clusters (see above) were reformatted to represent each pair found in the cluster outside of self-matches and syntenic blocks were generated between all combinations of genomes as described in ref..

48Local admixture analyses using an enhanced ADMIXTURE algorithm was used to assign each of the 62 strains to clusters representing these ancestral states.

49Detailed methods for the above sections can be found in the supplementary methods.

50Accession codes: Assembled and annotated sequences for the genomes of the T. gondii reference strains have been deposited in GenBank at NCBI under accession codes GT1 (GCA_000149715.2), VEG (GCA_000150015.2), ME49 (GCA_000006565.2), MAS (GCA_000224865.2), RUB (GCA_000224805.2), CAST (GCA_000256705.1), TgCtBr5 (GCA_000259835.1), P89 (GCA_000224885.2), VAND (GCA_000224845.2), COUG (GCA_000338675.1), ARI (GCA_000250965.1), TgCtPRC2 (GCA_000256725.1), GAB2-2007-GAL-DOM2 (GCA_000325525.2), TgCtCo5 (GCA_000278365.1), FOU (GCA_000224905.2) and TgCtBr9 (GCA_000224825.1). Assembled and annotated sequences for the genomes of the H. hammondi reference strains have been deposited in GenBank at NCBI under accession code GCA_000258005.2. The annotation of reference genomes are also available at ToxoDB (http://toxodb.org/common/downloads/WhitePaperProvisionalGenomes/): RNA-Seq reads for T. gondii strains have been deposited in the Sequence Read Archive at NCBI under accession codes GT1 (SRX099798) and ME49 (SRX099799). Sequence reads from the remaining T. gondii strains have been deposited in the Sequence Read Archive at NCBI under the following accession codes RH-88 (SRX160127, SRX160126), RH-JSR (SRX159850, SRX159851), TgCkCr1 (SRX160807), SOU (SRX160119), TgCkGy2 (SRX099785), CASTELLS (SRX099789), TgCatBr1 (SRX099790), TgCatBr18 (SRX099794), TgCatBr34 (SRX160051), BRC TgH 18002 GUY-KOE (SRX099796), BRC TgH 18003 GUY-MAT (SRX099783), GUY-2003-MEL (SRX160131), GUY-2004-ABE (SRX160132), TgRsCr1 (SRX160143), TgCkBr141 (SRX160124), TgCatBr25 (SRX160134), TgCatBr10 (SRX099791), TgCatBr64 (SRX160743), BRC TgH 18001 GUY-DOS (SRX099782), TgCkCr10 (SRX099784), B41 (SRX099774), RAY (SRX099793), TgCatBr44 (SRX160141), GAB1-2007-GAL-DOM10 (SRX159841, SRX159839), GAB3-2007-GAL-DOM9 (SRX160125), GAB5-2007-GAL-DOM6 (SRX159842, SRX159840), TgCatPRC3 (SRX09979), BRC TgH 18009 (SRX171132), BRC TgH 18021 (SRX160041), BRC TgH 20005 (SRX157499, SRX157465), BRC TgH 21016 (SRX099775), GUY-2004-JAG1 (SRX099776), B73 (SRX159844, SRX159843), PRU (SRX099792), M7741 (SRX159890, SRX159849), ROD (SRX160129, SRX160128), TgShUS28 (SRX160130), BOF (SRX099774), GAB5-2007-GAL-DOM1 (SRX160069), TgCatBr26 (SRX099780), TgCatBr72 (SRX160049), TgDogCo17 (SRX099787), TgCatBr15 (SRX099779), G662M (SRX160052), GAB3-2007-GAL-DOM2 (SRX160123), TgCatBr3 (SRX160142) and TgH26044 (SRX160050).

51How to cite this article: Lorenzi, H. et al. Local admixture of amplified and diversified secreted pathogenesis determinants shapes mosaic Toxoplasma gondii genomes. Nat. Commun. 7:10147 doi: 10.1038/ncomms10147 (2016).