Polymorphism of Antifolate Drug Resistance in Plasmodium vivax From Local Residents and Migrant Workers Returned From the China-Myanmar Border

Drug-resistant Plasmodium vivax malaria impedes efforts to control, eliminate, and ultimately eradicate malaria in Southeast Asia. P. vivax resistance to antifolate drugs derives from point mutations in specific parasite genes, including the dihydropteroate synthase (pvdhps), dihydrofolate reductase (pvdhfr), and GTP cyclohydrolase I (pvgch1) genes. This study aims to investigate the prevalence and spread of drug resistance markers in P. vivax populating the China-Myanmar border. Blood samples were collected from symptomatic patients with acute P. vivax infection. Samples with single-clone P. vivax infections were sequenced for pvdhps and pvdhfr genes and genotyped for 6 flanking microsatellite markers. Copy number variation in the pvgch1 gene was also examined. Polymorphisms were observed in six different codons of the pvdhps gene (382, 383, 512, 549, 553, and 571) and six different codons of the pvdhfr gene (13, 57, 58, 61, 99, 117) in two study sites. The quadruple mutant haplotypes 57I/L/58R/61M/117T of pvdhfr gene were the most common (comprising 76% of cases in Myitsone and 43.7% of case in Laiza). The double mutant haplotype 383G/553G of pvdhps gene was also prevalent at each site (40.8% and 31%). Microsatellites flanking the pvdhfr gene differentiated clinical samples from wild type and quadruple mutant genotypes (F ST= 0.259-0.3036), as would be expected for a locus undergoing positive selection. The lack of copy number variation of pvgch1 suggests that SP-resistant P. vivax may harbor alternative mechanisms to secure sufficient folate.


