Genome-wide association study of global Plasmodium vivax populations provides insights into the evolution of drug resistance

doi:10.21203/rs.3.rs-4235417/v1

Download PDF

Article

Genome-wide association study of global Plasmodium vivax populations provides insights into the evolution of drug resistance

https://doi.org/10.21203/rs.3.rs-4235417/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Increasing reports of chloroquine resistance (CQR) in Plasmodium vivax endemic regions has led to several countries, including Indonesia, to adopt dihydroarteminsin-piperaquine instead. Evidence for the major candidate, pvmdr1, as a putative determinant for CQR is conflicting. Using a genome-wide approach, we perform genomic analysis of 1,534 P. vivax isolates across 29 endemic countries, detailing population structure, patterns of relatedness, selection, and resistance profiling, providing insight into putative drivers of CQR. Differential selection metrics applied between isolates from low-grade and high-grade CQR regions revealed sweeps in a locus proximal to pvmdr1 and in transcriptional regulation genes. Our investigation of the temporal dynamics of selective sweeps in 106 isolates from Indonesian Papua, the epicentre of CQR, revealed pvmrp1 as an emerging candidate for piperaquine resistance. Overall, our work provides novel markers for resistance surveillance in candidate loci, supported by evidence of regions under recent directional selection in this continually evolving parasite.

Health sciences/Diseases/Infectious diseases/Malaria

Biological sciences/Genetics/Genomics

Plasmodium vivax is the most geographically widespread human malaria parasite and the leading cause of malaria outside of sub-Saharan Africa. In 2022, there were 6.9 million cases across 49 endemic countries across Central and South America, East Africa, Asia, and Oceania¹. Intensive efforts to combat the deadlier P. falciparum, particularly in areas that are co-endemic with P. vivax, has distributed resources away from P. vivax control programs, leading to its emergence as the dominant species, particularly in the Greater Mekong Subregion¹.

An added intricacy of P. vivax malaria is its complex biology, characterised by relapsing infections caused by dormant hypnozoite stages², extravascular hidden parasite reservoir^3,4, low parasitaemia infections that complicate diagnosis⁵, and the early circulation of transmissive stages⁶. The exclusive tropism for preferentially invading immature red blood cells has prevented the development of a long term continuous in vitro culture system for P. vivax⁷. Our understanding of the parasite’s life cycle, transmission, and biology has therefore been limited.

Chloroquine is the front-line treatment for P. vivax malaria in most endemic countries. First documented in the late 1980s, chloroquine resistant (CQR) P. vivax has led to the adoption of dihydroarteminsin-piperaquine (DHA-PPQ) as the front-line treatment in Indonesia, Papua New Guinea (PNG), and Cambodia^8,9. Indonesia and PNG have been the epicentre of high-grade CQR, but beyond these countries, the last two decades have seen increasing reports of CQR characterised by parasite persistence 28 days post-treatment, calling for a strategy change for P. vivax control.

Evidence of the molecular drivers of CQR in P. vivax is both weak and conflicting. The major candidate gene, pvmdr1, was initially posited as a mediator of CQR because of its high sequence homology with the orthologous pfmdr1, which is involved in CQR in P. falciparum¹⁰, despite contrasting reports of no correlation of pvmdr1 as a CQR determinant¹¹. Profiling of polymorphisms within the pvmdr1 gene led to the association of the Y976F mutation with CQR, due to increased chloroquine IC₅₀ in samples from Indonesian Papua and Thailand¹², cementing its use as a marker for CQR in subsequent studies across vivax endemic countries^13–15. However, the Y976F mutation has been reported in regions with no relationship to CQR^16–18. At least 50 pvmdr1 SNPs have been documented globally, yet none have emerged as a definitive CQR marker, questioning the extent and relevance of pvmdr1 in modulating CQR. Although profiling the prevalence of pvmdr1 polymorphisms has increased our understanding of its evolution across different drug pressure backgrounds^19,20, most studies focus solely on the pvmdr1 gene itself, biasing our understanding of the acquisition of CQR, which is hypothesised to be a multifaceted process involving numerous loci.

Extensive use of antimalarial drugs over several decades has resulted in high-resolution characterisations of recent selection events associated with drug resistance in Plasmodium spp, using genome-wide sequence data. Genome-wide analyses have shown evidence for recent positive selection in P. falciparum endemic regions with artemisinin resistance^21–23 and revealed selective sweeps around the pfcrt, pfmdr1, and pfaat1 genes, which are explicitly linked to CQR^24,25. Similar selection metrics applied to P. vivax genomes have revealed that the orthologous loci associated with antifolate resistance in P. falciparum have been subject to selective sweeps^19,20,26. More recent work has found evidence of selective sweeps around pvmrp1 in East Asian isolates, a gene associated with chloroquine and mefloquine resistance in P. falciparum, necessitating further investigation²⁰.

Understanding the population genetics and dynamics of P. vivax malaria, particularly in the context of drug resistance, is essential for successful control and elimination planning. With the expanding repertoire of whole genome sequence data for P. vivax, we can now investigate the temporal dynamics of populations pre- and post-chloroquine contraindication, to provide increased insight into the markers of CQR. Here, we leverage publicly available whole genome sequences to present the largest population genomics study to date of P. vivax. We provide an expanded insight into population structure, global ancestry, relatedness, and genomic diversity. Using a genome-wide approach, we perform intra- and inter-population analyses between isolates from low-grade and high-grade CQR regions to make inferences about novel loci that could be mediating CQR and the evolutionary forces that shape P. vivax populations.

P. vivax whole genome sequence data and clonality

A total of 499,206 high-quality bi-allelic SNPs were identified in the non-hypervariable regions of the P. vivax genome after filtering, comprising 1,534 isolates from 29 countries. In keeping with previous work²⁰, we divided these countries into 7 sub-regional populations based on the degree of geographic and genetic proximity: East Africa (N = 173), West Africa (1), South America (364), South Asia (156), South East Asia (SEA) (550), Maritime SEA (66), and Oceania (224) (Table S1).

Within-sample diversity, as a metric of multiclonality, was measured using the F_WS fixation index, where an F_WS ≥ 0.95 indicates an infection predominated by a single genotype. Here, a large proportion of isolates (71.7%) were monoclonal (F_WS ≥ 0.95). Regionally, mean F_WS was greatest in Maritime SEA (0.96), compared with South Asia (0.92), South America (0.92), East Africa (0.91), SEA (0.91), and Oceania (0.86) (Fig. S1, Table S1). These observations likely reflect trends in transmission intensity, where higher transmission intensity is often marked by higher clone multiplicity.

Global P. vivax isolates form distinct sub-populations to sub-continent level

An analysis of population structure using a principal component analysis (PCA) approach applied to 499,206 bi-allelic SNPs revealed P. vivax isolates form distinct and independent sub-populations, reflective of their continental or sub-continental origins (Fig. 1A). The first principal component separated East Asian and Oceanian populations from South American, South Asian, and African populations, producing three major geographical population centres, with South American populations being the most distinct. The population structure was supported by the corresponding SNP-based neighbour joining tree (Fig. 1B). The number of highly differentiating SNPs (F_ST ≥ 0.99) was positively correlated with geographic distance (Spearman’s Coefficient = 0.64) (Table S2). These findings suggest gene flow between neighbouring territories is contributing to this phylogeographic distribution and is supported by known waves of human migration²⁷.

