Temporal trends in prevalence of Plasmodium falciparum drug resistance alleles over two decades of changing antimalarial policy in coastal Kenya.

Molecular surveillance of drug resistance markers through time provides crucial information on genomic adaptations, especially in parasite populations exposed to changing drug pressures. To assess temporal trends of established genotypes associated with tolerance to clinically important antimalarials used in Kenya over the last two decades, we sequenced a region of the pfcrt locus encompassing codons 72-76 of the Plasmodium falciparum chloroquine resistance transporter, full-length pfmdr1 - encoding multi-drug resistance protein, P-glycoprotein homolog (Pgh1) and pfdhfr encoding dihydrofolate reductase, in 485 archived Plasmodium falciparum positive blood samples collected in coastal Kenya at four different time points between 1995 and 2013. Microsatellite loci were also analyzed to compare the genetic backgrounds of parasite populations circulating before and after the withdrawal of chloroquine and sulfadoxine/pyrimethamine. Our results reveal a significant increase in the prevalence of the pfcrt K76 wild-type allele between 1995 and 2013 from 38% to 81.7% (p < 0.0001). In contrast, we noted a significant decline in wild-type pfdhfr S108 allele (p < 0.0001) culminating in complete absence of this allele in 2013. We also observed a significant increase in the prevalence of the wild-type pfmdr1 N86/Y184/D1246 haplotype from 14.6% in 1995 to 66.0% in 2013 (p < 0.0001) and a corresponding decline of the mutant pfmdr1 86Y/184Y/1246Y allele from 36.4% to 0% in 19 years (p < 0.0001). We also show extensive genetic heterogeneity among the chloroquine-sensitive parasites before and after the withdrawal of the drug in contrast to a selective sweep around the triple mutant pfdhfr allele, leading to a mono-allelic population at this locus. These findings highlight the importance of continual surveillance and characterization of parasite genotypes as indicators of the therapeutic efficacy of antimalarials, particularly in the context of changes in malaria treatment policy.

Introduction

Understanding the evolution of resistance-associated genes is crucial in evaluating drug efficacy. Molecular trends underlying such phenotypes as tolerance or susceptibility can be effectively monitored by exploring loci selectively influenced by antimalarial pressure. Consequently, a temporal molecular map can be constructed from the adaptive changes observed in these markers over time, particularly in populations exposed to changing drug pressures. Extensive use of chloroquine (CQ) as a monotherapy led to significant increase in levels of resistance across many malaria-endemic countries prompting policy changes. In Africa, Malawi (in 1993) was the first to replace CQ with sulfadoxine/pyrimethamine (SP) as the first-line treatment for uncomplicated malaria, shortly followed by Kenya (in 1998) and a number of other countries (Shretta et al., 2000; Kamya et al., 2002; Eriksen et al., 2005). However, widespread reports of declining SP efficacy at the coast (Nzila et al., 2000) and other parts of Kenya (van Dillen et al., 1999; Omar et al., 2001) soon emerged prompting another first-line antimalarial policy change in 2004 (Amin et al., 2007) to the currently preferred Coartem™, an artemether–lumefantrine (AL) combination rolled out in government clinics since 2006.

Clinical resistance to CQ has been strongly associated with genetic replacements in the Plasmodium falciparum chloroquine resistance transporter, Pfcrt (PF3D7_0709000), with the lysine to threonine replacement at codon 76 (K76T) considered most critical (Fidock et al., 2000). However, the existence of chloroquine-sensitive (CQS) strains associated with K76T mutation suggests that other genes could also be involved in CQ resistance (Sa et al., 2009). Indeed, there is persuasive evidence that mutations in pfmdr1 (PF3D7_0523000), encoding the P. falciparum homolog of the human P-glycoprotein, are also involved in modulating CQ sensitivity as parasites bearing pfmdr1 86Y, 1034C, 1042D and 1246Y alleles have been shown to exhibit impaired transportation and accumulation of CQ into the food vacuole hence reduced CQ sensitivity (Koenderink et al., 2010). On the other hand, the molecular basis of resistance to SP in vitro has been linked to point mutations in the parasite’s dihydrofolate reductase, pfdhfr (PF3D7_0417200) and dihydropteroate synthase, pfdhps (PF3D7_0810800) genes (Peterson et al., 1988; Triglia et al., 1997). Alterations in pfdhfr proceed stepwise, with the gatekeeper mutation from serine to asparagine at codon 108 (S108N) preceding subsequent changes at codons 50, 51, 59 and 164 that further compound the extent of resistance. Treatment failure with SP occurs when one or more mutations are also present in pfdhps (Wang et al., 1997; Hallett et al., 2006).

