Prevalence of pvmrp1 Polymorphisms and Its Contribution to Antimalarial Response

As more sporadic cases of chloroquine resistance occur (CQR) in Plasmodium vivax (P. vivax) malaria, molecular markers have become an important tool to monitor the introduction and spread of drug resistance. P. vivax multidrug resistance-associated protein 1 (PvMRP1), as one of the members of the ATP-binding cassette (ABC) transporters, may modulate this phenotype. In this study, we investigated the gene mutations and copy number variations (CNVs) in the pvmrp1 in 102 P. vivax isolates from China, the Republic of Korea (ROK), Myanmar, Papua New Guinea (PNG), Pakistan, the Democratic People’s Republic of Korea (PRK), and Cambodia. And we also obtained 72 available global pvmrp1 sequences deposited in the PlasmoDB database to investigate the genetic diversity, haplotype diversity, natural selection, and population structure of pvmrp1. In total, 29 single nucleotide polymorphisms reflected in 23 non-synonymous, five synonymous mutations and one gene deletion were identified, and CNVs were found in 2.9% of the isolates. Combined with the antimalarial drug susceptibility observed in the previous in vitro assays, except the prevalence of S354N between the two CQ sensitivity categories revealed a significant difference, no genetic mutations or CNVs associated with drug sensitivity were found. The genetic polymorphism analysis of 166 isolates worldwide found that the overall nucleotide diversity (π) of pvmrp1 was 0.0011, with 46 haplotypes identified (Hd = 0.9290). The ratio of non-synonymous to synonymous mutations (dn/ds = 0.5536) and the neutrality tests statistic Fu and Li’s D* test (Fu and Li’s D* = −3.9871, p < 0.02) suggests that pvmrp1 had evolved under a purifying selection. Due to geographical differences, genetic differentiation levels of pvmrp1 in different regions were different to some extent. Overall, this study provides a new idea for finding CQR molecular monitoring of P. vivax and provides more sequences of pvmrp1 in Asia for subsequent research. However, further validation is still needed through laboratory and epidemiological field studies of P. vivax samples from more regions.


Introduction
Plasmodium vivax was responsible for up to 4.5 million cases of malaria in 2020 [1]. Although it is rarely fatal, P. vivax malaria is the leading cause of malaria-related deaths outside of Africa [2]. In most vivax-endemic areas, a combination of chloroquine (CQ) is the first-line treatment for uncomplicated vivax malaria. Despite inexpensive, well tolerated, and widely available chloroquine resistance (CQR), vivax malaria was first reported from Papua New Guinea (PNG) in 1989 [3], and has been documented in more than ten countries [4], especially in multiple regions of Myanmar [5][6][7]. P. vivax resistance to other antimalarial drugs, such as mefloquine (MQ), sulfadoxine-pyrimethamine (SP) and primaquine (PQ), has also been reported widely [8,9]. The rise of these drug-resistant parasites threatens the global efforts to control malaria. However, molecular markers of drug resistance in P. vivax remain elusive [10].
Generally, genetic variation in the expression of transporter proteins could contribute to evading antimalarial action [11]. ATP-binding cassette (ABC) transporters are transmembrane proteins that can carry various substrate types, such as drugs and metabolic products [12]. The overexpression or mutation of many ABC transporters can lead to drug resistance [13]. As one of the members of the ABC subfamily C, the multidrug resistanceassociated proteins (MRPs) are associated with antimalarial resistance. The increased expression of pfmrp1 has been associated with the resistance of MQ and CQ, and gene polymorphisms in pfmrp1 with in vivo selection after SP and artemether/lumefantrine treatment [14][15][16]. Furthermore, the deletion of this gene in the CQR P. falciparum strain results in increased sensitivity to CQ, Quinine (QN), artemisinin, piperaquine (PPQ) and PQ [17]. Since P. vivax cannot be cultured continuously in vitro [18], most research on the molecular mechanism of drug resistance in P. vivax has focused on the homologous genes related to the drug resistance of P. falciparum [10,19]. Hence, we chose to study the association between pvmrp1 and the development of P. vivax resistance to antimalarial drugs.
In this study, we genotyped 94 P. vivax isolates from China, the Republic of Korea (ROK), Myanmar, PNG, Pakistan, the Democratic People's Republic of Korea (PRK), and Cambodia, and combined the 72 available global pvmrp1 sequences deposited in the PlasmoDB database to analyze the characterization of genomic variation and population genomics methods, including genetic differentiation, haplotype network, linkage disequilibrium (LD) and the phylogenetic tree. Meanwhile, the correlation between pvmrp1 and antimalarial drug susceptibility observed in the previous in vitro assays was further analyzed.