An ADMIXTURE ancestral analysis inferred that global P. vivax isolates descend from ten ancestral populations (K1-K10), comprising two in East Africa (K1 and K5), one in South Asia (K1), three in South America (K2, K3, K9), three in SEA (K1, K8, K10), two in Maritime SEA (K4 and K6), and one in Oceania (K7) (Figs. 1C-1D). There was some concordance between ADMIXTURE populations and country of origin, where several ancestral populations were made up of predominantly (> 95%) one country, including K2 (Panama), K3 (Brazil), and K4 and K6 (Malaysia) (Table S1, Fig. S2). Except for K1, which was characterised by African, South Asian, and SEA isolates, all ancestral populations comprised isolates from the same sub-regional grouping. This population structure was supported by a PCA plot coloured by dominant ancestry (Fig. S3).

Pairwise relatedness supports model of isolation by distance

To further dissect the global P. vivax population structure, pairwise relatedness of isolates was inferred by calculating identity-by-descent (IBD), which describes shared evolutionary history. Countries which presented the highest median pairwise IBD include Panama (0.97), Mexico (0.23), and Malaysia (0.19), suggesting reduced outcrossing within these populations (Table S3). While the fractional IBD values of Malaysian and Mexican populations approach the expected value (0.25) of half-siblings in an outbred population, the value for Panama is inflated due to the persistence of a single clone in the region for a decade²⁸. Within the remaining subpopulations, IBD sharing was low (< 6%), revealing that samples were predominantly weakly related, with a minority of very highly related samples (Fig. S4). Regionally, median pairwise IBD was highest in Maritime SEA (0.15) and lowest in SEA (0.01). In agreement with F_ST observations, there was a moderate correlation between inter-regional median pairwise IBD and geographic distance (Spearman’s Coefficient = -0.39) (Table S4), consistent with the model of isolation by distance²⁹ which posits that populations in closer geographic proximity have increased genetic similarity.

IBD sharing reveals putative selective sweeps at P. vivax drug resistance loci

To investigate patterns of shared ancestry intrachromosomally, we analysed genome-wide IBD fractions calculated across 10 kb sliding windows (Fig. S5). We focused on windows with the top 5% of mean IBD fractions, within or around drug resistance loci, to identify candidates that have most likely experienced selective sweeps. Due to the high genetic relatedness of isolates from Malaysia, Mexico, and Panama (Figs. 1A-1B) and their inflated pairwise IBD values (Fig. S6, Table S3), we excluded them from further analysis. Overall, regions with the highest IBD fractions spanned antigenic loci (pvmsp1, PVP01_0728900; pvmsp5, PVP01_0418400; pvdbp, PVP01_0623800), genes involved in life-cycle specific processes (pvlisp2, PVP01_0304700), and candidate drug resistance loci (pvmrp1, PVP01_0203000; pvdhfr, PVP01_0526600; pvmdr1, PVP01_1010900; pvdhps, PVP01_1429500) (Table S5).

The strongest signals overall (0.35) were observed in Indonesia on chromosome 10 between 320–330 kb (Table S5), a region encompassing two conserved Plasmodium proteins of unknown function (PVP01_1007200, PVP01_1007250), with the former being a pseudogene. This region is of particular interest as it is ~ 140 kb downstream of pvmdr1, the gene putatively responsible for CQR in P. vivax. This peak of IBD sharing at 320–330 kb is also found within a large tract of significant IBD values spanning 200–500 kb on chromosome 10 (Fig. S7, Fig. S8, Table S6). Although peaking at 320–330 kb, the median fractional IBD value of the region was 0.12, and the IBD value at pvmdr1 was 0.11. Interestingly, we found high proportional IBD sharing at this downstream pvmdr1 locus (320–330 kb; PVP01_1007200, PVP01_1007250) in isolates from PNG (0.19), Sudan (0.13), and Brazil (0.11). As with Indonesian isolates, in the Sudanese population, this was found within a tract spanning 280–500 kb, with a median IBD value of 0.17. At pvmdr1 itself, countries with significant IBD values were Sudan (0.26), Peru (0.19), Ethiopia (0.19), and Brazil (0.13).

The pvdhps and pvdhfr genes on chromosomes 14 and 5 respectively are associated with resistance to the combination antimalarial therapy sulfadoxine-pyrimethamine (SP). Although not used to treat vivax malaria, it is used as chemoprevention in regions co-endemic with P. falciparum. We observed a trend of differential IBD values surrounding the genes, where isolates had high fractional IBD for one of the SP-resistance genes, but not the other. Isolates from Eritrea, Indonesia, and Myanmar had much greater fractional IBD values at pvdhps (0.25, 0.22, and 0.21, respectively) than at pvdhfr (0.11, 0.09, < 0.01, respectively) (Table S5).

Finally, we observed significant IBD values in pvmrp1, a gene on chromosome 2 associated with resistance to chloroquine, primaquine, and mefloquine in P. falciparum³⁰ in Vietnamese (0.22), Thai (0.21), Cambodian (0.18), Burmese (0.16), and Colombian (0.14) isolates. We found no evidence of significant IBD sharing in any population at the resistance candidates pvcrt-o (PVP01_0109300), pvpm4 (PVP01_1340900), pvk13 (PVP01_121100), or in the sister genes to pvmrp1 and pvmdr1, pvmrp2 (PVP01_1447300) and pvmdr2 (PVP01_1259100). However, there was a signal from pvaat1 (PVP01_1120000), a gene newly implicated in CQR in P. falciparum²⁵, in isolates from PNG (0.06).

Evidence of selective sweeps in candidate drug resistance loci

A genome-wide scan for the top 1% of genes under recent positive selection in monoclonal isolates was performed using the iHS metric (Table S7). Multiple strong signals (N = 125) were found in the genes encoding the surface proteins pvmsp1 and pvmsp5 (P < 1 x 10^− 24), which are under selection pressure due to their interactions with the host immune response. The most significant values were in Afghanistan (P < 1 x 10^− 23), Thailand (P < 1 x 10^− 14), Cambodia (P < 1 x 10^− 11), and Vietnam (P < 1 x 10^− 10). Similarly, other hotspots of selection pressure were found in the cytoadherence linked asexual protein (pvclag, Ethiopia (P < 1 x 10^− 12), Thailand (P < 1 x 10^− 10)), the liver-specific protein 2 (pvlisp2, Afghanistan (P < 1 x 10^− 14), Indonesia (P < 1 x 10^− 6), Cambodia (P < 1 x 10^− 6)), ApiAP2 transcription factors (PVP01_1418100, Cambodia (P < 1 x 10^− 7); PVP01_1440600, Cambodia and Vietnam (both P < 1 x 10^− 5)), and the invasion proteins pvrbp1a (Cambodia, Indonesia, Pakistan; P < 1 x 10^− 10) and pvrbp2b (Afghanistan and India; P < 1 x 10^− 5). In agreement with our IBD findings, a series of 18 SNPs (P < 1 x 10^− 8) were under selection in Indonesian isolates downstream of pvmdr1 on chromosome 10 (243480–319423), with peaks at a SNP upstream a translation initiation factor (PVP01_1005600, P < 1 x 10^− 8) and a SNP in the PVP01_1007200 gene (P < 1 x 10^− 7). Isolates from Thailand and Pakistan also had SNPs under selection in PVP01_1007200 (P < 1 x 10^− 5).