INTRODUCTION
Plasmodium vivax, the most broadly distributed source of severe malarial morbidity and mortality, puts almost 40% of the world population at risk (Naing et al., 2014;Rahimi et al., 2014;WHO, 2020), resulting in cases that can be life threatening (Gupta et al., 2015;Gupta et al., 2016). Data from 2017 shows that about 3.3 billion people live in the P. vivax transmission areas, and about 1.5 billion people are at risk of exposure in foci of stable transmission (Battle et al., 2019). Southeast Asian countries shoulder the highest proportion of global burden, with about 7.3 million people infected with P. vivax in 2017, alone (Battle et al., 2019).
Asia is the focus of P. vivax transmission, and Myanmar is among the most severely affected in the region. As one of the most severe malaria burden countries from this area, Myanmar is a source of cases exported to neighboring countries such as mainland China, where low endemic transmission is otherwise leading to successful eradication. From 2010 to 2014, there were 1893 P. vivax malaria cases imported from Myanmar into China, far exceeding importation from any other Asian country . Most imported cases are adult male Chinese laborers working in farming, construction, and mining . Reported cases in Myanmar have risen steadily since 2012 (WHO, 2016).
Most cases imported to China occur in Yunnan province, the southwestern frontier which borders Myanmar to the west, and Laos and Vietnam to the south (Lai et al., 2019). In particular, Tengchong County (in Yunnan, China) reported the highest number of imported malaria cases from Myitsone, Myanmar from 2010 to 2014 (Wang et al., 2015). Previous studies identified boundary townships, where roads span the China-Myanmar Border, facilitate the spread of P. vivax Xu et al., 2016).Therefore, the area merits significance for global malaria control.
The recommended first-line medication for vivax malaria patients in Myanmar is a radical cure with chloroquineprimaquine (CQ-PQ) (Htun et al., 2017). However, most local cases of malaria derive from Plasmodium falciparum (or from mixed infections with P. falciparum and P. vivax). Although the prevalence of mixed-species malaria infections was previously estimated at <2% in some surveys of Asian locales, therapeutic studies identify coinfection in over 30% of cases (Mayxay et al., 2004). The incidence of P. vivax infection after falciparum malaria treatment exceeds that which would be expected from entomological inoculation rates in Southeast Asia (Douglas et al., 2011). A survey from 2008 to 2012 in Yunnan found 17% of malaria cases were co-infections with P. vivax and P. falciparum, underscoring the difficulty in diagnosing low-level and coinfection cases (Zhou et al., 2014).
Given such high rates of coinfection, treatments meant to target chloroquine-resistant P. falciparum likely impose substantial selective pressure for resistance to sulfadoxinepyrimethamine (SP) in the local P. vivax population, and may induce high-grade antifolate resistance in P. vivax. SP has been the recommended treatment for falciparum and other non-vivax malaria since the 1980s, replacing chloroquine (CQ) as the front-line anti-malaria treatment when large-scale CQ resistance developed in many countries; but growing drug resistance soon necessitated replacement of SP with artemisinin-based combination therapy (ACT) (Jiang et al., 2020). In addition to typical treatment proscribed by WHO guidelines (WHO, 2013), SP is also used for intermittent preventive treatment in infants (IPTi) and pregnant women (IPTp) in areas endemic for malaria.
Meanwhile, China developed its own strategy of malaria control: seasonal mass drug administration with pyrimethamine and primaquine. This approach attempts to treat carriers of P. vivax before seasonal transmission begins (Hsiang et al., 2013). These several uses of such antimalarials near the China-Myanmar Border necessitate continuous monitoring of antifolate resistance in P. vivax in this area. Given growing concerns of drug resistance, it is essential to monitor antifolate resistance in P. vivax in this area, which has undergone changes in transmission intensity and drug therapy.
Antifolate resistance derives, in part, from point mutations of enzymes used by the parasite to synthesize folate. These include dihydropteroate synthase (DHPS) and dihydrofolate reductase (DHFR). Barriers to continuous culture of P. vivax severely limit efforts to link pvdhfr and pvdhps genotypes with parasite responses to SP in vitro. Instead, investigators assay drug sensitivity by expressing P. vivax DHPS and DHFR in E. colior yeast (Auliff et al., 2010). In P. falciparum, SP resistance has also been attributed to amplification of GTP-cyclohydrolase (gch1), the first gene in the folate biosynthesis pathway (Nair et al., 2008); this phenomenon has not yet been demonstrated in P. vivax. Drug resistance was last surveyed from the vicinity (Nabang Town, Yingjiang County, Yunnan Province, and Kachin, Myanmar were surveyed in 2015 and 2017) Zhao et al., 2020); sampling from Myitsone has not yet been reported.
Despite their likely contribution to emerging regional drug resistance, the presence and prevalence of mutations in genes of the folate biosynthesis pathways in P. vivax remain uncharacterized. This study aims to investigate the extent and spread of SP drug resistance alleles in P. vivax from local people in Laiza, Myanmar and from migrant Chinese laborers returning from Myitsone, Myanmar between 2012 and 2015.To do so, we sought signatures of selection by sequencing SP resistance-associated molecular markers and genotyped microsatellite markers around target genes.