Study Sites and Participants
Clinical blood samples with P. vivax infections (n = 102) were obtained from seven countries, including China (n = 46), ROK (n = 27), Myanmar (n = 21), PNG (n = 3), Pakistan (n = 3), PRK (n = 1), and Cambodia (n = 1). The Chinese samples were collected from local hospitals or centers for disease control and prevention in central China from 2005 to 2008 [20]; The South Korean samples were from local hospitals in endemic areas, such as the ROK, from 2007 to 2009 [21]; The samples of Myanmar were collected from Wet-Won Station Hospital, Yangon, Myanmar, in 1999 [21]; Other samples of P. vivax isolates were imported malaria cases in China. The protocol was reviewed and approved by the National Institute of Parasitic Diseases, the Chinese Center for Disease Control and Prevention, the Kangwon National University Hospital Human Ethics Committee, and the Myanmar Department of Health. The sequence of pvmrp1 from 72 P. vivax isolates deposited in the PlasmoDB database were downloaded and analyzed, and originated from 10 countries, including South America: Columbia (n = 21), Peru (n = 16), Mexico (n = 13), and Brazil (n = 3); Southeast Asia: Myanmar (n = 4) and Thailand (n = 7); Oceania: PNG (n = 4); East Asia: PRK (n = 1); South Asia: India (n = 2) and Africa: Mauritania (n = 1). All sample information was listed in Spreadsheet S1, and all sequences have been uploaded to Genbank with accession numbers from ON933478 to ON933571.

Single Nucleotide Polymorphisms (SNPs) Identification in pvmrp1 Gene
According to the manufacturer's instructions, genomic DNAs from the 102 whole blood samples were individually extracted using a QIAamp DNA blood kit (Qiagen, Valencia, CA, USA) and were stored at −20 • C in the previous studies. The pvmrp1 gene was amplified by nested or semi-nested PCR using specific primers (Table 1). All reactions were performed in 20 µL containing 4 µL of 5× Phusion HF Buffer (7.5 mM Mg 2+ plus), 0.2 mM of each dNTP, 0.25 µM of each outer primer, 0.4 U Phusion High Fidelity DNA Polymerase (New England Biolabs, Ipswich, MA, USA) and 1 µL of genomic DNA or the amplicon from the first PCR. The PCR was performed with initial denature at 98 • C for 30 s, followed by 35 cycles of 98 • C for 10 s, 56-61 • C for 30 s (Table 1), and 72 • C for 1.5 min, and a final extension period at 72 • C for 10 min. The second round of PCR amplification products was purified and sequenced by GenScript (Nanjing, China). The nucleotide sequences were compared and spliced using Lasergene software (DNASTAR, Madison, WI, USA).

Determination of pvmrp1 Copy Number (CN)
The copy number variations (CNVs) of pvmrp1 were measured by TaqMan-BHQ1 probe quantitative PCR assays performed on the Roche LC480 thermal cycler. A reference plasmid was constructed with pvmrp1 (nt, 1762-1872) and pvtubulin (nt, 1644-1765) fragments in a ratio of 1:1 similar as described previously [21]. Probes and primers used for amplification of both genes were listed in Table 1. The probe target pvmdr1 was labeled at the 5 end with the FAM reporter dye, and pvtubulin was the HEX reporter dye, both of which were labeled at the 3 end with the quencher dye BHQ1. A real-time PCR was conducted in 10 µL volumes containing 5 µL of Lightcycler ® 480 Probes Master 2 × (Roche Applied Science, Penzberg, Germany), 400 nM of each forward and reverse primer, 250 nM of each probe, and 1 µL of template DNA. Amplifications were performed in triplicate and the cycling parameters were as follows: 95 • C for 10 min, then 40 cycles of 95 • C for 15 s and 58 • C for 30 s. The single copy of the pvtubulin gene served as an internal control. The relative CN of pvmrp1 was calculated using a relative standard curve method as normal, and the amplifications were repeated according to specifications [20,21].