To identify signals of directional selection, the cross-population metric, Rsb, was used at both a country and regional level. The most common SNPs under differential selection encompassed the cluster of SNPs with significant iHS scores in the Indonesian population downstream of pvmdr1 (N = 248), differentiating Indonesian isolates from those in SEA (P < 1 x 10^− 13), East Africa (Eritrea and Ethiopia; P < 1 x 10^− 7), South Asia (Afghanistan, India, Pakistan; P < 1 x 10^− 10) and South America (Peru; P < 1 x 10^− 6) (Fig. 2, Table S8). Beyond SNPs in merozoite surface antigens, we found evidence of selective sweeps via extended haplotype homozygosity at pvmrp1 in isolates from SEA (Cambodia, China, Vietnam, Thailand) when compared to isolates from South Asia (Afghanistan and India, P < 1 x 10^− 6). Most of these SNPs were in the promoter region, except for L1361F and 1232I present exclusively in Thai isolates. We also found SNPs within pvmrp2 under differential selection between isolates from Vietnam and China, Cambodia, and India (P < 1 x 10^− 6). Similarly, we observed differential signals in pvdhps in SEA isolates (Cambodia, Thailand, Vietnam; P < 1 x 10^− 11), and in pvdhfr in South Asian isolates (Afghanistan, India; P < 1 x 10^− 7). The SNPs S513R and L845F in MDR1 were under differential selection between Thai isolates compared with Vietnamese (P < 1 x 10^− 9) and Cambodian isolates (P < 1 x 10^− 6).

Sub-regional differences in the frequency of drug resistance mutations in P. vivax populations

Identifying mutations in putative drug resistance loci can reflect the epidemiology of transmission in both local and global contexts. We therefore evaluated the prevalence of resistance haplotypes (non-synonymous mutations) in all isolates in genes with prior evidence of contributing to antimalarial resistance in P. vivax: pvmdr1 (33), pvmrp1 (53), pvdhfr (17), and pvdhps (20) (Table 1) and several additional genes of interest (Table S9). Multiple pvmdr1 mutations have been previously associated with CQR or reduced chloroquine susceptibility^10,14,18, including S513R, G698S, M908L, Y976F, and F1076L, and all, except S513R, are present in the P. vivax PvP01 reference. We observed the frequencies of the putative resistance alleles at 70.9% (698S), 97.9% (908L), 100% (976F), and 73.4% (1076L). The F1076L mutation was observed in all isolates from East Africa and Oceania and had lowest frequencies in South American (81.6%) and SEA (19.4%) isolates. Similarly, the G698S mutation was fixed in isolates from Oceania and Maritime SEA and had lowest frequency in South American (14.0%) isolates. The S513R mutation was found exclusively in isolates from SEA (14.0%), South Asia (37.2%), and East Africa (74.0%). There were lower frequency SNPs with regional specificity, including D500N (South America, 24.2%), V221L (South America, 21.1%), and L845F (Across Asia and Oceania, 9.5%). Only one isolate from Indonesia had a non-synonymous SNP in pvmdr1 (L845F).

Table 1

Prevalence of SNPs in genes that have been associated with drug resistance in *P. vivax*
Gene	Position	Mutation	East Africa (N = 173)	South Asia (N = 156)	SEA (N = 550)	Maritime SEA (N = 66)	Oceania (N = 224)	South America (N = 364)
pvdhfr	1077534	R58K					0.5	5.5
pvdhfr	1077535	R58S	37.6	62.8	2.55	59.1	6.3	61.3
pvdhfr	1077711	N117T	6.9	50.0	30.7	90.9	69.2	17.0
pvdhfr	1078180	N273K	1.7	10.3	0.2
pvdhfr	1077878	I173L				56.1		8.52
pvdhfr	1077530	F57I		0.64	78.0	89.4	71.0	0.5
pvdhfr	1077543	T61M				31.8	41.5	27.5
pvdhps	1270119	A647V	24.9	4.5	0.2
pvdhps	1270401	A553G		10.3	29.8	86.3	10.7
pvdhps	1270911	G383A*	71.1	83.3	9.6	4.5	16.1	39.8
pvdhps	1270914	S382C			0.4			12.9
pvdhps	1270915	S382A			10.9
pvdhps	1271444	M205I	71.7	0.6	96.2	6.1	1.8	58.5
pvdhps	1271634	E142G	64.7			59.1	10.7
pvmdr1	479908	L1076F*		1.9	19.5	1.5		81.6
pvmdr1	480412	L908M*	0.58					8.52
pvmdr1	480552	A861E		5.8	3.6			0.82
pvmdr1	480601	L845F		14.7	12.0	7.6	0.4
pvmdr1	481042	S698G	27.2	51.9				86.0
pvmdr1	481595	S513R	74.0	37.2	14.0
pvmdr1	481636	D500N						23.6
pvmdr1	482473	V221L						21.2
pvmrp1	154107/8	A1606V						5.2
pvmrp1	154168	H1586Y		5.1	2.6			19.8
pvmrp1	154215/6	D1570F					14.7
pvmrp1	154992	I1478V	24.3	17.9	0.2			9.9
pvmrp1	154668	G1419A	13.4	16.7				24.7
pvmrp1	154831	L1365F					17.9
pvmrp1	155305	L1207I			81.5
pvmrp1	156208	E906Q
pvmrp1	156563	E787D	0.6	3.8	3.5	1.5	1.8
pvmrp1	158148	R259I	35.3	27.6	0.6			11.3
pvmrp1	158223	T234M			78.4
pvmrp1	158272	Y218D	82.7	64.1	92.4	98.5	87.5	60.4
pvmrp1	158545	V127I	82.7	60.3	92.4	98.5	87.9	65.4
Allele frequencies were calculated in isolates if, at a country level, their frequency was ≥ 10. Allele frequency is bolded if ≥ 50%. East Africa (Eritrea, Ethiopia, Sudan); South Asia (Afghanistan, India, Pakistan); SEA (Cambodia, China, Myanmar, Thailand, Vietnam); Maritime SEA (Malaysia); Oceania (Indonesia, Papua New Guinea); South America (Brazil, Colombia, Mexico, Panama, Peru). Asterisked SNPs are those where the reference has the putative resistance allele. SNP prevalence was calculated solely for isolates with homozygous alternate calls.

Mutations in the pvmrp1 orthologue, pfmrp1, decrease susceptibility to chloroquine³⁰. Although there were several SNPs approaching fixation globally (V127I; 81.2%, Y218N, 80.4%), most SNPs showed regional specificity. South American isolates had the greatest number of unique pvmrp1 SNPs including A1606V (5.22%), T1525I (4.4%), and C1018Y (3.3%). Other SNPs with regional specificity of varying frequencies included T234M (SEA, 78.4%), L1365F (Oceania, 17.9%), F560I (East Africa (8.1%) and South Asia (1.9%)). We observed a novel mutation, D1570F, exclusive to Oceanian isolates (Indonesia (16.5%), PNG (10.0%)) (Table 1).

