Partner-Drug Resistance and Population Substructuring of Artemisinin-Resistant Plasmodium falciparum in Cambodia

Abstract

Plasmodium falciparum in western Cambodia has developed resistance to artemisinin and its partner drugs, causing frequent treatment failure. Understanding this evolution can inform the deployment of new therapies. We investigated the genetic architecture of 78 falciparum isolates using whole-genome sequencing, correlating results to in vivo and ex vivo drug resistance and exploring the relationship between population structure, demographic history, and partner drug resistance. Principle component analysis, network analysis and demographic inference identified a diverse central population with three clusters of clonally expanding parasite populations, each associated with specific K13 artemisinin resistance alleles and partner drug resistance profiles which were consistent with the sequential deployment of artemisinin combination therapies in the region. One cluster displayed ex vivo piperaquine resistance and mefloquine sensitivity with a high rate of in vivo failure of dihydroartemisinin-piperaquine. Another cluster displayed ex vivo mefloquine resistance and piperaquine sensitivity with high in vivo efficacy of dihydroartemisinin-piperaquine. The final cluster was clonal and displayed intermediate sensitivity to both drugs. Variations in recently described piperaquine resistance markers did not explain the difference in mean IC₉₀ or clinical failures between the high and intermediate piperaquine resistance groups, suggesting additional loci may be involved in resistance. The results highlight an important role for partner drug resistance in shaping the P. falciparum genetic landscape in Southeast Asia and suggest that further work is needed to evaluate for other mutations that drive piperaquine resistance.

plasmodium, partner drug, piperaquine, mefloquine, drug resistance, ACT, malaria, population genetics, pfmdr1, kelch, ex vivo susceptibility

Introduction

Mortality due to malaria nearly halved over the last 15 years, due in part to the adoption of artemisinin-based combination therapies (ACT). ACT is comprised of a potent, short-acting artemisinin derivative paired with a longer-acting partner drug. This allows ACT to be administered over just 3 days, but leads to a long terminal phase of partner-drug monotherapy (Fairhurst 2015; Taylor and Juliano 2014). Thus, as artemisinin resistance increases and parasites take longer to clear, the probability of developing resistance to the partner drug also increases, leading to ACT failure. Resistance to ACT represents a major threat to malaria-control efforts worldwide and has now been reported in five countries in Southeast Asia (Ashley et al. 2014). Efforts are underway to stem the rising tide of resistance in the region, and new antimalarial regimens are desperately needed.

Intensive research into the mechanisms of resistance to artemisinin, mefloquine (MQ) and piperaquine (PPQ) have revealed molecular markers of resistance to the last two ACTs used in the region, artesunate-mefloquine (AS-MQ) and dihydroartemisinin-piperaquine (DHA-PPQ) (Leang et al. 2015; Spring et al. 2015). Artemisinin resistance is associated with mutations in kelch (K13) (Ariey et al. 2014; Miotto et al. 2015; Takala-Harrison et al. 2012; MalariaGEN Plasmodium falciparum Community Project 2016). MQ resistance has for some time been associated with increases in Plasmodium falciparum multidrug resistance gene 1 (pfmdr1) copy number. More recently, two mutations have been associated with resistance to PPQ, a copy-number variation encompassing plasmepsin II and III and a Glu415Gly substitution within an exonuclease (Amato et al. 2017; Witkowski et al. 2017). However, it is unknown whether these mutations explain all ACT failures in the region (Mukherjee et al. 2017).

The emergence of drug resistance in the region has led to strong population substructuring, which has been ascribed to artemisinin resistance (Ariey et al. 2014; Miotto et al. 2015; Takala-Harrison et al. 2012; MalariaGEN Plasmodium falciparum Community Project 2016). However, previous population genomics analyses have not explored the potential role of ACT partner drugs, in conjunction with artemisinin derivatives, in shaping these subpopulations. Increasing PPQ resistance has been accompanied by increasing MQ sensitivity (Leang et al. 2015; Chaorattanakawee et al. 2015, 2016) which has led to interest in the use of triple ACT (TACT), combining an artemisinin derivative with both MQ and PPQ (Maxmen 2016). TACT is predicated upon the premise that parasites with resistance to both partner drugs have not yet developed, and that PPQ and MQ may exert opposing selection pressures if used simultaneously. Understanding the relationship between subpopulations in regard to partner drug resistance may provide evidence to support the use of TACT in the region.

To better understand how resistance to artemisinin and these two partner drugs (MQ and PPQ) has evolved in Cambodia, we investigated the genetic architecture of 78 P. falciparum parasites drawn from three provinces in western Cambodia between 2011–2014, including 60 from a cohort with a >50% rate of DHA-PPQ failure (Chaorattanakawee et al. 2015; Spring et al. 2015). Whole-genome sequencing revealed stark genetic structuring among the samples. We integrated the genetic findings with ex vivo PPQ and MQ resistance profiles as well as in vivo responses to DHA-PPQ therapy, finding resistance to ACT partner drug correlated with different parasite subpopulations. This differentiation of piperaquine versus mefloquine-resistant parasite subpopulations supports a triple ACT strategy. A single parasite subpopulation in particular was associated with high PPQ IC₉₀s and DHA-PPQ treatment failure. Although recently described candidate markers for PPQ resistance are predominantly mutant in this highly resistant subpopulation, they were also present in a subpopulation with intermediate piperaquine sensitivity and a lower rate of DHA-PPQ failure. Thus, further studies to identify other loci involved in PPQ resistance are warranted. Further molecular surveillance and in vitro studies of these subpopulations could help identify additional PPQ resistance loci and help to determine the likelihood that P. falciparum resistant to TACT regimens currently being deployed will emerge in the near future.