Samples
Between 2012 and 2015, finger-prick blood samples were collected from symptomatic patients with acute P. vivax infection. Infections were diagnosed by microscopy using Giemsa-stained thick and thin blood films at 13 local clinics around the Laiza township (N 24°45' 24", E 97°33' 02"), Myanmar, and in returning migrant Chinese laborers from the Myitsone area (N 25°41' 23", E 97°31' 04"), Myanmar ( Figure 1). Patients who provided blood were well-informed by our study team and signed a consent form before sampling. The Institutional Review Board of Kunming Medical University, China, and the local Bureau of Health, Kachin State, Myanmar approved the sampling procedures.
Parasite DNA was extracted from each sample using the High Pure PCR Template Preparation Kit (Roche, Germany). Multiclonal infections were identified by genotyping the singlecopy antigenic genes msp-3a and msp-3b as described previously (Khan et al., 2014); these were excluded from the ensuing analysis, leaving a total of 308 single-clone P. vivax isolates (101 from Myitsone, 207 from Laiza). A sample of fewer than 100 would suffice to achieve 80% power (with a type I error <.05) given reported mutation frequencies of in pvdhfr and pvdhps.
Pvdhfr and pvdhps Gene Sequencing, Pvgch1 Copy Number Quantification, and Microsatellite Genotyping Parasite genes pvdhfr and pvdhps were amplified by nested PCR using Premix Taq Version 2.0 plus dye (TaKaRa Bio, Japan), according to the manufacturer's instructions. Primers and PCR conditions were as described in Supplementary Table S1. PCR products were purified and sequenced by Sanger at SinoGenoMax (Kunming, China). For accuracy, sequences from the forward and reverse were examined by alignment using DNASTAR software (v.7.1). To avoid technical contamination, all PCR and sequencing assays include positive control (MRA-41 Sal I, gDNA, MR4/Bei resources) and negative controls (uninfected samples and water). Mutations were scored only when confirmed by each bidirectional sequencing read with high quality BioEdit (v.7.2.5) (Sarmah et al., 2017) was used to compare individual se quences with the r eference sequences pvdhfr (PVX_089950) and pvdhps (PVX_123230). Nucleotide sequences were translated into amino acid sequences to examine mutant codons.
To determine copy number variation in pvgch1, PCR primers were designed (Supplementary Table S1) and quantitative realtime (qRT-PCR) performed using SuperRealPreMix Plus (SYBR Green) (TIANGEN, China) on the Stratagene Mx3000P PCR machine (Stratagene, CA, USA), according to the manufacturer's instructions. The threshold cycle value (Ct value) of pvgch1 was estimated using pvaldolase the reference gene. In particular, we assessed, via real-time PCR, the proportion of mutant alleles of pvgch1 compared with the reference allele. This was calculated using Pfaffl method for each patient sample as amplification efficiency of the target and reference was not similar. Results were recorded as a ratio calculated from the mean Ct values for each gene in each sample. The pUC57 plasmid inserted with a single copy of pvgch1 or pvaldolase was purchased from Sangon Biotech (Shanghai, China) and used for amplification efficiency calibration.

Data Analysis
Population genetic parameters, pairwise linkage disequilibrium (LD), and tests of neutrality were estimated using DNASP software (v.5.0.) (Mardani et al., 2012). The copy number of pvgch1 in each sample was estimated by the Pfaffl method using the Ct value, rounded to the nearest integer (Raju et al., 2019).
The Chi-square test and Fisher's exact test were employed to compare the independence of categorical variables, using GraphPad Prism 9.0.0 software. A statistical significance was set at the P-value of < 0.05. The length polymorphism of microsatellite loci and allele frequencies were calculated using GenAlEx software (v.6.5.) (He et al., 2018). Expected heterozygosity (H e ) was estimated for each locus based on allele frequencies following previously reported methods (Anderson et al., 2000). The sampling variance for H e was estimated using ARLEQUIN software (v.3.5.), allowing estimation of genetic diversity indices and assessment of population differentiation between haplotypes. The Bottleneck program (Cornuet and Luikart, 1996) was used to investigate the influence of selection on allele distributions under the infinite alleles (IMM) and step wise mutation models (SMM). Wilcoxon tests were used to determine the significance of departures from H e .
Nucleotide diversity in pvdhfr was almost twice as great in Laiza than in Myitsone (0.006496 v.s.0.003381); diversity was lower, and more regionally uniform, in pvdhps (0.002263 v.s.0.002037). For the Laiza population, neutrality tests showed significant positive values for Tajima's D in the pvdhfr gene, indicating a dearth of singletons and consistent with balancing selection. However, for the pvdhps gene, we observed negative values for Fu & Li's F and Tajima's D, which indicates an excess of low-frequency alleles more consistent with population expansions or positive selection.
In the pvdhfr gene, Myitsone samples and Laiza samples contained mutations at six codons (13,57,58,61,99,117) ( Table 2). Compared to the pvdhfr reference sequence in Salvador I (PVX_089950), more than ninety percent of samples from each sampling locale had at least one codon change. Mutations in codons F57I/L, S58R, T61M, and S117T/N were highly prevalent in both populations. In the Myitsone samples, the F57I/L, S58R, T61M and S117T/N mutation in pvdhfr was significantly more frequent than in the Laiza samples; the H99R/S mutation in the Laiza samples was significantly more frequent than in the Myitsone samples. The quadruple mutant L 57 R 58 M 61 T 117 was observed significantly more frequent in Myitsone (44.2%) than in Laiza (12.6%) ( Table 3). The quadruple mutant I 57 R 58 M 61 T 117 was observed in both areas (32.6% in Myitsone and 31.1% in Laiza). Genotypes R 99 and S 99 were only in Laiza (37.8%) ( Table 2). Rare genotypes with triple   Table 3).
For the pvdhps gene, nucleotide polymorphisms were observed at five different codons (382, 383, 512, 549, and 553) in samples from Myitsone and five codons (382, 383, 512, 553, and 571) from Laiza. As compared to the reference dhps sequence in Salvador I (PVX_123230), more than ninety percent and seventy percent of the samples had at least one codon change in Myitsone and Laiza samples, respectively. The most prevalent mutations were at codon A383G and A553G (Table 2). Similarly, these two codons also were more prevalent in the Myitsone population (84.7 and 65.3%, respectively) than in the Laiza population (74 and 38.1%, respectively). Mutation S382A also reached a high frequency (20.4%) in the Myitsone population. In the Myitsone samples the S382A/C and A553G mutation in pvdhps was significantly more frequent than in the Laiza samples; the E571Q mutation in the Laiza samples was significantly more frequent than in the Myitsone samples.
The parasites in Mytsome and Laiza each possessed certain unique pvdhps haplotypes (Table 3). For instance, parasites carrying the double mutations G 383 Q 571 were identified only in Laiza, where they reached a sample prevalence of 6.3%. The double mutation A 382 G 383 in pvdhps were 8.2% in Myitsone but not found in Laiza. Also, genotypes with triple mutations G 383 M 512 G 553 (in 2 cases, 2%), G 383 T 512 G 553 (in 3 cases, 3.1%), and G 383 D 549 G 553 (in 1 case, 1%) were only found in Myitsone. However, the genotype with G 383 E 512 G 553 were only found in 1 case (0.7%) in Laiza. Multiple mutations 382C/383G/512E/553G were found in 4 cases (3%) from Laiza. The distribution of major haplotypes of the pvdhfr and pvdhps genes in P. vivax from the two study areas are listed in Table 3.