While the discontinuation of CQ use was expected to at least disrupt the selective pressure on pfcrt and pfmdr1, artemisinin partner drugs have been documented to exert opposing pressure on these loci in East Africa (Dokomajilar et al., 2006; Humphreys et al., 2007; Mwai et al., 2009a; Sisowath et al., 2009; Conrad et al., 2014). In fact, studies in Tanzania suggest that AL selects for lumefantrine (LM)-tolerant parasites (Martensson et al., 2005; Sisowath et al., 2005; Malmberg et al., 2013a). Interestingly, these putatively LM-tolerant parasites have wild-type pfmdr1 (asparagine at codon 86) and, in some cases, wild-type pfcrt (lysine at position 76) alleles, both associated with CQ susceptibility. Mutations that render an organism resistant to drugs may be associated with loss of fitness and consequently, parasite populations with these mutations would be outgrown by their drug-sensitive counterparts when drug pressure is withdrawn (Levy, 1994). CQ has now been out of clinical use for 15 years in Kenya while SP, for nearly half the time – though still effective for intermittent preventive treatment in pregnancy (IPTp) with a nation-wide coverage of 30–39% as at 2011 (van Eijk et al., 2013). This is an index of the proportion of pregnant women protected by IPTp, computed as the total number of protected births divided by the total number of malaria-exposed births. Complete or partial reversion to CQS alleles has been reported in Malawi (Kublin et al., 2003; Frosch et al., 2014), Tanzania (Temu et al., 2006), western Kenya (Eyase et al., 2013), and the Kenyan coast (Mwai et al., 2009b; Mang’era et al., 2012), among other sites. On the other hand, antifolate-resistant genotypes has remained high along the coast (Kiara et al., 2009), presenting a threat to the long-term future of IPTp. However, in Kilifi – a malaria endemic area along coastal Kenya, the overall temporal structure of drug resistant alleles especially with the introduction of AL pressure and intermittent deployment of SP is yet to be determined. On the backdrop of such changing antimalarial pressures on the parasite population since 1998, it would be instructive to also characterize the genetic background flanking the aforementioned loci. This has been previously employed in profiling the spatial origins and dissemination of resistant alleles (Wootton et al., 2002; Roper et al., 2004) and more recently in determining if the parasite populations between different time points are genetically comparable (Laufer et al., 2010; Nwakanma et al., 2014). In this study, we sought to assess the frequency of alleles of the drug resistance genes pfmdr1, pfcrt, and pfdhfr during a 19-year period of changing antimalarial policy and compare parasites’ genetic backgrounds. Our results provide crucial insights into the parasites’ genomic adaptations as they adjust to a landscape of changing drug pressure and underline the need for comprehensive genotypic data that can be used to audit the therapeutic efficacy of drugs in clinical use and those previously withdrawn.

Materials and methods

Sample population and ethics statement

Isolates were selected from a database of frozen blood samples by identifying malaria-positive samples collected before administration of treatment from patients presenting to Kilifi District Hospital with malaria. Samples clustering within 4 time points spanning 19 years of changing drug policy i.e. 1995, 1999/2000, 2006/2007 and 2012/2013 were randomly chosen for analysis (Fig. 1). The extraction and use of these samples was reviewed and approved by the Ethics Review Committee of Kenya Medical Research Institute under protocol number SSC 2533.

Fig. 1

Flow chart showing the four P. falciparum populations spread out through a 19-year time scale and punctuated with changing drug policies. Seminal reports on various milestones in the epidemiology of antimalarial resistance in Africa are also highlighted.

DNA preparation and PCR

Parasite genomic DNA was extracted from frozen erythrocytes using the automated QIAxtractor system (Qiagen, UK) according to the manufacturer’s instructions and eluted DNA frozen at −20°C. A segment of pfcrt exon 2 encompassing codons 72–76 was amplified using primers described elsewhere (Chan et al., 2012). To determine the presence of any additional mutations (presumably due to drug pressure), we amplified full-length pfmdr1 and pfdhfr genes using High Fidelity Taq polymerase (Roche). Details of PCR conditions and amplification primers sequences are available in Supplementary Table 1. The generated PCR products were visualized on 1% agarose gels under ultraviolet illumination.

