Plasmodium falciparum phenotypic and genotypic resistance profile during the emergence of Piperaquine resistance in Northeastern Thailand

Malaria remains a public health problem in Thailand, especially along its borders where highly mobile populations can contribute to persistent transmission. This study aimed to determine resistant genotypes and phenotypes of 112 Plasmodium falciparum isolates from patients along the Thai-Cambodia border during 2013–2015. The majority of parasites harbored a pfmdr1-Y184F mutation. A single pfmdr1 copy number had CVIET haplotype of amino acids 72–76 of pfcrt and no pfcytb mutations. All isolates had a single pfk13 point mutation (R539T, R539I, or C580Y), and increased % survival in the ring-stage survival assay (except for R539I). Multiple copies of pfpm2 and pfcrt-F145I were detected in 2014 (12.8%) and increased to 30.4% in 2015. Parasites containing either multiple pfpm2 copies with and without pfcrt-F145I or a single pfpm2 copy with pfcrt-F145I exhibited elevated IC90 values of piperaquine. Collectively, the emergence of these resistance patterns in Thailand near Cambodia border mirrored the reports of dihydroartemisinin-piperaquine treatment failures in the adjacent province of Cambodia, Oddar Meanchey, suggesting a migration of parasites across the border. As malaria elimination efforts ramp up in Southeast Asia, host nations militaries and other groups in border regions need to coordinate the proposed interventions.

a population at risk for succumbing to multi-drug resistant malaria as well as facilitating transmission, necessitating targeted control strategies if the Greater Mekong Subregion (GMS) is to achieve malaria elimination.
Efforts in mapping parasite population structure and gene flow can assist in understanding and predicting the spreading pattern of resistance. Military populations at border areas who are deployed to endemic areas should be integrated in surveillance and monitoring efforts as they are on the front-lines of malaria transmission. Here we report resistance characteristics of P. falciparum clinical isolates collected from Thai soldiers and civilian patients between 2013 and 2015 and show the trends in drug resistance. At that time, the national treatment guidelines in Thailand recommended AS-MQ for P. falciparum malaria. During this same period, we reported dramatic loss of efficacy of DHA-PIP (54% treatment failure) in a study that was conducted by AFRIMS in Anlong Veng, 12 km from the Thai border of Sisaket Province 40 . Intensive malaria surveillance to track drug resistance is required in high risk military and other border populations to achieve malaria control.

