Human malaria is transmitted by anopheline mosquitoes and caused by four species in the genus Plasmodium. Of these, Plasmodium vivax is the major malaria agent outside Africa, annually causing 80 to 100 million cases¹. Often mistakenly regarded as benign and self-limiting, P. vivax treatment and control present challenges distinct from those of the more virulent P. falciparum. Biological traits including a dormant (hypnozoite) liver stage responsible for recurrent infections (relapses), early infective sexual stages (gametocytes), and transmission from low parasite densities in the blood², coupled with emerging antimalarial drug resistance³, render P. vivax resilient to modern control strategies. Recent evidence indicates that P. falciparum derives from parasites of great apes in Africa⁴, while P. vivax is more closely related to parasites of Asian Old World monkeys (OWM)^5–7, though not itself infective to OWM.

P. vivax cannot be cultured in vitro, and the small New World monkeys capable of hosting it are rare and do not provide an ideal model system. P. knowlesi, an Asian OWM parasite recently recognized as a significant zoonosis for humans⁸, offers a sequenced genome⁹, but the species is distantly related to P. vivax and phenotypically dissimilar. In contrast, P. cynomolgi, a simian parasite which can infect humans experimentally¹⁰, is the closest living relative (i.e., a sister taxon) to P. vivax and possesses most of its genetic, phenotypic, and biologic characteristics: importantly, periodic relapses caused by dormant hypnozoites, early infectious gametocyte formation, and invasion of Duffy blood group-positive reticulocytes. P. cynomolgi thus offers a robust model for P. vivax in a readily-available laboratory host, the Rhesus monkey, whose genome was recently sequenced¹¹. Here we report draft genome sequences of three P. cynomolgi strains, and comparative genomic analyses of P. cynomolgi, P. vivax¹² and P. knowlesi⁹, three members of the monkey malaria clade.

We sequenced the genome of P. cynomolgi ‘B’ strain, isolated from a monkey in Malaysia and grown in splenectomized monkeys (Online Methods). A combination of Sanger, Roche 454 and Illumina chemistries was deployed to generate a high quality reference assembly at 161-fold coverage, consisting of 14 supercontigs (corresponding to the 14 parasite chromosomes) and ~1,649 unassigned contigs, comprising a total length of ~26.2 Mb (Supplementary Table 1). Comparing genomic features of P. cynomolgi, P. knowlesi and P. vivax reveals many similarities, including GC content (mean GC 40.5%), fourteen conserved centromeres, and presence of intrachromosomal telomeric sequences (ITSs; GGGTT[T/C]A) discovered in the P. knowlesi genome⁹ but absent in P. vivax (Table 1, Supplementary Table 2 and Fig. 1).

We annotated the P. cynomolgi strain B genome using a combination of ab initio gene prediction programs trained on high quality datasets and sequence similarity searches against the annotated P. vivax and P. knowlesi genomes. Unsurprisingly for species from the same monkey malaria clade, gene synteny along the 14 chromosomes is highly conserved, though numerous microsyntenic breaks occur in regions containing multigene families (Fig. 2 and Table 2). The P. cynomolgi genome contains 5,722 genes, of which about half encode conserved hypothetical proteins of unknown function, as is the case in all Plasmodium genomes sequenced to date. A maximum likelihood (ML) phylogenetic tree using 192 conserved ribosomal and translation/transcription-related genes (Supplementary Figure 1) confirms the close relationship of P. cynomolgi to P. vivax compared to five other Plasmodium species. Approximately 90% of genes (4,613) have reciprocal best match orthologs in all three species (Fig. 3), enabling refinement of the existing P. vivax and P. knowlesi annotations (Supplementary Table 3). The high degree of gene orthology enabled us to identify specific examples of gene duplication (an important vehicle for genome evolution), including a P. cynomolgi duplicated homolog of P. vivax Pvs28, a sexual stage surface antigen that is a transmission-blocking vaccine candidate¹³ (Supplementary Table 4). Genes common only to P. cynomolgi and P. vivax (n = 214) outnumber those exclusive to P. cynomolgi and P. knowlesi (n = 100) or P. vivax and P. knowlesi (n = 17). Such figures establish the utility of P. cynomolgi as a model species for studying the more intractable P. vivax.