Antifolate (SP) resistance is attributed to mutations in DHFR and DHPS. In pvdhfr and pvdhps, numerous SNPs structurally correspond to resistance conferring mutations in P. falciparum. The pvdhfr S58R and S117T mutations and the pvdhps A383G mutation were present in the PvP01 reference, so we report the frequencies of the wild-type (reference) allele (Table 1). The triple pvdhfr mutants F57L-S58R-S117N (LRN) coupled with the pvdhps double mutant A383G-A553G (GG) are associated with clinical SP failure³¹. We found 3 instances of the pvdhfr LRN haplotype in isolates from Myanmar, PNG, and Peru (Table S10). The pvdhps GG haplotype was globally more prevalent, found in India (33.3%), Indonesia (1.5%), Malaysia (3.2%), and Thailand (1.9%). Another pvdhfr haplotype of concern is the pvdhfr L/IRMT(F57L/I-S58R-T61M-S117N/T) quadruple mutant, found at moderate prevalence in East Asian and Oceanian isolates: China (16.7%), Malaysia (30.6%), Thailand (65.4%), Myanmar (28.6%), Indonesia (64.9%), and PNG (10.0%). Although we did not observe the pvdhps SAKAV mutant (S382A-A383G-K512E-A553G-V585A)³², we observed the alternate variant, MSAAK (M205I-S382A-A383G-A553G-K512M), exclusive to Thai isolates (4.5%). Finally, although not structurally linked to sulfadoxine resistance, we observed the Q142G mutation in various haplotypes specific to isolates from East Africa (Ethiopia, Eritrea, and Sudan), and the S382C-M205I double mutant specific to isolates from South America (Brazil, Guyana, Panama, and Peru).

Moderate degree of sequence conservation in global mdr1 haplotypes

Haplotype networks were constructed to visualise the global diversity of pvmdr1 and explore haplotypes in low-grade vs. high-grade CQR regions. We observed 169 distinct haplotypes, of which 93 (55.1%) were singletons and 151 (89.3%) had less than 10 observations. A median-joining haplotype network was estimated for the remaining 18 (10.7%) haplotypes identified across 1,238 samples (Fig. 3). Haplotype IV, representing the 598P synonymous mutation, was the most prevalent, present in all geographic subregions (41.1%). Most Oceanian isolates (80.8%) were represented by this major haplotype, indicating that isolates from low and high-grade CQR regions have high pvmdr1 sequence similarity. Globally, 44.0% of isolates had pvmdr1 haplotypes comprised solely of synonymous mutations, suggesting a substantial degree of conservation at pvmdr1. Intra-region estimates of haplotype diversities (h) ranged from 0.33–0.93, with Oceanian isolates having the lowest genetic diversity (h = 0.33), and South Asian isolates with the highest (h = 0.93) (Table S11).

Temporal trends in markers under selection in Indonesian isolates

An ex vivo susceptibility study of P. vivax isolates from Indonesian Papua revealed the persistence CQR between 2005-2018³³. With conflicting evidence of pvmdr1’s role in mediating CQR and maintained phenotypic CQR, we sought to investigate temporal trends in genome-wide signals of differential selection within the Indonesian population. We divided all monoclonal isolates into two sub-populations: pre-2014 (N = 75) and post-2014 (N = 29), with 2014 being 10 years post-contraindication of chloroquine in Indonesia, and a timeframe we perceive adequate for genotypic changes to occur. Comparing iHS data revealed that the hotspot of SNPs downstream pvmdr1 was present solely in the pre-2014 population (Figs. 4A-B). A scan for the top 1% of genes under differential selection (Rsb) between the pre- and post-2014 populations revealed the novel pvmrp1 mutation D1570F had a significant hit (P < 1 x 10^− 8) (Figs. 4A-4C). Pairwise F_ST revealed that this SNP was the most highly differentiating between the two populations (F_ST = 0.34; Median F_ST = 0.027) (Table S12, Fig S6). This SNP first occurred in 2 isolates from 2011, and was present in 58.6% of post-2014 samples, compared with only 20.0% of pre-2014 isolates.

Taking this further, we divided all Indonesian isolates into 5 groups based on year of sample collection (2008-09 (N = 33), 2010-11 (N = 49), 2012-13 (N = 62), 2014-15 (N = 30), and 2016-17 (N = 18)) and applied the same intra- and inter-population selection metrics. Across all groups, hotspots of selection were observed on chromosomes 4 and 7, corresponding to pvmsp5 and pvsmp1, and on chromosome 14, corresponding to proteins expressed in gametocytes (G377, PVP01_1467200) and histone methylation machinery (CARM1, PVP01_142800; JmJC1, PVP01_1430400) (Fig. 4D-F, Table S13). Significant iHS scores at the downstream pvmdr1 locus were observed only in populations 2012-13 and 2014-15 (all P < 1 x 10^− 6), indicative of transient directional selection at this locus. Differential selection metrics applied between the divergent 2008-09 and 2016-17 groups (4 years vs. 12 years post-contraindication of chloroquine) revealed SNPs in the sporozoite-invasion associated protein 2 (pvsiap2, P < 1 x 10^− 8) and the cysteine-rich protective antigen (pvcyrpa, P < 1 x 10^− 7) under differential selection (Table S14). Contrastingly, the F_ST metric revealed the pvmrp1 SNP L1365F was the most highly differentiating (F_ST = 0.54) between the two populations, with prevalence of the SNP increasing from 0% in 2008-09 to 38.5% in 2016-17 (Fig. 5).

Understanding the epidemiological and evolutionary dynamics of malaria using genome-wide data can provide essential biological insights that can be harnessed in the global malaria control and elimination agenda. This is of particular importance in the context of P. vivax malaria, where the biological determinants of resistance to the front-line treatment, chloroquine, are still unknown. Here, using genome-wide SNPs, we were able to show distinct subpopulations of vivax isolates to sub-continent level, with classic signs of population differentiation by geographic separation. Through profiling SNPs at candidate resistance genes and across the vivax genome, we describe new loci with population differentiation effects between isolates from low-grade and high-grade CQR regions. We describe the temporal dynamics of SNPs proximal to pvmdr1, the putative determinant of CQR, in Papua Indonesian isolates due to long-standing high-grade CQR, presenting hypotheses for the absence of selection pressure within the gene itself over a decade post-chloroquine contraindication. Finally, we describe pvmrp1 as an emerging putative marker for primaquine resistance in P. vivax.

