No Evidence of Plasmodium falciparum k13 Artemisinin Resistance-Conferring Mutations over a 24-Year Analysis in Coastal Kenya but a Near Complete Reversion to Chloroquine-Sensitive Parasites

Antimalarial drug resistance is a substantial impediment to malaria control. The spread of resistance has been described using genetic markers, which are important epidemiological tools. We carried out a temporal analysis of changes in allele frequencies of 12 drug resistance markers over 2 decades of changing antimalarial drug policy in Kenya.

well as of multidrug-resistant (resistance to both artemisinin and the partner drug, piperaquine) P. falciparum parasites in Cambodia and neighboring Vietnam (10,11), are now a major concern. Though isolated, the recent possible identification of parasites with high artemisinin survival rates in Equatorial Guinea and Uganda (12,13) calls for continued surveillance of artemisinin resistance markers in Africa.
The identification of P. falciparum genetic mutations associated with antimalarial drug resistance has provided molecular markers for surveillance of resistance, both in real time and retrospectively, to assess geographic origins and migration patterns of drug-resistant parasites (14). The use of such markers provided information on the origin and spread of CQ (1) and sulfadoxine-pyrimethamine (SP) (2,3). Likewise, mutations in the kelch13 gene (k13), including F446I, N458Y, M476I, Y493H, R539T, I543T, P553L, R561H, and C580Y, are enabling the detection and tracking of artemisinin resistance in SE Asia (9,15). Some of these mutations (Y493H, P553L, R561H, and C580Y) have been observed at low frequencies (Ͻ0.1%) in Africa (16). Whole-genome studies of artemisinin-resistant parasites in SE Asia were shown to accumulate artemisinin resistance predisposing mutations at other genetic loci (17). These mutations were found in apicoplast ribosomal protein S10 precursor gene (arps10) codon V127M, chloroquine resistance transporter (crt) codon I356T, ferredoxin (fd) codon D193Y, and multidrug resistance (mdr) protein 2 codon T484I. While these mutations have been identified in parasites from Africa, they occur at very low frequencies (arps-10-V127M, 0%; fd-D193Y, 0.1%; mdr2-T484I, 0.1%; crt-N326S, 0.8%) (16). Nonetheless, it is crucial for regions outside SE Asia to monitor the emergence of artemisinin resistance signatures, including the selection of markers associated with changes in antimalarial drug policy, drug trials, and experimental analyses. Already, recent studies have shown that artemether-lumefantrine (AL) selects for polymorphisms in the chloroquine resistance transporter (crt) (K76 allele) and multidrug resistance 1 (mdr1) (N86/184F/D1246-NFD haplotype) genes (18)(19)(20)(21). Furthermore, the mdr1-NFD haplotype and an increase in mdr1 copy number have been linked with reduced susceptibility to lumefantrine (LM) and mefloquine (MQ), respectively (22)(23)(24). Moreover, several recent studies have shown the selection of additional polymorphisms associated with ACT pressure and warrant further investigation. A mutation (S69Stop) in the cysteine proteinase falcipain-2a gene has previously been selected for after artemisinin in vitro selection pressure, using a parasite isolate from Africa (15). A clinical trial conducted in western Kenya revealed that 2 to 3 days of ACT treatment selected for parasites with either the 160N or 160T allele in the AP-2 complex subunit mu gene (ap2-mu, S160N/T) and the allele 1528D in the ubiquitin carboxyl-terminal hydrolase 1 gene (ubp-1, E1528D) (21). In The Gambia, a temporal increase in the frequency of the K65 allele (K65Q) in the cysteine desulfurase (nfs) gene was observed 6 years after the introduction of ACTs. Moreover, the 50% inhibitory concentration (IC 50 ) values for LM were significantly higher in nfs-K65 wildtype field isolates than in the mutant 65Q isolates (25).
Mutations in the dihydrofolate reductase gene (dhfr 51I, 59R, and 108N triple mutant) combined with mutations in the dihydropteroate synthetase gene (dhps 437G and 540E double mutant) are associated with clinical and parasitological SP treatment failure in East Africa (26)(27)(28)(29). Although SP is the drug of choice for intermittent preventive treatment in pregnant women (IPTp) (30), regional differences in drug resistance critically impact the success of this intervention (31). An additional dhps mutation (A581G) has been associated with reduced IPTp efficacy when the prevalence of sextuple-mutant (dhfr 51I/59R/108N and dhps 437G/540E/581G) parasites exceeds 37% (32). Also, although its impact has yet to be elucidated, a novel dhps I431V allele has been detected in West Africa, Cameroon, and Nigeria (33,34).
These studies highlight the need for continued surveillance of known drug resistance markers and novel markers to maintain the gains made in using ACTs and IPTp in controlling malaria. crt, mdr1, and dhfr have well-described selection patterns in response to the withdrawal of CQ and SP. Based on previous findings, there was evidence of selection toward a predominance of wild-type alleles in crt and mdr1 and toward the fixation of the mutant dhfr alleles, over a 20-year period of changing 2-fold increases and decreases in frequency: the KNE repeat at codon 1514 (5.5% to 11.6% and to 5.41%, respectively) and N1518N (5.6% to 11.7% and to 5.5%, respectively). Of the 23 haplotypes, the 3D7 haplotype dominated throughout the sampling period, and no significant temporal trends were observed (Table S11).
In nfs, only codon K65Q was found to have a significant, albeit marginal, trend preand post-ACT introduction ( 2 ϭ 4.4, P ϭ 0.04). It was also in high linkage disequilibrium with codons S62N and E67G. The 2005/2006 period appears to be the point at which the mutant allele 65Q and wild-type 65K diverge in frequency in opposite directions, with the wild-type 65K decreasing and the mutant allele 65Q increasing in frequency ( Fig. 1D and Table S7). The 3D7 haplotype was the seventh most dominant of all the 24 haplotypes, and there were no significant temporal trends in haplotype frequencies (Table S13). Additionally, we examined for spatial variation for crt, mdr1, and nfs and did not see any variation in alleles according to geographical area (Table  S15).
serine-tRNA ligase, putative gene. Serine-tRNA ligase, a marker not associated with resistance or drug selection, had only one polymorphic codon observed across all time points (L84V), with frequencies ranging between 2 and 5%, while the rest were rare (Ͻ5%) ( Table S8). A total of 15 haplotypes were observed, and only the 3D7 haplotype occurred across all time points and with frequencies of Ͼ80%, with the rest of the haplotypes being rare (Ͻ5%) (Table S14).