Genetic Hitchhiking of Microsatellites Flanking the pvdhfr Gene
Polymorphisms of six microsatellite loci flanking pvdhfr (located at -93 kb, -38 kb, -2.6 kb, +4.9 kb, +37 kb, and +94 kb from the pvdhfr gene on chromosome 5) were examined in 100 samples from Laiza and 101 from Myitsone. Genotypes were obtained from 81 to 99 of the Myitone samples at these six loci, and from 65 to 97 of the Laiza samples ( Table 4).
The mean number of observed alleles was higher in Laiza (10.2 ± 4.2) than in Myitsone (6.3 ± 2.3). A significant reduction in variability (compared to wild-type isolates) occurred across at the +4.9 kb region in Laiza samples and over a 7.5 kb region (spanning -2.6 to +4.9 kb) in Myitsone samples (P<0.01, Tukey's multiple comparisons test) (Figures 2 and 3).
Tests using Bottleneck software, assuming the stepwise mutation model (SMM), indicated that parasites with the haplotype 57I/L-58R-61M-117T in the Myitsone population had  Table S4). However, significant heterozygosity excess (P<0.05) was found in parasites with genotype S58R-S117N in the Myitsone population, indicating that they could be suffering a bottleneck. There was no bottleneck effect detected in parasites with wild-type genotypes in either population when either an infinite allele (IMM) or stepwise mutation model (SMM) was assumed, suggesting parasite population under mutation-drift equilibrium. Significant genetic differentiation was observed between the wild type and the quadruple mutant genotypes 57I/L/58R/61M/117T (F ST =0.25492 for 57I/L/58R/61M/117T, P<0.00, F ST =0.30369 for 57I/L/58R/ 61M/117T, P<0.001) (as shown in Table 5).