Materials and methods
Study setting, protocol and subjects. This minimal risk surveillance study was open between July 2013-September 2015, enrolling Thai adults (aged 18 years and over) presenting with uncomplicated P. falciparum or P. falciparum/P. vivax mixed-infections at military health facilities. The enrollment population included soldiers, border police, or their family members and other villagers located near the Royal Thai Army (RTA) health clinics in Sisaket and Surin Provinces, located in northeastern Thailand on the Thai-Cambodia border. Inclusion criteria were asexual parasitemia per rapid diagnostic test (RDT) or blood film, and no malaria infection or antimalarials taken within the past seven days. The protocol was approved by the Walter Reed Army Institute of Research Institutional Board, Institute for Development of Human Research Protection, and RTA Institutional Review Board. All research was performed in accordance with relevant guidelines and regulations. All study subjects provided written informed consent prior to participation. Goal enrollment was 50 malaria cases per year, balancing the number of isolates needed to characterize resistant genotypes/phenotypes with the local population size and expected incidence of cases that could be enrolled. While treatment was uniformly provided per national guidelines at the time, volunteers were not followed up for treatment efficacy.
Sample collection and preparation. Patients diagnosed with P. falciparum infection at RTA clinics were subjected to peripheral venipuncture prior to treatment. Two microscopists examined Giemsa-stained peripheral blood smears for each volunteer to determine malaria species infection and parasite densities for blood stages. Venous blood samples were collected in EDTA tubes for DNA extraction and molecular characterization and in sodium-heparin tubes for ex vivo bioassay and in vitro drug sensitivity assay. Plasma was separated from blood, frozen at − 80 °C, and infected packed red blood cells were cryopreserved. All blood and processed blood samples were stored at − 80 °C and transported in dry ice to AFRIMS for molecular characterization, ex vivo bioassay, and in vitro culture adaptation and drug sensitivity testing.
Molecular markers of malaria drug resistance. Parasite genomic DNA was extracted from 200 µL of EDTA whole blood by using EZ1 DNA blood kits with an automated EZ1 Advanced XL purification system as per the manufacturer's instructions. The DNA was stored at − 20 °C. T100TM Thermal Cycler (Bio-Rad Laboratories) was employed to evaluate resistance makers including, P. falciparum kelch13 propeller domain (pfk13) (PF3D7_1343700), pfcrt SNP F145I (KM288867.1), and pfcytb SNPs (AY282930.1) 10,41 . ABI 7500 Real-time PCR system (Applied Biosystems) was employed to characterize pfcrt SNPs (72)(73)(74)(75)(76)  www.nature.com/scientificreports/ Copy number variation assay. To determine copy numbers of pfmdr1 and pfpm2 genes, real-time quantitative PCR (qPCR) was performed on genomic DNA as previously described 15,29,42 with some modifications. For pfmdr1, amplification reactions were performed according to the TaqMan Real-time PCR methods using ABI 7500 Real-time PCR system (Applied Biosystems) with 200 nM of each forward and reverse primer (Table S2) and 2 ng of DNA template while Rotor-Gene Q (QIAGEN, Valencia, CA) using Type-it® HRM™ kit was employed for pfpm2 29 . The primers and probes (Table S2) used were as previously described to amplify the following loci: pfmdr1 (PF3D7_0523000) and pfpm2 (PF3D7_1408000), respectively 29 . For the housekeeping gene, β-tubulin (PF3D7_1008700), β-tubulin forward and reverse primers were designed and used as a reference control for all experiments with the same validated PCR conditions as target primers. P. falciparum 3D7 and Dd2 were used as references for single and multiple copy numbers of pfmdr1, respectively. All samples including the references clones were performed in duplicate. The average copy number values for each genes were calculated using 2 −ΔΔCt method. Parasites with copy number greater than 1.5 copies for pfmdr1 15 and 1.6 copies for pfpm2 29 were interpreted to contain multiple copies of each gene.
In vitro culture adaptation and maintenance of parasites. Of 112 collected samples, 86 samples were attempted for in vitro culture adaptation. The exclusion criteria of the in vitro culture adaptation were P. falciparum/P. vivax mixed infections, % parasitemia < 0.05, and ex vivo bioassay activity > 250 nM (DHA equivalent). Culture adaptation of the parasites was performed using the modification method of Trager and Jensen 43 . Parasites were maintained in fresh human erythrocytes (O + ) in RPMI-1640 medium (Sigma), containing 15% AB + human serum (heat inactivated and pooled) and supplemented with 25 mM HEPES, 25 mM sodium bicarbonate, and 0.1 mg/mL gentamicin. Human blood products (erythrocytes and serum) were obtained from the Thai Red Cross. Culture flasks were gassed with 5% CO 2 , 5% O 2 , 90% N 2 and incubated at 37 °C. After culture adaptation the cultured parasites were cultivated for 3 cycles until enough material was obtained and they were synchronized with 5% D-sorbitol (Sigma) to obtain at least 70% ring forms before in vitro assays were run.
In vitro 72 h drug susceptibility by Histidine rich protein 2 (HRP2). Drug susceptibility test using HRP2-ELISA to measure 50% or 90% inhibitory concentration (IC 50 and IC 90 ) was performed following previously published methods [44][45][46] . Dried drug-coated plates containing antimalarial drugs as described in Chaorattanakawee et al. 36,45 were used and in vitro drug susceptibility testing was carried out for control reference clones (W2, D6, C2B) (Malaria Research & Reference Reagent Resource, Manassas, Vermont, USA) as described previously 36 . IC 50 s and IC 90 s were estimated by nonlinear regression analysis using GraphPad Prism version 6.0.

