Introduction

An estimate of 229 million malaria cases occurred worldwide in 20191 and artemisinin-based combination therapies (ACT) have been the first-line drug treatment for uncomplicated Plasmodium falciparum infection for most malaria endemic areas; however, resistance to both artemisinin (ART) and its partner drugs has increased over the years at a pace requiring intensive surveillance and monitoring. The first case reports of ART resistance emerged in Thailand in the late 1990s2 and in 2003–2005. In 2008, Thai Ministry of Public Health (MoPH) and Armed Forces Research Institute of Medical Sciences (AFRIMS) studies reported that first-line ACT regimens (2-day artesunate (AS) plus mefloquine (MQ)) were failing along both sides of the Thai-Cambodia border near Trat, Thailand and Pailin, Cambodia3. Militaries and other mobile forest goers in these border areas constitute 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.

Following the case reports of ACT failure, the Artemisinin Resistance: Confirmation, Characterization and Planning for Containment Project (ARC) in 2006 showed a delayed clearance resistance phenotype in Cambodia along the Thai border which continued to spread across mainland Southeast Asia4,5,6,7,8,9. It was not until 2014 that an ART resistance molecular marker was identified with mutations in the propeller domain of P. falciparum Kelch-13 gene (pfk13) (PF3D7_1343700)10,11,12,13, which can be confirmed by the ring-stage survival assay (RSA)10,14 and/or delayed parasite clearance times in clinical trials4,5. Multiple copies of P. falciparum multi drug resistance 1 (pfmdr1) (PF3DF_0523000) is a well-established marker for MQ resistance15,16, while the pfmdr1 single-nucleotide polymorphisms (SNPs) have been associated with modulation of parasite tolerance or susceptibility to several antimalarial drugs including quinine (QN), amodiaquine (AQ), chloroquine (CQ), MQ and lumefantrine (LUM)17. The main pfmdr1 SNPs associated with drug resistance include N86Y, Y184F, S1034C, and N1024D18,19,20,21,22,23. Mutation N86Y in pfmdr1 is associated with increased CQ resistance but increased sensitivity to MQ24. The pfmdr1 184F mutation alone was not associated with susceptibility response of CQ, MQ, LUM, QN, monodesethylamodiaquine (MDAQ), and dihydroartemisinin (DHA)25. However, the parasites carrying pfmdr1 86Y and Y184 haplotype was associated with increased susceptibility to MDAQ, LUM, and MQ. The possible role of pfmdr1 N86, 184F and 1246D alleles and pfmdr1 copy number associated with artemether-lumefantrine resistance is not confrmed26. An ACT alternative to AS-MQ is dihydroartemisinin-piperaquine (DHA-PIP) and the first PIP resistance markers identified were multiple copies of P. falciparum plasmepsin 2 (pfpm2) (PF3D7_1408000)27,28,29,30, followed by mutations in the P. falciparum chloroquine resistance transporter (pfcrt) (PF3D7_0709000) downstream of the 4-aminoquinoline resistance locus (positions 72–76 with K76T)31,32,33,34, and the E415G mutation in P. falciparum exonuclease (pfexo) (PF3D7_1362500)27,28,35. PIP resistance can be characterized by the elevated IC90 values36 and PIP survival assay (PSA) with a survival rate of more than 10%, defining a PIP resistant phenotype31. Atovaquone-proguanil (ATQ-PG) is not an ACT but ATQ resistance subsequently emerged, which was linked to specific mutations in P. falciparum mitochondrial cytochrome B (pfcytb) (AY282930.1) in particular the mutation at positions 258 and 26837,38,39.

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 Province40. 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–76) and pfmdr1 SNPs (Positions 86, 184, 1034, 1042). Primers used to identify pfkelch13, pfcytb, pfcrt and pfmdr1 SNPs are shown in Tables S1 and S2.

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 described15,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 pfpm229. The primers and probes (Table S2) used were as previously described to amplify the following loci: pfmdr1 (PF3D7_0523000) and pfpm2 (PF3D7_1408000), respectively29. 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 pfmdr115 and 1.6 copies for pfpm229 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 Jensen43. 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% CO2, 5% O2, 90% N2 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 (IC50 and IC90) was performed following previously published methods44,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 previously36. IC50s and IC90s were estimated by nonlinear regression analysis using GraphPad Prism version 6.0.

