Population genomic evidence of Plasmodium vivax Southeast Asian origin

Abstract

Plasmodium vivax is the most common and widespread human malaria parasite. It was recently proposed that P. vivax originates from sub-Saharan Africa based on the circulation of its closest genetic relatives (P. vivax-like) among African great apes. However, the limited number of genetic markers and samples investigated questions the robustness of this hypothesis. Here, we extensively characterized the genomic variations of 447 human P. vivax strains and 19 ape P. vivax-like strains collected worldwide. Phylogenetic relationships between human and ape Plasmodium strains revealed that P. vivax is a sister clade of P. vivax-like, not included within the radiation of P. vivax-like By investigating various aspects of P. vivax genetic variation, we identified several notable geographical patterns in summary statistics in function of the increasing geographic distance from Southeast Asia, suggesting that P. vivax may have derived from a single area in Asia through serial founder effects.

Fig. 1. Geographical origin of Plasmodium isolates and patterns of genomic variation.

(A) Country of origin of the 447 P. vivax and 19 P. vivax-like isolates used in this study. Within each country, isolates were collected at different locations. The chimpanzee pictogram represents African P. vivax-like isolates. (B) Local variation of genetic relatedness along the genome of P. vivax isolates visualized by multidimensional scaling (MDS) based on the local PCA approach (39). Each point represents a non-overlapping genomic window of 100 SNPs. On the basis of the variation of the MDS-1 coordinate, each window was classified in two groups: subtelomeric hypervariable region (orange) and core region (purple). (C) Genome scan of the SNPs on chromosome 5 (left y axis) identified for both P. vivax and P. vivax-like. The gray area in the background shows the total number of SNPs identified (left y axis), and the black line represents the number of SNPs shared by the two species. The right y axis represents the MDS-1 coordinates against the middle point of each window. Windows were classified in the subtelomeric hypervariable region (orange) or the core region (purple).

Fig. 2. P. vivax-like strains are structured in two distinct clades that form a sister monophyletic lineage to the human P. vivax.

(A) ML phylogenetic tree illustrating the relationships between P. vivax (PV) and P. vivax-like (PVL) strains, rooted using P. cynomolgi (PC, M strain). Two distinct P. vivax-like groups were identified: PVL.grp1 (in orange) and PVL.grp2 (in brown). The animal pictograms on each leaf indicate the primate host, gorilla, chimpanzee, or mosquito (unknown primate host), colored according to the country of origin (Gabon, Cameroon, and Ivory Coast). Open and closed circles indicate nodes with >80 and >90% bootstrap support, respectively. (B) PCA displaying the two first eigen vectors (EVs) and the proportion of genetic variance they explain. (C) Genome-wide visualization of recent recombination events between P. vivax-like individuals. The horizontal black rectangles represent the 19 P. vivax-like genome sequences, and the vertical colored lines represent the recombining genomic segments on the recipient individual. Colors indicate the lineage membership of the donor individual. (D) Differences in nucleotide diversity (π) between P. vivax-like lineages and P. vivax strains. (E) Multiple sequentially Markovian coalescent (MSMC) estimates of the effective population size (N_e) in the two P. vivax-like groups (PVL.grp1 in orange and PVL.grp2 in brown). Lines in lighter color represent 50 bootstrap resampling replicates of randomly sampled segregating sites.

Fig. 3. The structure of P. vivax core region of the genome is mostly due to genetic isolation of natural populations.

(A) PCA using the SNPs from the core region of the genome of the 447 P. vivax isolates showing the top two EVs and the accounted proportion of genetic variance. (B) ML phylogenetic tree built using SNP data from the core region of the genome. The star indicates the root of the tree, leading to the P. vivax-like outgroup. Specific internal nodes are highlighted with open circles and represent >80% bootstrap support. (C) The genetic ancestry proportion for individual strains, depicted as a vertical bar, was estimated using ADMIXTURE for each of the K = 6 inferred ancestral populations. (D) These same ancestral proportions for each population are displayed as pie charts on the world map. (E) Isolation by distance among populations. Pairwise estimates of F_ST/1 − F_ST were plotted against the corresponding geographical distances between countries. The Spearman correlation coefficient (R²) and the P value estimated using a Mantel test with 1000 permutations are shown.

Fig. 4. Relationships and gene flow between P. vivax populations.

(A) TreeMix ML tree of P. vivax populations with four migration edges (arrows) rooted with P. vivax-like (PVL), including bootstrap node support. (B) TreeMix residual matrix from the tree in (A). PNG, Papua New Guinea.

Fig. 5. The Southeast Asian origin of P. vivax is supported by patterns of nucleotide diversity, LD, and AAF.

(A) Genetic diversity of P. vivax regressed on the great circle distance across the world. The value at each pixel of the map corresponds to the Spearman correlation coefficient (R) between the expected genetic diversity in each population and the geographic distance between this population and the pixel. Black dots represent the sampling sites used in the regression analysis (where n ≥ 5 individuals). The landmark on the map represents the most negative correlation coefficient indicative of the most probable source of the range expansion. (B) Correlation between nucleotide diversity (π) and geographic distance, measured as the distance between the milestone indicated in (A) and each population. (C) Correlation between LD at 100 bp (measured as the normalized R²) and the geographic distance. (D) Histograms of AAF proportions in P. vivax populations with an increasing sample size (N = 5, 7, 10, and 20).

Fig. 6. Coalescent-based inference of the demographic history in each P. vivax population.

(A) Effective population size (N_e) variations back to the TMRCA inferred using MSMC and (B) inferred TMRCA in each population. The analysis was performed using five individuals from each population, assuming a mutation rate per generation (μ·g) of 1.158 × 10⁻⁹ and a generation time (g) of 0.18. The gray area represents the first 1000 generations, at which low resolution of the recent past effective population size is observed, likely due to the presence of sites with high probability of being called inaccurately.