Data Analysis
The multiple sequence alignments of worldwide isolates containing the wild reference sequence of pvmrp1 were obtained using the MUSCLE in the MEGA v7.0.18 program to obtain SNPs [25,26]. Moreover, the average nucleotide diversity (π), haplotype diversity (Hd), and the neutrality tests (Tajima's D test and Fu and Li's D* test) were further analyzed to identify pvmrp1 gene polymorphism and determined whether it is under the neutral evolution model [26,27], in order to evaluate the evolutionary relationship of the pvmrp1 gene. The estimation of genetic differentiation (F ST ) of the pvmrp1 was analyzed. Haplotype network was also implemented to identify the genetic association of the pvmrp1 haplotypes by means of the Median-Joining method in the NETWORK v5.0 program [28]. The phylogenetic tree of the aligned sequences was constructed using the Neighbor-Joining method in the MEGA v7.0.18 program [29]. In addition, pairwise LD of pvmrp1 gene at different polymorphic sites was calculated using the DNASP v5.10.01 program [30]. Compared to the previous in vitro drug sensitivity test [20,21], Fisher's exact test and chi-square analysis were performed to analyze whether it was related to the polymorphism of pvmrp1. As the sample size is too small to carry out statistical analysis, the data would not be displayed. A value of p < 0.05 was considered significant.

Identification of Gene Mutations and CN in pvmrp1 among Collected Blood Samples in Asia
Of the 102 P. vivax isolates, 94 samples (92.1%) were sequenced successfully for the pvmrp1 gene, including China (n = 43), ROK (n = 27), Myanmar (n = 18), PNG (n = 2), Pakistan (n = 2), PRK (n = 1) and Cambodia (n = 1). Compared with Sal-I as the wild reference type, 29 polymorphic sites were observed, 23 (79.3%) of which resulted in nonsynonymous mutations, one site E533 (GAA) deletion, and five (17.2%) synonymous mutations were identified. Compared with the previously reported SNPs [10], nine nonsynonymous mutations were repeatable, and the present analysis revealed 20 different SNPs ( Figure 1B). Three non-synonymous substitutions, including T259R (97.87%), Y1393D (97.87%) and V1478I (95.74%), were high prevalence and approached fixation. The wild type T259 was found in two isolates from Oceania, and the other two wild types were found in isolates from South Asia and East Asia ( Table 2). Of the 29 mutations found in this study, 20 were shown to be region-specific, such as mutations R281K, S354N, E787D, A853, G949D, and V1360 which were only observed in East Asia (range: 12.7-38.0%), V879 and L1207I were unique to the Southeast Asia at high frequency (89.5%). Three mutations (T234M, Q906E and I1232) were more frequent in isolates from Southeast Asia compared to the other areas. Although only 2 isolates were from South Asia, four mutations (F271Y, T282M, F560I and G1419A) were identified exclusively. Also, two isolates were confirmed, I1620T specifically in Oceania (Table 2). In this study, the amplification of pvmrp1 was determined by the relative CN, which was calculated by the standard curve method. Except for one isolate from China, the CN of pvmrp1 was assessed successfully from the other 101 isolates. The estimates of pvmrp1 CN for these isolates ranged from 0.68 to 2.55. Most of the isolates carried one copy of the gene, and only three isolates had double CNs, two from the ROK and one from Pakistan ( Figure 2).