DISCUSSION
The stable frequencies of the common k13 A578S African mutation and K189T (16,40,41) suggest that these mutations are not under selective pressure from artemisinin.  Notably, a recent study confirmed that the A578S mutation does not confer resistance to artemisinin (41). The codon 137 Asn repeats, found at lower frequencies than in SE Asia, have been associated with day 3 positivity (36) and cooccur with k13 propeller mutants (37). However, it is not known how these polymorphisms modulate artemisinin resistance, and further studies are needed to investigate this. Similar to parasites from Africa, the N-terminal region of k13 had higher K a /K s ratios than did the C-terminal region, contrary to k13 data from SE Asia (16), implying that parasites from Africa acquire more changes in the N-terminal region than the C-terminal region. Moreover, as has been observed in parasites from Africa (16), there was no artemisinin resistancepredisposing mutations, and the mdr2 I492V mutation showed no evidence of selection over time. This mutation (I492V) has also been identified at 100% frequency (n ϭ 38) in Suriname (42). The withdrawal of CQ (43) has resulted in the rapid decline in CQR alleles crt 76T and mdr1 86Y and mdr1 1246Y in Kilifi. Nationally, a decline in CQR resistant alleles has previously been observed on the South Coast of Kenya, with the crt 76T and mdr1 86Y alleles decreasing from 88% to 63% and 75% to 54% from 1998 to 2008, respectively (44). In western Kenya, where malaria transmission is holoendemic, the crt 76T, mdr1 86Y, and mdr1 1246Y alleles dropped from 86% to 2%, 92% to 1%, and 67% to 6% between the years 2003 and 2014, respectively (45). In the same region, CQ median IC 50 values decreased significantly, from a median of 92 to 22 from 2008 to 2011 (P Ͻ 0.001) (46). Therefore, a national rerollout of CQ in the next few years is a possibility. However, it is also likely that resistant parasites may remain in the population below detectable levels, and reemergence from these parasites or the reintroduction from surrounding areas could be rapid. Regarding mdr1, the disappearance of the triple-mutant mdr1-YYY is also likely to be the result of CQ withdrawal. Likewise, the oscillation of the mdr1-NFD with the NYD (wild type) at a high frequency post-ACT introduction suggests that the NYD haplotype is likely to have risen in frequency due to the absence of CQ, restoring the wild-type parasite population, while the increase in the NFD haplotype frequency may be attributable to artemether-lumefantrine (AL) pressure (19).
The high frequency of SP resistance markers in Kenya (35) may be attributable to the continued distribution of SP for malaria case management in the private sector (47). Though SP maintains its utility in IPTp, the loss of IPTp efficacy has been noted when the prevalence of the sextuple-mutant (dhfr 51I/59R/108N and dhps 437G/540E/581G) parasites exceeds 37% (32), as seen in some sub-Saharan African settings. Consequently, the dhps 581G mutation needs to be monitored, given that it first occurs in our population in 2015/2016 at 3%. The dhps I431V allele was not detected, and its distribution could be restricted to West Africa (33,34).
We noted a decline in the nfs K65 wild-type allele in the Kilifi population since the introduction of AL in 2004, contrary to recent findings from The Gambia in West Africa (25), and it appears that there is an opposite trend of the K65 allele in this East African population. The Gambian parasites have shown increasing tolerance to lumefantrine (LM) in a study conducted between 2013 and 2015 (48); however, drug trials with AL in western Africa, including The Gambia, show that AL is still highly efficacious (49). Perhaps these discordant observations are due to differences in allele frequencies of other loci, such as pfmdr1, between East and West Africa (50), or differences in drug policy, such as the extensive use of amodiaquine (AQ) in West Africa compared to East Africa (51). AQ shows activity against CQR parasites, like in West Africa, where the prevalence of CQR parasites is still high (52). However, in East Africa, there are reports of AQ resistance (53,54) and AQ-resistant parasites are sensitive to LM (19). Thus, this geographical difference in drug use, the inverse resistance profile between AQ and LM, and the predominance of CQS parasites in East Africa suggest that AQ may drive LM resistance in West Africa.
The distinct change in frequency of resistance-associated alleles observed with crt, mdr1, dhps, and nfs were not seen in the other genes we evaluated, ap2-mu, falcipain-2a, and ubp-1. Notably, ap2-mu 160N showed no evidence of selection and has not been observed in parasites in SE Asia, where artemisinin resistance is common (16). Therefore, it appears to be restricted to Africa, as it has also been observed in Ghana (55). Recent evidence points to a different mutation, I592T (not detected in this population), that showed increased ring-stage survival following a dihydroartemisinin pulse (56). Further work is required to evaluate the potential role of ap2-mu mutations on the ACT response.
falcipain-2a was the most polymorphic gene in this study, and this can be attributed to drug pressure, as it is associated with the parasite's digestive vacuole, the target of many antimalarials (57). Nonetheless, The S69Stop mutation that results is artemisinin resistance in vitro (15) was not identified in our population. However, we found a significant temporal increase in the S59 allele post-ACT introduction, potentially due to AL pressure. It remains to be seen what the role of falcipain-2a is in modulating artemisinin resistance.
We did not identify the ubp-1 1528D mutation, contrary to a study by Henriques et al. (21), which showed that its frequency increased post-ACT treatment, although it has been found at a prevalence of 7.4% in Ghana (55). A different ubp-1 SNP (R3138H), outside the region we genotyped, has also been associated with artemisinin resistance on the Thai-Myanmar border (Cerqueira et al. [58]). Like ap2-mu and falcipain-2a, this calls for further work to understand the role of ubp-1 as a marker of ACT resistance.
We observed 2-fold increases or decreases in the frequencies of alleles and haplotypes of several genes, followed by a shift back to the original frequency in the durations spanning 2012/2013, 2015/2016, and through to 2017/2018. This could be due to the abrupt change in drug policy from SP to ACTs and hence the initial selection of long haplotypes followed by recombination breaking down the haplotype structure as the parasite adapts to changes in drug pressure. The lack of spatial variation in allele frequency is consistent with previous reports at this geographical scale (59); however, we had limited power to do a more detailed analysis beyond north and south Kilifi.
Following the introduction of ACTs in 2004, there has been a rapid increase in the CQ-sensitive population to near fixation, and this reignites the debate on the use of CQ for malaria treatment, such as in combination therapy. On the other hand, there is still a need for careful monitoring of the dhps A581G locus, since SP has proved useful in IPTp, significantly reducing morbidity in pregnant women (32). The decline in the novel marker (nfs), which potentially confers resistance to LM, contrary to the observations made in The Gambia, calls for additional studies to determine its role as a potential drug target. The artemisinin resistance-conferring SE Asian mutations in k13, such as C580Y, have not been identified in Kilifi, and many of the SNPs occurred in the N-terminal region of the gene with no evidence of drug selection. Consequently, due to a lack of the validated molecular markers of artemisinin and lumefantrine resistance, there appears to be no problem of resistance in the population; however, continued surveillance remains a requirement.