DISCUSSION
An increasing number of drug resistance markers in P. vivax have spread throughout Southeast Asia (Lu et al., 2010;Thongdee et al., 2013;Asih et al., 2015). These include polymorphisms at codon 373, 380, 382, 383, 384, 512, 553, 585, and 601 in the pvdhps gene, which are associated with resistance to sulfadoxine, and mutations at codon 15, 33, 50, 57, 58, 61, 64, 117, and 173 in the pvdhfr gene, contributing to treatment failure using pyrimethamine. We found polymorphisms in five different codons of the pvdhfr gene (13,57,58,61,117) in each locale; mutations in a sixth codon (99) were restricted to Laiza. These locales each had polymorphisms in codons 382, 383, 512, 553 of the pvdhps gene; observed polymorphisms in codon 549 were restricted to Myitsone, and in codon 571 were restricted to Laiza. These results were consistent with those reported previously in Southeast Asia (Lu et al., 2010;Thongdee et al., 2013;Asih et al., 2015). Though most of the pvdhfr and pvdhps haplotypes identified in this study have been reported in previous studies, some new haplotypes were also detected. For example, 2 cases each from Laiza and Myitsone bore quintuple-mutant 13L/57I/58R/61M/ 117T pvdhfr allele. To the best of our knowledge, this is the first report of this haplotype in the China-Myanmar border region. The S117N mutation has been suggested to represent the first mutational step conferring drug-resistance (Brega et al., 2004). However, as shown in the Table 6, only 5 cases in Myitsone bore the single S117N mutation, and the quadruple pvdhfr 57L/58R/61M/117T allele was the most prevalent; this differs from Yangon, Myanmar, where the double 58R/117N allele was most common (Lu et al., 2010).
Studies in China and Myanmar have found the prevalence of pvdhfr and pvdhps drug-resistant mutations in P. vivax in regions of Yunnan province, China, including the Xishuangbanna prefecture and areas along the Nu River, and these studies identified parasites highly resistant to SP in subtropical China and Yangon, Myanmar from 2010 to 2014 (Lu et al., 2010;Ding et al., 2013;Huang et al., 2014). The quadruple-mutant 57L/58R/61M/117T allele was the most prevalent pvdhfr allele in two study sites in the Yunnan province. However, there were no mutations at positions 13 or 173 in the pvdhfr gene reported in parasites from the Xishuangbanna. An altered 173 allele was, however, found in 3.18% of samples analyzed from areas along the Nu River. In addition to this, Yunnan province also reported haplotypes with quintuple allelic variants. Our study profiled mutations in marker genes associated with SP resistance among the P. vivax from the China-Myanmar border. The quadruple-mutant 57I/58R/61M/ 117T allele was also prevalent in Myitsone, Myanmar, and in Yunnan, China (Ding et al., 2013;Huang et al., 2014) but absent from our survey of Laiza and from a previous study in Yangon, Myanmar (Lu et al., 2010).   Beneficial mutations reduce variation in physically proximate regions of the chromosome through "genetic hitchhiking" (McCollum et al., 2008;Artimovich et al., 2015). The polymorphism in six microsatellite loci flanking mutant pvdhfr (S58R-S117N and the quadruple mutant genotypes 57I/L/58R/ 61M/117T) was significantly less than that would be expected in wild type genetic backgrounds. The valley of reduced variation spanning pvdhfr is resembles that found for samples from Yunnan and central China (Ding et al., 2013). In each, an approximately 100 kb region suggests local genomic variability has been suppressed by a recent selective sweep. These results strengthen the evidence for a relationship between the specific mutant alleles and selection for drug drug-resistant phenotypes.
Furthermore, favorable alleles experience increased linkage disequilibrium (LD), at least until sufficient time has elapsed to enable recombination to distribute alleles randomly among haplotypes in the haploid blood stages of Plasmodium parasites. Therefore, the multiple mutations with S117T can serve as a marker for drug resistance at the China-Myanmar Border. Interestingly, 36.4% of Laiza samples bore the single H99S mutation, which has been reported from other studies (Lu et al., 2012;Huang et al., 2014;Wang et al., 2020).
Our results also elucidate the drug pressure exerted by the antifolatesulfadoxine. Previous study has focused on the A383G and A553G in pvdhps gene, which are directly related to sulfadoxine resistance (Korsinczky et al., 2004). All of the mutant cases from our study carried at least one of these two mutations. The most prevalent pvdhps haplotype we observed was the double mutant 383G/553G allele. There are additional mutant alleles found in Yunnan, including one non-synonymous mutation at codon 512T/ M and triple mutant alleles at codons 372, 495, 561. We observed the triple mutant 382A/383G/553G allele and the double mutant 382A/ 383G allele in Myitsone, Myanmar, and Yunnan, China, but not in Laiza. Nor were these previously observed in Yangon, Myanmar (Lu et al., 2010). We observed the double mutant 383G/571Q only in nine Laiza cases. Besides the common triple mutant 382A/383G/ 553G allele described earlier, there are other types of triple mutant pvdhps alleles in Laiza and Myitsone. They included one case from Laiza with a 383G/512E/553G allele, one from Myitsone with a 383G/549D/553G allele, two from Myitsone with a 383G/  512M/553G allele, and three from Myitsone with a 383G/512T/ 553G allele. Moreover, we observed the quadruple-mutant 382A/ 383G/512E/553G allele in four Myitsone cases. SP combination therapy targets two enzymes in the folate-synthesis pathway, dihydropteroate synthase (DHPS) and dihydrofolate reductase (DHFR). The first enzyme in this pathway, GTPcyclohydrolase 1 (GCH-1), might also respond to drug pressure. Previous work on P. falciparum (Kümpornsin et al., 2014) suggested a compensating effect, by copy-number variation (CNV) of gch-1, in parasites bearing resistance mutations in dhps and dhfr. However, such a relationship has not been established in P. vivax.Our study, the first to explore such a relationship in P. vivax, identified little evidence of copy number variation in pvgch1, suggesting that SP-resistant P. vivax may harbor alternative mechanisms to secure sufficient folate. Although microsatellite analysis showed that drug resistance genes experience selective pressure, pvgch1 copy number variation was only observed only in a small number of cases, suggesting that P. vivax population might not undergo pvgch1 compensation.
The first-line medication for vivax malaria patients in Myanmar is CQ-PQ. Currently, the molecular markers of primaquine resistance are still unknown. And no molecular markers are confirmed for chloroquine resistance in P. vivax, while pvmdr1 and pvcrt-o have often been queried since these two genes are involved in chloroquine resistance in P. falciparum, but the possible role of Pvcrt-o and pvmdr1 in chloroquine resistance is controversial. At this area, the pvmdr1 and pvcrt-o genes were reported .The limit of detection of the PCR and sequencing assays used in the study may rule out the possibility of low-frequency mutation isolates. If we want to get the results, deepsequence needed to be involved in further study.