Sequencing

PCR products were purified using ethanol precipitation and directly sequenced using the PCR and additional sets of internal primers, BIG DYE terminator chemistry v3.1 (Applied Biosystems, UK) and an ABI 3130xl capillary sequencer (Applied Biosystems, UK). Nucleotide positions which displayed a peak within a peak in the electropherogram were noted as a “mixed” but excluded from further analysis. Sequences were assembled, edited and aligned using SeqMan and MegAlign (DNASTAR, Madison, WI). SNPs were identified and using their corresponding amino acids, haplotypes were defined. The sequencing primers are also listed in Supplementary Table 1.

Microsatellite analysis

We employed 8 microsatellite markers to compare CQS samples collected during CQ use (1995) and after withdrawal (2013). These comprised loci flanking pfcrt at −45.1 kb, −17.7 kb, −4.8 kb, −4.5 kb, 1.5 kb, 3.9 kb, 18.8 kb and 45.3 kb. We also interrogated the genetic relatedness of parasites bearing the triple mutant pfdhfr allele, before SP introduction (1995) and in 2013 by genotyping microsatellite loci flanking the gene at −7.5 kb, −4.4 kb, −3.8 kb, −0.06 kb, 0.1 kb, 0.45 kb, 1.3 kb, and 5.8 kb. In addition, we further analyzed 8 putatively neutral microsatellite loci selected from a set of 12 previously described (Anderson et al., 1999). The pfcrt and pfdhfr microsatellite positions, primers and cycling conditions were adopted as elsewhere (Alam et al., 2011) with slight modifications as detailed in Supplementary Table 2. Microsatellite allele scoring was done using the GeneMapper software, version 3.7 (Applied Biosystems), with samples presenting multiple alleles at any of the loci omitted from downstream analyses. Summary indices including allelic diversity and allelic richness were calculated using FSTAT Version 2.9.3.2. Allelic diversity was calculated for all microsatellite loci based on the allele frequencies, using the formula for ‘expected heterozygosity’ H = [n/(n − 1)][1 − ∑p2], where n is the number of isolates analyzed and p represents the frequency of each different allele at a locus. H has a potential range from 0 (no allele diversity) to 1 (all sampled alleles differ).

Statistical analysis

All statistical analyses were conducted using STATA version 11 (Stata, College Station, TX). Changes in the prevalence of alleles over time were evaluated for statistical significance using χ2 statistics for trend. For haplotype analysis, we excluded minority alleles (<5% frequency) as it is difficult to make meaningful statements about rare alleles. Logistic regression was used to assess temporal changes in allele prevalence and statistical significance confirmed using trend analysis for proportions (ptrend test). The odds ratio (OR) with corresponding 95% confidence interval (CI) represent the relative change between 2 years. To assess the extent of genetic diversity between neutral alleles and those under selection, differences between mean H values were compared using Student’s t-test. The significance level was assessed at 5% for all analyses.

Results

Prevalence of drug-resistant alleles

We evaluated 485 samples from microscopically-confirmed falciparum malaria cases, clustering within 4 time points spanning 19 years of changing drug policy (1995; n = 96, 1999/2000; n = 131, 2006/2007; n = 139 and 2012/2013; n = 119), to determine the prevalence of pfcrt, pfdhfr and pfmdr1 alleles in Kilifi. Of the 485 samples, 366 (75.5%) yielded single-genotype pfcrt sequences, 246 (50.7%) for pfdhfr and 231 (47.6%) for pfmdr1 as shown in Supplementary Table 3. The rest of the samples in each group either had multiple alleles (mixed genotype) or poor sequence data (1995; n = 62 [21.5%], 1999/2000; n = 136 [34.6%], 2006/2007; n = 209 [50.1%] and 2012/2013; n = 165 [46.2%]). Sequences were submitted to GenBank and are available under the accession codes KJ689814–KJ690044 for pfmdr1 and KJ715966–KJ716212 for pfdhfr. Note that only a short fragment of pfcrt was sequenced and as such did not meet the length criteria for submission to GenBank. Our data shows a resurgence in the proportion of the wild-type pfcrt alleles over time, with ∼82% of the isolates bearing the pfcrt K76 allele (and by extension the C72/V73/M74/N75/K76 haplotype) in 2013 compared to 38% in 1995. The CQR alleles were most prevalent in 1999/2000 (pfcrt C72/V73/74I/75E/76T = 93.2%). We did not observe any other polymorphisms along the entire 1827 bp pfdhfr sequence apart from the known N51I, C59R and S108N mutations (Supplementary Table 3). All the mutations existed as 6 different haplotypes with ∼89% of all the isolates at least polymorphic at one locus. Overall, the triple-mutant pfdhfr 51I/59R/108N haplotype was the most predominant at ∼55% followed by the double mutant pfdhfr 51I/C59/108N at 19.1%. The frequency of the mutant parasites was already high (in 1995) even before introduction of SP, as evidenced by a 76.6% prevalence among double and triple mutants. No novel mutations were observed on the full-length pfmdr1 either, with all the samples polymorphic at only codons 86, 184 and 1246. We also observed repeat sequence variation in the poly-asparaginated linker region as shown in Supplementary Table 4. Overall, the double-mutant 86Y/Y184/1246Y (32.9%) pfmdr1 haplotype was predominant, followed by the wild-type haplotype, N86/Y184/D1246 (23.8%). 86Y/184F/D1246 and 86Y/184F/1246Y haplotypes had the lowest frequencies at 0.87% and 0.43%, respectively. The pfmdr1 N86 SNP was observed in linkage with pfmdr1 D1246 (χ2 = 64.02; p < 0.0001) and pfcrt K76 (χ2 = 33.38; p < 0.0001) throughout the study period.