Ring-stage survival assay (RSA).
In vitro RSA 0-3 h was performed on 0-3 h post-invasion rings obtained from selected culture-adapted clinical isolates following published methods with slight modifications 14,46 . Parasite cultures were synchronized using 5% D-sorbitol and 75% Percoll to obtain 0 to 3-h post-invasion rings which were adjusted to 0.5-1% starting parasitemia with a 2% hematocrit in culture media (0.5% Albumax RPMI 1640 with 2.5% AB serum), and cultured in a 48-well microplates with 700 nM DHA and 0.1% DMSO in separate wells. The culture plates were then incubated for 6 h at 37 °C in modular incubator chambers and gassed with 5% CO 2 , 5% O 2 and 90% N 2 . Cells were then washed once, resuspended in a drug-free medium, and cultured further for 66 h. Susceptibility to DHA was assessed microscopically on thin films by estimating the percentage of viable parasites, relative to control (% survival rate). For the controls, the RSA 0-3 h was also performed on P. falciparum reference clones W2, IPC-4884 and IPC-5202 (BEI Resources, NIAID, NIH, USA). A survival rate > 1% was deemed resistant for RSA.
Ex vivo bioassay. P. falciparum-based bioassay was carried out to measure the antimalarial activity of patient plasma identifying if study volunteers were likely to have recently taken antimalarial drugs. Plasma was prepared from blood collected on the screening day for evaluation according to previously published methods 47,48 . In brief, the complement-inactivated samples were serially diluted and applied in one column to 96-well microplate at 50 µL/well. Two columns with serial dilutions of spiked plasma were added to each plate as controls. In addition to the plates with unknown samples, one plate was dosed with six serial dilutions in duplicate of DHA (from 100 to 2.5 ng/mL). A suspension of 200 µL of malaria parasite-infected red blood cells (W2 clone; 0.05% parasitemia with > 80% rings at a 1.7% hematocrit) was added to all wells. The microplates were placed into a chamber, flushed with a mixture of gas consisting of 5% CO 2 , 5% O 2 and 90% N 2 , and plated into an incubator at 37 °C for 72 h. HRP2-ELISA was performed as described above.

Liquid chromatography-tandem mass spectrometry (LC-MS/MS) analysis. To detect baseline
antimalarials prior to treatment in the study population, plasma samples were extracted by using protein precipitation twofold volume of acetonitrile containing internal standard, then 1-min vortex mixing, 10 min-centrifugation, filtration of supernatant with a 0. 22