CONCLUSION
We found that most P. vivax sampled from symptomatic patients at the China-Myanmar boarder harbor mutations in pvdhfr and pvdhps genes associated with drug resistance. In this region of intense transmission, which is the source of many cases of vivax malaria introduced to China, pvdhfr, a single mutation (S117T) serves as a marker for the presence of multiple mutations. All of the mutant cases from our study carried mutation A383G or A553G in this gene. Parasites throughout the studied region show genetic evidence of undergoing strong drug pressure. The lack of copy number variation of pvgch1 suggests that SP-resistant P. vivax may harbor alternative mechanisms to secure sufficient folate.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GenBank (pvdhps gene accession numbers: MZ234760-MZ234999 and pvdhfr gene accession numbers: MZ235000-MZ235236).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Kunming Medical University. The patients/ participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
SF, SW, and WZ performed the statistical analysis and wrote the first draft of the manuscript. DZ, YH, YB, YR, YS, HZ, QY, XL, XC, YZ, CL, ZX, YW, FC, PS, and BR wrote sections of the manuscript. ZY contributed to conception and design of the study. ZY organized the database. All authors contributed to the article and approved the submitted version.

FUNDING
This study was under the support of the National Science Foundation of China (31860604 and U1802286), Major science and technology projects of Yunnan Province (2018ZF0081) and International Science and Technology Cooperation-Yunnan International Science and Technology Cooperation Base (202003AE140004. Furthermore, XC and YZ were under sponsoring from the Yunnan Applied Basic Research Projects-Union Foundation (2018FE001-190, 2015FB034, respectively). Also, WZ was under the support of the Education Department Fund of Yunnan Province (2019J1184). BR was supported by USDA Project Detection and Control of Foodborne Parasites for Food Safety (8042-32420-007-00D).