Population genomics studies identify signatures of global dispersal and drug resistance in Plasmodium vivax

Abstract

Plasmodium vivax is a major public health burden, responsible for the majority of malaria infections outside Africa. We explored the impact of demographic history and selective pressures on the P. vivax genome by sequencing 182 clinical isolates sampled from 11 countries across the globe, using hybrid selection to overcome human DNA contamination. We confirmed previous reports of high genomic diversity in P. vivax relative to the more virulent Plasmodium falciparum species; regional populations of P. vivax exhibited greater diversity than the global P. falciparum population, indicating a large and/or stable population. Signals of natural selection suggest that P. vivax is evolving in response to antimalarial drugs and is adapting to regional differences in the human host and the mosquito vector. These findings underline the variable epidemiology of this parasite species and highlight the breadth of approaches that may be required to eliminate P. vivax globally.

Figure 1

Country of origin for 195 P. vivax patient isolates indicated on a map of P. vivax malaria endemicity. Asterisks mark the five published monkey-adapted reference genomes Brazil I, Salvador I, Mauritania I, North Korean and India VII and the published Belem and Chesson strains. The four additional monkey-adapted strains, Vietnam ONG, Nicaragua I, Panama and Thailand Pakchong, sequenced as part of this study are also indicated. Numbers of clinical isolates are given in parentheses. Lines do not represent accurate within-country locations; locations are listed in supplementary table 1. The map was adapted from the Malaria MAP Project

Figure 2

Enrichment of P. vivax DNA extracted from clinical samples using the hybrid selection method. (a) Enrichment of parasite DNA measured as the percentage of total sequence reads alignable to the P. vivax reference genome before versus after hybrid selection. Enrichment was variable but nearly universal. (b) Fold enrichment of parasite DNA, as compared to the percentage of alignable reads before hybrid selection. Hybrid selection resulted in the greatest fold enrichment for samples containing the lowest amount of parasite DNA before hybrid selection.

Figure 3

Global population structure of P. vivax. (a) Three-eigenvector PCA using variation data from all P. vivax isolates. Each isolate is colored according to its origin, and the percentage contribution of each eigenvector is indicated. (b) A maximum-likelihood phylogenetic tree with 50 bootstraps computed through alignment to the Salvador 1 reference genome. Colors and labels correspond to those in a, and the sister taxon to P. vivax, P. cynomolgi, was included as an outgroup. Select internal nodes are highlighted with open and closed circles, representing >40% and >70% bootstrap support, respectively. (c) Admixture analysis of the full set of variation data from 195 isolates using 100 bootstraps. Ancestry for each isolate was assessed according to an optimized cluster value of K = 5 (see supplementary Fig. 7 for admixture plots under six different K cluster values). Colors correspond to the subpopulations seen in a and b and were assigned to the five most deeply sampled subpopulations of P. vivax.

Figure 4

Diversity in orthologs across six P. vivax subpopulations and P. falciparum. Box-and-whisker plots of nucleotide diversity (π) calculated from 574 one-to-one orthologs in single-infection, high-quality isolates from Mexico, Colombia, Peru, Thailand, Myanmar and Papua New Guinea subpopulations are shown. The center and ends of the boxes represent the 25th, 50th and 75th percentiles of π values for each subpopulation, and the whiskers on each plot represent the range of π within 1.5 times the interquartile range of the lower and upper quartiles. Data points beyond the range of the whiskers are plotted as black dots. For P. falciparum, π was calculated from orthologs in a population of 12 African isolates previously sequenced at the Broad Institute (see URLs). P. falciparum was significantly less diverse than all plotted P. vivax populations (P < 0.001) using both a paired Student's t test and a Wilcoxon signed-rank test.

Figure 5

Signals of selection in P. vivax subpopulations. (a) Genome-wide view of the divergence (F_ST) between Old World and New World populations calculated from 1-kb windows. Plotted is a moving average across each of the 14 chromosomes, with a histogram of F_ST values on the right. The 99th percentile of the distribution of 1-kb windows forms a cutoff for outliers, which are labeled with the nearest gene annotation to the peak F_ST. (b) F_ST comparison as in a between the Thailand and Papua New Guinea populations. (c–e) Detailed plots for DHPS (c), DHFR-TS (d) and Pvs47 (e), overlaying F_ST values from a and b, the collapsed estimate of mean linkage disequilibrium (LD) calculated across a 50-kb window for each nucleotide, a moving average of π in 1-kb windows for isolates from Southeast Asia and Tajima's D in 1-kb windows. The location of each gene is shown (pink bar) above a schematic of the gene. Homozygous variants within coding sequences are indicated, including nonsynonymous changes, synonymous changes and deletions. Polymorphism columns are scaled by their allele frequency in the global population. (f) Haplotype map for the P. vivax genes DHFR-TS and Pvs47. Red blocks denote nonsynonymous changes, gray blocks indicate the presence of the reference allele and white blocks show regions of missing data. Nonsynonymous sites in DHFR-TS are limited to derived alleles.