Results
Study population, demographic, and parasitological parameters. In total, 117 individuals were enrolled but 5 individuals were excluded from analysis due to lack of P. falciparum on a blood smear in 4 patients, and P. vivax monoinfection in 1 patient. Therefore, 112 individuals with uncomplicated P. falciparum were included in the analysis ( www.nature.com/scientificreports/ Haplotype and copy number variation (CNV) of P. falciparum isolates. The parasite isolates were categorized into nine groups according to their genotypes of pfmdr1, pfk13, and pfcrt in combination with their CNV of pfmdr1 and pfpm2 (Table S3). Overall, the most prevalent parasites were those in group III, containing pfk13-580Y alleles, multiple pfpm2 copies, pfcrt-F145 alleles and pfmdr1 184F alleles with a single copy number. This was followed by parasites in group II which was the same as Group III except with pfpm2 single copy number.
The number of parasites in group II decreased from 55.6% to 13%, while those in group III increased from 7.4 to 56.5% over time. Parasites in group VI, containing the pfcrt-145I and multiple pfpm2 copies, were also found to be increasing over the study period. It was noted that parasites in group VII to IX, harboring the pfk13-539I/T alleles with no novel mutation on pfcrt, and pfmdr1-184F only held a single pfpm2 copy.

Prevalence of molecular markers for ART-, MQ-and PIP-Resistance. With limited number of
ACT options, we assessed if any parasites had all mutations for ART-, MQ-, and PIP-resistance (pfk13 SNPs, pfmdr1 CNV, pfpm2 CNV, and pfcrt-F145I mutation). No tested parasites in this study carried all the aforementioned markers, although this is largely driven by having only 5 isolates with multiple copies of pfmdr1. If that marker is excluded, only 11.6% of the isolates (13/112) harbored pfk13-C580Y, multiple pfpm2 copies, and pfcrt-F145I (3 isolates in 2014 and 10 isolates in 2015). All of the 16 parasites with pfk13-R539T carried a single pfpm2 copy number with no pfcrt-F145I mutation but four of them had multiple pfmdr1 copies. Figure 1 shows the prevalence of pfpm2 copy number variation, pfmdr1 copy number variation, pfcrt-F145I, pfk13-C580Y and pfk13-R539T mutations over time. The prevalence of parasites with multiple pfpm2 copies associated with PIP resistance increased from 2013 to 2015, similar to the prevalence of parasites harboring pfcrt-F145I and pfk13-C580Y mutations. In contrast, the prevalence of parasites with multiple pfmdr1 copies associated with MQ resistance and pfk13-R539T mutation decreased after 2013 through 2015.   In vitro ring stage survival assay (RSA) and pfk13 mutations. To better understand ART-resistance, 40 isolates from 2013 to 2015 were tested in in vitro RSA 0-3 h to measure % survival rate against DHA, and the association with pfk13 mutations assessed ( Fig. 2A). One isolate was excluded due to the growth rate between 0 and 72 h being less than 1.  www.nature.com/scientificreports/ In the light of an association between PIP resistance and haplotypes, Table S4 shows PIP-IC 90 s of the parasite isolates in different haplotype subgroups. Parasite isolates containing multiple pfpm2 copies alone (group III), the pfcrt-F145I alone (group V), or the combination of both markers exhibited noticeably high IC 90 values, indicative of PIP resistance phenotype.
Preexisting antimalarial treatment. During the enrollment process, participants were queried about recent malaria infection or antimalarial drug use and also past medical history of malaria. Of 112 participants, 27 (24.1%) had reported prior malaria infection: 15 (55.6%), 7 (25.9%), 3 (11.1%), and 2 (7.4) had 1, 2, 3, and 4 episodes of malaria in the last 12 months. Data on antimalarial treatment prior to enrollment in this study revealed that 78 participants (69.6%) had not taken any antimalarial drugs. The most commonly used antimalarials were primaquine (PQ) (26.8%), followed by AS and MQ (16.1%) and CQ (10.7%).  Table 4). Sixteen of the 21 had MQ, the first-line therapy for P. falciparum at the time of the study and with a half-life of 8 to 20 days. Six had PQ or cPQ; in 2013-2015 a single dose of PQ had not yet been widely used for P. falciparum anti-gametocyte activity but it was first-line for P. vivax infections. Interestingly four people had evidence of DHA or AS which has a very short half-life (2-3 h). When comparing both methods, 12 patient plasma samples were positive for antimalarial drugs in both the ex vivo bioassay and LC-MS/ MS methods while another 9 samples were negative in the ex vivo bioassay but were positive by LC-MS/MS. It is important to note that only 8 samples (38.1%) from the positive samples could be in vitro cultured and evaluated for their drug susceptibility, which showed no significant difference to those without the baseline medicines.