Temporal trends in SNPs and haplotype prevalence

Pfcrt

There was a statistically significant increase in pfcrt K76 between 1995 and 2012/2013 from 38% to 81.7% (Odds Ratio = 7.3; [95% CI 3.55–15.0]; p < 0.0001) as shown in Table 1. To probe the potential influence of AL on the prevalence of pfcrt alleles, we assessed our data on the basis of pre- and post-introduction of AL in 2006, and report a significant increase in pfcrt K76 from 24% pre-AL (1995 and 1999/2000 group) to 76% post-AL (Odds Ratio = 6.8; [95% CI 4.3–10.9]; p < 0.0001).

Table 1

Temporal trends in the prevalence of resistance-related haplotypes in Kilifi between 1995 and 2013.

Haplotype	1995 (freq %)	1999/00 (freq %)	2006/07 (freq %)	2012/13 (freq %)	Parametric trend test slope	Parametric trend test p-value
pfcrt _CVMNK∗	34.9 (30)	6.5 (7)	47.6 (50)	77.0 (67)	0.03	<0.0001
pfcrt _CVIET	57.0 (49)	88.9 (96)	49.5 (52)	17.2 (15)	−0.03	<0.0001
	N = 79	N = 103	N = 102	N = 82
pfmdr1 _NFD	29.1 (16)	7.5 (6)	14.0 (6)	31.9 (15)	0.00498	0.2228
pfmdr1 _NYD∗	14.6 (8)	12.5 (10)	14.0 (6)	66.0 (31)	0.02749	<0.0001
pfmdr1 _YYD	20.0(11)	16.3 (13)	34.9 (15)	2.1 (1)	−0.00574	0.1486
pfmdr1 _YYY	36.4 (20)	57.5 (46)	23.3 (10)	0.0 (0)	−0.02673	<0.0001
	N = 55	N = 75	N = 41	N = 47
pfdhfr _NCS∗	21.3 (16)	17.7 (11)	0 (0)	0 (0)	−0.013	<0.0001
pfdhfr _IRN	53.3 (40)	37.1 (23)	74.5 (38)	67.3 (35)	0.013	0.0041
pfdhfr _ICN	16.0 (12)	27.4 (17)	15.7 (8)	19.2 (10)	−0.0005	0.8921
pfdhfr _NRN	9.3 (7)	17.7 (11)	9.8 (5)	13.5 (7)	0.0006	0.8451
	N = 75	N = 62	N = 51	N = 52
pfmdr1 _NFD +p fcrt _K76	17.0 (8)	1.5 (1)	17.9 (5)	34.2 (13)	0.013	0.0013
pfmdr1 _NFD +p fcrt _76T	14.9 (7)	7.7 (5)	3.6 (1)	0.0 (0)	−0.008	0.0091
pfmdr1 _NYD +p fcrt _K76∗	12.8 (6)	3.1 (2)	7.1 (2)	52.6 (20)	0.023	<0.0001
pfmdr1 _NYD +p fcrt _76T	4.3 (2)	9.2 (6)	14.3 (4)	13.2 (5)	0.005	0.1421
pfmdr1 _YYD +p fcrt _76T	14.9 (7)	15.4 (10)	25.0 (7)	0.0 (0)	−0.007	0.0886
pfmdr1 _YYY +p fcrt _K76	14.9 (7)	1.5 (1)	17.9 (5)	0.0 (0)	−0.004	0.1616
pfmdr1 _YYY +p fcrt _76T	21.3 (10)	61.5 (40)	14.3 (4)	0.0 (0)	−0.022	<0.0001
	N = 47	N = 65	N = 28	N = 38