Materials and Methods

Study Population

We studied parasites collected from 78 adults with uncomplicated P. falciparum infection from three provinces of Western Cambodia: Oddar Meanchey in the North (n = 61), Battambang (n = 8) in the West and Kampot in the South (n = 9). These samples were collected between June 2011 and December 2013 and represent a subset of samples sequenced for a previous study (Parobek et al. 2016). Sixty of the samples in Oddar Meanchey came from an in vivo efficacy study of DHA-PPQ in which 54% of patients developed recrudescent parasitemia requiring rescue therapy with AS-MQ (clinical trial protocol WR1877 [NCT01849640]) (Spring et al. 2015). An additional 18 P. falciparum clinical isolates from inside and outside of Oddar Meanchey Province were selected at random from an in vitro surveillance study (protocol WR1576), for which clinical follow-up data was not collected (Chaorattanakawee et al. 2015). Both studies were approved by the Cambodian National Ethics Committee for Health Research and the Walter Reed Army Institute of Research Institutional Review Board. The molecular studies described herein were approved by the University of North Carolina Institutional Review Board (14-0419). All subjects provided informed consent.

Drug Susceptibility Testing

Immediate ex vivo drug susceptibilities to PPQ and MQ were measured using a histidine-rich protein 2 (HRP-2) enzyme-linked immunosorbent assay (Chaorattanakawee et al. 2015; Rutvisuttinunt et al. 2012; Saunders et al. 2015).

Whole-Genome Sequencing and CNV Determination

Leukodepleted DNA from 91 P. falciparum infections was used for whole-genome sequencing. The ratio of parasite DNA to host DNA was determined using a quantitative PCR assay, and isolates with ≥20% P. falciparum DNA were sequenced on the HiSeq2000 Sequencing System (Illumina, San Diego, CA) (Beshir et al. 2010). Sequence reads were aligned to the P. falciparum 3D7 (v3) genome using the bwa mem algorithm. Isolates with 5-fold coverage at ≥60% of the genome were considered for variant calling, resulting in calls for 78 isolates. Coverage over coding regions was not used as a criteria for isolate inclusion, though all included samples did have robust coverage over coding regions (see supplementary Results, Supplementary Material online). Variant calling was performed using the GATK HaplotypeCaller utility. Further details are described in the supplementary Methods, Supplementary Material online. WGS data was used for calling all variation, both SNP and CNV, as previously described (Parobek et al. 2016). CNVs at select loci were also evaluated by qPCR: pfmdr1 copy-number calls were integrated from previously published data (Chaorattanakawee et al. 2015; Lim et al. 2015; Eastman et al. 2011), and plasmepsin-II and -III CNs were ascertained experimentally (see supplementary Methods and table S3, Supplementary Material online). In both cases, isolates with with CN ≥1.5 were considered to represent increased copy number. Quantitative PCR data was used for all analysis.

Assessment of Population Structure and Demography

Whole-genome variant data were used to determine the principal components of genetic variation (Jombart and Ahmed 2011; Paradis 2010). Cluster assignments were determined using a nonparametric k-means approach, and clusters were assigned a unique “Cambodian Parasite” (CP) identifier. To create the network tree, variant calls were read into R using pegas (Paradis 2010) and a neighbor-joining tree was constructed using ape (Paradis et al. 2004) and bootstrapped to 100 replicates. The clonality of each population cluster was rigorously tested by comparing point estimates of linkage disequilibrium decay with bootstrap replicates drawn randomly from all sequenced isolates and by assessing pairwise total number of single nucleotide polymorphism (SNP) differences between each pair of isolates within a subpopulation. Observed minor allele frequencies were calculated and compared with the frequencies expected from an idealized population, modeled using ms (Hudson 2002). The folded site-frequency spectrum was used to infer demographic parameters and likelihoods under a number of two-population models (Gutenkunst et al. 2009). The tested scenarios were no split (Static, the two clusters are drawn from the same population), split with isolation (Iso, the two populations split at some past time T), and split with migration (IM, the two populations split at some past time T with subsequent migration between the split populations). Simulations have found that the 7000+ variants we report here are an order of magnitude more numerous than needed for robust demographic inference (Robinson et al. 2014). 100 bootstrapped replicates were constructed, and each of the above models was fitted to each of these replicates 100 times. Genewise F_ST between subpopulations was calculated and rendered in R (R CoreTeam, Vienna, Austria) (Pfeifer et al. 2014; Gu et al. 2014). For further details on population-genetic methods, see the supplementary Materials, Supplementary Material online.

Statistical Analysis

Median PPQ IC₉₀ and MQ IC₅₀ values were compared between PCA clusters using the Kruskal–Wallis test. The Wilcoxon–Mann–Whitney test was used to compare the median PPQ IC₉₀ values of isolates with and without select candidate resistance markers, as well as the median MQ IC₅₀ values of recrudescent versus nonrecrudescent parasites. Fisher’s exact test was used to compare the frequency of increased pfmdr1 copy number amongst recrudescent parasites and those that did not recrudesce and, in a post-hoc analysis, the frequency of increased pfmdr1 copy number, the exonuclease E415G mutation, and increased plasmepsin-II and -III copy number between PCA groups, with Bonferonni corrected thresholds for statistical significance to account for multiple comparisons. We analyzed 42-day cumulative risk of recrudescent malaria via survival-analysis with Kaplan–Meier curves and compared the risk by cluster and by PPQ resistance markers using Cox proportional hazards models. Statistical analyses were completed using SAS 9.4 (SAS Institute, Cary, NC) and figures were generated using R.