Exploration of the population structure of global P. vivax isolates revealed the genetic structure largely reflects geographic structure, with significantly higher F_ST values between more spatially disparate populations, consistent with previous findings^{19,20,34–36}. This observation also agrees with a model of isolation by distance, which indicates increased gene flow and genetic similarity between isolates that are more geographically proximal. At a local scale, we observed that intra-subcontinental populations were less structurally defined, except for South American and, to a lesser extent, SEA populations. These observations are not only a reflection of the complex geographic and ecological niches within these subcontinents, but in a South American context, could be reflective of independent introductions of P. vivax in the region, as seen in P. falciparum³⁷. ADMIXTURE analysis of nuclear SNP data revealed moderate shared South Asian ancestry in East African isolates from Eritrea, Uganda, Sudan, and Madagascar, which has been previously described^20,38, and coincides with human migration of South Asians to East African regions outside the Horn of Africa³⁹. However, this did not reveal an Ethiopian-specific ancestral population, which agrees with prior observations of Ethiopian populations being less structured and more diversified³⁶. In the SEA subgrouping, Cambodian and Vietnamese isolates were genetically and ancestrally distinct from isolates from Thailand, Myanmar, and China, which has recently been documented⁴⁰. This is likely a result of the intimate historical links between the two countries throughout the 20th century.

Genome-wide IBD analyses provided helpful insights into the genomic architecture of global P. vivax populations, especially in the East African context. We observed, across all three major populations, high fractional IBD at pvdhps (Eritrea) and pvdhfr (Ethiopia, Sudan). This is most likely a relic of the intensive use of SP as a first-line drug in the treatment of P. falciparum malaria, which is co-endemic with P. vivax in the region, and due to inadvertent SP drug pressure on P. vivax isolates by seasonal malaria chemoprevention regimes⁴¹. When we paired these data with our SNP profiling results, we observed low frequencies of single and double mutants in pvdhfr (S58R, S117N/T) that structurally corresponded with linked to pyrimethamine resistance in P. falciparum, directly contrasting prior findings of pvdhfr double mutants approaching fixation in Ethiopian vivax populations⁴². Similarly, in pvdhps, we observed an abundance of the wildtype alleles (553A/383G) or other pvdhps mutants that have not yet been associated with sulfadoxine resistance (M205I/E142G). The apparent heterogeneity in pvdhfr/pvdhps mutants, coupled with clear evidence of P. falciparum treatment strategies shaping the vivax resistance landscape necessitates further investigation in the East African context, especially due to the moderate rise in P. vivax cases in the region between 2019–2022¹.

While in vitro susceptibility studies have shown a role of the Y976F mutation in modulating CQR in isolates from Indonesian Papua¹², we found this mutation fixed in global populations (100%). Our observation agrees with no significant change found in chloroquine sensitivity in P. cynomolgi introduced with the Y976F mutation⁴³. Similarly, we observed the M908L mutation fixed globally (97.9%). With ex vivo measurements associating M908L with reduced susceptibility to chloroquine, DHA, and mefloquine¹⁷, if M908L and Y976F do mediate CQR, it implies that CQR is intrinsic in vivax populations. The F1076L mutation is approaching fixation in P. vivax isolates from Sabah, Malaysia^44,45 and Indonesian Papua¹⁰, yet again, as with the Y976F and M908L mutations, we found it at fixation in several countries with high clinical efficacy of chloroquine, such as Afghanistan⁴⁶. The presence of geographically diverse pvmdr1 haplotypes with representation of isolates from high-grade and low-grade CQR regions supports the finding of multiple independent local diversification events at pvmdr1⁴⁷, and strengthens the evidence against these mutations being reliable markers for CQR.

Operating on the assumption that these SNPs are indeed reliable markers, we found no evidence of selective sweeps encompassing the pvmdr1 gene itself within or between any population, despite our genome wide IBD analyses inferring putative selection events at pvmdr1 in Sudanese, Ethiopian and Brazilian isolates. Moreover, differential selection of pvmdr1 SNPs was within the SEA population only. Considering the trend of declining chloroquine sensitivity in regions such as Malaysia⁴⁴, Vietnam⁴⁸, and the China-Myanmar border¹⁷, we hypothesised that there should be signals of selection in pvmdr1 that could indicate this. One of the strongest selective sweeps in the Indonesian Papua population was at a locus proximal to pvmdr1, encompassing 18 SNPs in a ~ 76kb region. Strong directional selection of this region specific to the Indonesian population at genes PVP01_1007200-250 has been identified previously⁴⁹, however, investigations into its role, if any, in mediating CQR have been limited. From IBD data, we could also infer putative selection events at this locus in isolates from PNG and Sudan but limited further investigations due to limited sample size. Here, our study of the temporal dynamics of this cluster of SNPs in the Indonesian context revealed that the selective sweep takes place between 2012–2015. Given the temporally transient nature of this genomic signature, we hypothesise that they are either neutral alleles, and therefore not linked to CQR, or instead, reflect changing local dynamics and drug pressures within the country, considering that chloroquine is still readily available in the private sector in Indonesian Papua⁵⁰. It is unclear of the role these SNPs play in PNG and Sudan.

Indonesian Papua isolates have maintained a CQR phenotype, despite the removal of drug pressure exerted by chloroquine since its contraindication in 2004³³. Assuming the Indonesian Papua isolates herein are thus genotypically CQR, we identified genes involved in cytoadherence (CLAG; PVP01_1401400), gametocyte development (G377; PVP01_1467200) and genes encoding histone methylation machinery (CARM1; PVP01_1428800, JmJC1; PVP01_1430400) at the core of a selective sweep on chromosome 14. Similarly, isolates from Cambodia and Vietnam, both countries with declining chloroquine sensitivity⁵¹, had strong selective sweeps at an Ap2 domain transcription factor (PVP01_1418100) on chromosome 14, also observed in prior studies⁵². Isolates from PNG had large, non-overlapping tracts of high fractional IBD across genes involved in transcription and translation on chromosome 14. Expanding the diversity of transcriptional regulatory elements could alter the expression profile of their downstream targets, which could include drug resistance-associated genes, and warrants further investigation.

Strong signals of positive selection at the pvmrp1 gene has led to its emergence as a candidate drug resistance gene^20,52,53. We found the previously described L1207I SNP under selection in Thai isolates and the pvmrp1 gene with high fractional IBD in Vietnamese, Thai, and Cambodian isolates. Although not under selection, we observed the novel SNP D1570F in Indonesian isolates. Investigating the temporal dynamics of this novel mutation revealed an intriguing pattern. The first appearance of the corresponding SNP was in an isolate from 2011 and it had strong population differentiation effects between isolates pre- and post-2014. Further investigation of pvmrp1 between isolates specifically from 2008-09 and 2016-17 revealed that along with the D1570F mutation, the L1365F mutation was the most differentiating SNP between the two populations. With the adoption of DHA-PPQ as a novel front-line antimalarial in Indonesia, we observed, as DHA-PPQ drug pressure increases with time, prevalence of these pvmrp1 SNPs increases. Between 2008–2018, there has been a decline in P. vivax susceptibility to piperaquine in Indonesian Papua³³. In P. falciparum, knockout of the orthologous pfmrp1 gene resulted in an increased susceptibility to piperaquine³⁰. Moreover, in a P. falciparum genetic cross, progeny with mutant pfmrp1 were only found in parasites harbouring SNPs in kelch13 that are known drivers of artemisinin resistance⁵⁴. Taken together, the increasing frequency of this novel pvmrp1 SNP in the Indonesian population and the selective sweeps at pvmdr1 in SEA populations are of concern and should be the subject of surveillance. Leveraging the P. cynomolgi or P. knowlesi⁵⁵ in vitro culture models to introduce the D1570F and L1365F mutations to investigate drug susceptibility will be essential to gain important insights in the role of pvmrp1 in mediating piperaquine resistance.