Wild-type alleles are indicated with an asterisk (∗) and significant p-values highlighted bold. The negative sign on the values of the slope of the trend denote a decrease in frequency over time.

Pfmdr1

After excluding the rare alleles, our pfmdr1 analysis was only restricted to 4 haplotypes. We observed an increase in the 86Y/Y184/1246Y haplotype from 36.4% in 1995 to 57.5% in 1999/2000 (Odds Ratio = 2.8; [95% CI 1.4–5.7]; p = 0.005). There was also a significant increase in the wild-type N86/Y184/D1246 allele from 14.6% in 1995 to 66.0% in 2012/2013 (Odds Ratio = 11.4; [95% CI 4.3–29.8]; p < 0.0001) (see Table 1), an observation also noted in 2006/2007 although the increase was not statistically significant (Odds Ratio = 1.1; [95% CI 0.4–3.6]; p = 0.827). We also investigated the prevalence of haplotypes with regard to introduction of AL in 2006. The prevalence of the wild-type pfmdr1 allele also increased post AL introduction, with 67.3% of the samples bearing the wild-type N86/Y184/D1246 allele in 2006/2007 and 2012/2013 compared to 32.7% in 1995 and 1999/2000 (Odds Ratio = 1.9; [95% CI 0.5–6.3]; p = 0.319).

Pfdhfr

There was a significant trend in the decline of wild-type pfdhfr N51/C59/S108 (ptrend; p < 0.0001) allele culminating in the ultimate fixation of the mutant variants in 2012/2013. Hence in contrast to pfcrt, we observed no decline in the frequency of pfdhfr mutant alleles after replacement of SP with Coartem™ (Table 1). There was however a significant increase in the triple mutant 51I/59R/108N haplotype from 37.1% at SP introduction in 1999 to 67.3% in 2013 (Odds Ratio = 3.5; [95% CI 1.6–7.6]; p = 0.002). We did not observe any pfdhfr alleles with the 164L mutation in our analysis.

Pfcrt in combination with Pfmdr1

There is in vitro evidence strongly associating combined pfcrt 76T and pfmdr1 86Y genotypes with high CQ IC50s but increased sensitivity to LM (Mwai et al., 2009a). We therefore examined the trends in the selection of different allelic combinations in these two genes. We observed a significant increase in the combined wild-type pfmdr1 N86/Y184/D1246 + pfcrt K76 alleles from 12.8% in 1995 to 52.6% in 2012/2013 (Odds Ratio = 7.6; [95% CI 2.6–22.1]; p < 0.0001). Though the pfmdr1 N86/Y184/D1246 + pfcrt 76T combination also rose during this period, the increase was not significant (Odds Ratio = 3.4; [95% CI 0.62–18.7]; p = 0.157). We noted a significant decrease in pfmdr1 86Y/Y184/1246Y + pfcrt 76T from 21.3% in 1995 to 0% in 2012/2013 (ptrend; p < 0.0001) but an increase of the same allele from 21.3% to 61.5% during the period around extensive CQ use (Odds Ratio = 5.9; [95% CI 2.5–14.0]; p < 0.0001), only later declining in 2006/2007 and 2012/2013 (Table 1)