Results

Whole-Genome Sequencing, Network Analysis and Demographic Inference Show Subpopulations Undergoing Epidemic Expansion

Whole-genome sequencing of 78 P. falciparum clinical isolates produced, on average, >10-fold coverage in 94.7% of coding regions and >5-fold coverage in 99.5% of coding regions. Stringent variant calling within these samples identified 7,228 high-quality, genome-wide SNPs. The allele-frequency spectrum for the population showed two notable features: 1) an overrepresentation of rare alleles, and 2) an overrepresentation of intermediate-frequency alleles relative to the neutral expectation (supplementary fig. S1, Supplementary Material online). An excess of rare alleles can be observed in expanding populations, as mutations that have occurred since the onset of expansion have not had sufficient time to spread through the population (Maruyama and Fuerst 1984). An excess of intermediate frequency alleles suggests the existence of two or more genetically unique subpopulations (Hurst 2009; Nielsen and Rasmus 2005; Kreitman 2000). To explore the latter finding, we performed principal components analysis (PCA) of genetic variation to summarize the underlying population structure among these samples. The first two principal components (PCs) explained 52% of the variance and revealed four distinct population subclusters (fig. 1A and B). Components beyond the third PC were less informative (supplementary fig. S2A–D, Supplementary Material online). K-means clustering using the top three principal components assigned individuals to four genetically distinct clusters, which we named CP1 (n = 17), CP2 (n = 18), CP3 (n = 20), and CP4 (n = 23) (supplementary fig. S3, Supplementary Material online). We found that genetic dissimilarity was not explained by geography, as all clusters except for one (CP4) contained parasites from more than one geographic location (supplementary fig. S4, Supplementary Material online). Therefore, this population substructuring is driven by a factor that supersedes genetic drift.

Fig. 1.

—Principal components analysis demonstrates distinct clusters with non-overlapping resistance profiles. In the upper panels, nodes are scaled according to piperaquine IC₉₀ and mefloquine IC₅₀ values, with the largest nodes representing IC₉₀ values of 1,000 nM or greater for piperaquine. Panel A displays PPQ IC₉₀ values, and panel B displays MQ IC₅₀ values. Blue samples comprise cluster CP1, red CP2, green CP3, and yellow CP4. In the lower panels, IC₉₀ values by cluster are depicted for PPQ in panel C and IC₅₀ values for MQ in panel D. Boxplots display the median (horizontal line), interquartile range (box), and the lowest and highest data point within 1.5 times the IQR (whiskers) of the values [values outside this range are plotted as individual points]. One PPQ IC₉₀ outliers were clipped from cluster CP4 in panel C and labeled with an asterisk (*). There was a statistically significant difference among the four groups for PPQ (Kruskal–Wallis P = 0.0005) IC₉₀ values and for MQ (P = 0.004) IC₅₀ values. Noise has been added to the PCA to reduce overplotting. Samples are plotted only if information is available for them, thus some positions may appear different. Abbreviations: PC, principal component; IC₅₀, half-maximal inhibitory concentrations.

Network analysis shows that CP2 has, on average, longer terminal branch lengths and deeper internal nodes, findings that suggest an older most recent common ancestor for CP2 members than for CP1, CP3, or CP4 members (fig. 2). In addition, the outlying subpopulations appear more clonal by network analysis. Therefore, we assessed the population diversity and linkage disequilibrium (LD) within each subpopulation. Analysis of LD in each of the four populations confirm a higher degree of linkage disequilibrium within the CP1 and CP4 subpopulations versus CP2 (supplementary fig. S5, Supplementary Material online). The extreme clonality of CP3 prevented calculation of r² for this subpopulation. The number of pairwise SNP differences within each cluster also revealed that the ACT-resistant subgroups—and particularly CP3 and CP4—are clonal or nearly clonal, with greatly reduced within-group diversity compared with the ancestral-like CP2 population (supplementary fig. S6, Supplementary Material online). Strikingly, the greatest number of SNP differences observed in pairwise comparisons among CP3 isolates was four. Lastly, we compared the genetic diversity in our clusters to the recently reported genetic architecture of artemisinin-resistant P. falciparum (Miotto et al. 2015, 2013). Based on 7,214 SNPs common to both data sets, we determined that our populations fall within the scope of diversity previously described (supplementary fig. S7A–C, Supplementary Material online), and contain less genetic diversity than their overlayed population (supplementary table S1, Supplementary Material online). These data suggest that CP1, CP3, and CP4 may represent epidemic expansion of parasite subpopulations, potentially from CP2 and potentially due to drug resistance (Daniels et al. 2013).

Fig. 2.

—Neighbor-joining network tree of P. falciparum isolates. A rooted and bootstrapped (100 replicates) neighbor-joining network tree of the Cambodian P. falciparum isolates is shown. The specific isolates are colored by CP group: CP1 blue, CP2 red, CP3 green, and CP4 yellow. The tree is rooted with 3D7, an isolate of African origin.