Correlation between Polymorphisms and In Vitro Drug Susceptibilities
Of the 102 sequenced isolates mentioned above, partial samples from China and ROK were tested for antimalarial drugs susceptibility in the previous study [20,21]. Combined, a total of 39, 34, 39, and 13 isolates were cultured for more than 24 h and assayed for the susceptibility to CQ, QN, MQ and pyrimethamine (PYR), respectively (Table S1, Figure S1). For CQ, the geometric mean IC 50 was 20.92 nM (95% CI:   G949D, K1219N, Y1393D, V1478I, and H1586Y), did not appear to be correlated with the IC 50 values of the four antimalarial drugs (Table S2). Using an IC 50 value ≤ 220 nM as the sensitivity standard, the chi-square analysis was performed to analyze the correlation between the CQ sensitive/insensitive isolates and the mutation sites, the results showed that only the prevalence of S354N between these two CQ sensitivity categories revealed a significant difference (p < 0.05; Table 3). With regard to the variation of pvmrp1 CN, increased pvmrp1 CN did not appear to significantly alter parasites' susceptibilities to CQ, MQ and QN ( Figure S2).

The Polymorphism of pvmrp1 from Different Regions
To estimate the degree of genetic differentiation of the pvmrp1 in global isolates, the sequences of pvmrp1 obtained from the studied regions were compared with homologous sequences from the PlasmoDB database. An additional 72 pvmrp1 sequences from Plas-moDB were downloaded for further analysis, including Myanmar (n = 4), PNG (n = 4), PRK (n = 1), Columbia (n = 21), Peru (n = 16), Mexico (n = 13), Thailand (n = 7), Brazil (n = 3), India (n = 2) and Mauritania (n = 1). All sequences were aligned and cut to 4095 bp (1087-5181 bp) by MEGA7.0. An analysis of the polymorphism of pvmrp1 within the 166 global isolates revealed low nucleotide diversity (π = 0.0011) and high haplotype diversity (Hd = 0.9290). We also found significant differences in the gene polymorphism of pvmrp1 in different regions. Get rid of Africa, the polymorphism of pvmrp1 gene was the highest in South Asia (π = 0.0015), and the lowest in Southeast Asia (π = 0.0006). These re-sults suggested that pvmrp1 in different regions was subjected to different natural selection pressure and showed different levels of gene polymorphism (Table 4).

Natural Selection of Polymorphic Region of pvmrp1 from Different P. vivax Isolates
To determine whether natural selection promoted the generation of pvmrp1 gene diversity in global isolates, we calculated the ratio of non-synonymous to synonymous mutations (dn/ds). The dn/ds for pvmrp1 of all the 166 isolates was 0.5536, indicating that the pvmrp1 gene was affected by purifying selection. The overall Tajima's D test value for pvmrp1 of all the isolates was negative (Tajima's D = −1.1863, p > 0.1), and the Fu and Li's D* test of all the isolates was −3.9871 (p < 0.02) ( Table 4). It suggested that a neutral model of polymorphism occurrence with values for pvmrp1 was due either to a recent population expansion or genetic hitchhiking. In addition, we found that the significant Fu and Li's D* test value < 0 (0.01 < p < 0.05) were found in isolates from South America and Southeast Asia.