Overall, the work described here has provided the largest study to date of global P. vivax isolates, providing insights on population structure, admixture, markers of IBD, differentiation, and selection signatures in the context of drug resistance. Although the markers of CQR remain elusive, our findings of directional selection and selective sweeps in the candidate resistance gene pvmrp1 highlights the need for broader genotypic and phenotypic surveillance of P. vivax to complement elimination efforts.

Sequence data and raw reads processing

Illumina whole genome sequencing data were obtained from the publicly available Pv4 dataset generated by the MalariaGEN Community Project⁵⁶ for P. vivax (N = 1,895) and newly sequenced isolates sent to the UKHSA Malaria Reference Laboratory from imported UK cases (N = 47), totalling 1,942 isolates available for analysis. After quality control, the final dataset included 1,534 (79.0%) isolates from 29 countries: East Africa (N = 173, Eritrea (13), Ethiopia (138), Madagascar (4), Sudan (13), Uganda (5)), South America (N = 364, Brazil (139), Colombia (58), Guyana (3), Mexico (20), Panama (47), Peru (97)), South Asia (N = 154, Afghanistan (63), Bangladesh (1), India (48), Iran (5), Pakistan (37)), South East Asia (SEA) (N = 550, Cambodia (218), China (12), Laos (2), Myanmar (28), North Korea (1), Thailand (179), Vietnam (110)), Maritime SEA (Malaysia (62), The Philippines (4)), Oceania (n = 224, Indonesia (194), PNG (30)), and West Africa (N = 1, Mauritania (1)).

Variant calling and quality control

