Genetic diversity and natural selection of Plasmodium vivax multi-drug resistant gene (pvmdr1) in Mesoamerica

Background The Plasmodium vivax multidrug resistant 1 gene (pvmdr1) codes for a transmembrane protein of the parasite’s digestive vacuole. It is likely that the pvmdr1 gene mutations occur at different sites by convergent evolution. In here, the genetic variation of pvmdr1 at three sites of the Mesoamerican region was studied. Since 1950s, malarious patients of those areas have been treated only with chloroquine and primaquine. Methods Blood samples from patients infected with P. vivax were obtained in southern Mexico (SMX), in the Northwest (NIC-NW) and in the northeast (NIC-NE) of Nicaragua. Genomic DNA was obtained and fragments of pvmdr1 were amplified and sequenced. The nucleotide and amino acid changes as well as the haplotype frequency in pvmdr1 were determined per strain and per geographic site. The sequences of pvmdr1 obtained from the studied regions were compared with homologous sequences from the GenBank database to explore the P. vivax genetic structure. Results In 141 parasites, eight nucleotide changes (two changes were synonymous and other six were nonsynonymous) were detected in 1536 bp. The PvMDR1 amino acid changes Y976F, F1076FL were predominant in endemic parasites from NIC-NE and outbreak parasites in NIC-NW but absent in SMX. Thirteen haplotypes were resolved, and found to be closely related, but their frequency at each geographic site was different (P = 0.0001). The pvmdr1 codons 925–1083 gene fragment showed higher genetic and haplotype diversity in parasites from NIC-NE than the other areas outside Latin America. The haplotype networks suggested local diversification of pvmdr1 and no significant departure from neutrality. The F ST values were low to moderate regionally, but high between NIC-NE or NIC-NW and other regions inside and outside Latin America. Conclusions The pvmdr1 gene might have diversified recently at regional level. In the absence of significant natural, genetic drift might have caused differential pvmdr1 haplotype frequencies at different geographic sites in Mesoamerica. A very recent expansion of divergent pvmdr1 haplotypes in NIC-NE/NIC-NW produced high differentiation between these and parasites from other sites including SMX. These data are useful to set a baseline for epidemiological surveillance.