Genetic Differentiation, Haplotype Network and LD Analysis of Polymorphic Region of pvmrp1
The level of genetic differentiation of pvmrp1 was estimated by F ST values. As there is only one isolate from Africa, it was not considered. A low level of genetic differentiation was found in the isolates from South America and South Asia (the value of F ST was 0.0809), while others were in a moderate and high level of genetic differentiation (0.1661-0.5498). This was especially true of the southeast Asian isolates, which showed a great genetic differentiation based on the values of F ST (Table 5). To further verify the differential selection between groups, the 11 non-synonymous SNPs (E787D, Q906E, G949D, C1018Y, L1207I, L1287I, Y1393D, G1419A, V1478I, T1525I, and H1586Y), which occurred twice without deletion, were selected to construct a haplotype network among all 166 samples. The haplotype analysis of pvmrp1 in this study revealed 22 distinct haplotypes (Figure 3). Eight of the haplotypes were singleton haplotypes, of which H_4 and H_5 were found in East Asian isolates exclusively, H_10 and H_11 existed only in South Asian isolates, H_12 was found in Southeast Asia, and H_17, H_18, and H_19 were found in South America specifically. The mutant types H_2: EEGCLLDGITY (18.7%) and H_7: EEGCLLDGITH (18.1%) were the most common. H_2, H_3 (DQDCLLDGITH), and H_13 (EEGCLLDAITH) were the dominant haplotype in East Asia (81.8%, 100.0%, and 69.6%, respectively), H_14 (DQGCLLDGVTH) is mainly distributed in South America (63.6%), and H_7 is distributed worldwide, which mainly includes East Asia (40%), South America (33.3%) and Oceania (16.6%) (Figure 4, Table S3). In particular, 14 haplotypes found in the isolates from South America, which indicated the highest haplotype diversity, and nine haplotypes in East Asia, which was consistent with the results of haplotype diversity in Table 4. In addition, by pairwise LD of the pvmrp1 gene at 11 non-synonymous SNPs using the DNASP v5.10.01 program, we observed that there was a strong LD between E787D and Q906E, G949D, and H1586Y (p < 0.0001) ( Figure 5). Meanwhile, the strong LD was present in pvmrp1 gene between G949D and H1586Y, Q906E (p < 0.0001), and existed in the following pairs also: Q906E/L1207I, V1478I/T1525I (p < 0.0001). We found that most of the genes with LD occurred around the ABC transporter domain, such as E787D and Q906E, G949D.
Furthermore, phylogenetic analysis of the 166 isolates indicated that the genetic evolution of pvmrp1 varies in different regions, which have obvious population genetic structures related to geographical isolates. The genetic differentiation of P. vivax isolates indicated that the geographically distant was high; for example, the genetic differences between pvmrp1 gene from East Asia and South America were the highest ( Figure 6).