Ring-stage survival assay (RSA)

In vitro RSA0-3 h was performed on 0–3 h post-invasion rings obtained from selected culture-adapted clinical isolates following published methods with slight modifications14,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% CO2, 5% O2 and 90% N2. 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 RSA0–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 methods47,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% CO2, 5% O2 and 90% N2, 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 µm PTFE filter and then transferred to an HPLC vial. The LC separation was performed on ACQUITY UPLC (Waters) coupled with tandem mass spectrometer (Xevo TQ-S, Waters) and eluted on Waters Acquity UPLC® BEH C18, 2.1 × 50 mm, 1.7 µm column at a flow rate of 0.5 mL/minute in 8 min run time. Mobile phase consisted of (A) 5 mM ammonium acetate pH 4.5 in water and (B) 5 mM ammonium acetate pH 4.5 in acetonitrile:methanol (50:50 v/v). The gradient started with 10% B, raised to 98% B in 6 min and held at this composition for 30 s, decreased to 10% B and re-equilibrated for 1 min. Column temperature was set at 40 °C. Selective mass to charge (m/z) transition for each compound was monitored as follow: AS (407.4 > 261.23), DHA (307.13 > 261.26), MQ (379.1 > 321.0), carboxy primaquine (cMQ) (310.22 > 290.05), PG (256.1 > 171.9), cycloguanil (CYC) (251.95 > 194.9), primaquine (PQ) (260.07 > 85.95), carboxy primaquine (cPQ) (275.34 > 175.11), orthoquinone (OQ) (260.2 > 147.08), stable isotopic labeled artesunate (SIL-AS) (411.3 > 261.18), and stable isotopic labeled primaquine (SIL-PQ) (264.35 > 86.15).

Statistical analysis

Statistical analysis was performed using GraphPad Prism version 6.0 (GraphPad Software, Inc., San Diego, CA, USA). The difference of data between groups was assessed by nonparametric Kruskal–Wallis, Mann–Whitney or Dunn’s multiple comparison tests, as appropriate. Statistical significance was defined as a P value < 0.05.

Ethics approval and consent to participate

All participants or guardians provided written consent and samples were collected under approval Royal Thai Army Institutional Review Board (RTA IRB) and Walter Reed Army Institute of Research Institutional Review board.

Disclaimer

Material has been reviewed by the Walter Reed Army Institute of Research. There is no objection to its presentation/publication. The opinions or assertions contained herein are the private views of the author, and are not to be construed as official, or as reflecting true views of the Department of the Army or the Department of Defense. The investigators have adhered to the policies for protection of human subjects as prescribed in AR 70–25.

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 (Table 1). The median age of the participants was 23 years (IQR: 22–30.75) and most were male (98.2%), in the military (83%), and from Sisaket Province (96.4%). Among rapid diagnostic test and/or microscopy positive malarial isolates, 108 samples (96.4%) were confirmed by PCR as P. falciparum infections and 4 samples (3.6%) as mixed P. falciparum and P. vivax infections. The geometric mean of parasitemia of the participants was 12,097 parasites/µL.

Table 1 Patient and parasitological characteristics of 112 participants.

Mutations in pfmdr1, pfcrt, pfk13, and pfcytb gene

112 samples were evaluated for SNPs of pfmdr1, pfcrt (positions 72–76), pfcrt (position 145), pfk13, and pfcytb as molecular markers of drug resistance associated with MQ, CQ, PIP, ART, and ATQ, respectively. Of four point mutations on pfmdr1, only the Y184F mutation (97.3%) was detected (Table 2), suggesting the presence of two pfmdr1 haplotypes in the samples: NYSN and NFSN (representing amino acids positions 86, 184, 1034, and 1042, respectively) with NFSN as the most prevalent genotype. Genotyping of pfcrt showed all mutant haplotypes at positions 72–76 (CVIET), which is associated with CQ resistance. Regarding the novel pfcrt mutation (F145I) associated with PIP resistance, no F145I mutation was identified in 2013, but the F145I mutation of pfcrt was detected in 5 (12.8%) samples and 14 (30.4%) in 2014 and 2015, respectively. Concerning pfk13 as associated with ART resistance, one of the three nonsynonymous mutations (R539T, R539I, or C580Y) in the k13 propeller domain was observed in 112 samples and C580Y was the predominant pfk13 mutation, accounted for 84.8% of all pfk13 mutations. The number of C580Y mutants increased from 2013 (63%) to 2015 (100%). There were no mutations detected in pfcytb for amino acid positions 258 and 268.