MATERIALS AND METHODS
Study design. We used parasite DNA extracted from frozen blood obtained in 1995/1996, 1999/ 2000, 2006/2007, and 2012/2013 from a previous study (35) and from additional samples from 2015/2016 and 2017/2018 from this new study. One hundred fifty samples were collected for each time point, except for 2017/2018, in which 109 samples were used. Blood samples were obtained from patients presenting to Kilifi County Hospital with malaria (1 month to 15 years of age). The samples span 24 years of changing drug policy (Fig. 2). Ethical approval for this study was obtained from the ethics review committee of the Kenya Medical Research Institute under protocol number SERU 3149.
Sample preparation, PCR, and capillary sequencing. P. falciparum genomic DNA was extracted from frozen blood using the Qiagen DNA blood minikit (Qiagen, UK), as per the manufacturer's instructions. Amplicons were generated from the following genes using primers from previous studies and primers designed in this study (Table S1): crt (PF3D7_0709000), mdr1 (PF3D7_0523000), dhps (PF3D7_0810800), nfs (PF3D7_0727200), k13 (PF3D7_1347700), ap2mu (PF3D7_1218300), falcipain-2a (PF3D7_1347700), ubp-1 (PF3D7_0104300), and serine-tRNA ligase, putative gene (PF3D7_0717700). Additionally, four artemisinin resistance-predisposing mutations in arps10 (PF3D7_1460900), crt (PF3D7_0709000), fd (PF3D7_1318100), and mdr2 (PF3D7_1447900) were genotyped. Given that these mutations in these genes have been identified at low frequencies (Ͻ1%) in Africa (16), genotyping was first done for the 1995/1996 and 2015/2016 time points, and then the remaining time points were included if SNPs with Ͼ10% frequency were identified. We used the Expand high-fidelity PCR system, 0.5 l of template DNA, and primers and PCR conditions indicated in Table S1. The final reaction volume was 10 l, and the PCR amplification products were visualized on 1% agarose gels stained with RedSafe nucleic acid staining solution (iNtRON Biotechnology DR). PCR-negative samples were taken through a second and final round of PCR with 0.75 l of template DNA. Positive PCR products were purified using ExoSAP-IT (Thermo Fisher Scientific) and directly sequenced using the PCR primers, additional internal primers (Table S1), and BigDye Terminator v3.1 cycle sequencing kit v3.1 (Applied Biosystems, UK). Capillary sequencing was done at the International Livestock Research Institute (ILRI, Kenya) using an ABI 3730xl capillary sequencer (Applied Biosystems).
Sequence analysis. Sequence assembly was performed in CLC Main Workbench v7.9.1 (Qiagen, UK), and SNPs were identified and called based on the respective 3D7 reference sequences. Nucleotide positions which displayed a peak within a peak in the sequence chromatograms were noted as "mixed." Consensus sequences were extracted from the sequence assemblies using CLC Genomics Workbench v9.5.3 and used to construct multiple-sequence alignments in Clustal Omega v1.2.1 (60). SNP frequencies were calculated per gene per time point, and singletons were confirmed by an additional round of PCR and sequencing.
Statistical analysis. Nucleotide sequences were translated into amino acid sequences in AliView v1.26 (61). Haplotypes were then generated based on the amino acid residues from all the polymorphic codons that cut across all sequences and time points after excluding sequences with mixed bases. Data for crt and mdr1 from the 1995 to 2013 time points were obtained from a previously published study (35). The difference in the prevalences of alleles and haplotypes in pre-ACT and post-ACT periods was evaluated using a chi-square test. For this analysis, 2005/2006 was used as the point to divide the data into the pre-and post-ACT periods, since we do not expect complete ACT coverage following its formal introduction in 2004; hence, data for 2005/2006 were excluded from this part of the analysis. Additionally, it was the crossover point for the shift in high frequency from the crt 76T to crt K76 allele. The chi-square test for alleles included data for all SNPs with Ͼ10% frequencies at any time point, while for haplotypes, we included only the sequences that had data across all loci and with frequencies of Ͼ10% at any time point. Additionally, the chi-square test for haplotypes was conducted only for the two dominant haplotypes. All plots were generated using ggplot2 v3.1.1 (62) and ggpubr v0.2 packages (63) in R v3.6.0 (R Core Team, 2014). The K a /K s ratio was calculated for k13 by dividing the number of nonsynonymous substitutions per nonsynonymous site (K a ) by the number of synonymous substitutions per synonymous site over time.

SUPPLEMENTAL MATERIAL
Supplemental material for this article may be found at https://doi.org/10.1128/AAC .01067-19.