Discussion and Conclusions
Plasmodium vivax is the most geographically widespread cause of human malaria. Due to the lack of an in vitro culture and transgenic system, molecular markers represent a more practical tool to monitor the introduction and spread of drug resistance in P. vivax [31]. In our study, as a member of ABC transporters, PvMRP1, localizes at the parasite plasma membrane, and the predicted primary protein structure includes 11 TM and 2 NBDs, which function primarily as drug transports [32]. Consistent with the potential involvement of PvMRP1, it has been observed as a transporter with a broad range of substrates, including important endogenous substances such as glutathione and a lot of drugs with diverse structures. CQ was identified as a substrate for this ABC transporter, and the mutations of pfmrp1 were found to be associated with reduced susceptibility to CQ and also to QN [13]. In addition, pvmrp1 transcription level is decreased in the trophozoite stage, which is consistent with the phenotype that P. vivax trophozoites insensitive to CQ [33]. All of these suggest that pvmrp1 may be a potential molecular marker of drug resistance.
In this study, we amplified and sequenced the pvmrp1 gene of 102 whole blood samples originating from Asia, and 29 genetic mutations were found out of the 94 successfully sequenced isolates. Among them, most mutations have not been hitherto reported, and most of these mutations are found only in East Asia, which is likely due to the low number of east Asian isolates used in previous studies. The mutations T259R (97.87%), Y1393D (97.87%) and V1478I (95.74%) were approaching fixation in the sequenced samples, and the mutations V879 (89.5%) and L1207I (89.5%) were highly prevalent and present in Southeast Asia exclusively. A small number of imported isolates, such as South Asia (n = 2) and Oceania (n = 2) were not considered. We found significant differences in the types of mutations prevalent in East and Southeast Asia except for the mutations approaching fixation. A sequence similarity analysis of pvmrp1 and pfmrp1 indicated that the Y1393D and G1419A mutations of pvmrp1 overlap with pfmrp1 locations residing between the ABC transmembrane and the second ABC transporter domains which were associated with drug resistance [10]. Furthermore, a recent study has provided evidence that G1419A and V1478I had a significant association with the IC 50 to CQ and artesunate, and G1419A was also associated with the decreased susceptibilities to PPQ, MQ, and QN [34]. However, in combination with our previous in vitro drug susceptibility studies [20,21], we could not find a correlation between the polymorphisms of pvmrp1 and antimalarial drug susceptibilities (p > 0.05). The isolates used in our previous drug susceptibility studies, Y1393D and V1478I, were fixed, and no G1419A was observed in the pvmrp1 gene, which may limit the correlation analysis. Moreover, we found that the prevalence of S354N between the two CQ sensitivity categories revealed a significant difference by chi-square analysis, although more clinical samples are needed to confirm this conclusion. Numerous studies have shown that gene CN polymorphism is related to genetic and phenotypic variation, which is as important as SNPs [35,36]. Three CNVs were also determined from two ROK isolates and one Pakistan isolate in this study, which indicates the variation in pvmrp1 gene amplification in Asian isolates. However, we did not find a statistically significant correlation between CNVs and drug sensitivity.
The 94 sequenced samples combined with the 72 known sequences from the PlasmoDB database formed pvmrp1 sequences worldwide. The result displayed a low genetic diversity at the pvmrp1 with an average π of 0.0011. Haplotype analysis of pvmrp1 showed that haplotype diversity varies in different regions, and haplotype polymorphism was higher in South America and Oceania than in the other areas. These results suggested that pvmrp1 of P. vivax from different areas was subjected to various natural selection pressures and showed different levels of gene polymorphism and haplotype polymorphism. In this study, the rates of non-synonymous mutation to synonymous mutation of pvmrp1 were < 1, and the neutral evolutionary test statistic Tajima's D test was negative (Tajima's D = −1.1863, p > 0.1), indicating that the pvmrp1 gene was affected by purifying selection [27]. Fu and Li's D* test further confirmed that the pvmrp1 gene was under purifying selection (Fu and Li's D* = −3.9871, p < 0.02), which suggested that the mutations and accumulate at silent sites, there were likely to be lots of segregating sites, but not much heterozygosity. This could explain that nucleotide diversity (π) was small and the average number of nucleotide differences (K) was high in Table 4.
Furthermore, genetic differentiation, the haplotype network, and phylogenetic analysis of the 166 isolates indicated that the genetic evolution of pvmrp1 varies in different regions, which has obvious population genetic structures related to geographical isolates. The genetic differentiation between P. vivax isolates that were geographically distant was high, for example, the genetic differences between the pvmrp1 gene from Southeast Asia and South America were the highest. Perhaps parasite movement could be controlled by different factors, including geographical barriers, distance, poor road infrastructure, cultural and language barriers, and the effectiveness of malaria control interventions.
Overall, PvMRP1 localizes on the P. vivax plasma membrane with 11 TM, the gene nucleotide polymorphism, and CNVs were analyzed with the field P. vivax isolates in our study. We found some unreported point mutations of pvmrp1 in the collected samples, and S354N substitution may lead to CQR in P. vivax. In combination with the worldwide isolates, the analysis showed that the pvmrp1 gene had a high haplotype diversity, low nucleotide diversity and was under purifying selection. This study showed the polymorphisms of the pvmrp1 gene from worldwide isolates, and pointed to the potential contribution of the pvmrp1 gene in the CQR of P. vivax. However, the lack of correlation between the pvmrp1 polymorphisms and the IC 50 values of the four antimalarial drugs (CQ, MQ, QN and PYR), highlights the need for more informative tools to function the role of pvmrp1 gene.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/microorganisms10081482/s1, Table S1: IC 50 values to four antimalarial drugs of P. vivax isolates; Figure S1: Dot plots of in vitro susceptibilities of P. vivax isolates to four antimalarial drugs; Table S2: Association of SNPs in pvmrp1 with in vitro susceptibilities to CQ, MQ, QN and PYR; Figure S2: Comparison of IC 50 values for three antimalarials between parasites with single copy and multicopies of pvmrp1 gene; Table S3: Distribution of the number of 22 pvmrp1 haplotypes in worldwide isolates; Spreadsheet S1: Sample information.