Table 2 Mutations of malaria resistance molecular markers of P. falciparum in Thailand between 2013 and 2015.

Copy number variation of pfmdr1 and pfpm2

Single pfmdr1 copies were present in nearly all of evaluable cases (107/112), only 5 samples had multiple copies of pfmdr1, found in samples collected in 2013 and 2014. Multiple pfpm2 copies were found in about half of evaluable cases (58/112), and the number of parasites harboring multiple pfpm2 copies steadily rose from 7.4 to 78.3% from 2013 to 2015. Of 5 samples containing multiple pfmdr1 copies, all were single copy pfpm2 with pfcrt F145 wild-type. All parasites harboring multiple pfpm2 copies contained the pfpk13-C580Y mutation and of those samples, 13 (22.4%) contained the novel pfcrt F145I mutation.

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.

Figure 1
figure 1

Prevalence of antimalarial drug resistance mutations in Thailand from 2013–2015. Bars indicate the prevalence of parasites with single and amplified pfpm2 copy number. The line graphs indicate the prevalence of parasites harboring the pfcrt-F145I mutation, amplified pfmdr1 copy number and pfk13-C580Y and pfk13-R539T mutations.

In vitro drug susceptibility of P. falciparum isolates

Of 86 P. falciparum monoinfection samples that met the parasitemia level criteria for in vitro culture, in vitro drug susceptibility of only 40 samples (11, 11 and 18 samples from 2013, 2014, and 2015, respectively) against an antimalarial drug panel could be obtained by HRP2-ELISA (Table 3). In viewing the geometric mean (GM) IC50 data for all drugs by year, CQ, PIP, CYC, MQ, and QN all exhibited statistically significant decreases in drug susceptibility. Compared to the ART- and MQ-sensitive W2 reference clone, isolates had greatly reduced sensitivity to PIP, CYC, and DOX, with PIP and CYC having a noticeable drop in effectiveness between 2014 and 2015. In contrast, while there was a steady decline in CQ and MQ sensitivity from 2013 to 2015 the GM-IC50 values remained lower than that of the CQ-resistance W2 reference (IC50-CQ = 234 nM) and MQ-resistance D6 reference (IC50-MQ = 118 nM). There was a modest increase in DHA and AS susceptibility, which did not change dramatically over the 3-year period. All parasite isolates had an IC50-ATQ lower than that of the ATQ-resistance C2B reference (IC50-ATQ = 13,240 nM), suggesting no ATQ-resistant phenotype. No significant changes in drug susceptibility of DHA, LUM and ATQ were observed during this study period.

Table 3 In vitro susceptibility of P. falciparum isolates to dihydroartemisinin (DHA), artesunate (AS), artemisone (ATM), chloroquine (CQ), piperaquine (PIP), cycloguanil (CYC), mefloquine (MQ), quinine (QN), lumefantrine (LUM), doxycycline (DOX), atovaquone (ATQ).

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 RSA0-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.5, leaving 39 isolates evaluable by RSA0-3 h. Survival rates for the 11 evaluable isolates from 2013 ranged from 0.84 to 44.5% (median 19.9, IQR of 3.4–30.3), for the 11 evaluable 2014 isolates the range was 3.7 to 82.8% (median 16.1, IQR of 9.9 -60.9) and for the 17 evaluable isolates in 2015 ranged from 1.4 to 22.8% w (median 8.3, IQR of 2.8 -22.8), whereas reference clone survival rate was as expected with medians of 0.4, 2.5, and 23.2% for ART-sensitive W2, ART-resistant IPC-4884 and IPC-5202, respectively. In 2013, 90.9% (10/11) were deemed ART-resistant, all of which contained either pfk13-R539T or -C580Y mutations, while the parasite harboring pfk13-R539I was sensitive to ART. In 2014 and 2015, all tested parasites containing either pfk13-R539T or –C580Y mutation were deemed ART resistant. No significant difference in the proportion of resistant isolates between 2013 and 2014 were found (P > 0.05) but there was a significant decrease in median survival rate for 2015 (8.3), compared to those from 2014 (16.1) and 2013 (20.9), Mann–Whitney test, P = 0.004 and P = 0.02, respectively. It was noted that % survive rate of the parasite harboring pfk13-R539I was 0.84%, closer to the cut-off value.