We characterized microsatellite polymorphisms at 8 loci flanking pfdhfr in all 74 evaluable triple mutant samples (n = 39 in 1995 and n = 35 in 2012/2013) and pfcrt in all 95 evaluable wild-type samples (n = 30 in 1995 and n = 65 in 2012/2013). The triple mutant pfdhfr and the wild-type pfcrt alleles were used for this temporal microsatellite analysis since these were the two forms showing evidence of significant positive selection over time. We also typed 8 neutral microsatellite markers in 141 samples (n = 47 in 1995 and n = 94 in 2012/2013) to illustrate the selection landscape and diversity around pfdhfr and pfcrt. Among the resistant pfdhfr parasites, our results reveal substantial allele-sharing before and after SP introduction (Fig. 2a). Markers distal to pfdhfr (−7.5 kb, −4.4 kb, 1.3 kb and 5.8 kb) exhibited greater diversity, consistent with the tenets of selective sweep (Nair et al., 2003). The mean expected heterozygosity (H ± SD) at the 8 loci around pfdhfr were low but comparable between 1995 (0.23 ± 0.1) and 2012/2013 (0.21 ± 0.08) as shown in Supplementary Table 5a. Compared to the neutral loci in 1995 and 2012/2013, these means were significantly lower (unpaired Student’s t-test; p < 0.0001), thus affirming the selective sweep around the pyrimethamine-resistant (PYR-R) alleles as shown in Fig. 3A. In contrast, there was high diversity among the CQS parasites (Fig. 2b) with mean H around the C72/V73/M74/N75/K76 alleles recorded in 1995 (0.70 ± 0.09) and 2013 (0.72 ± 0.09) as shown in Supplementary Table 5b. This was comparable to the means around the neutral markers around the same period, 0.77 ± 0.05 in1995 and 0.73 ± 0.12 in 2013 (Fig. 3B).

Fig. 2a

Microsatellite haplotypes around a 13.3 kb region flanking pfdhfr in parasites collected in 1995 (n = 39) and 2013 (n = 35) bearing the triple mutant allele. The figure shows extensive allele-sharing among the samples and similarities in genetic backgrounds between Kenyan samples and Southeast Asian strains. Microsatellite sizes are indicated in nucleotide base pairs and alleles identical to triple mutant P. falciparum K1 strain are shown in gray shading.

Fig. 3

(a) and (b) Change in diversity in microsatellite loci around triple mutant pfdhfr and wild-type pfcrt alleles in 1995 and 2013. Panel A shows the variation in expected heterozygosity (H) (y-axis) around the triple mutant pfdhfr (51I/59R/108N) in 1995 and 2013. The dashed (1995) and solid (2013) horizontal lines represent the estimates of the mean He of neutral loci examined at both times and visually depict the low diversity in mutant pfdhfr relative to neutral loci. Panel B illustrates variation around CQS parasites and is juxtaposed to (A) to show difference in diversity. Diversity around wild-type pfcrt alleles was notably comparable to that of neutral alleles as evidenced by proximity of the plot to the mean He around neutral microsatellite (horizontal lines) in contrast to the plot of diversity around pfdhfr that lie much lower.

Fig. 2b

Microsatellite haplotypes around a 90.4 kb region flanking pfcrt in parasites collected in 1995 (n = 30) and 2013 (n = 65) bearing the wild-type allele. Alleles identical to the wild-type P. falciparum 3D7 strain are shown in gray shading. The high diversity among the wild-type samples is clearly evidenced by the number of unique alleles at each locus.

Discussion

Following widespread CQ resistance, Kenya switched to SP as the first-line antimalarial against uncomplicated malaria in 1998 (Shretta et al., 2000). However, clinical resistance to SP soon prompted the adoption of artemisinin-based combination therapy (ACT) with Coartem™ as the first-line regimen and SP relegated to intermittent use during pregnancy (Amin et al., 2007). Our results confirm the progressive resurgence of CQS parasite populations in Kilifi, and suggest that the mutant pfdhfr alleles are maintained at high frequencies a decade after withdrawal of SP. We have further demonstrated extensive genetic heterogeneity in CQS parasites before and after CQ withdrawal, in contrast to the near-clonal triple mutant pfdhfr population during the same period.