To confirm this finding, and to determine the timing and direction of the subpopulation splits, we fit various two-population demographic scenarios to the allele frequency spectra of CP subgroup pairs (see supplementary Materials, Supplementary Material online). To test a broad range of demographic models and parameter space, we used a diffusion-based inference approach based on the multi-population allele frequency spectrum (supplementary fig. S8, Supplementary Material online), rather than a coalescent-theory based Monte Carlo approach (Gutenkunst et al. 2009). This approach has been used for efficient inference of recent demographic events in populations within nonmodel species including Plasmodium (McCoy et al. 2014; Chang et al. 2012). These two-population scenarios assume an ancestral coalescent state and explicitly test which subpopulation is descended from a founder event. Fitting two-population models to the CP2 versus CP1/CP3/CP4 cluster pairs demonstrated that CP1 (n = 17) and CP4 (n = 23) likely derived from a CP2 ancestral-like population (n = 18), and that CP2–CP1 split occurred prior to CP2–CP4 split (supplementary table S2, Supplementary Material online). However, in both cases, the best fit models suggest split with back migration resulting in potential genetic mixing of the populations.

Distinct Subpopulations Correlate with Ex Vivo Partner-Drug Resistance

Our isolates follow a pattern of alternate resistance profiles, in which samples with high PPQ IC₅₀ values exhibit low MQ IC₅₀ values, and vice versa (fig. 3). We found that these differing PPQ- and MQ-resistance profiles correlate with the genetic substructure found by whole genome sequencing. Mean ex vivo resistance to PPQ was highest in one cluster, CP4, whereas mean ex vivo resistance to MQ was highest in another cluster, CP1. This is depicted in figure 1A and B, where points in the PCA represent individual clinical isolates and are sized according to their respective PPQ IC₉₀ and MQ IC₅₀ values. PPQ IC₉₀ values were analyzed because they are more predictive of resistance than IC₅₀ values for PPQ (Chaorattanakawee et al. 2016). Comparative analysis of median PPQ IC₉₀ and MQ IC₅₀ values by cluster shows different drug-resistance profiles among the four groups (fig. 1C and D). The median PPQ IC₉₀ values were: 62.4 nM (interquartile range 45.5–70.7) for CP1, 201.4 nM (97.8–360.4) for CP2, 288.6 nM (85.4–678.4) for CP3, and 667.4 nM (136.3–2523.7) for CP4 (Kruskal–Wallis P = 0.0005). The median MQ IC₅₀ values were: 143.0 nM (IQR 53.9–173.9) for CP1, 9.1 nM (1.2–63.7) for CP2, 47.8 nM (23.0–62.7) for CP3, and 33.1 nM (4.8–42.7) for CP4 (Kruskal–Wallis P = 0.004).

Fig. 3.

—Piperaquine and mefloquine IC₅₀ values associated with parasite subpopulations in Cambodia. Subjects who had both piperaquine and MQ IC₅₀ values are displayed. Those who later had an episode of recrudescence are displayed using solid orange triangles, whereas those who did not have an episode of recrudescence or did not have follow-up data available are displayed using solid navy circles. Two isolates that had recrudescence with very high piperaquine IC₅₀ values (307.0 and 397.1 nM) are represented with hollow symbols.

A Single Parasite Subpopulation Is Associated with DHA-PPQ Treatment Failures

The in vivo phenotypes for parasite clearance half-life and treatment efficacy were available for 60 of the isolates. Patients who failed DHA-PPQ therapy were disproportionately infected with parasites from the PPQ-resistant CP4 group. Among the 58 subjects followed for 7 weeks after DHA-PPQ treatment, there were no recrudescence events in CP1 (n = 11, 0%), two in CP2 (n = 8, 25%), four in CP3 (n = 16, 25%), and 18 in CP4 (n = 23, 78%). Survival analysis demonstrated that subjects infected with parasites from cluster CP4 were more likely to suffer recrudescence after treatment with DHA-PPQ than those from other clusters (hazard ratio 7.9, 95% CI 3.0–20, P < 0.001) (fig. 4). Though determination of recrudescence by standard WHO genotyping procedures is potentially complicated in populations that are highly clonal, the fact that initial and recurrent parasitemias in these patients shared a genotype would suggest a reinfection by the same clonal line if it was not a true recrudescence. CP4 parasites also harbored the C580Y mutation and seemed to display the longest parasite clearance half-lives, consistent with artemisinin resistance. However, C580Y enrichment and parasite clearance half-lives were not significantly different in the CP3 group (supplementary fig. S9, Supplementary Material online), despite lower rates of DHA-PPQ failure.

Fig. 4.

—Recrudescence after DHA-PPQ by cluster. Kaplan–Meier curves demonstrate the risk of recrudescence by cluster (Wald P < 0.001 for CP4 compared with other clusters). Red crosses represent patients censored for withdrawal or development of new infection.

Candidate Partner Drug Resistance Markers Explain Ex Vivo MQ Resistance but Do Not Differentiate Intermediate and High Ex Vivo PPQ Resistance

We used the sequence data to investigate single nucleotide polymorphism markers of artemisinin, PPQ and MQ resistance in the isolates. In addition, we evaluated copy number variation (CNV) of pfmdr1 and plasmepsin-II and -III using both whole-genome sequence data and real-time PCR. Relative to real-time PCR, sequence data had a sensitivity and specificity for CNV detection of 82.4% and 100% for pfmdr1 and 93.3% and 92.5% for plasmepsin-II and -III, respectively. Given inconsistent genome coverage characteristic of P. falciparum short-read sequencing projects, we chose to use real-time PCR data rather than whole genome sequencing data for CNV analyses.