Discussion
Even though the Greater Mekong Subregion (GMS) has long been the epicenter of antimalarial drug resistance, these countries are aiming to achieve malaria elimination by 2030 49 . Military populations along the Thailand-Cambodia border often patrol in forested areas where the majority of multidrug resistant malaria in this region has been identified. We strategically located our lab at the nearest clinic. Similar approaches of forward deployed laboratories can be applied to migrants and occupations where follow-up post treatment is difficult, including workers in forestry, agriculture or animal husbandry and refugee populations. Despite the limitations of the study, where clinical treatment outcomes were not readily available, we were able to characterize the changes in drug resistance markers overtime in this difficult to reach population.
Chloroquine resistance was first reported in the GMS the 1950s 50 , and several groups 51-55 have characterized CQ resistant parasite lineages as one of four mutant pfcrt genotypes at positions 72-76 (CVIET, SVMNT, CVMNT, and CVMET; mutation underlined). The present results show that all collected parasite isolates harbored the pfcrt 72-76 CVIET haplotype similar to previous reports [56][57][58] . Even though CQ sensitivity was still slightly higher than that of the CQ-resistance W2 reference clone, when the geometric mean IC 50 s of the CQsensitive isolates of 30.1 nM 59 was applied, all the parasites collected under this study are deemed CQ-resistant. Decreased CQ sensitivity was observed from 2013 to 2015 that can likely be attributed to continued CQ use in Thailand for treatment of P. vivax.
Parasite isolates carrying the artemisinin resistant gene have been reported in Thailand 58,60 . The current study showed that the pfk13-C580Y mutation was predominant, rising from 63% in 2013 to 100% in 2015. These findings are in good agreement with those reported by others 41,58 in that several pfk13 mutations developed in Table 4. Comparison of preexisting antimalarial activity in the samples obtained by ex vivo bioassay and liquid chromatography-mass spectrometry (LC-MS/MS). www.nature.com/scientificreports/ Sisaket (62% 580Y) 41 but by 2015 the C580Y mutation became the predominant mutation found 58 . In contrast, in southern Thailand (Yala Province) very few isolates with pfk13 mutations have been detected 58,61 . The pfk13-R539I found in this study has not yet been reported in Thailand. The pfk13-R539I mutation was previously observed at 0.9% in Ghanaian isolates 62 . As assessed by the RSA 0-3 h , the R539I mutation does not appear to be associated with ART resistance. All the collected isolates had a wild-type cytochrome b gene (pfcytb) at both codons 258 and 268 (known to be associated with atovaquone/proguanil treatment failure 37,39,63 ) which is consistent with other isolates collected along Thai-Myanmar and Thai-Cambodia borders during 1998-2005 64 . Of all the pfmdr1 SNPs analyzed, allelic variation was only observed in pfmdr1 position 184. This is in good agreement with the previous report by Thita et al. 65 , in which 89% of P. falciparum isolates from the Thai-Cambodia border in Chanthaburi and Trat province from 1988 to 2016 had the pfmdr1 184F allele. Previous studies in Uganda and Bioko Island suggest that this allele may play a certain role in mediating resistance to some antimalarials 26,66 , and while in the current study, decreased susceptibility to QN, CQ, and MQ was observed, causality cannot be established. Different results were observed from the P. falciparum isolates from the southern part of Thailand, with the pfmdr1 86Y allele significantly more common 67 .
It was previously reported that parasites containing a pair of point mutations in P. falciparum dihydrofolate reductase (pfdhfr) (A16V and S108T) are resistant to cycloguanil but not to pyrimethamine, while those with I164L in conjunction with S108N show high levels of resistance to both cycloguanil and pyrimethamine 68 . We did not sequence pfdhfr and P. falciparum dihydropteroate synthase (pfdhps) genes in our study. However, based on previously published data demonstrating fixed proguanil resistance in this region 69,70 , it would be no surprise that parasites in the present study would carry significant antifolate resistance on a background of pfdhfr-pfdhps mutations. This is consistent with the study results demonstrating high IC 50 values for CYC.
In addition to SNPs in established malaria resistance markers, pfpm2 and pfmdr1 increased copy numbers are associated with PIP and MQ resistance, respectively. The increased trend in multiple pfpm2 copy numbers were clearly seen with a decline in multiple pfmdr1 copy numbers. Novel mutations of the pfcrt gene downstream of the 4-aminoquinoline resistance locus (positions 72 to 76), including T93S, H97Y, F145I, I218F, M343L, or G353V, were recently shown to be associated with PIP resistance 31,33,34 . From this study, there was a relatively rapid rise in the F145I mutation, from 0% in 2013 to 30% in 2015. In this study, IC 90 values of PIP were employed to investigate the PIP resistant phenotypes associated with the tested genetic markers. Parasite isolates holding either multiple pfpm2 copies with pfcrt-F145 or the combination of single or multiple pfpm2 copy numbers with pfcrt-F145I mutation exhibited elevated IC 90 values of PIP. The high IC 90 values of parasite isolates with only multiple pfpm2 copy numbers could be stemmed from pfpm2 amplification or other novel pfcrt mutations that were not detected in this study. Previous studies have shown that the overexpression of pfpm2 and pfpm3 in the 3D7 genetic background did not alter PIP sensitivity, suggesting that the increase in pfpm2 copy number alone is not the sole modulator of PIP resistance 46,71 .
In 2013, an AFRIMS study in Anlong Veng, Cambodia, 12 km from the Cambodia-Thai border at Sisaket province, showed a similar rate of 65% C580Y pfk13 mutant 40 . When the 54% failure rate of DHA-PIP was observed in Anlong Veng in 2013 40 , the molecular markers for neither ART nor PIP were available so it is only retrospectively that the association of molecular markers and in vitro data can be seen. Similarly in Sisaket, the 87% failure rate of DHA-PIP from the TRACII study (2015-2018) 60 may have been predicted from this surveillance data.
The use of previous antimalarial drugs or preexisting antimalarial activity in patients can have an effect on malaria treatment 45 . A bioassay method was developed for the measurement of antimalarial activities of ART and its derivatives in either the plasma or sera of patients 72,73 , while bioanalytical LC-MS/MS has widely been used to detect several drugs 74,75 . For example, 3 samples found positively in the LC-MS/MS method contained PQ, the erythrocytic activity of which is little, hence not detected in ex vivo bioassay. Another 6 samples, on the other hand, were found negative in ex vivo bioassay but positive in LC-MS/MS method with the presence of MQ and cMQ at the concentration lower than 150 ng/ml. These results suggest that the bioassay was able to screen blood stage antimalarial activity in certain levels but could not specify which drugs or metabolites are present in the patient plasma samples whereas LC-MS/MS method could. In this study, medical records of the P. vivax infected patient were not included in our data set; however, P. vivax coinfection is common in this region 70 . It should be noted that PQ would not be expected to persist long in the serum, and the far more likely explanation is that PQ was given as a treatment to prevent transmission of P. falciparum. This is now common practice in Thailand.
The difficulty of choosing a national first-line therapy in the GMS is due to differing patterns of antimalarial resistance as we illustrated. High rates of AS-MQ failures along the Thai-Myanmar border 16,76 resulted in the first-line treatment being changed to DHA-PIP in 2015. However, the parasites from this high risk mobile military population had patterns developed over three years which were more similar to isolates in Cambodia, and then manifested a similar drop in DHA-PIP efficacy. The Thai government responded to that regional challenge in 2019, by changing the first-line treatment for P. falciparum to artesunate-pyronaridine (AS-PND) in Sisaket and Ubon Ratchathani Provinces. Recent reports suggest there are two lineages within the GMS. The C580Y mutation predominate in the eastern GMS and the F466I mutation has spread along the Myanmar border 58 . While there may be predominant lineages which have emerged, continued local surveillance will be important to determine if and how these patterns shift. The emergence of the novel pfcrt-F145I mutation in 2014 from this study followed the emergence and DHA-PIP treatment failures in Cambodia 40,45,77 . Shrestha, et al. 78 showed that more than 98% of the parasites collected from northern Cambodia carried newly emerged pfcrt mutations and after 2014, the prevalence of parasites with pfcrt-F145I mutation started to decrease, being out-competed by other less resistant, but more fit pfcrt variants such as pfcrt-T93S and I218F 34 . Targeting surveillance studies on more mobile human populations along border areas such migrants and militaries, can help predict drug efficacy or first detect any introduction of novel mutations which may interrupt the path to malaria elimination in the GMS. www.nature.com/scientificreports/ A few key limitations of this study should be considered. This includes the fact that the study samples were collected more than 5 years ago from one geographical area. Therefore, drug resistance profiles reported here should not be construed for all of Thailand. Due to inability to culture adapt all the samples, only 33% of the collected samples could be analyzed for in vitro drug susceptibility. The majority of the samples were genotyped for the presence of drug resistance SNPs; however, it is realized that other SNPs not included in this study may contribute to drug resistance. Finally, multiplicity of infections was not examined. However, future studies will employ whole genome sequencing to address the complexity of multiple infections.

Data availability
All data generated or analyzed during this study are included in this published article and its supplementary information files. www.nature.com/scientificreports/