Remarkably, most of the genes specific to a particular species belong to multigene families (excluding hypothetical genes; Table 2 and Supplementary Table 5). This suggests repeated lineage-specific gene duplication and/or gene deletion in multigene families within the three monkey malaria clade species. Moreover, copy number of multigene families was generally greater in the P. cynomolgi/P. vivax lineage than in P. knowlesi, suggesting repeated gene duplication in the ancestral lineage of P. cynomolgi/P. vivax (or repeated gene deletion in the P. knowlesi lineage). Thus the genomes of P. cynomolgi, P. vivax and P. knowlesi can almost be completely characterized by variations in the copy number of multigene families. Examples of such families include those that encode proteins involved in evasion of the human immune system (vir, kir and SICAvar), and invasion of host red blood cells (dbp and rbp), described below.

In malaria parasites, invasion of host erythrocytes, mediated by specific interactions between parasite ligands and erythrocyte receptors, is a crucial component of the parasite life-cycle. Of great interest are the ebl and rbl families, which encode parasite ligands required for the recognition of host erythrocytes. The ebl genes encode erythrocyte binding-like (EBL) ligands such as the Duffy binding proteins (DBP) that bind to Duffy Antigen/Receptor for Chemokines (DARC) on human and monkey erythrocytes. The rbl genes encode the reticulocyte binding-like (RBL) protein family, including reticulocyte binding proteins (RBPs) in P. cynomolgi and P. vivax, and normocyte binding proteins (NBPs) in P. knowlesi, which bind to unknown erythrocyte receptors¹⁴. We confirmed the presence of two dbp genes in P. cynomolgi¹⁵ (Supplementary Table 6), in contrast to the one dbp and three dbp genes identified in P. vivax and P. knowlesi, respectively. This raises an intriguing hypothesis -- that P. vivax lost one dbp gene, and hence its infectivity to OWM erythrocytes, after divergence from a common P. vivax/P. cynomolgi ancestor. This hypothesis is also supported by our identification of only single copies of dbp genes in two other closely related OWM malaria parasites P. fieldi and P. simiovale, which are incapable of infecting humans¹⁶. These latter OWM species lost one or more dbp genes during divergence that conferred infectivity to humans, while P. cynomolgi and P. knowlesi retained dbp genes that allow invasion of human erythrocytes (Supplementary Figure 2).

We found multiple rbp genes, some truncated or present as pseudogenes, in the P. cynomolgi genome (Fig. 1 and Table 2). Phylogenetic analysis shows that rbls from P. cynomolgi, P. vivax and P. knowlesi can be classified into three distinct groups, RBP/NBP-1, RBP/NBP-2 and RBP/NBP-3 (Supplementary Figure 3), and suggests that these groups existed before the three species diverged. All three groups of RBP/NBP are represented in P. cynomolgi, whereas P. vivax and P. knowlesi lack functional genes of the RBP/NBP-3 and RBP/NBP-1 groups, respectively. Thus rbl gene family expansion appears to have occurred after speciation, indicating that the three species have multiple species-specific erythrocyte invasion mechanisms. Intriguingly, we found an ortholog of P. vivax rbp1b in some strains of P. cynomolgi, but not in others (Supplementary Table 6). This is the first example of an rbp copy number variant between strains of a single Plasmodium species that we are aware of, and highlights how repeated birth and death of rbl genes, a signature of adaptive evolution, may have enabled species of the monkey malaria clade to expand or switch hosts between monkeys and humans.

The largest gene family in P. cynomolgi, consisting of 256 cyir (cynomolgi interspersed repeat) genes, is part of the pir (plasmodium interspersed repeat) superfamily that include P. vivax virs (n = 319) and P. knowlesi kirs (n = 70; Table 2). Pir proteins reside on the surface of infected erythrocytes and play an important role in immune evasion¹⁷. Most cyir genes have sequence similarity to P. vivax vir genes (n = 254; Supplementary Table 7) and are found in subtelomeric regions (Fig. 1), but interestingly 11 cyirs have sequence similarity to P. knowlesi kirs (Supplementary Table 7) and occur internally as the kirs do in P. knowlesi. As with the ‘molecular mimicry’ in P. knowlesi⁹, one CYIR protein (PCYB_032250) has a 56-amino acid region highly similar to the extracellular domain of primate CD99 (Supplementary Figure 4), a molecule involved in the regulation of T-cell function. A novel finding is that P. cynomolgi has two genes whose sequences are similar to P. knowlesi SICAvar genes (Supplementary Table 7) that are expressed on the surfaces of schizont-infected macaque erythrocytes and are involved in antigenic variation¹⁸.