Investigation of molecular markers of partner drug resistance showed an expected increase in pfmdr1 CNV in the MQ-resistant CP1 parasites (fig. 5) (Chaorattanakawee et al. 2015). Newly described candidate markers of PPQ resistance (Witkowski et al. 2017; Amato et al. 2017), including the plasmepsin-II and -III copy number variant and the exonuclease E415G mutation, were associated with treatment failure; subjects with mutant alleles at both markers were more likely to suffer recrudescence than those who lacked either or both markers (hazard ratio 16.8, 95% CI 2.2–127, P = 0.006, fig. 6). There was a range of one to three plasmepsin-II and -III copies detected in our isolates. However, these alleles, which discriminated between the populations with lower PPQ IC₉₀ values (CP1 and CP2) and those with higher PPQ IC₉₀ values (CP3 and CP4), did not discriminate between the intermediate PPQ resistance (CP3) and higher PPQ resistance (CP4) populations (fig. 5). This suggests that other genetic loci may drive this difference. Lastly, all parasites in the study were wild type at Plasmodium falciparum chloroquine resistance transporter (pfcrt) loci, at residues 97, 343, and 353, previously described by Duru et al. (2015).

Fig. 5.

—Proportion of isolates bearing partner drug resistance markers by cluster. The proportion of subjects with increased pfmdr1 copy number (A), increased plasmepsin-II and -III copy number (B) and the exonuclease E415G mutation (C) are shown. Comparison of CP3 and CP4 did not reveal statistically significant differences in the proportion of isolates bearing plasmepsin-II and -III copy number or exonuclease E415G mutation, with P values of 0.06 and 1.0, respectively (Bonferroni corrected threshold for significance of 0.0083).

Fig. 6.

—Recrudescence after DHA-PPQ by PPQ resistance markers. Kaplan–Meier curves demonstrate the risk of recrudescence by PPQ resistance markers (Wald P = 0.006 for subjects with both a plasmepsin-II and -III copy number variant and the exonuclease E451G mutation compared with others). Red crosses represent patients censored for withdrawal or development of new infection.

Clinical and Molecular Data Support Alternating Sensitivities to Partner Drugs

Available clinical and molecular data support the hypothesis that DHA-PPQ treatment failures in our cohort are MQ-sensitive. In the original clinical study, 46 subjects who failed DHA-PPQ treatment and developed recurrent malaria were treated with AS-MQ, per national guidelines (Spring et al. 2015; Chaorattanakawee et al. 2015). After this rescue treatment, only one subject (2.2%) had a second recurrent malaria infection during limited follow-up (median 15 days, range 4–33 days). In that subject, PCR genotyping demonstrated that the initial recurrence was a new infection, not a recrudescence due to DHA-PPQ failure. Thus, no early treatment failures occurred in 46 patients when AS-MQ was given for DHA-PPQ failure, though follow-up times were insufficient to rule out late treatment failures. Additionally, data from the original study showed only one of 35 (2.9%) recrudescent isolates among DHA-PPQ treatment failures had increased pfmdr1 CN (associated with MQ resistance), whereas 10 of 59 (16.9%) isolates from patients who successfully cleared parasites after DHA-PPQ treatment had increased copy number (P = 0.05) (Spring et al. 2015; Chaorattanakawee et al. 2015). Accordingly, lower MQ IC₅₀ values were measured in DHA-PPQ treatment failures (median 32.3 nM, interquartile range 9.3–51.6) compared with cures (48.0 nM, 21.0–108.1) (P = 0.03).

Extended K13 and Partner Drug Resistance Marker Haplotypes Are Shared among Parasite Subpopulations

In addition to differing partner-drug resistance profiles, our subpopulations were also correlated with specific K13 (PF3D7_1343700) genotypes. C580Y, the most common K13 mutation in Cambodia, was found in all CP3 and CP4 isolates (supplementary fig. S10 and S1, Supplementary Material online). The second most common K13 mutation, R539T, occurred exclusively in CP1, but this group also contained C580Y mutations. Lastly, CP2 contained a few different K13 mutations, including many isolates with C580Y. This pattern of subpopulations associated with specific K13 subtypes is similar to previous reports of “KH” subpopulations in Cambodia defined by K13 alleles (Miotto et al. 2015), and suggests that population fracturing may have occurred along artemisinin resistance subtypes. We investigated SNPs occurring 50 kb upstream and downstream of C580Y. Although we did find multiple C580Y haplotypes in the overall population consistent with the C580Y mutation arising on multiple genetic backgrounds, a single identical extended C580Y haplotype predominated and was shared across all clusters (fig. 7A). This suggests that despite their substructuring, parasites from all clusters contain C580Y mutations from a single origin, consistent with recent evidence of a dominant haplotype sweeping the region (Imwong et al. 2017). Thus, the separation of CP3 and CP4 cannot be explained by K13-associated artemisinin resistance. A similar analysis of extended haplotypes in CP1, the MQ-resistant population, revealed a single R539T haplotype (fig. 7B) accompanied by two C580Y haplotypes that are both shared with the diverse CP2 population (fig. 7A).

Fig. 7.