The significant increase in pfcrt C72/V73/M74/N75/K76 allele in 2006/2007 and 2012/2013, coincides with the period after CQ withdrawal. This fits the expectations of a fitness cost-associated selection model where the survival advantage conferred to CQ-resistant parasites in the presence of the drug is lost on withdrawal of CQ pressure as observed in The Gambia (Ord et al., 2007), and is consistent with recent reports from coastal Kenya (Mwai et al., 2009b; Mang’era et al., 2012) and other parts of Africa (Laufer et al., 2010; Ndiaye et al., 2012; Malmberg et al., 2013b). This increase is also attributable to AL use which has been demonstrated to select for LM-tolerant parasites, which coincidentally harbor the wild-type pfcrt K76 allele (Martensson et al., 2005; Sisowath et al., 2005; Henriques et al., 2014). Though these trends reveal a recovery in the frequency of CQS parasites from 38% to ∼82% in 19 years (1995–2013) compared to 5–40% in 13 years (1993–2006) from the same population (Mwai et al., 2009b), this rate is however still lower compared to changes in some parts of Africa (Table 2). This could be due to extensive use of the CQ analog, amodiaquine (AQ), in parts of Kenya (including Kilifi) as second-line antimalarial even before CQ withdrawal and long after SP introduction (Amin et al., 2007), maintaining selective pressure on CQR parasites. Also, CQ was still widely retailed for self-medication even 4 years after its official withdrawal (Amin et al., 2007) hence also maintaining pressure on the resistant variants, thus highlighting the implications of unsynchronized cross-over in treatment policies. Recent reports on the re-emergence of CQS parasites have prompted debate on the possible re-introduction of CQ (albeit in combination with another drug) in the event of widespread resistance to LM. Our analysis of a different set of microsatellite markers to those used in Malawi (Laufer et al., 2010) showed that high genetic diversity is maintained in CQS populations between 1995 and 2013, similar to the observations in Malawi. It therefore seems CQS diversity may not have been entirely extinguished under decades of drug pressure, as indicated by the high mean expected heterozygosity values comparable to the neutral loci (Fig. 3b). These findings also corroborate observations from Ghana where a similar degree of diversity was noted among the CQS parasites (Alam et al., 2011) (see Table 3).

Table 2

Comparative pfcrt K76 and pfmdr1 N86 allele frequency changes in various malaria-endemic African countries relative to withdrawal and introduction of CQ and ACTs, respectively.

Country	CQ Withdrawal/ACT	Year of Study	% Frequency Change		Reference
	Introduction		Pfcrt _K76	Pfmdr1 _N86
Malawi	1993/2008	1992–2000	15.0–87.0	69.0–75.0	Kublin et al. (2003)
Mozambique	2002/2008	2006–2010	3.90–67.6	25.3–69.1	Raman et al. (2011)
Zanzibar	2001/2003	2003–2010	4.00–37.0	25.0–48.0	Froberg et al. (2012)
Mozambique	2002/2008	2009–2010	43.9–66.4	64.7–84.1	Thomsen et al. (2013)
Tanzania	2001/2006	2006–2011	49.0–85.0	14.0–61.0	Malmberg et al. (2013b)
Uganda	2000/2004	2003–2012	0.00–17.0	10.0–51.0	Mbogo et al. (2014)
Senegal	2003/2006	2000–2009	27.6–40.5	67.0–78.0	Ly et al. (2012)
The Gambia	2004/2008	2000–2008⁎	23.7–40.7	21.7–74.2	Nwakanma et al. (2014)

This Gambian study was conducted between 1984 and 2008. Over subsequent survey time points, proportions of isolates with resistant pfcrt 76 and pfmdr 86 alleles increased progressively to peak in 2000. This, therefore, is the point from which we begin to analyze the frequency change from mutant to wild-type alleles.

Table 3

Allelic diversity (expected heterozygosity, He) and allelic richness (Rs) at 8 neutral microsatellite loci in various chromosomes within the genome in samples collected at two different time points.

Sample population – 1995 (n = 47)			Sample population – 2013 (n = 94)
Microsatellite locus	Allelic richness (Rs)	Expected heterozygosity (He)	Allelic richness (Rs)	Expected heterozygosity (He)
Population sampled (n = number of individual isolates)
Poly-α	12.0	0.819	13.0	0.839
PfPK2	9.0	0.739	12.0	0.861
ARA2	11.0	0.791	11.0	0.688
TA87	9.0	0.702	18.0	0.847
TA42	14.0	0.760	12.0	0.649
2490	11.0	0.715	13.0	0.661
TA60	11.0	0.785	9.0	0.510
TA109	14.0	0.849	18.0	0.780
Mean ± SD	11.4 ± 1.9	0.77 ± 0.05	13.3 ± 3.2	0.73 ± 0.12

The sample population represents the evaluable genotypes in the two time points. Though the original total samples available for genotyping was 96 and 119 in 1995 and 2012/2013 respectively, samples presenting >1 allele at any of the 8 loci were excluded leading to the loss of a substantial number of samples (ultimately n = 47 and n = 94 in 1995 and 2012/2013, respectively). This sampling variance, however, did not occasion any significant difference between the mean He in 1995 and 2013.