The ability to form a dormant hypnozoite stage is common to both P. cynomolgi and P. vivax, and was first shown in laboratory infections of monkeys by mosquito-transmitted P. cynomolgi¹⁹. In a search for candidate hypnozoite genes, we identified nine coding for ‘dormancy-related’ proteins and possessing ApiAP2 motifs²⁰ necessary for stage-specific transcriptional regulation at the sporozoite (pre-hypnozoite) stage (Supplementary Table 8). The candidates include kinases that are involved in cell-cycle transition; hypnozoite formation may be regulated by phosphorylation of proteins specifically expressed at the pre-hypnozoite stage. Our list of P. cynomolgi candidate genes represents an informed starting point for experimental studies on this elusive stage.

We sequenced P. cynomolgi strains Berok (from Malaysia) and Cambodian (from Cambodia) to 26× and 17 coverage, respectively, to characterize P. cynomolgi genome-wide diversity through analysis of SNPs, CNVs and microsatellites. A comparison of the three P. cynomolgi strains identified 178,732 SNPs (Supplementary Table 9), a frequency of one SNP per 151 bp, a polymorphism level somewhat similar to that found when P. falciparum genomes are compared^21,22. We calculated the pairwise nucleotide diversity (π) as 5.41 × 10⁻³ across the genome, which varies little between the chromosomes. We assessed genome-wide copy number variants (CNVs) between P. cynomolgi B and Berok strains, using a robust statistical model in the program CNV-seq²³, identifying 1,570 CNVs (1 per 17 kb), including one containing the rbp1b gene on chromosome 7 described above (Supplementary Figure 5). Finally, mining of P. cynomolgi B and Berok strains identified 182 polymorphic intergenic microsatellites (Supplementary Table 10), the first set of genetic markers developed for this species. These provide a toolkit for genetic diversity and population structure studies of laboratory stocks or natural infections of P. cynomolgi, many of which have recently been isolated from screening hundreds of wild monkeys for the zoonosis P. knowlesi²⁴.

We estimated the difference between the number of synonymous changes per synonymous site (Ds) and the number of non-synonymous changes per non-synonymous site (Dn) over 4,605 pairs of orthologs within P. cynomolgi strain B and Berok, and between these two P. cynomolgi genomes and P. vivax Salvador I, using a simple Nei-Gojobori model²⁵. We found 196 genes with Dn > Ds within the two P. cynomolgi strains, and 1,742 genes with Ds > Dn (Supplementary Table 11). Genes with relatively high Dn > Ds include transmembrane proteins such as antigens and transporters, among which is a transmission blocking target antigen Pcyn230 (PCYB_042090). Interestingly, the P. vivax ortholog (PVX_003905) does not show evidence for positive selection²⁶, suggesting species-specific positive selection. We explored the degree to which evolution of orthologs has been constrained between P. cynomolgi and P. vivax, and found 154 genes under possible accelerated evolution but 1,613 genes under possible purifying selection (Supplementary Table 12). This conservative estimate indicates that at least 35% of loci have diverged under strong constraint, compared with 3.3% of loci under less constraint or positive selection (Figure 1), indicating that overall the genome of P. cynomolgi is highly conserved in single locus genes when compared to P. vivax, and emphasizing the value of P. cynomolgi as a biomedical and evolutionary model for studying P. vivax.

Our generation of the first P. cynomolgi genome sequences is a critical component in the development of a robust model system for the intractable and neglected P. vivax species²⁷. Comparative genome analysis of P. vivax and the OWM malaria parasites P. cynomolgi and P. knowlesi presented here provides the foundation for further insights into traits such as host specificity that will enhance the prospects for the eventual elimination of vivax malaria and global malaria eradication.