—Extended haplotypes surrounding artemisinin, piperaquine, and mefloquine drug resistance loci. The extended haplotypes of all SNPs extending 50kb upstream and downstream around K13 in C580Y bearing parasites (A), around K13 in R539T bearing parasites (B), around plasmepsin-II and -III (C) and around pfmdr1 (D) are shown. In Panel A, one haplotype is shared among all four subpopulations. A second haplotype is shared between CP1 and CP2. In Panel B, all K13 R539T mutations occurred within CP1 on a single haplotype. In Panel C, one dominant haplotype from CP3 and CP4 is shared among all four subpopulations. Parasites with an elevated plasmepsin-II and -III copy number are marked with an orange bar within their subpopulation. Parasites with no evident plasmepsin-II and -III duplication are marked with a purple bar within their subpopulation. In Panel D, a single extended haplotype dominates in CP1 parasites, irrespective if there is pfmdr1 duplication or not, though a highly similar haplotype is found in some nonduplicated CP4 parasites. Parasites with an elevated pfmdr1 copy number are marked with an orange bar within their subpopulation. Parasites with no evident pfmdr1 duplication are marked with a purple bar within their subpopulation. SNPs are numbered relative to the center of the gene and dark blue boxes represent nonreferent alleles. The gene is marked by the red shaded box.

Extended haplotypes surrounding partner drug resistance genes show different profiles. For the plasmepsin-II and -III copy number variant, a pattern similar to that at K13 emerges. All plasmepsin CNVs occur on a single haplotype, which predominates in CP3 and CP4 and which also occurs in the other two subpopulations (fig. 7C). Two CP2 isolates bear the plasmepsin CNV and are associated this extended haplotype, though additional CP2 isolates bear the plasmepsin haplotype without the CNV itself. Conversely, the extended haplotypes surrounding pfmdr1 within CP1 are not shared with any other subpopulation and are relatively conserved within CP1 whether or not a parasite had a detected duplication (fig. 7D).

Genomic Regions Differentiating CP3 and CP4 May Be Linked to High DHA-PPQ Clinical Failure Rate in CP4

Though this study was underpowered to detect statistically significant associations between individual genetic variants and clinical phenotypes, we nonetheless used population genetic approaches to explore genomic differences between CP3 and CP4. We chose to compare these two populations as there was little genetic diversity between CP3 and CP4, but a large clinical phenotype difference between subpopulations. Pairwise F_ST of CP3 versus CP4 revealed large genomic regions with genes that were either identical or highly divergent (supplementary fig. S12, Supplementary Material online). Importantly, as suggested by supplementary figures S10 and S11, Supplementary Material online and figure 7, the genomic regions containing candidate PPQ resistance loci are identical between the two populations, suggesting other loci may drive the clinical phenotype difference between these groups. Pairwise F_ST calculated between CP4–CP2 and CP3–CP2 subpopulations showed a less dichotomous pattern, with a more diverse range of values (supplementary fig. S12, Supplementary Material online). In addition, differences between CP1 and CP2 showed less striking divergence between the groups (supplementary fig. S12, Supplementary Material online) as suggested in our network analysis (fig. 2).

To explore which variants in the CP3–CP4 divergent regions contribute to the failure phenotype, we investigated the contribution of each SNP (determined by loading values) to principal component two (PC2), as PC2 explains nearly all the genetic difference between these two subgroups (supplementary fig. S13, Supplementary Material online). This can be seen by comparing the PC1–PC2 PCA (fig. 1A and B;supplementary fig. S2B, Supplementary Material online), where the populations split along PC2, to the PC1–PC3 PCA (supplementary fig. S2C, Supplementary Material online), where the populations collapse on each other. We found that the majority of variation between CP3 and CP4 is driven by 1,006 variants, including 330 missense or nonsense mutations exclusive to CP4, which are located in 206 genes (supplementary Appendix S2 and S3, Supplementary Material online). After Bonferroni correction, gene-ontology analysis showed no enrichment of specific biological processes among these genes.

Discussion

The P. falciparum population in Southeast Asia has been the focus of intense malaria control and elimination efforts over the past two decades. This has shaped a complex population structure marked by fractured subpopulations. Most work describing this genetic architecture has focused on artemisinin resistance as a selective force, showing that subpopulations associate with specific K13 mutations (Miotto et al. 2015, 2013). However, until recently, little attention has been given to the role of partner drug resistance in shaping these subpopulations. Rising ACT resistance in Southeast Asia leading to the failure of the last two ACT regimens in Cambodia makes it critical to understand how partner drugs contribute to changes in the parasite population. Here, we show that the P. falciparum population substructure in Cambodia correlates with partner-drug resistance and represents epidemic expansion of nearly clonal parasite clades. We also show that current candidate markers for PPQ resistance do not fully explain resistance and further studies of one subpopulation, CP4, will likely shed light on additional loci involved in clinical resistance.