Our results also indicate a steady increase in the prevalence of pfmdr1 N86 and D1246 alleles while pfmdr1_184F only slightly increased in 2006 and 2013. This wane in the pfmdr1 86Y and 1246Y mutant alleles, coupled with the rise of the Y184F mutation, alludes to disparate selective pressure on this locus, eliminating some mutations while driving others to high prevalence. Indeed, there is compelling evidence implicating AL in these trends. In Zanzibar, a 2.7 fold increase in frequency of pfmdr1 N86 was observed after 42 days following treatment with AL (Sisowath et al., 2005) while pfmdr1 N86 and 184F alleles have recently been associated with in vivo selection by AL in east (Dokomajilar et al., 2006; Gadalla et al., 2011) and west Africa (Lekana-Douki et al., 2011; Dahlstrom et al., 2014). In addition, pfmdr1 184F has been found to be under selection among parasite populations in Cambodia (Vinayak et al., 2010), where artemisinin delayed parasite clearance has been described. These have implications for the useful therapeutic life of Coartem™ since the increase of parasites harboring combined wild-type pfmdr1 N84/Y186/D1246 and pfcrt K76 alleles in the population could be the first step in the selection of LM-tolerant parasites which would consequently form the backdrop for developing Coartem™ resistance, perhaps mediated by changes at other loci.

The high prevalence of PYR-R parasites in our population mirrors results from other studies using samples from this location (Kiara et al., 2009; Mwai et al., 2009b) and could be partly due to SP use in IPTp as the PYR component of the drug selects for fitter drug-tolerant variants. However, the high parasite proportions already bearing the resistant genotypes before its introduction absolve intermittent SP use alone as primary driver for the high mutant frequencies. Selection pressure could possibly have been enhanced by similar-acting antifolate combination drugs, notably cotrimoxazole. This drug possesses only mild antimalarial potency but is a common prophylactic prescription against opportunistic respiratory tract infections among HIV patients (White, 2004), hence may also have perpetuated the mutant populations. Despite reports of the pfdhfr 164L mutation in western (McCollum et al., 2006; Hamel et al., 2008) and coastal Kenya (Kiara et al., 2009), this allele was absent in our analysis. However, there is need to continually monitor pregnant women and pediatric cases which are potential sources of amplification and dissemination of parasites bearing this allele due to their predisposition to IPT. The reduced mean heterozygosity in the loci flanking pfdhfr relative to the neutral loci indicates that the pfdhfr 51I/59R/108N haplotype has undergone rapid expansion in coastal Kenya. Most samples with this allele bore microsatellite profiles identical to those of Southeast Asian strains, supporting earlier assertions of a Southeast Asian origin of PYR-R east African parasites (Roper et al., 2004). Nonetheless, we also observed few unique profiles specific to Kilifi, which could either be pfdhfr 51I/59R/108N indigenes or recombinant hybrids of the Southeast Asian and local parasites. Despite the high proportion of parasites harboring resistance-associated mutations, SP-IPT has been effective in preventing the adverse consequences of malaria on maternal and fetal outcomes in Africa (World Health Organization, 2012). However, recent reports on alarming rates of recrudescence following SP-IPTp (Mutabingwa et al., 2009; Moussiliou et al., 2013) coupled with our microsatellite data revealing clonality in pfdhfr parasite genotypes that can endure SP pressure, raise concern about the continued use of SP in IPT strategies.

Conclusion

We have shown increases in the pfcrt C72/V73/M74/N75/K76 and pfmdr1 N84/Y186/D1246 alleles over time in Kilifi after withdrawal of CQ and introduction of AL. The temporal selection of CQS alleles which are also putatively LM-tolerant raises concern on the effectiveness of LM as a partner drug since it could potentially form the starting point for AL resistance. We have also captured the early events in the dynamics of the resistant pfdhfr alleles through to their fixation in the population. The significance of such retrospective surveillance brings into focus the need for temporal monitoring of the recently identified artemisinin resistance marker (Ariey et al., 2014) to track its progression in populations. We concede that the study would have been even more comprehensive had it been powered and designed to also explore adaptive copy number evolution in pfmdr1 and pfdhfr over time. This phenomenon, in pfmdr1, has been associated with reduced sensitivity to LM (Price et al., 2006) while GTP-cyclohydrolase 1 (encoding the first enzyme in the folate pathway) has been shown to exhibit antifolate-selected copy number polymorphism (Nair et al., 2008). Nonetheless, this report reiterates the need for continued surveillance while seeking more suitable alternative drugs or a vaccine.