Plasmodium falciparum have been responsible for a fluctuating number of malaria cases during the last decade. From 2008 to 2011, there were less than 1000 cases reported, but the following years registered an increase in this number; in 2014, a total of 1163 cases were confirmed by microscopy. More than 70% of the cases reported in this country are produced in the Northeast or North Caribbean Coast Autonomous Region, while transmission is sporadic (with periodic outbreaks) in the Northwest side of the country. The most affected area of countries Nicaragua and Honduras is the indigenous Miskito region, which lies on both sides of the border (Northeast). Also in 2014, the number of reported malaria cases was 4931 in Honduras and 3380 in Guatemala [2], the latter country between Honduras and Mexico. Since 1999, only P. vivax has been transmitted in Mexico, with 696 cases reported in 2014. This country is in the pre-elimination phase, similar to El Salvador and Belize [2].
In Mesoamerica, two drugs have been used to treat malaria since the 1950s: chloroquine (CQ) and primaquine (PQ) [3]. CQ interferes with haemoglobin digestion by the parasite and causes death of asexual and sexual blood parasites. However, PQ is needed to eliminate dormant stages of P. vivax and young gametocytes of P. falciparum [4]. The first cases of P. vivax resistance to CQ were reported in Oceania and Asia, and later in South America [5]. The multidrug resistance 1 gene (mdr1) codes for a transport protein, which is comprised of two transmembrane domains and two ATP binding sites [6]. In Southeast Asia, a low in vitro sensitivity of P. vivax to CQ was associated with the PvMDR1 amino acid substitution Y976F [7,8]. The gene fragment containing this and other polymorphisms (e.g. F1076L, also a transmembrane residue), have been found in distinct geographic sites and at different frequencies [9]. In South America, four haplotypes comprising residues (976/1076; Y/F, Y/L, F/F and F/L) were reported as well as one polymorphism (T958M) at a marginal transmembrane position [6,7,[10][11][12]. In P. vivax from Honduras, the PvMDR1 haplotypes 976/1076; Y/L, F/F and F/L were described, but the wild type (Y/F; Sal I strain) was the most frequent (>70%) [13]. It is remarkable that in clinical studies at different geographic sites no association between PvMDR1 polymorphism and CQ resistance has been found [12,14,15]. Besides, a recent study suggested that Pvmdr1 mutations emerged at different geographic sites possibly by convergent evolution [11].
Seeking to extend the genetic characterization of current P. vivax in Mesoamerica, parasites were obtained from three different geographic sites. Then, the occurrence of natural selection and the genetic relationships of P. vivax studying mdr1 polymorphisms and haplotype frequency were investigated.

Blood samples and geographic origin
All infected blood samples were obtained from symptomatic patients seeking malaria diagnosis in Mexico and Nicaragua (in each case, after informed consent was given). From 2008 to 2010 in Mexico, 93 samples were obtained in the laboratory facility at the Regional Research Center for Public Health (CRISP-INSP) from patients living in the Tapachula municipality and its surroundings in southern Chiapas (SMX), the southernmost point of Mexico that borders with Guatemala [16]. From 2011 to 2012 in Nicaragua, 107 samples were selected from the sentinel laboratory network established by the Health Ministry at head municipalities, 83 from the Northeast (NIC-NE; RACCN, North Caribbean Coast Autonomous Region) and 24 from the Northwest (NIC-NW; North Pacific Coast). Another 14 samples were included that had been obtained in NIC-NW during 2006-2007 [17]. Distances between these sites are the following: ≈611 km from SMX to NIC-NW, ≈853 km from SMX to NIC-NE, and ≈331 km from NIC-NE to NIC-NW. The NIC-NE in comprised by the Miskito and the mining regions.
The diagnosis of P. vivax was carried out by microscopic examination of stained thick blood smears, coming from capillary blood samples (before treatment) used to impregnate filter paper (Whatman #2). Then the participants were administered the CQ-PQ combination treatment, which is in accordance with health standards in Mexico [18,19] and in Nicaragua [20]. Genomic DNA was extracted from P. vivax-infected bloods impregnated on filter paper, using a commercially available QIAmp DNA blood Minikit (Qiagen, USA), and following the manufacturer's instructions.
All amplified products were resolved in agarose gels at 1%, and stained with 0.2 µg/mL ethidium bromide using an electrophoresis chamber Midicell primo (Thermo EC330, New York, USA). A molecular marker of 100 bp ladder was used (Invitrogen Corporation, Carlsbad, CA, USA). PCR products were then purified using a MiniElute PCR Purification Kit (Qiagen, Valencia, CA, USA), according to the manufacturer's instructions. The purified products were Sanger sequenced by using forward and reverse primers (at the High Throughput Genomics Unit, Department of Genome Sciences, University of Washington, Seattle, WA, USA).
The quality of pherograms with the forward and reverse nucleotide sequences was verified manually and by using BioEdit v7.1.3 software. DNA sequences with rare polymorphisms were re-amplified and sequenced. All pherograms showed a single genotype pattern. The consensus sequences obtained for each gene fragment were submitted to the NCBI-Gen Bank [Accession Numbers: KX180164-KX180638].

Data analysis
The Salvador I strain (Sal I) mdr1sequence was used as reference: XM_001613678. Sequences forward and reverse were aligned using ClustalW Multiple Alignment of BioEdit v7.0 [21] and revised manually. Nonsynonymous and synonymous changes were identified, and the frequency of each nucleotide change was calculated per geographic site in Mesoamerica (NIC-NE, NIC-NW and SMX), as well the haplotypes were constructed to background mutations at codons 976/1076 with the dnaSP program v5.1 [22]. Statistical analysis was carried out with STATA v12.1.
In order to explore the regional and global parasite relationships based on pvmdr1, haplotype networks were constructed using TCS 1.21 [23]. Mutational steps represent the connections between haplotypes, and empty squares showed the non-sampled or extinct haplotypes. The colour of the circles represents the geographic origins of each haplotype, while the size of the circle represents the frequency of each haplotype.
Indexes of nucleotide diversity were calculated in dnaSP v5.1. To test whether natural selection and gene flow shaped the evolution of pvmdr1 gene in parasite populations of Mesoamerica, the number of synonymous (s) and nonsynonymous (ns) nucleotide changes and the difference in the rate of nonsynonymous versus synonymous changes (dN˗dS) were determined by using the Nei Gojobori proportion method with 1000 bootstrap replicates in MEGA v6.0 software [24]. Tajima's D test and the minimal number of recombination events were calculated using dnaSP. To estimate the P. . Additionally, other DNA sequences were reconstructed, using the aforementioned sequence reference of the Sal I strain as a template (AY618622, with codon 958 atg/M; AY571984 and XM001613678 with codon 958 acg/T), and only if the length and coordinates of the pvmdr1 gene fragment, synonymous and nonsynonymous nucleotide changes, and haplotype frequency were clearly described: IND, n = 25 [28]. Honduras (HON), n = 37 [13]. Nepal (NEP), n = 39 [29]. Ecuador (ECU), n = 17; Sri Lanka (SLK), n = 119; Pakistan (PK), n = 24; Sudan (SUD), n = 4; Sao Tomé (SAT), n = 3 [11].

Pvmdr1 gene polymorphism
Upon analysing 1536 bp of pvmdr1, two synonymous and six nonsynonymous nucleotide changes were detected (Table 1). From 163 parasites (66, 23 and 74 from NIC-NE, NIC-NW and SMX, respectively), two changes (one nonsynonymous and one synonymous) were detected in the gene fragment of 411 bp (nucleotides 1330-1740). The nonsynonymous change was observed (only in NIC-NW) at codon 500 D → N). While the synonymous change at codon 529, was highly frequent in SMX and NIC-NE, but not detected in parasites from NIC-NW. Interestingly, from 155 parasites (65, 26 and 64 from NIC-NE, NIC-NW and SMX, respectively) a second fragment of 549 bp (nucleotides 2767-3315) was obtained, and six nucleotide changes were detected. The nonsynonymous changes at codons Y976F and F1076L were highly frequent in NIC-NE and NIC-NW, but not detected in SMX. Similarly, the synonymous change at codon 1021 was highly frequent in NIC-NE and NIC-NW, but rare in SMX. The change at codon A927T exclusive to NIC-NE, it was detected in parasites from different municipalities of the Miskito and the mining regions. The substitution T958M was highly frequent at the three sites, while the F1070L change was detected in only two isolates from NIC-NE (Rosita and Prinzapolka municipalities) and two isolates from SMX (from the year 2010). The third gene fragment of 576 bp (nucleotides 2032-2607), obtained from 148 parasites (63, 17 and 72 from NIC-NE, NIC-NW and SMX, respectively), was identical to the Sal I strain sequence.

Pvmdr1 haplotypes and its geographic frequency
There were 141 parasite isolates (59 from NIC-NE, 59 from SMX and 23 from NIC-NW) with the nucleotide information at variant codons 500, 529-927, 958, 976, 1021, 1070 and 1076 resolved 13 pvmdr1 haplotypes (m1 → m13; ordered by their frequency) (Fig. 1). There were differences in haplotype frequencies between the three sites (P = 0.000): haplotype m1 GG-GTACTT was predominant (93.2%) in SMX; m3 AA-GTTTTC was highly frequent (78.2%) in NIC-NW; m2 GA-GTTTTC and m4 GG-ATATTC were at 62.7 and 22%, respectively, in NIC-NE. Other haplotypes were detected in two or all sites, although at a very low frequency. For instance, haplotype m1 was also detected in one isolate from NIC-NE and in another from NIC-NW, m6 and m10 were observed at SMX and NIC-NE, and m5 was found in NIC-NW and NIC-NE. One haplotype (m8) was exclusive to SMX and resembled the Sal I strain sequence. Another three haplotypes (m9, m11 and m13) were exclusive to NIC-NE and two (m7 and m12) to NIC-NW (Fig. 1). There were four haplotypes having nucleotide changes at codons 976F (a → t) and 1076L (t → c); m2, m3, m7: GG-GTTCTC, m12: GG-GTTTTC; those haplotypes differed among them by one to three mutations. Two other haplotypes showed only 1076F/L; m4 and m11: GA-ATATTC. The haplotypes m2 and m4, highly frequent at NIC-NE were also frequent in all sampled municipalities: Rosita, Bonanza, Waspam and Puerto Cabezas. The higher number of exclusive haplotypes were detected in Waspam (n = 3), the Miskito region that borders Honduras.

Genetic relationship of pvmdr1 haplotypes from SMX, NIC-NW and NIC-NE
The haplotype network of P. vivax mdr1 comprising mutations at codons 500, 529-927, 958, 976, 1021, 1070 and 1076 for parasites from the study sites is shown in  Fig. 2. All 13 haplotypes were found closely related; most haplotypes were at one mutational step from each other, and at 1-9 mutational steps among them. Haplotypes from SMX were only separated by 1-2 mutational among them. The most frequent haplotype (m1) in SMX was at one mutational step from the Sal I strain sequence. While haplotypes m2 and m3 were at about 6-9 mutational steps from m1; they were highly frequent in NIC-NE and NIC-NW, respectively.

Genetic diversity
Parasites from NIC-NE and mostly from the Mesoamerican region showed a nucleotide diversity (π = 0.0036) similar to that found at other geographic origins as well as the global level (π = 0.0028). Moreover, the minimal number of recombination events were from very low to cero (the value for Mesoamerica was Rm = 2) ( Table 2).

Natural selection
Tajima's D test statistical value was negative for SMX and NIC-NW, while being positive for NIC-NE and other sites, but no significant. The lowest value was for Brazil (−1.681) and the highest for Ecuador (0.949), but in either case, it was found to be no significant. Although a higher number of nonsynonymous than synonymous changes were detected, the dN-dS values showed no significant departure from neutrality. Parasites from Nepal had the highest dN-dS value with marginal significance (1.845; P = 0.067) ( Table 2).

Haplotype relationships
When using global sequences of P. vivax mdr1, the haplotype network evidenced the relationship of 19 global haplotypes (g1 → g19) identified in 610 sequences (Fig. 3). The most frequent haplotype g1 and probably the ancestral haplotype, it was predominant in parasites from Sri Lanka, India, Nepal, Pakistan and Madagascar and existed in few parasites from Latin America. In contrast, the haplotype g2 was highly frequent in parasites from Latin America (SMX, Brazil, Honduras and Ecuador).
Haplotype g2 was at one mutational step from g5, the latter represent Sal I strain (the only haplotype with codon 958 → acg) which was found in parasites from Latin America and India. Moreover, from g2 emerged a branch of four descendant haplotypes (restricted to American parasites), one present in SMX and NIC-NE (connected to another haplotype present in parasites from SMX, NIC-NE and Brazil), and two exclusive to NIC-NE (Fig. 3). Whereas, haplotype g3 was included in one parasite from Brazil and another from NIC-NW, is at one mutational step from g4, the latter of which showed high frequency in NIC-NE and NIC-NW parasites and was found in one parasite from Brazil (Fig. 3

Discussion
The Mesoamerican region has ecological and demographic conditions that favor malaria transmission. Chloroquine (CQ) and primaquine (PQ) have been the only drugs used in this region to treat patients infected with P. vivax and/or P. falciparum [19,20,30] and, so far, both Plasmodium species seem susceptible to CQ. According to this, mutations at pvmdr1 might not be associated to CQ resistance. Therefore, PvMDR1 976F and CQ resistance in vitro seemed to be a counterfeit association [8]. In fact, clinical studies have shown that polymorphisms at pvmdr1 e.g. 976F and others adjacent are not associated with CQ resistance [12,14,28], and there are no functional studies exposing pvmdr1 mutations participating in the development of CQ resistance. Schousboe et al. [11] suggested that pvmdr1 mutations 958M, 976F and 1076L arose independently and several times at different geographic sites. In Mesoamerica, the selective force of natural selection on pvmdr1 differed spatially, with lack of statistical significance. Besides, recent genomic studies have not revealed signature of selection on pvmdr1 in parasites from different geographic sites [31][32][33], including parasites from the Southeast Asia [33] where high levels of CQ resistance have been reported [8,34]. Malaria has been historically hypo endemic in Mesoamerica. From 2000 to 2010 a reduction in more than 75% was recorded in several countries including Nicaragua and Mexico. Contrarily from 2012 to 2014, Nicaragua was the only country in Mesoamerica that reported a significant increase in malaria cases [2]. In this country, while intense and effective control strategies caused reduction in the number of malaria cases, climate changes and delay in the delivery of control measures have probably caused the increase in malaria transmission. Alike southern Mexico, in NIC-NW and NIC-NE several haplotypes exposing changes 976F and/or 1076L were exclusive and highly frequent. The haplotype networks suggest that these mutations resulted from a recent diversification in Mesoamerica. It is remarkable that haplotypes expressing 1076L, present in parasites from NIC-NE (haplotype g6) seemed the result of different evolutionary pathway than in other geographic sites (haplotype g1).
The presence of some pvmdr1 haplotypes in all sites (SMX, NE-NIC and NW-NIC), may be attributed to ancient parasite flow across the Mesoamerican region. It was notorious a predominance of blood infections with a single pvmdr1 haplotype, presumably monoclonal infections. A clonal transmission of P. vivax might be occurring in the region. It is supported by the lack of mixed genotype infections for pvmsp1 42 a gene fragment, of moderate diversity, reported in NIC-NE/NW [35]. In SMX, the rate of mixed genotype infections has been negligible during the last decade [16,36]. In Mesoamerica, parasite flow might be controlled by different factors, including separation by distance, geographic barriers, cultural and language barriers, poor road infrastructure, the effectiveness of the malaria control interventions. There is a volcanic arc in between NIC-NE and NIC-NW. At either side of the arc exclusive pvmdr1 haplotypes were detected. Similarly, for the blood stage antigen pvmsp1 42 a divergent haplotype was exclusive in Chinandega municipality (NIC-NW) [35]. In NIC-NE, the most frequent pvmdr1 haplotypes in the region, so were at the municipality level. A similar spatial pattern of haplotypes was reported for pvmsp1 42 [35]. Yet, the presence of exclusive haplotypes suggests partial parasite flow between the Miskito and mining regions in NIC-NE. By using microsatellites, the analysis of P. falciparum exposed that in the mining region there was a genetic population distinct from that detected in the Miskito area [37]. The Miskito ethnic group (Waspam and Puerto Cabezas) live in multiple rural communities and it is partly inaccessible due to cultural and language barriers. Alike P. falciparum, P. vivax cause relapse episodes of long latency, an additional risk of parasite dispersion [38], and microsatellite markers provided higher resolution than single nucleotide polymorphism, and subtle genetic changes can be detected [39]. It is probable that some pvmdr1 mutations might have appeared at different sites in Mesoamerica as formerly suggested for other geographic sites [11].  7). Each circle indicates a different haplotype (n = 19) connected by lines. The empty squares indicate haplotypes not sampled or extinct. Each colour indicates a different country. Three frequent haplotypes were detected at different geographic sites, respectively. One haplotype was frequent in Latin America. Signs of gene diversification were found at the regional and at local level. Sal I haplotype: XM001613678. A gene fragment of 479 bp (nt 2772-3250) was used. NIC-NW, Northwest, Nicaragua. SMX, Southern Mexico The results from this study and the epidemiological and demographic scenario suggest that a strong genetic drift might have influenced the spatial differences of the pvmdr1 haplotype frequencies in Mesoamerica. The Sal I type was only detected in SMX and at low frequency. In Honduras, parasites isolated from 2004 to 2009 exhibited haplotypes with PvMDR1 976F and/or 1076L changes at a frequency of <30% [13]; while in the Northeast Nicaragua area the frequency of those changes was almost 90%. These changes were also present in strains causing an outbreak in the Northwest Nicaragua in 2012, it corresponded to an exclusive haplotype. In South America, mutations 976F and 1076L were reported to have low frequency [12,31].
Plasmodium vivax arriving at the American continent faced adaptive processes triggered by different local factors e.g. the specificity of the Anopheline species, the host immune pressure, the anti-malarial drugs employed, and the ecological conditions to mention some reasons [40][41][42]. Genomic studies have exposed that P. vivax exhibits low diversity and high differentiation in Latin America [42,43]. Analysis of the mitochondrial genome exposed low F ST values between South America and Central America, the last ones corresponded to ten strains obtained in the 1970s [44]. Evidence based on the entire pvmdr1 gene also supports a high genetic differentiation between old and new world parasites [32]. In the present study, the F ST value using the pvmdr1 codons925-1083 matched partially the previous results, also suggesting a strong differentiation globally but low-moderate differentiation at regional level. Excepting parasites from NIC-NE/NW, these were highly differentiated from all other parasites including those from SMX and Honduras, possibly given by the recent expansion of pvmdr1 divergent haplotypes.
Regardless of its origin, P. vivax has been at some extent circulating in the American continent. Southern Mexico is a transit area for human migrations coming from South and Central America [45], and still a proportion of humans migrants entering North America can be potentially infected with P. vivax [46,47]. Therefore, further genomic studies are necessary to uncover the evolution of P. vivax in the Mesoamerican region comparing to other sites of America and worldwide. Furthermore, molecular monitoring of P. vivax mdr1 at different sites in Mesoamerica would help to understand transmission dynamics. Because of its variability and spatial patterns, pvmdr1 is an important marker for epidemiological surveillance to support the current initiatives to control and eliminate malaria from Mesoamerica.

Conclusions
These results suggest a recent diversification of pvmdr-1which likely occurred at a regional level, natural selection varied spatially and was mostly in a weakly to-mild manner. In this hypo-endemic region, a strong genetic drift might have caused the great differences in pvmdr1 haplotype frequency among the studied sites, and hence the high differentiation of parasites between NIC-NW/ NIC-NE and SMX or other sites. These data are relevant to set a baseline level for epidemiological surveillance in the region.