At our study site in northwestern Cambodia, the first where high-grade DHA-PPQ resistance was reported, these subpopulations correlate with different ex vivo and in vivo patterns of antimalarial resistance. The most concerning finding was that roughly four-fifths of all of the parasites in CP4 subpopulation recrudesced after DHA-PPQ therapy. This group also had the highest PPQ IC₉₀ values. Surprisingly, parasites within this group carried mutations in the newly described PPQ resistance markers (Witkowski et al. 2017; Amato et al. 2017) at a similar frequency as CP3, a group with lower PPQ IC₉₀ values and only 25% treatment failure after DHA-PPQ (equivalent to CP2, which had even fewer resistant parasites by these molecular markers). Differences between CP3 and CP4 for other antimalarial IC₅₀ values were not significant other than for quinine (supplementary fig. S14, Supplementary Material online). F_ST comparisons of SNPs between CP3 and CP4 populations showed large genomic regions of near perfect homology contrasting with regions of high divergence (supplementary fig. S12, Supplementary Material online). These highly divergent genomic regions suggest—but cannot pinpoint—genetic drivers that may confer increased fitness to CP4 parasites in the face of DHA-PPQ therapy (supplementary Appendices S2 and S3, Supplementary Material online). Though we saw no enrichment for biological processes among the genes with mutations unique to CP4 parasites, there are nonetheless intriguing polymorphisms exclusive to this highly resistant lineage. For example, parasites of the CP4 group bore three missense polymorphisms in PF3D7_0419900, a phosphatidylinositol 4-kinase, which has been investigated as a potential drug target for P. falciparum (Ghidelli-Disse et al. 2014; McNamara et al. 2013) and been associated with drug resistance in Asia (Cerqueira et al. 2017). In addition, among the SNPs that differentiate CP3 and CP4 (supplementary S2, Supplementary Material online), we identified SNPs in proximity, but not identical to other SNPs identified in multiple GWAS studies of PPQ resistance (MAL8:280180, MAL12:1787279, MAL12:1959796, MAL14:2385550, MAL14:2395752, and MAL14:2411942) (Wang et al. 2016; Amato et al. 2017). Additional studies involving larger numbers of isolates with elevated PPQ IC₉₀ values are needed to further understand genetic determinants of the phenotypic differences between CP3 and CP4 parasites and will help us understand how these subpopulations continue to evolve.

Though the genetic determinants of DHA-PPQ failure are not fully elucidated, understanding how these resistant subpopulations arose is of key interest. We hypothesized that the populations CP1, CP3, and CP4 were derived from the centrally positioned CP2 population. Past work has posited the opposite hypothesis—that the KHA population (which shares significant overlap with CP2) is an admixture of other populations circulating in Cambodia (Miotto et al. 2013). Because these contrasting hypotheses have never been tested, we used modeling to elucidate the demographic history of these subpopulations. The two-population demographic comparisons presented here suggest that CP2 is an ancestral-like population from which the other populations arose. Consistent with this, recently derived CP4 (and likely CP3) subpopulations have undergone rapid expansion despite malaria control efforts. As evidence, the network tree (fig. 2), diversity (supplementary fig. S6, Supplementary Material online), and LD (supplementary fig. S5, Supplementary Material online) analyses all suggest that CP1, CP3, and CP4 have lower diversity than CP2, with varying degrees of clonality. In addition, our isolates, which were collected relatively recently, show less genetic diversity than the overlaying subpopulations described by Miotto et al. (supplementary fig. S7, Supplementary Material online), suggesting a recent acceleration of subpopulation growth.

Beyond the demographic models, extended haplotypes surrounding drug resistance genes can also provide insights to the origins of mutations. Similar to previous work, our subpopulations associate with specific K13 mutations, with CP3 and CP4 isolates exclusively C580Y, CP1 primarily R539T, and the putative ancestral-like population CP2 a mix of C580Y, other infrequent kelch mutations, and wild-type isolates (supplementary S1, Supplementary Material online). In order to better understand the origins of the C580Y and R539T mutations, we evaluated the extended SNP haplotypes surrounding K13 (fig. 7A). Isolates with the R539T mutation exhibited a single haplotype within a solitary subpopulation, suggesting a single origin of the mutation amongst these isolates. For C580Y, a single extended haplotype was contained in CP3 and CP4, as well as being present in CP2. Although K13 mutations historically have arisen through multiple origins, consistent with “soft selective sweeps” of drug resistance, the more contemporary nature of these samples suggests a dominant haplotype is now spreading through in a “hard selective sweep.” This is consistent with other findings recently reported in Cambodia (Imwong et al. 2017) and with the split with migration results from our demographic models.

The pattern is similar for extended haplotypes surrounding plasmepsin-II and -III copy number variants. All plasmepsin CNVs occur on a single dominant haplotype, which predominates in CP3 and CP4 and which also occurs in the other two subpopulations (fig. 7C). Two CP2 isolates bear the plasmepsin-II and -III CNV and its associated haplotype, though additional CP2 isolates bear the plasmepsin-II and -III haplotype without the CNV itself. Thus, it is possible that the plasmepsin-II and -III CNV arose once in a CP2 isolate and was associated with a subsequent differentiation and expansion of CP3 and CP4 isolates. Consistent with this, the predicted breakpoints for the duplication in our isolates were highly conserved [Chr14: 289558 (±34) to 298819 (±16)] and align to a similar position to those previously reported (Amato et al. 2017). Irrespective of the demographic history of these loci, the fact that K13 mutations, the exonuclease E415G, and the plasmepsin-II and -III CNV exist at similar frequencies in CP3 and CP4 groups and at different frequencies in CP2 and CP3 groups—yet CP2 and CP3 groups have similar clinical failure rates while CP3 and CP4 groups have dramatically different clinical failures rates—suggests that these two loci do not fully explain DHA-PPQ clinical failure. Additional loci remain to be found.

