|Home | About | Journals | Submit | Contact Us | Français|
Plasmodium knowlesi is an intracellular malaria parasite whose natural vertebrate host is Macaca fascicularis (the ‘kra’ monkey); however, it is now increasingly recognized as a significant cause of human malaria, particularly in southeast Asia1,2. Plasmodium knowlesi was the first malaria parasite species in which antigenic variation was demonstrated3, and it has a close phylogenetic relationship to Plasmodium vivax4, the second most important species of human malaria parasite (reviewed in ref. 4). Despite their relatedness, there are important phenotypic differences between them, such as host blood cell preference, absence of a dormant liver stage or ‘hypnozoite’ in P. knowlesi, and length of the asexual cycle (reviewed in ref. 4). Here we present an analysis of the P. knowlesi (H strain, Pk1(A+) clone5) nuclear genome sequence. This is the first monkey malaria parasite genome to be described, and it provides an opportunity for comparison with the recently completed P. vivax genome4 and other sequenced Plasmodium genomes6-8. In contrast to other Plasmodium genomes, putative variant antigen families are dispersed throughout the genome and are associated with intrachromosomal telomere repeats. One of these families, the KIRs9, contains sequences that collectively match over one-half of the host CD99 extracellular domain, which may represent an unusual form of molecular mimicry.
The P. knowlesi genome sequence was produced by whole-genome shotgun sequencing to eightfold coverage, with targeted gap closure and finishing (Supplementary Table 1). The 23.5-megabase (Mb) nuclear genome is composed of 14 chromosomes and contains the expected complement of non-coding RNA (ncRNA) genes with known function (Supplementary Table 2) and a large number of novel structured ncRNA candidate genes (Supplementary Figs 1-5 and Supplementary Tables 3 and 4). The presumed centromeres are similar to those found in other Plasmodium species4,6, and are positionally conserved within regions sharing synteny with P. vivax (see Fig. 1 of ref. 4). The overall G+C base composition is 37.5%. A total of 5,188 protein-encoding genes were identified, which is slightly lower than the predicted proteome size of P. falciparum and P. vivax4,6.
Unusually for Plasmodium species, (G+C)-rich repeat regions containing intrachromosomal telomeric sequences (ITSs, containing the heptad sequence GGGTT[T/C]A) are found at multiple internal sites in the P. knowlesi chromosomes, arrayed tandemly or as components of larger repeat units (Fig. 1). These sequences appear infrequently in P. vivax and P. falciparum at internal chromosome sites (Supplementary Figs 6 and 7). In the protozoan parasite Trypanosoma brucei10, ITSs may be the templates for recombination events that result in gene conversion among variant antigen VSG genes11. In mammalian genomes12, ITSs are common and may represent the ‘scars’ of double-stranded DNA break repair12. Alternatively, ITSs may have a role in transcriptional control.
For approximately 80% (4,156 out of 5,185) of predicted genes in P. knowlesi, orthologues could be identified in both P. falciparum and P. vivax (for details, see ref. 4). The P. knowlesi-specific variant antigen gene families, SICAvar genes13 and kir genes9, form the largest groups of P. knowlesi-specific expansions (Supplementary Tables 5 and 6). Five distinct gene families of unknown function, with 4-15 paralogous members, are unique to P. knowlesi (referred to as Pk-fam-a to Pk-fam-e in Supplementary Table 7). Pk-fam-a and Pk-fam-b each have more than nine paralogous members (Supplementary Fig. 8), which have a two-exon gene structure with a signal peptide, a carboxy-terminal transmembrane region, but lack typical export motifs14,15. Members of the protein family Pk-fam-c and Pk-fam-e represent two new families with putative protein export signals (Supplementary Fig. 8 and Supplementary Table 8).
A comparison of Pfam domains16 between the predicted proteomes of P. knowlesi, P. vivax and P. falciparum (Supplementary Table 9, Supplementary Information) revealed major differences in domains that distinguish species-specific protein families involved in antigenic variation. The remainder of the proteome was relatively conserved albeit with some interesting copy number variations of a few key housekeeping enzymes (Supplementary Fig. 9 and Supplementary Table 9).
In other Plasmodium genomes sequenced so far, variant gene families involved in antigenic variation (Supplementary Figs 6 and 7) are typically arranged in the subtelomeres, and only a few members of these families have hitherto been found at intrachromosomal sites. Notably, the P. knowlesi genome sequence has revealed that the major variant gene families (that is, SICAvar13 and kir9) are randomly distributed across all 14 chromosomes (Fig. 1) and often co-localize with ITS-containing repeats (Supplementary Information). Although all of the telomeres were not fully assembled, we know that in the case of chromosome 7, P. knowlesi and P. vivax have atypical gene content—the subtelomere encodes proteins associated with merozoite invasion (for example, MAEBL and members of the reticulocyte-binding-like (RBL) family) (Supplementary Fig. 10).
Variant SICA (schizont-infected cell agglutination) antigens on the surface of infected red blood cells5 are associated with parasite virulence17 and are encoded by the SICAvar gene family13—the largest variant antigen gene family in P. knowlesi. Switching of variant types underlies the establishment of a chronic infection in the vertebrate host, a process that is essential in all species, to ensure mosquito transmission and the completion of the life cycle. Full-length SICAvar genes have 3-14 exons (Supplementary Table 5 and Supplementary Fig. 11), resulting in a range of sizes for the predicted proteins of 53-247 kDa. Although many of the SICAvar genes are present only as fragments, we estimate that there are up to 107 members in the H strain of P. knowlesi based on the number of conserved final exons.
Twenty-nine predicted SICAvar genes have complete gene structures and were divided into two subtypes (Fig. 2). The type I SICAvar genes with 7-14 exons predominate, with a few containing unusually long introns (Fig. 2). The type II subgroup represents small SICAvar genes with 3-4 exon structures. Unusually large introns (5.8-13.6 kb) are a unique feature of SICAvar genes and have not previously been seen in any other sequenced apicomplexan gene (Fig. 2).
SICA antigens have a modular structure (Fig. 3, Supplementary Fig. 12) comprising a variable number of highly diverged cysteine-rich domains (CRDs) encoded by multiple exons, a transmembrane domain and a cytoplasmic domain. A high level of sequence diversity was observed, with the exception of the 3′ terminal exon13.We investigated the domain organization of the CRDs using profile hidden Markov models (HMMs; Fig. 3 and Supplementary Fig. 13). The full-length SICA proteins contain a distinct five-cysteine CRD (termed SICA-α) at the amino terminus, which occurs once or twice and may have a stabilizing role analogous to the cysteine-rich N-terminal capping motifs of extracellular leucine-rich repeat proteins18. There are 1-8 CRDs (referred to as SICA-β) with 7-10 conserved cysteine residues. The transmembrane domain and a conserved domain follow at the C terminus (termed SICA_C in Supplementary Figs 12 and 13).
Although P. knowlesi and P. falciparum are phylogenetically distant, the SICA and P. falciparum erythrocyte membrane protein 1 (PfEMP1) variant antigens share many fundamental biological characteristics (reviewed in ref. 19). Common regulatory mechanisms involving post-transcriptional gene silencing have been proposed between the var gene family in P. falciparum and the SICAvar family in P. knowlesi19. We have identified conserved sequence motifs between the single var intron and SICAvar introns (Supplementary Figs 14-18) in the region thought to be the origin of a ncRNA transcript involved in the silencing of var genes20, indicating possible commonality in regulatory mechanisms.
We searched for evidence of gene conversion within the SICAvar family, using the predicted sequences of 20 type I full-length SICAvar genes (Supplementary Information). It is clear that exon shuffling has an important role in SICAvar evolution13. The low-complexity repeat regions found within introns might facilitate recombination through misalignment during mitosis; this could explain the presence of SICAvar fragments found throughout the genome and/or SICAvar gene models with partial intron/exon structures. These comprise whole, and apparently intact, exons that might provide a reservoir for diversification analogous to that seen with VSG genes in Trypanosoma brucei11 (Supplementary Information).
Kirs represent the second largest variant gene family. They encode predicted proteins of 36-97 kDa that are hypothesized to be expressed at the surface of infected erythrocytes and undergo antigenic variation9. There are 68 predicted kir genes, 4 of which have incomplete structures (Supplementary Table 6). They were divided into four types depending on the number of exons (Supplementary Fig. 19). Most (58 out of 64) kir genes belong to types I and II. The domain organization of all predicted KIR proteins was also determined using profile HMMs (Fig. 3 and Supplementary Fig. 20). They contain 1-3 domains, followed by a transmembrane domain at the C terminus (referred to as KIR TM in Supplementary Fig. 20). A BLAST analysis of KIR proteins revealed stretches of up to 36 amino acids within the predicted extracellular domain that have 100% identity to host proteins, the most striking of which is to CD99. These matches were evident in several KIR proteins. Interestingly, different family members contain matches to different regions of CD99, such that together, they represent over one-half of the CD99 extracellular domain (Fig. 4). Tests were performed to assess the possibility that such matches could occur by chance (Supplementary Table 10). We have compared the sequences to Macaca mulatta, African green monkey and human. The matches exclude conserved cysteine regions and the degree of sequence identity decreases noticeably as the evolutionary distance to the natural host increases (Fig. 4 and Supplementary Table 10). CD99 has a critical role as a immunoregulatory molecule in T-cell function (see http://www.ncbi.nlm.nih.gov/omim/). These exact matches may interfere with recognition of parasitized erythrocytes by the host immune system or act as CD99 analogues that interfere by competing with T cells for CD99 partner molecules.
We undertook a more systematic search for other such instances of parasite proteins containing extensive stretches of identical host sequences, using the PMATCH algorithm (Supplementary Information). Unsurprisingly, a large number of matches to highly conserved housekeeping genes were observed, but in addition regions of perfect identity to another host protein (known as AHNAK, see http://www.ncbi.nlm.nih.gov/omim/) were detected in two KIRs and one SICA-like protein (Supplementary Fig. 21 and Supplementary Table 10). Analogous searches using the predicted exported protein repertoires (exportome) of P. vivax and P. falciparum found no such matches to host proteins (Supplementary Table 11). The identity to host proteins is maintained at the amino acid sequence rather than DNA sequence level (data not shown).
Acquisition of host proteins, and thus the ability to mimic their function, has been observed in many bacterial and viral pathogens21. In parasitic protozoa there are known cases where stretches of amino acids present on a parasite-encoded cell surface protein match perfectly to regions of host proteins22. However, in all such cases, the matches correspond to a common amino acid repeat that is shared between them22-24. Malaria parasites are known to have a potential immunomodulatory role either by secreting functional homologues of host molecules or by binding to host antigen-presenting cells25,26. This is the first observation of its kind in a malaria protein that shows acquisition of host peptide sequences that are likely to be on the infected cell surface and thus may interact with the host. The mechanism by which these host sequences have arisen remains to be clarified. Possible explanations include convergent evolution or horizontal transfer followed by gene degeneration events.
During the intraerythrocytic life cycle, malaria parasites significantly remodel the erythrocyte by exporting numerous proteins14,15. This depends on a short motif, termed the plasmodium export element (PEXEL) or vacuolar transport signal (VTS), which is present in over 300 P. falciparum proteins and is common to all Plasmodium species sequenced so far27. In addition to the members of the PHIST family27, an additional 100 proteins in P. knowlesi have typical PEXEL-like motifs (Supplementary Table 8 and Supplementary Fig. 22).
Like the PfEMP1 protein in P. falciparum, the SICAs and KIRs lack a signal peptide and a typical PEXEL-motif. We have identified a novel motif in the N-terminal region of SICA-α domains with a positionally conserved tryptophan residue surrounded by hydrophilic residues (Supplementary Fig. 22) that may be the export signal. Similarly, 75% of KIR proteins have a conserved Z-L-P-S motif (where Z denotes a hydrophilic residue) at the beginning of the KIR domain that may also facilitate export (Supplementary Fig. 22). In summary, approximately 280 predicted P. knowlesi proteins may be exported to the infected erythrocyte surface via the PEXEL-dependent or PEXEL-independent pathways. By comparison, the exportome of P. vivax is considerably larger than that of P. knowlesi and seems to be much bigger than previously thought27. About 145 P. vivax proteins contain typical PEXEL motifs including the members of the PHIST family and a small subgroup of 12 VIRs.
Genome sequencing of P. knowlesi and its comparison with other malaria genomes has highlighted several novel features of this emerging and potentially life-threatening human malaria parasite, and underscores the importance of full genome sequencing of new Plasmodium species. Major differences in both content and organization of its genome were revealed that involve the host-parasite interface, reinforcing the notion that malaria species have evolved specific mechanisms for enhancing their survival within their respective hosts. The P. knowlesi genome will also greatly enhance the utility of this human-infective species as a model for addressing questions pertinent to all Plasmodium species.
The random shotgun approach was used to obtain roughly eightfold coverage of the whole nuclear genome sequence from the erythrocyte stage of the Pk1(A+) clone of the H strain of P. knowlesi5. Sequence reads were assembled (as described in the Supplementary Information) and positional information from sequenced read pairs were used to resolve the orientation and position of the contigs. The assembled P. knowlesi contigs were iteratively ordered and oriented by alignment to P. vivax assembled sequences (described in ref. 4) and by manual checking. Automated predictions from the gene finding algorithms were manually reviewed by comparison to orthologues in other Plasmodium species. Artemis and Artemis Comparison Tool (ACT) were used (as described previously28) for annotation and curation and viewing the TBLASTX comparisons of regions with conserved synteny between P. knowlesi, P. vivax and P. falciparum. This also allowed us to curate gene models and identify local interruptions of synteny. Functional annotations were based on standard protocols as described previously6.
Supplementary Information is linked to the online version of the paper at www.nature.com/nature.
We acknowledge the support of the Wellcome Trust Sanger Institute core sequencing and informatics groups. The study was funded by the Wellcome Trust through its support to the Pathogen Sequencing Unit at the Wellcome Trust Sanger Institute. We thank J. Barnwell for providing the Pk1(A+) clone of the H strain of the parasite for the generation of genomic DNA by A. Thomas. We thank A. Voorberg-vd Wel (BPRC, Rijswijk) for technical assistance. We thank D. Fergusson for providing us with the electron micrograph image of the erythrocyte, used in Fig. 2. Part of this work was supported by the Netherlands Organization for Scientific Research, NIH, BioMalPar and the Virimal contract. This work is dedicated to the memory of Marie-Adele Rajandream.
Genomic DNA was isolated from blood drawn from an infected rhesus monkey at 10% ring stage parasitaemia. Blood was Plasmodipur-filtered five times to remove white blood cells and erythrocytes were lysed in 0.1% saponin. Total parasite DNA was isolated using the PUREGENE DNA isolation kit (Gentra Systems), according to the manufacturer’s instructions. All experimental animal work in these studies was carried out under protocols approved by the independent Institutional Animal Care and Use Committee and performed according to Dutch and European laws.
We sequenced the P. knowlesi genome from plasmid clones containing small fragments of up to 4 kb inserted into pUC19 vector. Problems associated with high G+C sequence were addressed by optimizing the sequence mixture. The quality of reads for the project was as follows: 97.6% of P. knowlesi reads had a quality score of (derived from the PHRED score generated by GAP429) >70 (P = 1 × 10-7). Regions containing repeat sequences or an unexpected read depth were manually inspected. In addition, a P. knowlesi fosmid library was constructed in pCC1FOS vector and end sequences were produced (10.5-fold clone coverage) to obtain paired-end information from 40-kb inserts. In particular, we re-examined regions with apparent breaks in synteny for potential misassembly errors and location of several intrachromosomal telomeric-repeat (GGGTT[T/C]A) sequences associated with SICAvar and kir genes. Sequence reads were assembled with PHRED/PHRAP on the basis of overlapping sequence and were manually edited in GAP4 database29. Information from oriented read pairs, together with additional sequencing from selected large-insert clones and synteny with P. vivax chromosomes, were used to resolve potential misassemblies. Using long-range sequence information from the fosmid end sequences, we were able to bridge 142 out of 190 total gaps (Supplementary Table 1).
Annotation (PK4 version of assembly) was performed using the Artemis30 and ACT software31. Genes were identified by manual curation of the output of the gene finding software SNAP32 and Annotaid (an extension of the comparative gene prediction program Projector33; I. M. Meyer, unpublished). A set of 100 manually curated Plasmodium knowlesi genes was used as the training set for SNAP predictions. Annotaid was optimized for genome-wide analysis by training its parameters with a manually curated training set of 180 orthologous gene pairs from P. knowlesi and P. falciparum.
Functional assignments were based on assessment of BLAST and FASTA similarity searches against public databases and searches in protein domain databases such as InterPro34. In addition, TMHMMv2.035, SignalPv3.036 and t-RNA scan37 were used to identify transmembrane domains, signal peptides and t-RNA genes.
To define the orthologous and paralogous relationships between the predicted proteomes of three Plasmodium species (P. falciparum, P. knowlesi, P. vivax), the OrthoMCL protein clustering algorithm38 was used with an inflation value of 1.5.
To search for parasite proteins containing stretches of perfectly matched host sequences, the PMATCH algorithm (R. Durbin, unpublished) was used to report exact matches of 15 amino acids or greater after screening out low complexity sequences (details are provided in Supplementary Information).
Sequence alignments and dotter39 analysis of SICA proteins revealed the presence of a distinct N-terminal cysteine-rich domain (termed SICA-α: in some cases there are two copies of this domain), multiple central cysteine-rich domains (SICA-β) and a C-terminal cytoplasmic encoding domain (SICA_C). For each domain, a profile HMM (using HMMer, http://hmmer.janelia.org/) was constructed and searched against the P. knowlesi genome to find all examples of the domain (significant matches had E-values <0.001). The HMMs were rebuilt, using alignments constructed using all significant hits, and re-searched until no additional examples of the domain were found.
The program Phobius40 was used to identify the putative transmembrane region located between the end of the last SICA-β domain and the SICA_C domain in all cases. An identical procedure was used to identify the domains in the KIR proteins. In this case, a single domain type was found on all KIR proteins, repeated between one and three times. Putative transmembrane proteins were identified as before, but only ~50% of KIR proteins had a predicted transmembrane region. Visual inspection of the corresponding C-terminal regions from sequences, both with and without predictions, showed the presence of a common hydrophobic patch. To investigate whether the Phobius40 software was insufficiently sensitive to identify all of the KIR transmembrane regions, the predicted transmembrane regions were aligned and used to build a HMM of the transmembrane region. This was then used to iteratively search the whole genome as before.
Full Methods and any associated references are available in the online version of the paper at www.nature.com/nature.