Evolution of the Plasmodium vivax multidrug resistance 1 gene in the Greater Mekong Subregion during malaria elimination
Parasites & Vectors volume 13, Article number: 67 (2020)
The malaria elimination plan of the Greater Mekong Subregion (GMS) is jeopardized by the increasing number of Plasmodium vivax infections and emergence of parasite strains with reduced susceptibility to the frontline drug treatment chloroquine/primaquine. This study aimed to determine the evolution of the P. vivax multidrug resistance 1 (Pvmdr1) gene in P. vivax parasites isolated from the China–Myanmar border area during the major phase of elimination.
Clinical isolates were collected from 275 P. vivax patients in 2008, 2012–2013 and 2015 in the China–Myanmar border area and from 55 patients in central China. Comparison was made with parasites from three border regions of Thailand.
Overall, genetic diversity of the Pvmdr1 was relatively high in all border regions, and over the seven years in the China–Myanmar border, though slight temporal fluctuation was observed. Single nucleotide polymorphisms previously implicated in reduced chloroquine sensitivity were detected. In particular, M908L approached fixation in the China–Myanmar border area. The Y976F mutation sharply decreased from 18.5% in 2008 to 1.5% in 2012–2013 and disappeared in 2015, whereas F1076L steadily increased from 33.3% in 2008 to 77.8% in 2015. While neutrality tests suggested the action of purifying selection on the pvmdr1 gene, several likelihood-based algorithms detected positive as well as purifying selections operating on specific amino acids including M908L, T958M and F1076L. Fixation and selection of the nonsynonymous mutations are differently distributed across the three border regions and central China. Comparison with the global P. vivax populations clearly indicated clustering of haplotypes according to geographic locations. It is noteworthy that the temperate-zone parasites from central China were completely separated from the parasites from other parts of the GMS.
This study showed that P. vivax populations in the China–Myanmar border has experienced major changes in the Pvmdr1 residues proposed to be associated with chloroquine resistance, suggesting that drug selection may play an important role in the evolution of this gene in the parasite populations.
Worldwide, around 100 million cases of Plasmodium vivax infections are registered annually with majority of them occurring in the Asian Pacific Region [1, 2]. Most countries in Southeast Asia are making steady progress in reducing the malaria burden; the six countries in the Greater Mekong Subregion (GMS) have set their goals to eliminate malaria by 2030 [1, 3]. This elimination plan is, however, challenged by the difficulties to eliminate P. vivax because of its several biological features such as very low blood parasitemia that is often missed by conventional detection methods, and formation of hypnozoites in the liver of an infected individual that are responsible for subsequent relapses . Despite increased control efforts in the GMS, P. vivax transmission along international borders remains high [5, 6].
Understanding the genetics of drug resistance in P. vivax is important for implementing an effective chemotherapeutic strategy and monitoring the progress of elimination . Whilst the mechanisms of drug resistance in Plasmodium falciparum are much better understood, those in P. vivax are largely unknown. Chloroquine (CQ) has been withdrawn from treating P. falciparum malaria in most endemic countries due to widespread resistance to this drug [8, 9], but CQ-primaquine (PQ) combination is still the first-line treatment for P. vivax infections in most endemic countries [10, 11]. Unfortunately, there is an increased number of reports of reduced susceptibility of P. vivax parasites to CQ from malaria-endemic areas, including the GMS countries [2, 8, 12,13,14,15,16,17,18,19,20,21]. Despite this, there is still a lack of a confirmed marker(s) for CQ resistance in P. vivax. Several studies have indicated that mutations in the multidrug resistant 1 gene (Pvmdr1) may be used as markers for CQ resistance surveillance [22, 23]. In vitro drug susceptibility assays identified an association between higher copy numbers of the Pvmdr1 and increased CQ IC50 values [24, 25], although the cut-off IC50 value for CQ resistance is uncertain. More recently, a connection has been made between the copy number of the Pvmdr1 harboring the Y976F/F1076L mutations and treatment failure in severe P. vivax malaria cases [26, 27]. In addition, the M908L and T958M mutations were shown to be associated with reduced in vitro CQ sensitivity . However, some studies failed to detect a link between the Pvmdr1 mutations and reduced CQ sensitivity, raising doubts on the suitability of the Pvmdr1 mutations as markers for CQ resistance [29, 30].
Population genomics studies revealed great diversity of the P. vivax parasites as compared to P. falciparum [31, 32], indicating more stable populations. Moreover, signals of natural selection have been detected in P. vivax, highlighting the ability of P. vivax to evolve in response to antimalarial drug pressure and changing environments in the human host as well as in the mosquito vector . For instance, dihydropteroate synthase and dihydrofolate reductase genes that are associated with resistance to antifolate drugs were found to be selected in P. vivax . In the GMS, P. vivax parasites were found to exhibit high levels of genetic diversity in Thailand , southern China, and Myanmar . In this study, we focused on the genetic diversity of the Pvmdr1 gene in the vivax-endemic area along the China–Myanmar border, hoping to understand the evolution of the parasites amid the falling CQ treatment efficacy  and increased proportions of vivax malaria in most areas of the GMS .
Study sites and samples
Clinical P. vivax samples were collected from 330 patients with acute P. vivax malaria attending different malaria clinics. Among them, 39 and 16 were from Anhui Province of central China in 2004 and 2006–2008, respectively. For the longitudinal samples from the China–Myanmar border, 27, 129 and 119 samples were collected in 2008, 2012–2013 and 2015, respectively, giving a total of 275 samples from this border region. Finger-prick blood samples of microscopy-confirmed P. vivax cases were spotted onto Whatman 3M filter papers.
Sequencing of the Pvmdr1 gene
Genomic DNA was extracted from dried blood spots on filter paper using the QIAamp DNA Mini kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. Genotyping of two polymorphic genes (msp3α and msp3β) by PCR/RFLP was done to distinguish single from mixed-strain infections [36, 37]. For PCR amplification of the Pvmdr1 gene, primary PCR was performed using primers P1F and P1R and two fragments were amplified by semi-nested PCR with primer pairs P1F × N-PR and N-PF × P1R, respectively (Additional file 1: Table S1). PCR was performed using the Advantage 2 polymerase mix (Takara Bio, Mountain View, USA) and PCR products were sequenced in both directions using the Sanger method on an ABI DNA analyzer. The Pvmdr1 sequences were assembled and edited using DNAStar (Lasergene, Madison, USA). The Pvmdr1 sequences generated from this study are available in GenBank with the accession numbers: MN891946–MN891972; MN891973–MN892091; MN892092–MN892220; MN892221–MN892236; and MN892237–MN892275. In addition, 98 Pvmdr1 sequences from parasites collected in western (Tak and Kanchanaburi provinces) and eastern Thailand (Ubon Rachathani Province) were also used for analysis . All sequences were aligned with the reference Pvmdr1 sequence from the Salvador I strain (PVX_080100) using Clustal Muscle 3.8  incorporated in the MEGA7 software .
Assessment of genetic diversity
All the Pvmdr1 sequences were scanned for the presence of single nucleotide polymorphisms (SNPs). The genetic diversity of the Pvmdr1 gene was assessed using DnaSP software v6.10 . Haplotype diversity (Hd) of the Pvmdr1 gene was estimated based on the number and frequency of the haplotypes, while nucleotide diversity was measured using two parameters: π, the average number of pairwise nucleotide differences per site  and θw, the number of segregating sites.
Tests for detecting selection
To determine whether natural selection played a role in the evolution of Pvmdr1, we first performed a series of frequency-based tests including Tajima’s D test , Fu and Li’s F test , and Fu and Li’s D test  using the DnaSP v6.10 software. The Tajima’s D statistic calculates the normalized differences between the two measures of nucleotide diversity θw and π . Both the Fu and Li’s D and F statistics rely on the difference between the number of polymorphic sites in external branches (polymorphisms unique to an extant sequence) and the number of polymorphic sites in internal phylogenetic branches (polymorphisms shared by extant sequences) . For all statistical analyses, a P-value of ≤ 0.05 was considered significant. We also used McDonald–Kreitman (MK) test to examine departure from neutrality using Plasmodium knowlesi mdr1 sequence as the outgroup . The MK test compares the ratio of nonsynonymous to synonymous polymorphism within a species (Pn/Ps) and the ratio of nonsynonymous to synonymous substitutions between closely related species (dN/dS). Fisher’s exact test was used to assess statistical significance.
We then determined the nucleotide substitutions and the ratio of nonsynonymous (dN) to synonymous (dS) substitutions per site (dN/dS), using the Nei-Gojobori method  after Jukes-Cantor correction for multiple substitutions. Under the neutral model of evolution, dS is expected to be equal to dN. An excess of nonsynonymous substitutions (dN > dS) can be interpreted as positive selection , indicating that replacement substitutions increase fitness of the parasite, whereas a rarity of replacement changes (dN < dS) specifies that purifying selection might be working to remove such substitutions from the gene pool . Statistical significance of the difference was estimated using the codon-based Z-test of selection in MEGA7 .
Finally, since selection is often directed at a few amino acids of a gene and sometimes can be camouflaged by purifying selection also acting on the gene , we conducted maximum likelihood tests in the HyPhy package implemented in the Data Monkey Web Server  to determine the specific amino acids targeted by selection . Significant recombination events were tested in the DnaSP program and by genetic algorithm for recombination detection  incorporated in the Data Monkey Web Server before running the tests of selection.
Prediction of possible effects of the Pvmdr1 mutations on protein function
To predict if any of the Pvmdr1 mutations might affect the protein structure and function, we mapped these residues on a modeled 3D structure using the Sal I reference sequence. The homology model of PvMDR1 was built based on the structures of the multidrug transporter P-glycoprotein (Pgp) from Caenorhabditis elegans (4F4C) and mouse (4M1M and 3G61) using the multiple threading alignment in I-TASSER . A confidence score (C-score) for estimating the quality of predicted models by I-TASSER was calculated . Web-based software PROVEAN and SIFT (Sorting Intolerant from Tolerant) were utilized to predict the effect of amino acid mutations in PvMDR1 . Mutations predicted to be deleterious according to both the software were mapped on the predicted 3D structure of PvMDR1.
Population differentiation and linkage disequilibrium (LD)
To determine the genetic interrelationships among all parasite isolates, a phylogenetic tree was constructed using the Maximum Likelihood algorithm with 1000 bootstraps as implemented in MEGA7. The Sal I reference strain was represented as the wild type. In addition to the 275 Pvmdr1 sequences obtained from this study, a total of 180 complete or nearly complete Pvmdr1 sequences retrieved from GenBank and PlasmoDB (plasmodb.org) representing parasite isolates from 11 countries were also analyzed: 6 from China; 98 from Thailand; 5 from Papua New Guinea (PNG); 7 from Madagascar; 14 from Mexico; 20 from Colombia; 24 from Peru; 3 from Brazil; and one each from India, North Korea, and Mauritania. Each sequence was trimmed to remove low-quality segments, yielding 4137 bp of the 4395 bp Pvmdr1 open reading frame. To estimate the proportion of genetic variance of the Pvmdr1 gene due to population subdivision, Wright’s fixation index of inter-population variance in allele frequencies (FST) was calculated. Pairwise linkage LD was used to determine the degree of random association between different mutations within this gene. The correlation coefficient (R2) between paired alleles was estimated using DnaSP v6.10 and the significance of each association was determined using the Fisher and Chi-square tests after Bonferroni correction.
Haplotype network analysis
To visualize the distribution of the Pvmdr1 polymorphisms across different P. vivax populations, haplotypes were constructed from nonsynonymous SNPs that were observed in more than two isolates. A minimum spanning tree was drawn using the median-joining algorithm in the PHYLOViZ software (https://www.phyloviz.net).
Genetic diversity of the Pvmdr1 gene
Mutations in the Pvmdr1 gene have been associated with CQ resistance in P. vivax in some endemic areas. In light of the deteriorating CQ clinical efficacy for treating vivax malaria at the China–Myanmar border , we followed the evolution of the Pvmdr1 gene in parasite populations from this region over a seven-year period and sequenced the full-length Pvmdr1 gene in 275 P. vivax clinical samples. To put this study in context with parasites from other regions in the GMS, we also analyzed 98 Pvmdr1 sequences from the western and eastern borders of Thailand . For the longitudinal P. vivax clinical samples from the China–Myanmar border, 27 isolates collected in 2008 harbored 22 SNPs, of which 20 are nonsynonymous. In the 129 samples collected in 2012–2013, eight synonymous and 15 nonsynonymous SNPs were identified. In the 119 samples collected in 2015, 34 synonymous and 59 nonsynonymous SNPs were found. For the 55 temperate-zone P. vivax isolates collected from central China in 2004–2008, 24 SNPs were found, of which 17 were nonsynonymous. Several of the nonsynonymous mutations in the Pvmdr1 gene had allele frequencies of at least 5% (Table 1); eight were common in the China–Myanmar border populations from the three time points: S513R, G698S, L845F, A861E, M908L, T958M, F1076L and K1393N (Table 1). Among them, G698S, M908L, and T958M reached or almost reached fixation (97–100%). Six of the eight mutations were also detected in the central China parasite population with allele frequencies of ≥ 5% (Table 1). For the three SNPs (T958M, Y976F and F1076L) proposed to be associated with CQ resistance [23, 52], F1076L was fixed in the central China population. In the China–Myanmar border parasite populations, the prevalence of F1076L continually increased over time, from 33.3% in 2008 to 41.7% in 2012–2013 and 77.8% in 2015. In contrast, the Y976F mutation was not present in the central China parasite population, and its frequency in the border parasite populations was moderately high at 18.5% in 2008, but sharply decreased to 1.5% in 2012–2013 and was completely absent in the 2015 samples (Table 1).
Overall, the genetic diversity of the Pvmdr1 gene in the China–Myanmar border parasite populations was relatively high (π = 0.0009–0.0012) with slight fluctuation over the years. Similarly, haplotype diversity was also high in the China–Myanmar border parasite populations: 21, 33 and 75 haplotypes were identified in the 2008, 2012–2013 and 2015 samples, respectively (Table 2). Compared with the China–Myanmar border parasite populations, the genetic diversity of the Pvmdr1 gene in parasites from central China was much lower (π = 0.0006). Yet, haplotype diversity of these temperate-zone parasites was high (0.914 ± 0.026). Except for the population of 2012–2013, the π value was lower than the θw value for the rest of populations, suggesting that most SNPs in the Pvmdr1 gene were rare alleles.
Compared with P. vivax parasite populations from other parts of the GMS, haplotype diversity in the China–Myanmar border in 2015 was similarly high (0.971 ± 0.008) as that from the Thailand–Myanmar border (0.974 ± 0.009) (Table 2). For all parasite populations from the three border areas (China–Myanmar, Thailand–Myanmar, and Thailand–Cambodia) in this study, the π value was lower than the θw value (Table 2), suggesting the prevalence of rare alleles as mentioned above.
Mutations within the putative 3D model of PvMDR1
PvMDR1 is a member of the ATP-binding cassette (ABC) protein superfamily with two symmetric domains. Each domain has a transmembrane domain (TMD), consisting of three external loops and two internal helices that link six TMDs followed by a nucleotide binding domain (NBD) . According to the protein alignment and domain mapping analyses, the two TMDs contain 5–6 transmembrane helices (at amino acids 62–84, 99–121, 171–193, 197–216, 281–303, 323–345, 825–847, 867–889, 940–962, 966–985 and 1062–1084), while the two NBDs, also referred to as the AAA domains, are located at amino acids 410–662 and 1191–1433, respectively (Fig. 1). The predicted domains in PvMDR1 show high sequence homology to the corresponding PfMDR1 functional domains . For the 71 mutations reported in this study, 38 are predicted to be deleterious by at least one of the prediction programs (Additional file 2: Table S2), whereas 19 are predicted to be detrimental according to both Provean and SIFT analysis (Additional file 3: Table S3). To predict the effect of these 19 mutations on protein structure, we mapped the mutated residues on the predicted tertiary structure. The homology model of PvMDR1 built using I-TASSER aligned well with the C. elegans multidrug transporter P-glycoprotein. Except I595, 18 of the 19 amino acids are conserved in PfMDR1 protein sequence, indicating functional conservation and significance.
Of these 19 predicted deleterious mutations, V324G is located in the first TMD, whereas three (L845F, I951K, and V959G) are in the second TMD (Fig. 1). Notably, I951K represents a drastic change from a hydrophobic to a positively charged residue, which may disrupt the integrity of the TMD. Five mutations (Y348D, Y359D, E911K, D932N and P1177T) are located in the predicted inter-domain regions. Interestingly, most of the mutations with predicted adverse effects lie in the first NBD/AAA domain (K456T, L470H, V562G, A593T, I595F, L610F, D611K, V618G, N623I), whereas only one mutation (K1393N) is present in the second AAA domain. All these mutations can potentially affect the structural integrity of the protein by altering the charges, hydrophobicity, or size of the amino acids. In addition, the residues corresponding to L470, L610, and D611 in PfMDR1 are predicted to be involved in NDB dimerization , and mutations at these positions may hinder dimerization of the protein. Insights into the ATP-binding pockets can help decipher if any of these residues are involved in direct binding to ATP or some antimalarial drugs, allowing the determination of their potential role in transport and resistance.
Departure from neutrality
Frequency-based neutrality tests were used to assess the evolution of the Pvmdr1 gene. All the neutrality tests yielded negative values for all the sample sets except the 2012–2013 parasite population from the China–Myanmar border (Table 2). However, only the 2015 parasites from the China–Myanmar border significantly deviated from neutral with an excess of low-frequency polymorphisms, suggesting that the Pvmdr1 gene in China–Myanmar border possibly experienced either a directional selection or population expansion during the seven years.
The dN-dS statistic generated by the Z-test was negative in all the populations, suggestive of purifying selection on Pvmdr1, albeit it was not statistically significant (Table 3). The MK test also indicated that the Pvmdr1 gene was evolving under purifying selection in the central China and the China–Myanmar border populations. The number of recombination events was found to be very high in the 2015 border parasite population (Rm = 15), corroborating the high diversity observed in this population (Table 2).
The likelihood-based algorithms (SLAC , FEL , and FUBAR  implemented in the Datamonkey webserver ) all revealed positive as well as purifying selection at specific codons (Table 4). The number of sites negatively selected increased with the years in the China–Myanmar border populations. Among the mutations presumptively associated with CQ resistance, F1076L seemed to be positively selected in the 2008 and 2012–2013 China–Myanmar border populations, whereas M908L was positively selected in the 2015 population (Table 4). In comparison, none of these mutations were selected in the central China population (Table 4). It is noteworthy that of the two mutations M908L and T958M associated with reduced in vitro CQ sensitivity  and also highly prevalent in Thailand , M908L was positively selected in the China–Myanmar border parasites but not in the Thailand–Myanmar or the Thailand–Cambodia border parasite populations (Table 4). Furthermore, most of mutations with frequencies of ≥ 5% such as A861E, L845F, and K1393N (Table 1) were also selected in the China–Myanmar parasites (Table 4). In the Thailand–Cambodia border population, only F1076L was found to be positively selected (Table 4), whereas in the Thailand–Myanmar border populations, S513R, G698S, A861E, F1076L and K1393N were positively selected (Table 4).
The Hudson and Kaplanʼs lower bound on the minimal number of recombination events in an infinite site model computed with DnaSP revealed 4, 6, 15, and 3 minimum recombination events in the Pvmdr1 gene from the China–Myanmar border populations in 2008, 2012–2013, and 2015 and the central China population, respectively (Table 2). In the Thailand–Myanmar and Thailand–Cambodia border regions, the minimum recombination events were 5 and 4, respectively (Table 2). On the other hand, analysis by using GARD (genetic algorithm for recombination detection) , a model-based approach that searches for putative breakpoints delimiting sequence regions having distinct phylogenies, found evidence of a breakpoint only in the Thailand–Myanmar border parasite population. The low number of recombination events in Thailand–Myanmar and Thailand–Cambodia border populations might be due to a recent expansion in effective parasite population. Focusing on the China–Myanmar border parasites, our study detected an increasing number of minimum recombination events from 2008 to 2015, highlighting a possible reduction of P. vivax effective population size during these years.
Global distribution of the Pvmdr1 haplotypes and geographical differentiation
Based on the Pvmdr1 amino acid sequences, a total of 188 haplotypes were found in 510 parasites isolates from the world (Fig. 2). Most geographical regions except Madagascar and Mauritania had more than one predominant haplotype. There were significant differences in the number of haplotypes and prevalence of individual haplotypes among all the countries considered. The China–Myanmar border parasites had the highest number of 79 unique haplotypes, followed by parasites from Thailand (46 haplotypes). A minimum spanning network clearly showed geographical clustering of the haplotypes (Fig. 2). A continental, followed by a country-wise and then region-wise difference were observed. In this regard, the majorities of the Asian and African parasites are separated from the New World parasites (Fig. 3). Also, parasites from the Thailand–Myanmar border were distinctive from those from the Thailand–Cambodia border (Fig. 2). Similarly, temperate-zone P. vivax parasites from central China were completely separated from those of the China–Myanmar border. Only 20 haplotypes out of the 188 were shared across all populations in the world (Fig. 2). Besides, some long branches are present within the torso of the network (Fig. 2), highlighting a local genetic difference of the parasites. Significant sharing of the haplotypes was detected between physically connected Asian countries, suggesting potential genetic exchanges in the past between these populations. Interestingly, the African cluster was linked to one of the predominant haplotypes of the Asian cluster.
Phylogenetic analysis using the maximum likelihood method further corroborated the high-degree genetic differentiation of parasite populations from the three GMS border regions, as well as parasites from the rest of world (Fig. 3). The result showed a clear population substructure, in particular, within the China–Myanmar populations (Fig. 3).
Population differentiation examined through the estimation of FST, the Wright’s fixation index of inter-population variance in allele frequencies, revealed large degrees of variation in population differentiation between countries (FST = 0.099–0.77), not considering India, Mauritania and North Korea due to very limited sample size. Overall, the FST estimate of the worldwide populations was 0.36, indicating that about 36% of the variation was apportioned between parasite populations. Great genetic differentiation was not only denoted between countries or regions, but also within some countries/regions such as the China–Myanmar border. Although there was a low degree of difference between parasites from the China–Myanmar border, the difference varied over the years (Table 5), with FST ranging from 0.025 to 0.105, suggesting extensive genetic inheritance. High degrees of genetic difference were detected among the South American countries (FST = 0.13–0.77, Table 5). Significant difference was also denoted between the Thailand–Myanmar border and the China–Myanmar border populations as well as between the Thailand–Myanmar border and the Thailand–Cambodia border populations (Table 5). Taken together, the result confirmed observations made in the phylogenetic and haplotype network analyses and is concordant with the hypothesis of the presence of genetic substructure.
High LD was detected in the 2015 parasite population from the China–Myanmar border as compared to 2008 and 2012–2013 (Fig. 4), suggesting an effective population size reduction with the years. In contrast, limited LD was detected in central China, the Thailand–Cambodia and Thailand–Myanmar border populations (Fig. 4), suggesting an effective population expansion and isolation.
The GMS countries are aiming to eliminate malaria by 2030, and intensified malaria control efforts have led to a continuous decline of malaria incidence . Malaria transmission is concentrated along international borders as these places are often mountainous, hard to reach, and inhabited by economically deprived populations [5, 6]. Some border areas had increased migrant populations as a result of political instability and military conflicts, leading to substantial changes in malaria epidemiology . The Kachin civil wars resulted in the establishment of camps for internally displaced people in the China–Myanmar border area to host human population migrated from other malaria endemic areas. Human migration, together with the poor public health infrastructure, has led to a rising trend of malaria incidence and even malaria outbreaks in recent years . Furthermore, cross-border human migration also poses another threat of malaria introduction to neighboring countries [6, 61]. Therefore, the China–Myanmar border region represents an interesting scenario for characterizing the parasite population structure and understanding their evolution during the course of malaria elimination.
This study focused on the molecular evolution of the Pvmdr1 gene, a potential marker for CQ resistance, in order to address the deteriorating CQ efficacy at the China–Myanmar border areas . The results showed an increase of the Pvmdr1 genetic diversity at the China–Myanmar border over the years in spite of the intensified control measures in place. The genetic diversity of the Pvmdr1 gene was high in most GMS border areas, though it fluctuated over the seven-year study period. As the frontline treatment for P. vivax malaria is CQ-PQ throughout the GMS, the genetic diversity of different parasite populations likely reflected the local differences in malaria epidemiology. Since most malaria endemic areas in the GMS also have sympatric P. falciparum transmission, different ACTs, especially the quinoline partner drugs, could have also exerted divergent selection pressures on the Pvmdr1 gene. Moreover, this high genetic diversity might exemplify the impact of parasite introduction as found in preceding studies . Three mutations (T958M, Y976F and F1076L) have been associated with P. vivax CQ resistance [23, 28, 30, 52]. The T958M mutation was fixed or approached fixation (98.3–100%) in the P. vivax populations from the China–Myanmar border as well as from the Thailand–Myanmar and Thailand–Cambodia borders, thus it is unlikely responsible for the reduced sensitivity of the parasites to CQ. However, the F1076L mutation in the longitudinal samples from the China–Myanmar border had an increasing trend in prevalence; its frequency reached 77.8% in 2015, almost doubled from that (41.7%) in 2012–2013. Interestingly, this mutation had a moderate prevalence in the Thailand–Myanmar border area (~30–62%) over the period of 2008–2016 and in Ubon Ratchathani (28%), but it almost reached fixation in Chanthaburi Province at the Thailand–Cambodia border [33, 62]. Conversely, the prevalence of the Y976F mutation in the China–Myanmar border populations progressively decreased and was not detected in the 2015 samples. Similarly, the Y976F prevalence also showed a decreasing trend in the Thailand–Myanmar and Thailand–Cambodia border samples collected in 2008 and 2014 [33, 62]. It is noteworthy that the Y976F mutation was associated with a low-level reduction of in vitro susceptibility to CQ [24, 25, 52]. It would be interesting to determine whether the reverse trends of the Y976F and F1076L mutations are associated with the decline of CQ efficacy in the GMS.
All the neutrality tests yielded negative values in most P. vivax populations, suggestive of the occurrence of rare alleles and parasites experiencing a directional selection or population expansion. Only the 2012–2013 parasites from the China–Myanmar border had positive values of those statistics, which may indicate a signature of balancing selection or population size decline [42, 43, 63, 64]. Inference of selection identified that the Pvmdr1 gene has evolved under purifying selection, reflecting Pvmdr1 as an essential gene. However, many tests that rely on differences between nonsynonymous and synonymous changes do not systematically take into account that positive selection often acts only on small regions of a gene product . In fact, zooming in specific regions of Pvmdr1 also identified individual codons to be under positive selection in the studied parasite populations. Whereas F1076L was the only position found to be positively selected in Thailand–Cambodia border, several loci, including G698S, M908L and F1076L were under positive selection in the Thailand–Myanmar border and China–Myanmar border populations. This finding further corroborated an earlier analysis of the publicly available P. vivax genomes gathered from various sources, which similarly revealed that T958M and M908L, F1076L, G698S and S513R were under directional selection . Of note, mapping the SNPs to the putative 3D model of PvMDR1 structure identified amino acid changes S513R, L845F, F1076L, K1393N and S1450L, all found to be under positive selection, could have impact on protein function [33, 52].
Both haplotype network and phylogenetic analysis revealed considerable clustering of the haplotypes relevant to the countries/continents of origin. This is intuitively understandable, as parasite populations in geographically separated continents or countries have been evolving under isolation, creating significant divergence among themselves. Under this scenario, parasites from Asia were more closely related among themselves than from those of the American origins, and parasites from the temperate region of central China were, to a large degree, separated from those of the tropical regions of the GMS. These comparisons also identified relatively little differentiation of parasite populations from the same geographical locations. Although the Pvmdr1 genes from the GMS displayed high diversity, there was extensive sharing of the major haplotypes among these border parasite populations, suggesting little differentiation of parasite populations within the GMS. This observation was further reinforced by the very low FST value obtained among these GMS populations (FST = 0.025–0.085). Studies on vaccine candidate genes such as PvAMA1 genes also showed high diversity and little differentiation of the P. vivax parasites from China–Myanmar border . Nevertheless, parasite populations from the GMS did fall into several distinctive clades, suggesting the presence of gene flow barriers or/and divergent selection on the Pvmdr1 protein. This is plausible, as intensified control efforts of the malaria elimination campaign may have led to separated pockets of transmission hotspots, and these isolated parasite populations may have been evolving independently as what has been observed for the P. falciparum populations in the China–Myanmar border region . This has also been the case of some genes such as the P. falciparum gamete surface protein gene Pfs48/45  and the P. vivax gamete surface protein gene Pvs48/45 .
A temporal increase in LD was denoted in parasite populations from the China–Myanmar border. Central China as well as Thailand–Myanmar and Thailand–Cambodia borders had limited LD. This might suggest high level of inbreeding  and a history of bottleneck risen by an effective population size decline on the China–Myanmar border , whereas in the other borders and central China, there might be an expansion of the effective population size. Interestingly, despite observation of high numbers of minimal recombination events in the history of the China–Myanmar P. vivax populations, no recombination breakpoint was found in all of them. This latter finding reinforces our hypothesis of reduced effective population size resulting in high-level inbreeding and consequently strong LD on the China–Myanmar border. Frequent inbreeding and recombination between parasite genotypes also play a role in contributing to high genetic diversity within populations .
Our study showed that the Pvmdr1 gene in P. vivax populations of the China–Myanmar border area have undergone a strong diversification process with evidence of purifying selection on the whole gene and positive selection on certain loci of the gene. Furthermore, there was a low level of genetic differentiation among the GMS parasite populations, suggesting extensive gene flow within the GMS. The increased diversity of P. vivax found parasites from the China–Myanmar border suggests parasite introduction, which might be associated with the migration of human population from other P. vivax endemic regions of Myanmar due to military conflicts. The increase in LD with the years indicated expansion of particular parasite genotypes associated with the recent outbreaks of P. vivax malaria in recent years. These findings emphasize once more that effective management of clinical vivax cases and monitoring of human migration are indispensable for malaria elimination in the GMS.
Availability of data and materials
The datasets supporting the conclusions of this article are available the in additional files.
Greater Mekong Subregion
nucleotide binding domain
Plasmodium vivax multidrug resistance 1
single nucleotide polymorphism
Baird JK. Asia-Pacific malaria is singular, pervasive, diverse and invisible. Int J Parasitol. 2017;47:371–7.
Bright AT, Manary MJ, Tewhey R, Arango EM, Wang T, Schork NJ, et al. A high resolution case study of a patient with recurrent Plasmodium vivax infections shows that relapses were caused by meiotic siblings. PLoS Negl Trop Dis. 2014;8:e2882.
WHO. World malaria report 2018. Geneva: World Health Organisation; 2018. https://www.who.int/malaria/publications/world-malaria-report-2018/en/. Accessed 29 Jan 2020.
Olliaro PL, Barnwell JW, Barry A, Mendis K, Mueller I, Reeder JC, et al. Implications of Plasmodium vivax biology for control, elimination, and research. Am J Trop Med Hyg. 2016;95:4–14.
Cui L, Yan G, Sattabongkot J, Cao Y, Chen B, Chen X, et al. Malaria in the Greater Mekong Subregion: heterogeneity and complexity. Acta Trop. 2012;121:227–39.
Delacollette C, D’Souza C, Christophel E, Thimasarn K, Abdur R, Bell D, et al. Malaria trends and challenges in the Greater Mekong Subregion. Southeast Asian J Trop Med Public Health. 2009;40:674–91.
Arnott A, Barry AE, Reeder JC. Understanding the population genetics of Plasmodium vivax is essential for malaria control and elimination. Malar J. 2012;11:14.
Baird JK. Resistance to therapies for infection by Plasmodium vivax. Clin Microbiol Rev. 2009;22:508–34.
Price RN, von Seidlein L, Valecha N, Nosten F, Baird JK, White NJ. Global extent of chloroquine-resistant Plasmodium vivax: a systematic review and meta-analysis. Lancet Infect Dis. 2014;14:982–91.
Chu CS, Phyo AP, Lwin KM, Win HH, San T, Aung AA, et al. Comparison of the cumulative efficacy and safety of chloroquine, artesunate, and chloroquine-primaquine in Plasmodium vivax malaria. Clin Infect Dis. 2018;67:1543–9.
Saravu K, Kumar R, Ashok H, Kundapura P, Kamath V, Kamath A, et al. Therapeutic assessment of chloroquine-primaquine combined regimen in adult cohort of Plasmodium vivax malaria from primary care centres in southwestern India. PLoS ONE. 2016;11:e0157666.
Cooper RD. Studies of a chloroquine-resistant strain of Plasmodium vivax from Papua New Guinea in Aotus and Anopheles farauti s.l. J Parasitol. 1994;80:789–95.
Baird JK. Chloroquine resistance in Plasmodium vivax. Antimicrob Agents Chemother. 2004;48:4075–83.
de Santana Filho FS, Arcanjo AR, Chehuan YM, Costa MR, Martinez-Espinosa FE, Vieira JL, et al. Chloroquine-resistant Plasmodium vivax, Brazilian Amazon. Emerg Infect Dis. 2007;13:1125–6.
Guthmann JP, Pittet A, Lesage A, Imwong M, Lindegardh N, Min Lwin M, et al. Plasmodium vivax resistance to chloroquine in Dawei, southern Myanmar. Trop Med Int Health. 2008;13:91–8.
Ketema T, Getahun K, Bacha K. Therapeutic efficacy of chloroquine for treatment of Plasmodium vivax malaria cases in Halaba district, South Ethiopia. Parasit Vectors. 2011;4:46.
Lee SW, Lee M, Lee DD, Kim C, Kim YJ, Kim JY, et al. Biological resistance of hydroxychloroquine for Plasmodium vivax malaria in the Republic of Korea. Am J Trop Med Hyg. 2009;81:600–4.
Ratcliff A, Siswantoro H, Kenangalem E, Wuwung M, Brockman A, Edstein MD, et al. Therapeutic response of multidrug-resistant Plasmodium falciparum and P. vivax to chloroquine and sulfadoxine-pyrimethamine in southern Papua, Indonesia. Trans R Soc Trop Med Hyg. 2007;101:351–9.
Singh RK. Emergence of chloroquine-resistant vivax malaria in south Bihar (India). Trans R Soc Trop Med Hyg. 2000;94:327.
Thanh PV, Hong NV, Van NV, Louisa M, Baird K, Xa NX, et al. Confirmed Plasmodium vivax resistance to chloroquine in Central Vietnam. Antimicrob Agents Chemother. 2015;59:7411–9.
Tjitra E, Anstey NM, Sugiarto P, Warikar N, Kenangalem E, Karyana M, et al. Multidrug-resistant Plasmodium vivax associated with severe and fatal malaria: a prospective study in Papua, Indonesia. PLoS Med. 2008;5:e128.
Chung DI, Jeong S, Dinzouna-Boutamba SD, Yang HW, Yeo SG, Hong Y, et al. Evaluation of single nucleotide polymorphisms of pvmdr1 and microsatellite genotype in Plasmodium vivax isolates from Republic of Korea military personnel. Malar J. 2015;14:336.
Huang B, Huang S, Su XZ, Tong X, Yan J, Li H, et al. Molecular surveillance of pvdhfr, pvdhps, and pvmdr-1 mutations in Plasmodium vivax isolates from Yunnan and Anhui provinces of China. Malar J. 2014;13:346.
Suwanarusk R, Chavchich M, Russell B, Jaidee A, Chalfein F, Barends M, et al. Amplification of pvmdr1 associated with multidrug-resistant Plasmodium vivax. J Infect Dis. 2008;198:1558–64.
Suwanarusk R, Russell B, Chavchich M, Chalfein F, Kenangalem E, Kosaisavee V, et al. Chloroquine resistant Plasmodium vivax: in vitro characterisation and association with molecular polymorphisms. PLoS ONE. 2007;2:e1089.
Fernandez-Becerra C, Pinazo MJ, Gonzalez A, Alonso PL, del Portillo HA, Gascon J. Increased expression levels of the pvcrt-o and pvmdr1 genes in a patient with severe Plasmodium vivax malaria. Malar J. 2009;8:55.
Melo GC, Monteiro WM, Siqueira AM, Silva SR, Magalhaes BM, Alencar AC, et al. Expression levels of pvcrt-o and pvmdr-1 are associated with chloroquine resistance and severe Plasmodium vivax malaria in patients of the Brazilian Amazon. PLoS ONE. 2014;9:e105922.
Chehuan YF, Costa MR, Costa JS, Alecrim MG, Nogueira F, Silveira H, et al. In vitro chloroquine resistance for Plasmodium vivax isolates from the western Brazilian Amazon. Malar J. 2013;12:226.
Hamedi Y, Sharifi-Sarasiabi K, Dehghan F, Safari R, To S, Handayuni I, et al. Molecular epidemiology of P. vivax in Iran: high diversity and complex sub-structure using neutral markers, but no evidence of Y976F mutation at pvmdr1. PLoS ONE. 2016;11:e0166124.
Schousboe ML, Ranjitkar S, Rajakaruna RS, Amerasinghe PH, Morales F, Pearce R, et al. Multiple origins of mutations in the mdr1 gene—a putative marker of chloroquine resistance in P vivax. PLoS Negl Trop Dis. 2015;9:e0004196.
Hupalo DN, Luo Z, Melnikov A, Sutton PL, Rogov P, Escalante A, et al. Population genomics studies identify signatures of global dispersal and drug resistance in Plasmodium vivax. Nat Genet. 2016;48:953–8.
Neafsey DE, Galinsky K, Jiang RH, Young L, Sykes SM, Saif S, et al. The malaria parasite Plasmodium vivax exhibits greater genetic diversity than Plasmodium falciparum. Nat Genet. 2012;44:1046–50.
Kittichai V, Nguitragool W, Ngassa Mbenda HG, Sattabongkot J, Cui L. Genetic diversity of the Plasmodium vivax multidrug resistance 1 gene in Thai parasite populations. Infect Genet Evol. 2018;64:168–77.
Ngassa Mbenda HG, Zeng W, Bai Y, Siddiqui FA, Yang Z, Cui L. Genetic diversity of the Plasmodium vivax phosphatidylinositol 3-kinase gene in two regions of the China–Myanmar border. Infect Genet Evol. 2018;61:45–52.
Yuan L, Wang Y, Parker DM, Gupta B, Yang Z, Liu H, et al. Therapeutic responses of Plasmodium vivax malaria to chloroquine and primaquine treatment in northeastern Myanmar. Antimicrob Agents Chemother. 2015;59:1230–5.
Bruce MC, Galinski MR, Barnwell JW, Snounou G, Day KP. Polymorphism at the merozoite surface protein-3α locus of Plasmodium vivax: global and local diversity. Am J Trop Med Hyg. 1999;61:518–25.
Yang Z, Miao J, Huang Y, Li X, Putaporntip C, Jongwutiwes S, et al. Genetic structures of geographically distinct Plasmodium vivax populations assessed by PCR/RFLP analysis of the merozoite surface protein 3β gene. Acta Trop. 2006;100:205–12.
Edgar RC. MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinformatics. 2004;5:113.
Kumar S, Stecher G, Tamura K. MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:1870–4.
Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.
Tajima F. Evolutionary relationship of DNA sequences in finite populations. Genetics. 1983;105:437–60.
Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95.
Fu YX, Li WH. Statistical tests of neutrality of mutations. Genetics. 1993;133:693–709.
Zhang LB, Ge S. Multilocus analysis of nucleotide variation and speciation in Oryza officinalis and its close relatives. Mol Biol Evol. 2007;24:769–83.
McDonald JH, Kreitman M. Adaptive protein evolution at the Adh locus in Drosophila. Nature. 1991;351:652–4.
Yang Z, Bielawski JP. Statistical methods for detecting molecular adaptation. Trends Ecol Evol. 2000;15:496–503.
Poon AF, Frost SD, Pond SL. Detecting signatures of selection from DNA sequences using Datamonkey. Methods Mol Biol. 2009;537:163–83.
Yang Z, Nielsen R. Codon-substitution models for detecting molecular adaptation at individual sites along specific lineages. Mol Biol Evol. 2002;19:908–17.
Pond SL, Frost SD, Grossman Z, Gravenor MB, Richman DD, Brown AJ. Adaptation to different human populations by HIV-1 revealed by codon-based analyses. PLoS Comput Biol. 2006;2:e62.
Yang J, Yan R, Roy A, Xu D, Poisson J, Zhang Y. The I-TASSER suite: protein structure and function prediction. Nat Methods. 2015;12:7–8.
Siddiqui FA, Cabrera M, Wang M, Brashear A, Kemirembe K, Wang Z, et al. Plasmodium falciparum falcipain-2a polymorphisms in southeast Asia and their association with artemisinin resistance. J Infect Dis. 2018;218:434–42.
Lu F, Lim CS, Nam DH, Kim K, Lin K, Kim TS, et al. Genetic polymorphism in pvmdr1 and pvcrt-o genes in relation to in vitro drug susceptibility of Plasmodium vivax isolates from malaria-endemic countries. Acta Trop. 2011;117:69–75.
Ferreira PE, Holmgren G, Veiga MI, Uhlen P, Kaneko A, Gil JP. PfMDR1: mechanisms of transport modulation by functional polymorphisms. PLoS ONE. 2011;6:e23875.
Patel SK, George LB, Prasanth Kumar S, Highland HN, Jasrai YT, Pandya HA, et al. A computational approach towards the understanding of Plasmodium falciparum multidrug resistance protein 1. ISRN Bioinform. 2013;2013:437168.
Pond SLK, Frost SD. Not so different after all: a comparison of methods for detecting amino acid sites under selection. Mol Biol Evol. 2005;22:1208–22.
Murrell B, Moola S, Mabona A, Weighill T, Sheward D, Kosakovsky Pond SL, et al. FUBAR: a fast, unconstrained Bayesian approximation for inferring selection. Mol Biol Evol. 2013;30:1196–205.
Delport W, Poon AF, Frost SD, Kosakovsky Pond SL. Datamonkey 2010: a suite of phylogenetic analysis tools for evolutionary biology. Bioinformatics. 2010;26:2455–7.
Hewitt S, Delacollette C, Chavez I. Malaria situation in the Greater Mekong Subregion. Southeast Asian J Trop Med Public Health. 2013;44(Suppl. 1):46–72 (discussion 306–7).
Cui L, Cao Y, Kaewkungwal J, Khamsiriwatchara A, Lawpoolsri S, Soe TN, et al. Malaria elimination in the Greater Mekong Subregion: challenges and prospects. In: Manguin S, Dev V, editors. Towards malaria elimination: a leap forward. London: IntechOpen; 2018. p. 179–200.
Geng J, Malla P, Zhang J, Xu S, Li C, Zhao Y, et al. Increasing trends of malaria in a border area of the Greater Mekong Subregion. Malar J. 2019;18:309.
Lo E, Lam N, Hemming-Schroeder E, Nguyen J, Zhou G, Lee MC, et al. Frequent spread of Plasmodium vivax malaria maintains high genetic diversity at the Myanmar-China border, without distance and landscape barriers. J Infect Dis. 2017;216:1254–63.
Tantiamornkul K, Pumpaibool T, Piriyapongsa J, Culleton R, Lek-Uthai U. The prevalence of molecular markers of drug resistance in Plasmodium vivax from the border regions of Thailand in 2008 and 2014. Int J Parasitol Drugs Drug Resist. 2018;8:229–37.
Das A, Mohanty S, Stephan W. Inferring the population structure and demography of Drosophila ananassae from multilocus data. Genetics. 2004;168:1975–85.
Fu YX. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147:915–25.
Wagner A. Rapid detection of positive selection in genes and genomes through variation clusters. Genetics. 2007;176:2451–63.
Benavente DE, Ward Z, Chan W, Mohareb FR, Sutherland CJ, Roper C, et al. Genomic variation in Plasmodium vivax malaria reveals regions under selective pressure. PLoS ONE. 2017;12:e0177134.
Zhu X, Zhao P, Wang S, Liu F, Liu J, Wang J, et al. Analysis of Pvama1 genes from China–Myanmar border reveals little regional genetic differentiation of Plasmodium vivax populations. Parasit Vectors. 2016;9:614.
Zeng W, Bai Y, Wang M, Wang Z, Deng S, Ruan Y, et al. Significant divergence in sensitivity to antimalarial drugs between neighboring Plasmodium falciparum populations along the eastern border of Myanmar. Antimicrob Agents Chemother. 2017;61:e01689–716.
Conway DJ, Machado RL, Singh B, Dessert P, Mikes ZS, Povoa MM, et al. Extreme geographical fixation of variation in the Plasmodium falciparum gamete surface protein gene Pfs48/45 compared with microsatellite loci. Mol Biochem Parasitol. 2001;115:145–56.
Feng H, Gupta B, Wang M, Zheng W, Zheng L, Zhu X, et al. Genetic diversity of transmission-blocking vaccine candidate Pvs48/45 in Plasmodium vivax populations in China. Parasit Vectors. 2015;8:615.
Rogers AR. How population growth affects linkage disequilibrium. Genetics. 2014;197:1329–41.
Koepfli C, Ross A, Kiniboro B, Smith TA, Zimmerman PA, Siba P, et al. Multiplicity and diversity of Plasmodium vivax infections in a highly endemic region in Papua New Guinea. PLoS Negl Trop Dis. 2011;5:e1424.
We would like to thank the medical staff at the border clinics for collecting the clinical P. vivax samples.
This study was supported by grants from the National Institute of Allergy and Infectious Diseases (U19 AI089672) and the Fogarty International Center (D43 TW006571), National Institutes of Health, USA.
Ethics approval and consent to participate
Ethical approvals were obtained from the institutional review boards of China Medical University, Kunming Medical University, Mahidol University and Pennsylvania State University. Written consent was obtained from all donors of the blood samples or their guardians.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Ngassa Mbenda, H.G., Wang, M., Guo, J. et al. Evolution of the Plasmodium vivax multidrug resistance 1 gene in the Greater Mekong Subregion during malaria elimination. Parasites Vectors 13, 67 (2020). https://doi.org/10.1186/s13071-020-3934-5