All raw Illumina sequencing reads were first trimmed using trimmomatic⁵⁷ (v0.39, parameters PE -phred33, LEADING:3, TRAILING:3, SLIDINGWINDOW:4:20 MINLEN:36). Trimmed reads were aligned to the reference genome, PvP01_v1⁵⁸ using bwa-mem (v0.7.17, default parameters)⁵⁹. The resultant BAM files were processed with samtools functions fixmate and markdup (v1.9, default parameters)⁶⁰. A training-set of high-quality P. vivax single nucleotide polymorphisms (SNPs) from previously published literature was used to calibrate variant calling²⁰. Using this set, the GATK BaseRecalibrator and ApplyBQSR functions⁶¹ were run to produce improved BAM files for all isolates^19,20. SNPs and indels were identified using the GATK HaplotypeCaller function (v.4.4.0.0, parameters: -ERC GVCF) to produce individual sample variant call format (VCF) files. The resultant VCFs were merged to create a multi-sample VCF using the GATK CombineGVCF function (default parameters). A total of 3,671,958 unfiltered SNPs were identified across 1,942 isolates. The merged VCF was filtered iteratively (described in full here²⁰) to produce the final dataset for subsequent analyses. Briefly, variants identified in the hypervariable subtelomeric regions were removed by mapping the 14 chromosomal sequences against the P. vivax PvP01_v1 reference, leaving variants only within the core P. vivax genome. Variants with a Variant Quality Score Log-Odds (VQSLOD) score < 0, representing variants that are most likely false, were excluded. Variants where > 40% of SNPs had missing genotype data were excluded. Monomorphic SNPs, heterozygous SNPs, and indels were also excluded. A set of 175 isolates from Indonesia that passed quality control, were taken from the original Pv4 release VCF (ftp://ngs.sanger.ac.uk/production/malaria/Resource/30) and merged into the VCF created within the present study. The final dataset encompassed 499,206 high-quality bi-allelic SNPs across 1,534 samples. SNPs were annotated with predictions of their downstream coding effects using the SnpEFF software⁶².

Population genetics analyses

A binary matrix of pairwise genetic distances was constructed from the filtered biallelic VCF in PLINK (v1.9, default parameters)⁶³. Using the distance matrix, population structure was assessed by conducting a Principal Component Analysis (PCA) and constructing a neighbour joining (NJ) tree using the R package ape (v5.7). The biallelic VCF of 499,206 SNPs was filtered using bcftools (v1.17) view function and used to select SNPs with a minor allele frequency > 1%, producing a multi-sample VCF of 115,022 high-quality bi-allelic SNPs. A maximum-likelihood (ML) tree was created using IQTREE (v2.2.3) on a nucleotide FASTA alignment of the MAF filtered VCF, with 1000 bootstraps to measure branch quality and robustness. Both NJ and ML trees were visualised in iTOL⁶⁴. To further investigate population structure, we used the ADMIXTURE software⁶⁵, a tool used to estimate individual ancestries and population allele frequencies in SNP genotype datasets. PLINK (v1.9)⁶³ was used to convert the MAF-filtered multi-sample VCF to a BED file. Using ADMIXTURE, the most likely number of subpopulations (K) was obtained using cross-validation error of 1–10 dimensions of eigenvalue decay. The output was visualised in R (v4.2.3). The population differentiation metric, fixation index (F_ST), was calculated pairwise to assess SNPs driving allele frequency differences between populations at a country and regional level using the vcftools (v0.1.17) function -weir-fst-pop.

Within-infection genomic diversity or multiplicity of infection (MOI), expressed as the variant of inbreeding coefficient F_WS, was calculated at an individual level using the R package moimix (v0.0.2.9001, https://github.com/bahlolab/moimix). The F_WS metric expresses the probability of the heterozygosity of parasites within an individual against the heterozygosity within a parasite population. Samples with an F_WS score of ≥ 0.95 are highly considered to be monoclonal, whereas samples with an F_WS <0.95 suggests mixed strain infections, and therefore, a poorly defined population sub-structure. In all selection analyses, only monoclonal isolates (F_WS ≥0.95) were considered.

IBD and directional selection analyses

Relatedness between all samples was explored in a pairwise manner through identity by descent (IBD) analysis, conducted by the hmmIBD⁶⁶ package (v2.0.4, default parameters). hmmIBD implements a hidden Markov model for inference of pairwise IBD between haploid genotypes, enabling detection of DNA segments with shared ancestry. Only populations of monoclonal samples (F_WS ≥ 0.95) at both a country and regional level with > 10 isolates and biallelic SNPs with a MAF > 1% were used for analysis. An additional filtering step to the binary SNP matrix replaced all missing calls to a reference call and all mixed calls to alternative calls. Using sliding windows of 10 kb, IBD was cumulatively calculated and plotted by chromosomal location in R (v4.2.3) using in-house scripts.

To detect signals of recent positive selection, all samples were screened in a pairwise manner at both a country and regional level using the R package rehh⁶⁷ on SNPs with MAF > 1%. We calculated the integrated haplotype score (iHS)⁶⁸ for within-population selection and the cross-population ratio of extended haplotype homozygosity (EHH) expressed as Rsb and XP-EHH, for between-population selection⁶⁹. For iHS analysis, we describe the ancestral and derived alleles as the reference and non-reference alleles respectively. A positive iHS score suggests that the reference allele has undergone selection, whereas a negative iHS score suggests selection of a non-reference allele. Critical loci were identified using 100 kb sliding windows, which included at least 5 SNPs with a p value < 1 x 10^− 4 for iHS and < 1 x 10^− 5 for Rsb and XPEHH. XP-EHH detects selective sweeps in which a selected allele has approached or achieved fixation in one population while remaining polymorphic in the other population by comparing the lengths of the haplotypes associated with the selected allele in both populations⁷⁰. A positive XP-EHH indicates selection occurring in population 1, whereas a negative XP-EHH indicates selection occurring in population 2. The Rsb metric is the ratio of EHH between two populations, normalised to 1.

Haplotype network estimation

The aligned FASTA files for pvmdr1 and its downstream locus were used to estimate haplotype networks using the Pegas (v0.11) package in R, along with nucleotide and haplotype diversity statistics⁷¹.

DATA AVAILABILITY

All raw data used in this work is publicly available from the European Nucleotide Archive (see Table S1 for accession numbers and sample provenance).

FUNDING

G.C.N-J is funded by a BBSRC LIDo PhD studentship. T.G.C. and S.C. are funded by UKRI MRC (IAA2129, MR/R026297/1, and MR/X005895/1) and EPSRC (EP/Y018842/1), and Wellcome iTPA Translational Accelerator Award (214227/Z/18/Z) grants. C.R.F.M was supported by the São Paulo Research Foundation-FAPESP (2020/06747-4 and 2022/13150-0) and the National Council for Scientific and Technological Development-CNPq (302917/2019-5). J.G.D was supported by fellowships from FAPESP (2019/12068-5 and 2022/02771-3). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

AUTHOR CONTRIBUTIONS

T.G.C and S.C. conceived and designed the study. G.C.N-J performed all bioinformatic analyses and interpreted results under the supervision of T.G.C and S.C. J.E.P and E.M contributed bioinformatic tools and insights. J.G.D, R.L.D.M, C.R.F.M, D.N, F.N, and C.J.S contributed samples and sequence data. S.C sequenced samples. G.C.N-J wrote the first draft of the manuscript with supervision from T.G.C and S.C. All authors provided comments on the versions of the manuscript, approving of the final manuscript.

ETHICS DECLARATIONS

The authors declare no competing interests.

World Health Organisation World Malaria Report 2023
Wells TNC, Burrows JN, Baird JK (2010) Targeting the hypnozoite reservoir of Plasmodium vivax: the hidden obstacle to malaria elimination. Trends Parasitol 26:145–151
Obaldia N et al (2018) Bone Marrow Is a Major Parasite Reservoir in Plasmodium vivax Infection. mBio 9:e00625–e00618
Angrisano F, Robinson LJ (2022) Plasmodium vivax – How hidden reservoirs hinder global malaria elimination. Parasitol Int 87:102526
Moreira CM, Abo-Shehada M, Price RN, Drakeley C (2015) J. A systematic review of sub-microscopic Plasmodium vivax infection. Malar J 14:360
Douglas NM et al (2013) Gametocyte Dynamics and the Role of Drugs in Reducing the Transmission Potential of Plasmodium vivax. J Infect Dis 208:801–812
Bermúdez M, Moreno-Pérez DA, Arévalo-Pinzón G, Curtidor H, Patarroyo MA (2018) Plasmodium vivax in vitro continuous culture: the spoke in the wheel. Malar J 17:301
Chu CS, White NJ (2016) Management of relapsing Plasmodium vivax malaria. Expert Rev Anti-infective Therapy 14:885–900
Price RN, Commons RJ, Battle KE, Thriemer K, Mendis K (2020) Plasmodium vivax in the Era of the Shrinking P. falciparum Map. Trends Parasitol 36:560–570
Brega S et al (2005) Identification of the Plasmodium vivax mdr-like gene (pvmdr1) and analysis of single-nucleotide polymorphisms among isolates from different areas of endemicity. J Infect Dis 191:272–277
Sá JM et al (2005) Plasmodium vivax: allele variants of the mdr1 gene do not associate with chloroquine resistance among isolates from Brazil, Papua, and monkey-adapted strains. Exp Parasitol 109:256–259
Suwanarusk R et al (2007) Chloroquine Resistant Plasmodium vivax: In Vitro Characterisation and Association with Molecular Polymorphisms. PLoS ONE 2:e1089
Barnadas C et al (2008) Plasmodium vivax Resistance to Chloroquine in Madagascar: Clinical Efficacy and Polymorphisms in pvmdr1 and pvcrt-o Genes. Antimicrob Agents Chemother 52:4233–4240
Orjuela-Sánchez P et al (2009) Analysis of Single-Nucleotide Polymorphisms in the crt-o and mdr1 Genes of Plasmodium vivax among Chloroquine-Resistant Isolates from the Brazilian Amazon Region. Antimicrob Agents Chemother 53:3561–3564
Mekonnen SK et al (2014) Return of chloroquine-sensitive Plasmodium falciparum parasites and emergence of chloroquine-resistant Plasmodium vivax in Ethiopia. Malar J 13:244
Faway E et al (2016) Plasmodium vivax multidrug resistance-1 gene polymorphism in French Guiana. Malar J 15:540
Li J et al (2020) Ex vivo susceptibilities of Plasmodium vivax isolates from the China-Myanmar border to antimalarial drugs and association with polymorphisms in Pvmdr1 and Pvcrt-o genes. PLoS Negl Trop Dis 14:e0008255
Spotin A et al (2020) Global assessment of genetic paradigms of Pvmdr1 mutations in chloroquine-resistant Plasmodium vivax isolates. Trans R Soc Trop Med Hyg 114:339–345
Benavente ED et al (2017) Genomic variation in Plasmodium vivax malaria reveals regions under selective pressure. PLoS ONE 12:e0177134
Benavente ED et al (2021) Distinctive genetic structure and selection patterns in Plasmodium vivax from South Asia and East Africa. Nat Commun 12:3160
Cheeseman IH et al (2012) A major genome region underlying artemisinin resistance in malaria. Science 336:79–82
Miotto O et al (2013) Multiple populations of artemisinin-resistant Plasmodium falciparum in Cambodia. Nat Genet 45:648–655
Takala-Harrison S et al (2013) Genetic loci associated with delayed clearance of Plasmodium falciparum following artemisinin treatment in Southeast Asia. Proc Natl Acad Sci U S A 110:240–245
Mu J et al (2010) Plasmodium falciparum genome-wide scans for positive selection, recombination hot spots and resistance to antimalarial drugs. Nat Genet 42:268–271
Amambua-Ngwa A et al (2023) Chloroquine resistance evolution in Plasmodium falciparum is mediated by the putative amino acid transporter AAT1. Nat Microbiol 8:1213–1226
Winter DJ et al (2015) Whole Genome Sequencing of Field Isolates Reveals Extensive Genetic Diversity in Plasmodium vivax from Colombia. PLoS Negl Trop Dis 9:e0004252
Rodrigues PT et al (2018) Human migration and the spread of malaria parasites to the New World. Sci Rep 8:1993
Buyon LE et al (2020) Population genomics of Plasmodium vivax in Panama to assess the risk of case importation on malaria elimination. PLoS Negl Trop Dis 14:e0008962
Slatkin M, ISOLATION BY DISTANCE, IN EQUILIBRIUM AND NON-EQUILIBRIUM POPULATIONS (1993) Evolution 47:264–279
Raj DK et al (2009) Disruption of a Plasmodium falciparum multidrug resistance-associated protein (PfMRP) alters its fitness and transport of antimalarial drugs and glutathione. J Biol Chem 284:7687–7696
Desai M et al (2016) Impact of Sulfadoxine-Pyrimethamine Resistance on Effectiveness of Intermittent Preventive Therapy for Malaria in Pregnancy at Clearing Infections and Preventing Low Birth Weight. Clin Infect Dis 62:323–333
Korsinczky M et al (2004) Sulfadoxine resistance in Plasmodium vivax is associated with a specific amino acid in dihydropteroate synthase at the putative sulfadoxine-binding site. Antimicrob Agents Chemother 48
Marfurt J et al (2021) Longitudinal ex vivo and molecular trends of chloroquine and piperaquine activity against Plasmodium falciparum and P. vivax before and after introduction of artemisinin-based combination therapy in Papua, Indonesia. Int J Parasitology: Drugs Drug Resist 17:46–56
Koepfli C et al (2015) Plasmodium vivax Diversity and Population Structure across Four Continents. PLoS Negl Trop Dis 9:e0003872
Hupalo DN et al (2016) Population genomics studies identify signatures of global dispersal and drug resistance in Plasmodium vivax. Nat Genet 48:953–958
Rougeron V et al (2020) Human Plasmodium vivax diversity, population structure and evolutionary origin. PLoS Negl Trop Dis 14:e0008072
Yalcindag E et al (2012) Multiple independent introductions of Plasmodium falciparum in South America. Proceedings of the National Academy of Sciences 109, 511–516
Gartner V et al (2024) Genomic insights into Plasmodium vivax population structure and diversity in central Africa. Malar J 23:27
Culleton R et al (2011) The Origins of African Plasmodium vivax; Insights from Mitochondrial Genome Sequencing. PLoS ONE 6:e29137
Liu Y et al (2023) Retrospective analysis of Plasmodium vivax genomes from a pre-elimination China inland population in the 2010s. Front Microbiol 14
Flegg JA et al (2022) Spatiotemporal spread of Plasmodium falciparum mutations for resistance to sulfadoxine-pyrimethamine across Africa, 1990–2020. PLoS Comput Biol 18:e1010317
Kebede AM et al (2023) Genomic analysis of Plasmodium vivax describes patterns of connectivity and putative drivers of adaptation in Ethiopia. Sci Rep 13:20788
Ward KE et al (2023) Integrative Genetic Manipulation of Plasmodium cynomolgi Reveals Multidrug Resistance-1 Y976F Associated With Increased In Vitro Susceptibility to Mefloquine. J Infect Dis 227:1121–1126
Auburn S et al (2018) Genomic analysis of a pre-elimination Malaysian Plasmodium vivax population reveals selective pressures and changing transmission dynamics. Nat Commun 9:2585
Cheong F-W, Dzul S, Fong M-Y, Lau Y-L, Ponnampalavanar S (2020) Plasmodium vivax drug resistance markers: Genetic polymorphisms and mutation patterns in isolates from Malaysia. Acta Trop 206:105454
Awab GR et al (2010) Dihydroartemisinin-piperaquine versus chloroquine to treat vivax malaria in Afghanistan: an open randomized, non-inferiority, trial. Malar J 9:105
Schousboe ML et al (2015) Multiple Origins of Mutations in the mdr1 Gene—A Putative Marker of Chloroquine Resistance in P. vivax. PLoS Negl Trop Dis 9:e0004196
Thanh PV et al (2015) Confirmed Plasmodium vivax Resistance to Chloroquine in Central Vietnam. Antimicrob Agents Chemother 59:7411–7419
Pearson RD et al (2016) Genomic analysis of local variation and recent evolution in Plasmodium vivax. Nat Genet 48:959–964
Kenangalem E et al (2019) Malaria morbidity and mortality following introduction of a universal policy of artemisinin-based treatment for malaria in Papua, Indonesia: A longitudinal surveillance study. PLoS Med 16:e1002815
Thuan PD et al (2016) A Randomized Comparison of Chloroquine versus Dihydroartemisinin–Piperaquine for the Treatment of Plasmodium vivax Infection in Vietnam. Am J Trop Med Hyg 94:879–885
Dharia NV et al (2010) Whole-genome sequencing and microarray analysis of ex vivo Plasmodium vivax reveal selective pressure on putative drug resistance genes. Proc Natl Acad Sci U S A 107:20045–20050
Flannery EL et al (2015) Next-Generation Sequencing of Plasmodium vivax Patient Samples Shows Evidence of Direct Evolution in Drug-Resistance Genes. ACS Infect Dis 1:367–379
Mok S et al (2023) Mapping the genomic landscape of multidrug resistance in Plasmodium falciparum and its impact on parasite fitness. Sci Adv
Mohring F et al (2019) Rapid and iterative genome editing in the malaria parasite Plasmodium knowlesi provides new tools for P. vivax research. Elife 8:e45829
MalariaGEN et al (2022) An open dataset of Plasmodium vivax genome variation in 1,895 worldwide samples. Wellcome Open Res 7:136
Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30:2114–2120
Auburn S et al (2016) A new Plasmodium vivax reference sequence with improved assembly of the subtelomeres reveals an abundance of pir genes. Wellcome Open Res 1:4
Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25:1754–1760
Li H (2011) A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics 27:2987–2993
McKenna A et al (2010) The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res 20:1297–1303
Cingolani P et al (2012) A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff. Fly (Austin) 6:80–92
Purcell S et al (2007) PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. Am J Hum Genet 81:559–575
Letunic I, Bork P (2021) Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic Acids Res 49:W293–W296
Alexander DH, Lange K (2011) Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinformatics 12:246
Schaffner SF, Taylor AR, Wong W, Wirth DF, Neafsey DE (2018) hmmIBD: software to infer pairwise identity by descent between haploid genotypes. Malar J 17:196
Gautier M, Vitalis R (2012) rehh: an R package to detect footprints of selection in genome-wide SNP data from haplotype structure. Bioinformatics 28:1176–1177
Voight BF, Kudaravalli S, Wen X, Pritchard JK (2006) A Map of Recent Positive Selection in the Human Genome. PLoS Biol 4:e72
Tang K, Thornton KR, Stoneking M (2007) A New Approach for Using Genome Scans to Detect Recent Positive Selection in the Human Genome. PLoS Biol 5:e171
Sabeti PC et al (2007) Genome-wide detection and characterization of positive selection in human populations. Nature 449:913–918
Paradis E (2010) pegas: an R package for population genetics with an integrated–modular approach. Bioinformatics 26:419–420

There is NO Competing Interest.

Download PDF

Version 1

posted

You are reading this latest preprint version

Genome-wide association study of global Plasmodium vivax populations provides insights into the evolution of drug resistance

Status:

Version 1

Abstract

Figures

INTRODUCTION

RESULTS

Pairwise relatedness supports model of isolation by distance

Evidence of selective sweeps in candidate drug resistance loci

Temporal trends in markers under selection in Indonesian isolates

DISCUSSION

METHODS

Sequence data and raw reads processing

Variant calling and quality control

Population genetics analyses

IBD and directional selection analyses

Haplotype network estimation

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1