Although Plasmodium vivax is responsible for the majority of malaria infections outside Africa, little is known about its evolution and pathway to humans. Its closest genetic relative, P. vivax-like, was discovered in African great apes and is hypothesized to have given rise to P. vivax in humans. To unravel the evolutionary history and adaptation of P. vivax to different host environments, we generated using long- and short-read sequence technologies 2 new P. vivax-like reference genomes and 9 additional P. vivax-like genotypes. Analyses show that the genomes of P. vivax and P. vivax-like are highly similar and colinear within the core regions. Phylogenetic analyses clearly show that P. vivax-like parasites form a genetically distinct clade from P. vivax. Concerning the relative divergence dating, we show that the evolution of P. vivax in humans did not occur at the same time as the other agents of human malaria, thus suggesting that the transfer of Plasmodium parasites to humans happened several times independently over the history of the Homo genus. We further identify several key genes that exhibit signatures of positive selection exclusively in the human P. vivax parasites. Two of these genes have been identified to also be under positive selection in the other main human malaria agent, P. falciparum, thus suggesting their key role in the evolution of the ability of these parasites to infect humans or their anthropophilic vectors. Finally, we demonstrate that some gene families important for red blood cell (RBC) invasion (a key step of the life cycle of these parasites) have undergone lineage-specific evolution in the human parasite (e. g. , reticulocyte-binding proteins [RBPs]).
Eleven P. vivax-like genotypes were obtained from 2 different kinds of samples: 10 infected chimpanzee blood samples collected during successive routine sanitary controls of chimpanzees living in the Park of La Lékédi (a sanctuary in Gabon) and 1 infected Anopheles mosquito (An. moucheti) collected during an entomological survey carried out in the same park (S1 Table) [11]. For blood samples, white blood cells were depleted using the CF11 method [12] to reduce the amount of host DNA. After DNA extraction, samples were subjected to whole-genome amplification (WGA) in order to obtain sufficient parasite DNA for library preparation. Sequencing was then performed using short-read Illumina technology. For one sample (Pvl06), long-read sequencing (Pacific Biosciences [PacBio] technology) was performed in order to get a better coverage of regions containing subtelomeric gene families. Among the 11 samples, 10 presented mixed infections with other Plasmodium species (S1 and S2 Tables). Four samples containing P. gaboni or P. malariae-like co-infections were used in other studies (see S1 Table) [13,14]. In order to obtain the P. vivax-like genotypes and to preclude errors due to co-infections with other Plasmodium species, sequencing reads were extracted based on their similarity to the reference genome sequence of P. vivax, PvP01 [15], and any reads mapping against a reference genome of a Laverania chimpanzee-infecting species (e. g. , PrG01, PbilcG01, and PGAB01; see Otto and colleagues [13]) were removed (S2 Table). Concerning multiple infections with several P. vivax-like strains, only 2 (Pvl09 and Pvl10) seemed to be multiply infected as suggested by the reference allele frequency (RAF) distributions (see Materials and methods section and S1 Fig). In order to avoid any bias in the genomic analysis due to multiple infections with several strains of P. vivax-like, only 1 variant was extracted per sample. This was done by considering the allele with the highest frequency in the calling and filtering analysis of Single Nucleotide Variants (SNVs) (see Materials and methods section). Sequencing reads from 2 samples, one obtained using Illumina sequencing, Pvl01, and another using PacBio technology, Pvl06, were used to perform de novo genome assemblies and were annotated to produce reference genomes for P. vivax-like (S1 Table). Of the 2 assemblies, Pvl01 is of considerably higher quality (4,570 one-to-one orthologues to the PvP01 reference genome compared with 2,925 for Pvl06 [Table 1]). Both assemblies consist of 14 supercontigs (corresponding to the 14 P. vivax chromosomes) —and 1,176 and 351 unassigned contigs—comprising a total of 27. 5 Mb and 18. 8 Mb in size, respectively, for Pvl01 and Pvl06, respectively. After annotation with Companion [16], these 2 genomes contained 5,532 and 4,953 annotated genes (Table 1). We obtained for the 9 remaining samples between 2. 9 and 86 million of reads that paired with the P. vivax PvP01 reference genome (mean ± SD = 2. 17 ± 25. 71), with a mean depth per high-quality position ranging from 13. 98 to 335. 1 (mean ± SD = 93. 43 ± 99. 51; see S3 Table). These genome sequences were used for SNV calling for population genetic and phylogenetic analyses.
Comparing the P. vivax-like reference genomes to those of P. vivax (PvP01 and SalI) [3,15], P. cynomolgi (B and M strains) [8,17], and P. knowlesi (H strain) [18] reveals several similarities, including a similar guanine–cytosine (GC) content and extensive collinearity and conservation of gene content and organization (Table 1). The P. vivax-like core genome sequences are completely syntenic to the P. vivax PvP01 reference genome sequence (S2 and S3 Figs). Because multigene families are known to evolve extremely rapidly in their genome structure, obtaining the full genomes of species closer to human P. vivax is fundamental for a better understanding of its evolution, adaptation, and emergence in different host species. For Plasmodium parasites, most species-specific genes are part of large gene families, such as var genes in P. falciparum or pir genes that are present in all Plasmodium genomes studied [19,20]. Table 2 provides a summarized view of gene content and copy number of the main multigene families in P. vivax-like in comparison with P. vivax, P. knowlesi, and P. cynomolgi. Even if certain subtelomeric regions of our reference genomes (S2 and S3 Figs) are not complete, at least one copy of each major gene family was detected (Table 2). In comparison with P. vivax (and expected because of the partial subtelomeric sequencing coverage), the number of copies in each family was generally lower or equal in P. vivax-like. For these families, all genes were functional except for the Cytoadherence-linked asexual gene (clag) family. For the clag family, all genes are functional except the one situated on chromosome 8 for P. vivax-like (confirmed for both Pvl01 and Pvl06) (S4 Fig). The clag family, strictly conserved in malaria parasites, is an essential gene family in host–parasite interactions, playing a role in merozoite invasion, parasitophorous vacuole formation, and in the uptake of ions and nutrients from the host plasma [21,22]. The pseudogenization of the clag gene on chromosome 8 for P. vivax-like suggests that this species lost the clag gene during its adaptation to the ape host. During the life cycle of Plasmodium parasites, host RBC invasion is mediated by specific interactions between parasite ligands and host erythrocyte receptors. Two major multigene families are involved in RBC invasion: the Duffy-binding proteins (dbp) and the reticulocyte-binding proteins (rbp) multigene families [23]. DBP is a protein secreted by the micronemes of the merozoite stage that binds to the Duffy Antigen Receptor for Chemokine (DARC) to invade RBCs. P. vivax is characterized in its genome by 2 dbp genes (dbp1 on chromosome 6 and dbp2 on chromosome 1) that seem to be essential for RBC invasion, as demonstrated by their inability to infect individuals not expressing the Duffy receptor on the surface of their RBCs (i. e. , Duffy-negative individuals) [24–26]. In the reference genomes of P. vivax-like (Pvl01 and Pvl06), we observe the dbp1 gene as in the P. vivax genome PvP01 (Table 2) and also in the other Plasmodium species (P. knowlesi and P. cynomolgi); however, we did not observe the dbp2 gene (no read obtained mapping to this region). This observation was confirmed in the other genotypes sequenced in this study. Knowing that gorillas and chimpanzees are all described as Duffy positive [10], we propose that P. vivax-like parasites infect only Duffy-positive hosts, which could be associated with the absence of the dbp2 gene. This would be in accordance with the fact that the only described transfer of P. vivax-like to humans was in a Caucasian Duffy-positive individual [9] and that no transfers of P. vivax-like have been recorded in Central African Duffy-negative populations despite the fact that they live in close proximity to infected ape populations [27]. rbp genes encode a merozoite surface protein family present across all Plasmodium species and known to be involved in RBC invasion and host specificity [23]. Among rbp family genes, 3 gene classes (rbp1, rbp2, and rbp3) exist and are associated with the ability of Plasmodium parasites to invade different maturation stages of RBCs. In this study, comparison of the organization and characteristics of the rbp gene family between P. vivax, P. vivax-like, P. knowlesi, and P. cynomolgi (Fig 1 and Table 2) first reveals that gene classes RBP2 and RBP3 are ancestral to the divergence of all these species except P. knowlesi. Second, an expansion of the rbp2 class is observed in the P. vivax/P. vivax-like/P. cynomolgi lineage (Fig 1A), suggesting that, in this lineage, specific expansion likely occurred during the evolution of these species. Finally, rbp3 genes, which are supposed to confer the ability to infect normocytes, are functional in all species except in P. vivax (for which the gene is pseudogenized in both SalI and PvP01 strains), suggesting that P. vivax lost the ability to infect normocytes or has developed an ability to infect specifically only reticulocytes during its adaptation to human RBCs (Fig 1A and 1B).
Conservation of the gene content between P. vivax-like with the other primate-infective Plasmodium species has enabled us to reconstruct with confidence the relationships between the different species and to estimate the relative age of the different speciation events. This analysis confirmed the position of P. vivax-like as the closest sister lineage of P. vivax (Fig 2). Regarding the estimation of divergence times using genomic information, different methods were recently used for Plasmodium, such as the one implemented in Generalized Phylogenetic Coalescent Sampler (G-PhoCS) [28] or the one developed by Silva and colleagues [29]. G-PhoCS uses a Bayesian Markov Chain Monte Carlo (MCMC) approach to infer, based on the information provided by multiple loci, the divergence time between species. This method has been applied in 2 recent studies for Plasmodium parasites—one aiming at estimating the relative split times between the 2 P. ovale subspecies and between P. malariae and P. malariae-like [14], the other to estimate the divergence time within the Laverania subgenus, a subgenus including P. falciparum and all its closest ape relatives [13]. The Silva method is based on the estimate of the sequence divergence in different proteins and comparison of this divergence measured between different lineages [29]. In this method, the regression slope of the divergence between the proteins in 2 lineages reflects their relative age. The advantage is that it does not rely on an estimate of mutation rate. Finally, it has already been used in a recent study estimating avian and primate Plasmodium species divergence times [30]. Here, without calibration points and a good estimation of the P. vivax and other Plasmodium species substitution rates, our aim was to evaluate the divergence time between P. vivax and P. vivax-like relative to the divergence time between the other primate–nonprimate Plasmodium species pairs. To evaluate the influence of the approach on the estimations, we used 2 strategies to estimate divergence time between P. vivax and P. vivax-like relative to the other divergence events within the tree: the Silva method [29] and the RelTime method (see Materials and methods section for a description of the two methods) [31,32]. We first used the Silva method [29], checking for the influence of the model of evolution and the reference species pair on the analysis (see Materials and methods section). Because neither influenced the results of the relative ages of species pairs (see Fig 2, S5 Fig and S4 Table), we only report the results for the Silva method using the JTT (Jones, Taylor, and Thornton) model of evolution and the P. ovale curtisi–P. ovale wallikeri pair as the species pair of reference as in Böhme and colleagues in Fig 2 [30]. While the 2 methods are mostly in agreement (with RelTime showing larger confidence intervals [CIs]), they show a discrepancy concerning the relative divergence time of the P. malariae and P. malariae-like parasites (see Fig 2 and S4 Table). However, both methods suggest that the evolution of P. vivax in humans did not occur at the same time as the other human malaria agents of the Plasmodium genus (i. e. , P. malariae and P. falciparum). The analyses show that the time of the split between P. vivax and P. vivax-like happened before the divergence between P. falciparum and P. praefalciparum (about 2. 3 times earlier). It is unclear from our analyses whether the divergence between the 2 P. malariae species occurred before or after the divergence of the other 2 human–nonhuman parasite pairs. Unlike us, a previous study estimated, using the G-PhoCS method, the relative dating of the split between P. reichenowi (a chimpanzee parasite) and P. falciparum (P. praefalciparum was not available at that time) and that of the split between P. malariae and P. malariae-like (note that similar estimates were recently obtained using another set of data using the Silva method) [14] to have occurred at the same time. This discrepancy between methods suggests that the same strict molecular clock (which is a hypothesis of the Silva method) may not apply over the entire tree (especially for the Laverania subgenus because of their extremely low GC content in comparison with other Plasmodium species). Whatever the method used, all these estimates nevertheless suggest that the evolution of P. vivax in humans did not occur at the same time as the other human malaria species and that the transfer of Plasmodium parasites to humans may have happened several times independently over the history of the Homo genus [33–37].
To analyze the relationship between our 11 P. vivax-like isolates and human P. vivax, we completed our dataset with 19 published human P. vivax genomes (S1 Table) [38]. All sequencing reads were aligned against the PvP01 reference genome [15], and SNVs were called and filtered as described in the Materials and methods section. Maximum likelihood phylogenetic trees were then produced based on 100,616 SNVs. Our results clearly demonstrate the presence of 2 significantly distinct genetic clades (with a bootstrap value of 100) composed of P. vivax-like strains on one side and human P. vivax isolates on the other side (Fig 3). This result differs from previous results suggesting that human strains formed a monophyletic clade within the radiation of ape P. vivax-like parasites [10]. One explanation for this difference with previous published results could be that it is due to a phenomenon called Incomplete Lineage Sorting (ILS) or to a lack of phylogenetic signal for phylogenies performed on a single or few genes. ILS is the discordance observed between some gene trees and the species or population tree due to the coalescence of gene copies in an ancestral species or population [39]. Such a phenomenon is often observed when species or population divergence is recent, which is the case for P. vivax/P. vivax-like [40,41]. ILS may thus result in the wrong conclusion of P. vivax and P. vivax-like populations being intermixed and P. vivax diversity being included in the diversity of P. vivax-like. A lack of phylogenetic signal, which occurs frequently when species diverged recently, would have similar consequences. In order to test the implication of ILS or lack of phylogenetic signal, we generated a phylogenetic tree and a reticulated network on partial mitochondrial genomes of P. vivax and P. vivax-like obtained in Liu and colleagues [10], Prugnolle and colleagues [9], and in the current study. These analyses show that partial mitochondrial genetic information is not enough to make a final conclusion on the origin of P. vivax parasites (S6 Fig). When considering the polymorphism level of this portion of mitochondrial genomes, over the 2,483 bp, only 127 positions showed variability, 20 of them being a position specific to the outgroup P. cynomolgi. Among the 107 variable positions found in the P. vivax and P. vivax-like samples, 87 were singletons, meaning that only 20 SNVs were shared by more than 2 individuals. This suggests that the discrepancy between our phylogeny and the one of Liu and colleagues [10] is probably because of a lack of phylogenetic signal. Indeed, in our study, the use of significantly more genetic information from throughout the genome, both in genic and intergenic regions, provides a more accurate picture of the genetic relationships between the different parasite species. Reducing our genetic data to single genes (as performed in previous studies) or a limited number of SNVs also generates phylogenies in which P. vivax is included within the diversity of P. vivax-like (see S7 Fig). Another explanation for the discrepancy between our results and the results from Liu and colleagues [10] could be that we are missing part of the diversity of P. vivax-like given that we obtained the genomes from a limited number of parasites isolated from a small population of apes in Gabon. Nevertheless, this hypothesis does not seem to hold in light of the tree and network produced using all the available mitochondrial sequences of P. vivax-like and our data, as isolates are distributed all over the currently known genetic diversity of P. vivax-like (S6 Fig). Our results show that P. vivax-like is composed of 2 distinct lineages: one including the 2 reference genomes (Pvl01 and Pvl06) and 7 other isolates that will hereafter be referred to as P. vivax-like 1, and another one including 2 isolates (Pvl09 and Pvl10) that will hereafter be referred to as P. vivax-like 2 (Fig 3). This sub-structuration of P. vivax-like is confirmed in the mitochondrial-reticulated network obtained using a larger number of isolates (see S6 Fig). These 2 lineages may thus reflect an ancient split within P. vivax-like or be the consequence of a recent introgression or hybridization event between P. vivax-like and P. vivax in Africa. A search of recent recombination events between lineages using SplitsTree (http: //www. splitstree. org/) [77] does not support this latter hypothesis (S6 Fig). Previous studies highlighted the high genetic diversity of P. vivax-like populations in comparison with P. vivax worldwide [9,10]. In this genome-wide analysis of nucleotide diversity π, we confirm that P. vivax-like populations are significantly more diverse than P. vivax populations (P < 0. 001; Wilcoxon test), with P. vivax-like samples showing nearly 10 times higher nucleotide diversity (πP. vivax = 0. 0012; πP. vivax-like = 0. 0096). Such genetic diversity in P. vivax-like strains in comparison with human P. vivax has already been described in other studies [9,10], suggesting that P. vivax-like strains probably display extremely high genetic diversity in Central African regions. This suggests that P. vivax-like parasites of African great apes are probably more ancient than the human P. vivax strains and that the human P. vivax species (as for human Plasmodium species like P. falciparum) went through a bottleneck and only recently underwent population expansion.
Before the discovery of the ape P. vivax-like, the main hypothetical scenario concerning the origin of P. vivax was that of an “Out of Asia. ” More specifically, it was considered that P. vivax emerged in humans following its transfer to humans from Asian monkeys, as has been recently described for P. knowlesi [42]. However, the recent discovery of P. vivax-like in African great apes [10] and the analysis of their genetic characteristics led researchers to propose an “Out of Africa” origin of P. vivax [10,36]. Based on phylogenetic analyses of partial mitochondrial genomes and nuclear sequences of P. vivax and P. vivax-like parasites isolated from great apes and humans, Liu and colleagues [10] suggested that all extant human P. vivax parasites derived from one single ancestor that was transferred from great apes to humans. Our results do not bring any new evidence in favor of one or the other scenario. We think, nevertheless, that the origin of P. vivax is more complex than recently proposed (a single host switch from apes to humans in Africa) and that some previous data—especially those regarding the paraphyly of P. vivax-like and the inclusion of P. vivax within the diversity of P. vivax-like (the arguments used in support of an African scenario) —have been overinterpreted. The inclusion of P. vivax into the P. vivax-like diversity based on only a couple of nuclear genes and partial mitochondrial genomes is not a definitive proof of transfer of P. vivax-like to humans in Africa because alternative explanations can be provided. Indeed, as discussed above, such a phylogenetic pattern can be obtained because of ILS or because of a lack of phylogenetic signal for the sequences used, 2 phenomena that are frequent when species diverged recently (which is the case for P. vivax and P. vivax-like). Our study indeed shows that incongruent phylogenies may be obtained when one limits the analyses to a single or a couple of genes (see section above and S6 and S7 Figs). In our mind, it is still impossible to conclude that P. vivax has an African origin and that it was transferred from African apes to humans, especially when several observations are more in favor of an Asian origin, such as the following: (i) P. vivax evolved in a clade of parasites infecting Asian monkeys, and (ii) the highest genetic diversity of P. vivax is observed in Asian populations, and its diversity decreases toward Africa. This pattern of diversity is the opposite for P. falciparum, which has a well-established African origin [43]. For this parasite, the highest genetic diversity is found in Africa and decreases toward Asia accompanying the human migration [43]. Because the data that are currently available are not easy to interpret and are somehow contradictory, we think that more complex scenarios regarding the origin of P. vivax should be envisaged in light of phylogenetic and population genetic evidences (as proposed in Prugnolle and colleagues [9]) and more genomic data need to be obtained in support of these scenarios.
Comparison of the P. vivax genome to its closest sister lineage (P. vivax-like) and to the other primate Plasmodium parasites provides a unique opportunity to identify P. vivax–specific adaptations to humans. We applied a branch-site test of positive selection to detect events of positive selection that exclusively occurred in the P. vivax lineage. Within the reference genome P. vivax-like (Pvl01), 418 genes exhibited significant signals of positive selection (S5 Table). In the human P. vivax genome PvP01, the test allowed the identification of 255 genes showing significant signals of positive selection (S6 Table). Among these genes presenting a significant dN/dS ratio, 71 were shared between P. vivax and P. vivax-like (genes indicated in orange in the S6 Table), including 56 encoding for proteins with unknown function and 15 encoding for proteins that are involved either in energy metabolism regulation (n = 9), chromatid segregation (n = 2), or cellular-based movement (n = 4). We then took into consideration those 255 genes detected under positive selection in P. vivax and compared them to those obtained in P. falciparum (172 genes under selection; see Supplementary Table 4 in [13]). We identified a subset of 10 genes under positive selection in both the human P. vivax and P. falciparum parasites (P < 0. 05) (S7 Table). Among these 10 genes, 5 code for conserved Plasmodium proteins with unknown function and 3 for proteins involved in either transcription or transduction. Interestingly, the 2 remaining genes under positive selection in these 2 human Plasmodium parasites code for the oocyst capsule protein, which is essential for malaria parasite survival in the Anopheles’ midgut, and for the rhoptry protein ROP14, involved in protein maturation and the host cell invasion. These results suggest that these proteins could be essential for infection of humans or their vectors, and future studies should focus on the involvement of these proteins in human parasite transmission and infection.
Through technical accomplishments, we produced and assembled the first P. vivax-like reference genomes, the closest sister clade to human P. vivax—an indispensable step for a better understanding of this enigmatic species. We established that P. vivax-like parasites form a genetically distinct clade from P. vivax. Concerning the relative divergence dating, we estimated that the divergence between both species occurred probably before the split between P. malariae species. This suggests that the transfer of Plasmodium parasites to humans happened several times independently over the history of the Homo genus. Our genome-wide analyses provided new insights into the adaptive evolution of P. vivax. Indeed, we identified several key genes that exhibit signatures of positive selection exclusively in the human P. vivax parasites and show that some gene families important for RBC invasion have undergone species-specific evolution in the human parasite, e. g. , rbp and dbp. Are these genes the keys to the emergence of P. vivax in the human populations? This pending question will need to be answered through functional studies associated with deeper whole-genome analyses. Among the genes identified to be under positive selection, 2 were also identified to be under positive selection in the other main human malaria agent, P. falciparum, thus suggesting their key role in the evolution of the parasites in their ability to infect humans or their anthropophilic vectors. To conclude, this study provides the foundation for further investigation into P. vivax traits that are of public health importance, such as features involved in host–parasite interactions, host specificity, and species-specific adaptations.
P. vivax-like samples were identified by molecular diagnostic testing during a continuous survey of great ape Plasmodium infections carried out in the Park of La Lékédi, in Gabon, by the Centre International de Recherches Médicales de Franceville (CIRMF) [9]. In parallel, a survey of Anopheles mosquitoes circulating in the same area was conducted in order to identify potential vectors of ape Plasmodium [11]. Specifically, mosquitoes were trapped with CDC light traps in the forest of the Park of La Lékédi. Anopheles specimens were retrieved and identified using a taxonomic key [44] before proceeding to dissection to isolate the abdomen. Samples were then stored at −20 °C until transportation to the CIRMF, Gabon, where they were stored at −80 °C until processed. Blood samples of great apes were treated using leukocyte depletion by CF11 cellulose column filtration [45]. P. vivax-like samples were identified either by amplifying and sequencing the Plasmodium Cytochrome b (Cytb) gene as described in Ollomo and colleagues [46] or directly from samples already used for studies of other Plasmodium species [13,14]. This allowed the detection of 11 P. vivax-like samples, 10 from chimpanzees and 1 from an An. moucheti mosquito. Most of these samples were co-infected with other Plasmodium species and/or probably with multiple P. vivax-like isolates (see below and S1 Table). The identification of intraspecific P. vivax-like co-infections was made by analyzing the distribution of the RAF [47]. The animals’ well-being was guaranteed by the veterinarians of the “Parc of La Lékédi” and the CIRMF, who were responsible for the proceeding sanitary procedures (including blood collection). All animal work was indeed conducted according to relevant national and international guidelines. These investigations were approved by the Government of the Republic of Gabon and by the Animal Life Administration of Libreville, Gabon (no. CITES 00956). It should be noted that our study did not involve randomization or blinding. To overcome host DNA contamination, within 6 hours after blood collection in the La Lékédi park, each 5 mL chimpanzee blood sample was diluted with 1 volume of PBS and passed through CF11 cellulose powder columns to remove host leukocytes (i. e. , leukocyte depletion) [45]. DNA was then extracted for the 11 samples using Qiagen Midi extraction kits (Qiagen) following the manufacturer’s recommendations. Next, DNA samples were enriched using a WGA using REPLI-g Mini Kit (Qiagen) following the strategy described in Oyola and colleagues [48] to optimize WGA of AT-rich genomes from low DNA quantities. The Illumina isolates were sequenced using Illumina Standard libraries of 200- to 300-bp fragments, and amplification-free libraries of 400- to 600-bp fragments were prepared and sequenced on the Illumina HiSeq 2500 and the MiSeq version 2 according to the manufacturer’s standard protocol (S1 Table). Because of its high DNA content, the Pvl06 isolate was sequenced using PacBio chemistry. After a greater than 8-kb–size selection of DNA fragments using the BluePippin Size-Selection System (Sage Science), the library for this Pvl06 sample was sequenced using P6 polymerase and PacBio chemistry version C3/P5 in 20 SMRT cells. Raw sequence data are deposited in the European Nucleotide Archive. The accession numbers can be found in S1 Table. For the P. vivax-like Pvl01 and Pvl06, P. vivax PvP01 and SalI, P. cynomolgi B strain, and P. knowlesi H strain genomes, gene variants were detected and counted using Geneious software [60]. This allowed the presence/absence of the variants of the different gene families across Plasmodium species to be evaluated. For each gene family, the number of variants identified in the 2 reference genomes Pvl01 and Pvl06 was confirmed by manually checking the presence/absence of these in the other P. vivax-like genotypes obtained using ACT and bamview [61,62]. Orthologous groups across (1) P. vivax PvP01 [15], P. vivax-like Pvl01, P. cynomolgi B strain [8], and P. knowlesi H strain [18] reference genomes and (2) the 13 Plasmodium reference genomes used for the phylogeny (the here-generated P. vivax-like Pvl01, P. vivax PvP01 [15], P. cynomolgi M version 2 [17], P. coatneyi strain Hackeri [63], P. knowlesi H strain [18], P. falciparum 3D7 [64], P. praefalciparum G01 [13], P. reichenowi PrCDC [65], P. gallinaceum 8A [30], P. ovale wallikeri PowCR01, P. ovale curtisi PocCR01, P. malariae PmUG01, and P. malariae-like PmlGA01 [14]) were identified using OrthoMCL version 2. 09 [66,67] after an all-against-all BLASTp (E-value threshold: 1 × 10−6). From these, we extracted different sets of one-to-one orthologous genes for the subsequent analyses: a set of 4,056 genes that included the one-to-one orthologues among the 4 restricted species—P. vivax, P. vivax-like, P. cynomolgi, and P. knowlesi (the Pv4sp set) —and a set of 2,943 among the 13 Plasmodium species considered here for the interspecies phylogenetic analysis (the Pl13sp set). The first set of orthologous groups identified was used for the detection of selection (see below); the second one was used to build the phylogeny and for the dating analyses. Amino acid sequences of the one-to-one orthologues were aligned using MUSCLE [68] or MAFFT [69], respectively, for the first and second set of orthologous groups. Prior to aligning codon sequences, we removed the low-complexity regions identified on the nucleotide level using dustmasker [70] and then in amino acid sequences using segmasker [71] from ncbi-blast, using default parameters. After MUSCLE/MAFFT alignments, we finally excluded poorly aligned codon regions using Gblocks (parameters: -t = p -b5 = h -p = n -b4 = 2) [72]. After going through all these alignment cleaning steps, we ended with a low number of missing data and no gap in our dataset of one-to-one orthologues (31 missing data over the 2,784 alignments of more than 50 amino acids). Most of the current methods available to estimate the timing of species splits make strong assumptions on the evolutionary models and often require an accurate mutation rate or calibration points, which are not always available. Here, we estimated the relative divergence times of Plasmodium species to free from temporal calibration of phylogenies and used 2 methods that do not depend on the estimation of a mutation rate: a method based on pairwise amino acid sequence divergences and Total Least Squares (TLS) regressions but assuming a constant rate of evolution across the Plasmodium lineages—the so-called Silva method [29]—and the RelTime method developed in Tamura et al. [31,32], which does not assume any specific model for the rate of evolution. We first built the phylogenetic tree of the 13 Plasmodium species considered here for the dating analyses using the software RAxML version 8. 2. 8 [73]. The tree was constructed using the concatenation of the 2,784 orthologous groups that have a final alignment of more than 50 amino acids (among the Pl13sp set of 2,943 groups). RAxML was called using the following options in order to automatically determine the substitution model that best fits the data model and so that 100 bootstraps were performed to assess the tree robustness: “-m PROTGAMMAAUTO -f a -# 100. ” The best amino acid model of substitution identified by RAxML for the dataset (the JTT model with empirical base frequencies JTT + F) was used in the subsequent analyses. The idea behind the Silva method [29] for relative age estimation is that the divergence between amino acid sequences in independent lineages is correlated and that the divergence regression slopes of the proteins of a species pair of reference with the divergence of the proteins of other species pairs reflect the relative age of those splits (see Silva and colleagues [29] for a detailed description of the method). Following this method, we obtained, for each of the 2,784 orthologous groups of the Pl13sp set that have a final alignment of more than 50 amino acids and for each species pair, the amino acid sequence divergence dAA through a pairwise comparison using PAML version 4. 7 [74], with the option “cleandata” set to 1 to remove sites with missing data from the dataset. The sequence divergences dAA were estimated using 4 different models of substitution (JTT, WAG, LG, and Dayhoff) to evaluate their influence on the estimates of the relative ages. An R script from the authors of this method [29] allowed the estimation of α, the slopes of the TLS regressions of the divergence of the proteins between every possible species pairs and the reference species pair, with the 95% CI by bootstrapping (n = 10,000 bootstraps). The slope α is an estimator of the relative age of the 2 considered species, relative to the reference species pair. To evaluate the influence of the choice of the reference species pair on the results, we used different species pairs to consider multiple divergence references: We considered the relative distance of the split of the speciation of (i) P. vivax and P. knowlesi (reference pair considered in the original paper), (ii) the 2 P. ovale species, and (iii) the P. malariae-like and P. malariae species. Because the model of Silva assumes a strict molecular clock [29]—which would not apply to all Plasmodium species, specifically P. falciparum because of its extreme GC content in comparison with other Plasmodium species—we used the clock calibration-free method RelTime [31,32] implemented in MEGA7 [75] to estimate divergence times when evolutionary rates vary. RelTime estimates relative rates for each branch of the tree without requiring a distribution of rate heterogeneity. The relative rates of 2 sister lineages are estimated assuming that the time to their most recent common ancestor is equal in the presence of contemporaneous sampling, and the rate of their ancestral branch is estimated as the average of their branch-specific relative rates. This framework is applied to subtrees of 3 or 4 taxa following a bottom-up strategy, from tips to the root to estimate relative rates and times, relative rates being scaled by setting the lineage relative rate of the ingroop root node to 1 (see Tamura and colleagues [32]). RelTime analysis was performed on the concatenation of the alignments of the same 2,784 orthologous groups we used to build the phylogenetic tree, removing sites with missing data (Complete-Deletion option) and using the same substitution model (JTT + F). As required by the RelTime analysis, we specified P. gallinaceum 8A [30] as the outgroup. The phylogenetic relationships of the RBPs in the P. cynomolgi–P. knowlesi–P. vivax–P. vivax-like group of species was reconstructed using a maximum likelihood analysis using RAxML [73] with 100 bootstrap replicates. The tree was visualized using Geneious software [60]. To investigate relationships between the P. vivax and P. vivax-like populations, we constructed a maximum likelihood tree using the filtered variant call set of SNVs limited to the higher allelic frequency genotypes identified within each sample using RAxML and PhyML (using general time-reversible [GTR] models) [73,76]. Trees were visualized using Geneious software [60]. All approaches showed the same final phylogenetic tree described in the “Relationships to worldwide human P. vivax isolates” section. We looked for signatures of introgression and/or ILS by performing a phylogenetic network using SplitsTree 4. 14. 6 [77] based on the alignment of a portion of the mitochondrial genome containing the cox1 and cytb genes, as in Liu and colleagues [10]. We included all but 2 (KF618538 and KF618534) samples from Liu and colleagues [10], all from Prugnolle and colleagues [9], and 3 of our samples with a high-quality–sequenced mitochondrial genome (Pvl08, Pvl09, and Pvl10), and we added P. cynomolgi as the outgroup (GenBank accession number: AB434919. 1). Gaps and monomorphic positions were excluded for the analysis, and the dataset resulted in an alignment of 85 sequences of 127 variable positions. We conducted a reticulate network using the RECOMB 2007 [78] and the NeighborNet distance transformation methods [79]. A phylogenetic tree based on the complete alignment, without exclusion of sites (i. e. , an alignment of 2,530 bp), was constructed using FastTree 2. 1. 5 (maximum likelihood) implemented in Geneious [60] for comparison with results obtained from Liu and colleagues [10]. To explore the influence of the amount of phylogenetic signal on phylogenetic reconstruction, we also performed phylogenetic trees as described in the previous section using a subset of SNVs: Trees were reconstructed using 100,200,300,400,500,600,800,1, 000, or 5,000 SNVs. The nucleotide diversity (π) is the average number of nucleotide differences per site between 2 sequences. This parameter is interesting to estimate because it gives valuable information about variation in prevalence and demographic histories of the parasites. For the P. vivax and P. vivax-like populations, we calculated the genome-wide nucleotide diversity (π) from VCF files using VCFTools [59]. The means of nucleotide diversities was compared between P. vivax and P. vivax-like species based on the Wilcoxon-Mann-Whitney nonparametric test. In order to identify genomic regions involved in parasite adaptation to the human host, i. e. , regions under positive selection, we performed branch-site tests. To search for genes that have been subjected to positive selection in the P. vivax lineage alone, after the divergence from P. vivax-like, we used the updated branch-site test of positive selection [80] implemented in the package PAML version 4. 4c [74]. This test detects sites that have undergone positive selection in a specific branch of the phylogenetic tree (foreground branch) using Likelihood Ratio Tests (LRTs) to compare models allowing positive selection or not. The selective pressure is defined as the ratio (ω) of nonsynonymous (i. e. , dN, keeping the amino acid) to synonymous (i. e. , dS, changing the amino acid) substitutions (dN/dS). Under neutrality, rates of synonymous dS and nonsynonymous dN substitutions are equivalent, so the ω ratio is expected to be equal to 1. Under purifying selection, the rate of nonsynonymous substitutions dN is expected to be lower than the rate of synonymous substitution dS because selection will prevent amino acid modifications: The ω ratio observed is then under 1. This process is acting in reducing the fixation of deleterious mutation. Finally, under positive selection, amino acid replacement will be chosen, with the dN rate expected to be superior to the dS rate, leading to an observed ω ratio superior to 1. This indicates the fixation of advantageous mutations. All coding sequences in the core genome with a one-to-one orthologous relationship among P. vivax Sal1, P. vivax-like Pvl01, P. knowlesi H strain, and P. cynomolgi B strain were used for this test (i. e. , the 4,056 gene Pv4sp set of orthologous genes). We obtained dN/dS ratio estimates per branch and gene for P. vivax and P. vivax-like lineages alone using Codeml (PAML version 4. 4c [74]) with a free-ratio model of evolution. Genes with a significant signal of positive selection in P. vivax only were compared to the ones obtained in P. falciparum from Otto and colleagues [13] (S4 Table in Otto and colleagues [13]), in order to identify, e. g. , essential proteins for human or vector infection. The data are deposited in the Dryad repository (doi: 10. 5061/dryad. 32tm1k4) [81].
Among the 5 species responsible for malaria in humans, Plasmodium vivax is the most prevalent outside Africa. It causes severe and incapacitating clinical symptoms with significant effects on human health. Yet little is known about its evolution, adaptation, and emergence in humans. The recent discovery in African great apes of its closest known relative—P. vivax-like—may help resolve the evolutionary history of P. vivax. This study aims to characterize the genome of P. vivax-like isolated from infections of the closest ape relative to humans to get a better understanding of the evolution of this parasite. A total of 11 P. vivax-like samples were obtained from infected chimpanzee blood samples and an infected mosquito collected in Gabon. We generated, and here present, the first 2 new genomes of P. vivax-like and 9 additional draft sequences. Genome-wide analyses provide new insights into the biology and adaptive evolution of P. vivax to different host species. Indeed, they highlight lineage-specific evolution of some gene families involved in key steps of the life cycle of P. vivax. Analyses also revealed that the divergence between P. vivax and P. vivax-like occurred before the one between P. falciparum and its sister species P. praefalciparum.
taxonomy
parasite groups
medicine and health sciences
methods and resources
plasmodium
parasite evolution
parasitic diseases
parasitic protozoans
parasitology
apicomplexa
phylogenetics
data management
protozoans
phylogenetic analysis
mammalian genomics
research and analysis methods
sequence analysis
computer and information sciences
sequence alignment
bioinformatics
malarial parasites
evolutionary systematics
animal genomics
eukaryota
database and informatics methods
genetics
biology and life sciences
genomics
evolutionary biology
organisms
Plasmodium vivax-like genome sequences shed new insights into Plasmodium vivax biology and evolution