The MQ-resistant parasites from CP1 exhibited relatively conserved extended haplotypes around pfmdr1 (fig. 7D), with a dominant haplotype and two additional highly related haplotypes. The dominant extended pfmdr1 haplotype exists in parasites with and without a duplication in pfmdr1 within this group. This dominant extended haplotype associated with ex vivo MQ resistance does not occur within other subpopulations. However, CP4 contains a highly similar haplotype (different at a single position). Similar to Nair et al. (2007), we saw multiple other genetic backgrounds surrounding pfmdr1 duplications among all four subgroups. The high prevalence of a single haplotype within the group exhibiting ex vivo MQ resistance is consistent with an artemisinin-resistant mutation (e.g., K13 R539T) arising on one of the pfmdr1 backgrounds and sweeping to a high frequency between the survey done by Nair and our study.

As many have remarked, ACT in its current form may no longer be a viable strategy in Southeast Asia (Shanks et al. 2015; Spring et al. 2015; Saunders and Lon 2016). Current trials of TACT regimens, which incorporate both PPQ and MQ as partner drugs (Maxmen 2016), are underway. These are predicated on the premise that parasites with resistance to both partner drugs have not yet developed and that the two drugs exert opposing selective pressure. Theoretically, opposing selective pressure between these drugs may occur. MQ resistance has been shown to sensitize parasites to chloroquine in vitro, consistent with historical observations of a reciprocal relationship between chloroquine and MQ resistance (Veiga et al. 2014; Barnes et al. 1992; van Es et al. 1994). PPQ and amodiaquine are 4-aminoquinolines that are related to chloroquine. It follows that PPQ and MQ may exert opposing selective pressures based on PPQ’s similarity to chloroquine. Opposing selective pressure is well described for other antimalarials but has yet to be documented for MQ and PPQ (Venkatesan et al. 2014; Conrad et al. 2014). Our findings provide some support for TACT by showing that resistance to PPQ and MQ tend to exist on distinct genetic backgrounds. However, this does not preclude the possibility that a single parasite could harbor resistance to both drugs. Indeed, one clinical isolate (from CP3) demonstrated moderate tolerance to both PPQ (IC₉₀ 1060.1 nM) and MQ (IC₅₀ 102.63 nM). This sample contained the C580Y K13 mutation, the exonuclease E415G mutation, and the plasmepsin-II and -III CNV, but lacked a duplication in pfmdr1. Additionally, the presence of increased pfmdr1 copy number in several parasites in the CP4 group suggests the possibility of growing MQ resistance. The importance of these findings to public health is unclear, but they could be a harbinger of TACT-resistant P. falciparum.

In conclusion, we identified a strong relationship between DHA-PPQ treatment failure and a specific parasite genetic background among Cambodian P. falciparum parasites and determined that parasites of different genetic backgrounds carry different partner drug resistance profiles. This suggests that the substructuring of parasite populations was likely driven by a combination of artemisinin and partner-drug resistance, rather than artemisinin alone. Our demographic models and network analysis suggest that these partner-drug resistant subpopulations split out of a central ancestral-like population and have been rapidly expanding; rather than admixture of multiple preexisting subpopulations into a central population. Importantly, we found that differential DHA-PPQ treatment failure rates are not fully explained by recently described genetic markers of PPQ resistance. Mutations in these markers failed to differentiate a subpopulation with high PPQ IC₉₀ and high in vivo failure (CP4) from one with modestly increased PPQ IC₉₀ and a similar in vivo failure rate to the least drug-resistant population (CP3). This strongly suggests that additional loci differentiate these two subpopulations may affect PPQ resistance. Future interventions and studies should be guided by an understanding of this evolutionary history and genetic landscape, which can be leveraged to refine our understanding of PPQ resistance and the spread of ACT resistance (Fairhurst 2015).

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Authors Contribution

C.M.P., J.B.P., J.T.L., C.L., D.L.S., S.R.M., and J.J.J. designed the study. C.M.P., J.B.P., L.N., J.T.L., J.J.J., C.A.L., C.L., M.D.S., D.L.S., S.C., and P.G. collected the data. C.M.P., J.B.P., N.F.B., L.N., S.C., P.G., E.B., C.A.L., C.L., D.L.S., J.T.L., and J.J.J. did experiments and analyzed the data. C.M.P., J.B.P., N.F.B., J.A.B., J.T.L., D.L.S. and J.J.J. interpreted the data and prepared the report. J.T.L., J.J.J., C.L., and D.L.S. oversaw the project.

Funding

This work was supported by the National Institutes of Health and the National Institute of Allergy and Infectious Diseases [5T32AI007151 to JBP, T32GM007092 and F30AI109979 to CMP, R01AI089819 to JJJ, and K08 AI110651 to JTL]; the Armed Forces Health Surveillance Center/Global Emerging Infections Surveillance and Response System; Military Infectious Disease Research Program; and the American Society of Tropical Medicine and Hygiene/Burroughs Wellcome Fund [fellowships to JTL and CMP]. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The corresponding author had full access to all the data in the study and had final responsibility for the decision to submit for publication.

Acknowledgments

We thank the Cambodian and Thai laboratory and clinical support staff, and all patients who volunteered for this study. PPQ, artesunate, and other antimalarials used in the in vitro drug susceptibility assays were provided by the Chemical Repository of the Walter Reed Army Institute of Research. The views expressed in this article are those of the authors and do not reflect the official policy of the Department of the Army, Department of Defense, or the U.S. Government. We would also like to thank Corbin Jones, Praveen Sethupathy, and Kristina De Paris for their insightful comments.

Author notes

Associate editor: Geoff McFadden

Data deposition: This project has been deposited at the NCBI Sequence Read Archive under the accession PRJNA292536.

†

Joint first authors.

‡

Joint senior authors.