Figure 2
figure 2

In vitro characterization. (A) RSA and pfk13 mutations. RSA0-3 h survival rate for standard laboratory-adapted clones (W2 for ART-sensitive control, IPC-4884 and IPC-5202 for ART-resistance control) and culture-adapted clinical isolates. The dashed line represents the 1% survival rate cut-off that differentiates ART-resistance (≥ 1% red-dashed line) from ART-sensitive (< 1%) parasites in RSAs. Median and interquartile ranges are shown. (B) IC90 to PIP in parasites with pfcrt-F145I mutation, with and without pfpm2 amplification. Copy number variations of pfpm2 are shown on the x-axis, and PIP-IC90 is shown on the y-axis (log10 scale).

PIP-susceptibility of parasites with PIP molecular markers

With the emergence of PIP resistance associated with DHA-PIP treatment failures identified in Cambodia from 2010–201341, the distribution of PIP-IC90 values for parasites from Thailand (with newly emerged pfcrt-F145I mutation with or without amplified pfpm2) is displayed in Fig. 2B. In 2013, there was only one isolate with multiple pfpm2 copies and it had an elevated IC90 value. In 2014, no significant differences in IC90 values in the isolates with single and multiple pfpm2 copies were observed (P > 0.05). In 2015 there was a significant difference in parasites with single (226 nM) and multiple pfpm2 copies (62,347 nM) without pfcrt-F145I (P = 0.036) whereas no significant differences in IC90 values in parasites with single (53,930 nM) and multiple pfpm2 copies (77,882 nM) with pfcrt-F145I (P > 0.05). Isolates with single pfpm2 and wild type pfcrt had more than tenfold lower IC90s than parasites with multiple pfpm2 and/or pfcrt-F145I. With the exception of 1 isolate, parasites harboring pfcrt-F145I mutation exhibited higher IC90 values than those without the mutation (PIP-IC90s for standard clone W2, D6 and C2B: 124, 157, and 232 nM, respectively).

In the light of an association between PIP resistance and haplotypes, Table S4 shows PIP-IC90s 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 IC90 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%).

Many antimalarials are long-acting and participants may not disclose use of antimalarial drugs; therefore, antimalarial drug levels in plasma samples were determined by two independent techniques, ex vivo bioassay and LC–MS/MS. In the ex vivo bioassay, 12 of 112 patient plasma samples (10.7%) evaluated against P. falciparum W2 clone had significant antimalarial activity (> 17.6 nM DHA equivalents) (Table 4), while 89.3% of evaluable samples were considered negative (≤ 17.6 nM). The positive results were measured in DHA activity equivalents ranging from 49.8 to 3,855 nM with a median activity of 234.2 nM (IQR of 94.3 to 751.4 nM). When LC–MS/MS was employed to detect antimalarial drugs in patient plasma samples, 21 samples (18.8%) were found to contain some antimalarial drugs (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.

Table 4 Comparison of preexisting antimalarial activity in the samples obtained by ex vivo bioassay and liquid chromatography-mass spectrometry (LC–MS/MS).

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 203049. 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 1950s50, and several groups51,52,53,54,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 reports56,57,58. Even though CQ sensitivity was still slightly higher than that of the CQ-resistance W2 reference clone, when the geometric mean IC50s of the CQ-sensitive isolates of 30.1 nM59 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 Thailand58,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 others41,58 in that several pfk13 mutations developed in Sisaket (62% 580Y)41 but by 2015 the C580Y mutation became the predominant mutation found58. In contrast, in southern Thailand (Yala Province) very few isolates with pfk13 mutations have been detected58,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 isolates62. As assessed by the RSA0-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–200564. 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 antimalarials26,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 common67.

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 pyrimethamine68. 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 region69,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 IC50 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 resistance31,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, IC90 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 IC90 values of PIP. The high IC90 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 resistance46,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 mutant40. When the 54% failure rate of DHA-PIP was observed in Anlong Veng in 201340, 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 treatment45. A bioassay method was developed for the measurement of antimalarial activities of ART and its derivatives in either the plasma or sera of patients72,73, while bioanalytical LC–MS/MS has widely been used to detect several drugs74,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 region70. 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 border16,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 border58. 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 Cambodia40,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 I218F34. 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.

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.