Which genotype is most likely to survive in a region of the world that contains malaria?

Introduction

Artemisinin combination therapy, the use of a short-acting artemisinin derivative and a long-acting partner drug, is recommended worldwide for the treatment of Plasmodium falciparum malaria.

Treatment with dihydroartemisinin–piperaquine, a current front-line artemisinin combination therapy in Cambodia, Vietnam, Thailand, Myanmar, China, and Indonesia is now failing in Cambodian provinces where artemisinin resistance

has emerged and spread.

, 

, 

, 

, 

, 

, 

9

• Leang R
• Barrette A
• Bouth DM
• et al.

Efficacy of dihydroartemisinin–piperaquine for treatment of uncomplicated Plasmodium falciparum and Plasmodium vivax in Cambodia, 2008 to 2010.

• Crossref
• PubMed
• Scopus (141)
• Google Scholar

, 

, 

, 

12

• Leang R
• Taylor WR
• Bouth DM
• et al.

Evidence of Plasmodium falciparum malaria multidrug resistance to artemisinin and piperaquine in western Cambodia: dihydroartemisinin-piperaquine open-label multicenter clinical assessment.

• Crossref
• PubMed
• Scopus (209)
• Google Scholar

This situation probably arose because the survival of parasites after artemisinin exposure increases the chance that they develop spontaneous genetic resistance to piperaquine, survive declining piperaquine plasma concentrations, and propagate via mosquitoes to other human beings. Through such processes, artemisinin resistance, which has evolved across the Greater Mekong subregion,

, 

, 

threatens to compromise the efficacy of dihydroartemisinin–piperaquine and other artemisinin combination therapies in the global treatment and elimination of malaria.

The natural selection of low-frequency, pre-existing resistance to piperaquine might also be occurring and further contributing to this problem.

Research in context

Evidence before this study

We searched PubMed using the terms “artemisinin”, “piperaquine”, “resistance”, and “Cambodia” without any date or language restrictions. This search identified 26 articles, five of which reported original clinical trials that documented dihydroartemisinin–piperaquine failures in the treatment of uncomplicated Plasmodium falciparum malaria in Cambodia. Collectively, these studies reported failure rates of 11–54% for dihydroartemisinin–piperaquine in Pailin, Pursat, Oddar Meancheay, and Preah Vihear provinces in 2010–13, and associated an increased risk of treatment failure with several parasite genotypes and phenotypes: the presence of kelch23 gene mutations that confer artemisinin resistance; the presence of two mutations (MAL10:688956 and MAL13:1718319) linked to slow parasite clearance after artesunate treatment; the absence of mdr1 gene amplification, a molecular marker of mefloquine resistance; and decreased piperaquine susceptibility in vitro. These studies did not identify molecular markers of both dihydroartemisinin–piperaquine failure and piperaquine resistance.

Added value of this study

This study, which investigated the whole-genome sequences and piperaquine susceptibilities of 297 P falciparum clinical isolates from Cambodia, provides two new molecular markers of piperaquine resistance and dihydroartemisinin–piperaquine failures. The first marker is a single-nucleotide polymorphism, which encodes a glutamic acid-to-glycine aminoacid substitution, in a putative exonuclease gene. This exonuclease marker is associated with decreased piperaquine susceptibility in vitro, and a low rate (0·38) of dihydroartemisinin–piperaquine treatment efficacy. This marker is found predominantly in Cambodia but not in other areas of the world where dihydroartemisinin–piperaquine has not been used. The second marker is an amplification of the plasmepsin 2 and plasmepsin 3 genes. This plasmepsin marker is also associated with decreased piperaquine susceptibility in vitro, and a low rate (0·41) of dihydroartemisinin–piperaquine treatment success. This marker is a strong candidate for causal mutation because plasmepsins encode aspartic proteases, which are involved in the parasite’s haemoglobin degradation pathway targeted by related antimalarial drugs such as chloroquine. These exonuclease and plasmepsin markers are commonly found together, along with kelch23 gene mutations, which confer artemisinin resistance, and a single copy of the mdr1 gene, which is associated with decreased piperaquine, but increased mefloquine, susceptibility in vitro.

Implications of all the available evidence

In Cambodia, emerging treatment failure with dihydroartemisinin–piperaquine is associated with the presence of plasmepsin and exonuclease markers. The plasmepsin marker is a strong candidate for a functional polymorphism and can be used to monitor the spread of dihydroartemisinin–piperaquine treatment failures in Cambodia and neighbouring countries in mainland southeast Asia, where dihydroartemisinin–piperaquine is a first-line treatment for malaria. Mapping the extent of this marker will enable national malaria control programmes to immediately recommend alternative therapies to improve the treatment of patients and drive the elimination of malaria.

Recent increases in treatment failures with dihydroartemisinin–piperaquine

, 

9

• Leang R
• Barrette A
• Bouth DM
• et al.

Efficacy of dihydroartemisinin–piperaquine for treatment of uncomplicated Plasmodium falciparum and Plasmodium vivax in Cambodia, 2008 to 2010.

• Crossref
• PubMed
• Scopus (141)
• Google Scholar

, 

, 

12

• Leang R
• Taylor WR
• Bouth DM
• et al.

Evidence of Plasmodium falciparum malaria multidrug resistance to artemisinin and piperaquine in western Cambodia: dihydroartemisinin-piperaquine open-label multicenter clinical assessment.

• Crossref
• PubMed
• Scopus (209)
• Google Scholar

and piperaquine 50% inhibitory concentrations (IC50s)

15

• Chaorattanakawee S
• Saunders DL
• Sea D
• et al.

Ex vivo drug susceptibility testing and molecular profiling of clinical Plasmodium falciparum isolates from Cambodia from 2008 to 2013 suggest emerging piperaquine resistance.

• Crossref
• PubMed
• Scopus (57)
• Google Scholar

suggest that piperaquine resistance has emerged in Cambodia. These findings, and the discovery that piperaquine-resistant parasites are sensitive to the former artemisinin combination therapy partner drug mefloquine,

have recently led Cambodia’s national malaria control programme and WHO to recommend artesunate–mefloquine as the first-line artemisinin combination therapy in ten Cambodian provinces, including Pursat and Preah Vihear. Molecular markers are urgently needed for large-scale surveillance programmes to predict treatment failures with dihydroartemisinin–piperaquine in Cambodia and other countries in the Greater Mekong subregion, and investigate the molecular mechanism of piperaquine resistance.

Several genetic variations have been associated with decreased piperaquine susceptibility: a single copy of the mdr1 gene has been associated with dihydroartemisinin–piperaquine treatment failures in Cambodian patients,

whereas amplification of a region downstream of mdr1 and the crt single-nucleotide polymorphism (SNP) C101F has been associated with raised piperaquine IC50s in vitro.

16

• Eastman RT
• Dharia NV
• Winzeler EA
• Fidock DA

Piperaquine resistance is associated with a copy number variation on chromosome 5 in drug-pressured Plasmodium falciparum parasites.

• Crossref
• PubMed
• Scopus (92)
• Google Scholar

These genetic variants are not useful as molecular markers, however, because the first is wild type, the second is very uncommon, and the third has not yet been observed in MalariaGEN P falciparum Community Project’s global catalogue of variation in clinical samples. Another study

found that a triple mutation (kelch23-C580Y, MAL10:688956 SNP, and MAL13:1718319 SNP) linked with slow parasite clearance is associated with a 5·4 times greater risk of dihydroartemisinin–piperaquine treatment failure in Cambodia’s Oddar Meancheay province; however, these SNPs were not linked to piperaquine resistance.

The purpose of this genome-wide association study (GWAS) was to discover molecular markers of dihydroartemisinin–piperaquine treatment failures in Cambodia. In designing it, we reasoned that such markers would associate with raised piperaquine IC50s, increase in prevalence over time in areas where artemisinin resistance is common, and be uncommon where dihydroartemisinin–piperaquine has not been used. Therefore, we sequenced the genomes and measured the piperaquine IC50s of parasites collected in 2010–13 from patients with P falciparum malaria in Pursat, Preah Vihear, and Ratanakiri, where the prevalences of kelch23 mutations, genetic markers for artemisinin resistance,

, 

, 

were 77%, 34%, and 11%,

and where the prevalences of dihydroartemisinin–piperaquine treatment failures were 46%, 16%, and 2%,

respectively, in 2012–13.

Methods

Study design and participants

To obtain samples for this GWAS study, we enrolled patients with uncomplicated P falciparum malaria into parasite clearance rate

, 

and drug efficacy

studies in 2010–13 in three provinces where piperaquine resistance is common (Pursat), emerging (Preah Vihear), or uncommon (Ratanakiri). Written informed consent was given by adult patients, or a parent or guardian of child patients. Protocols were approved by the Cambodian National Ethics Committee for Health Research and the National Institute of Allergy and Infectious Diseases Institutional Review Board and are registered with ClinicalTrials.gov, numbers NCT00341003, NCT01240603, and NCT01736319.

Using parasitised blood samples from these studies, we measured piperaquine IC50s ex vivo or after short-term culture in vitro, and obtained whole-genome parasite sequence data, whenever possible. The GWAS was designed to identify genetic markers of raised piperaquine IC50s, whereas a dihydroartemisinin–piperaquine efficacy study

was used to test for association between GWAS candidate markers and parasite recrudescence. PCR genotyping using msp1, msp2, and glurp microsatellites as genetic markers distinguished recrudescences from newly acquired infections.

A summary of samples, according to province of origin and year of collection, is shown in the appendix.

GWAS and survival analyses

The preparation, sequencing, genotyping, and phenotyping of samples are described in the appendix. The GWAS and correction for population structure were done using a linear mixed model, implemented in FaST-LMM

version 2.06.

We tested 11 630 SNPs with minor allele frequency greater than 0·033, using genotypes encoded as the number of non-reference alleles (0 or 1), and excluding heterozygous calls to minimise confounding effects of mixed infections. Piperaquine IC50 was used as the continuous dependent variable. When whole-genome sequence or phenotype data were available for both initial and recrudescent samples, we only analysed data from the recrudescent samples. Different approaches (using initial sample data only, using the phenotype and genotype of the same sample, averaging initial and recrudescent IC50s) did not alter results. A relationship matrix was calculated using a subset of 6678 SNPs with minor allele frequency greater than 3% (options: –maf 0·03), missing data rate below 25% (options: –geno 0·25), and unlinked (in windows of 100 SNPs, shifted forward by ten SNPs each time, removing one from each pair of SNPs with linkage disequilibrium >0·3; options: –indep-pairwise 100 10 0·3) and extracted using PLINK

version 1.07. In estimating the relationship matrix, we found that excluding proximal SNPs (within 10 kb or 100 kb from the tested variant) did not substantially affect results. Given the number of independent SNPs used, we applied Bonferroni correction to define a significance threshold of p of 8·6 × 10−7 or less for GWAS analyses. We also defined a suggestive threshold of p of 10−4 or less to identify relatively high-ranking loci.

To adjust for potential confounder effects, we treated the geographical origin of the sample, the presence of mdr1 amplification, and the presence of kelch23 resistance alleles collectively as covariates in the linear mixed model. These covariates reduced the genomic inflation factor λGC from 3·455 to 1·964 (appendix). We also treated genetic similarity across samples as a random effect to correct for confounding effects of population structure, which further reduced λGC to 1·106. Although this value is still slightly above 1·0, this is unlikely due to un-accounted population stratification. When we confined the estimation of λGC to only unlinked SNPs, its value dropped to 1·013, suggesting that residual inflation is due to extended homozygosity haplotypes in the samples and thus has little effect on the GWAS.

We did a survival analysis using the R package survival. We fitted a Cox proportional hazard regression model, in which treatment success (recrudescent vs non-recrudescent infection) represented the survival status, and we added the age of the patient (in years), parasitaemia on day 0 (log-scaled), and dose of piperaquine given (mg/kg, in five-unit increments starting at 35 mg/kg) as covariates.

21

WorldWide Antimalarial Resistance Network (WWARN) DP Study Group
The effect of dosing regimens on the antimalarial efficacy of dihydroartemisinin-piperaquine: a pooled analysis of individual patient data.

• Crossref
• PubMed
• Google Scholar

We then included the two markers (exo-E415G and plasmepsin 2–3) as covariates, and estimated adjusted hazard ratios (aHRs) and adjusted survival curves.

Copy number variations

We subsequently tested the association of copy number variations (CNVs) with piperaquine IC50s using the same method and parameters just described. We tested 43 CNVs present in at least five samples; genotypes represented the presence or absence (encoded as 0 or 1) of the CNV. In this analysis, mdr1 amplification was not included as a covariate in the linear mixed model. To call CNVs across the genome, we modified a procedure used previously.

, 

Briefly, we first divided the genome into 300 bp non-overlapping bins and calculated for each sample the number of reads whose alignment started within each bin. We then normalised these binned read counts by dividing by the median read count of the core regions of chromosome 9. We excluded bins where guanine–cytosine (GC) content was less than 20% due to coverage bias in most samples. Copy number state for each bin was predicted in each sample by fitting a Gaussian hidden Markov model to the normalised coverage data.

After determining the most likely state for each window (Viterbi algorithm), in each sample independently we merged adjacent windows having the same state prediction and allowed a maximum gap of 5000 bp between them. We then further removed windows less than 700 bp or not supported by face-away reads, and data for 43 samples with excessive variation in read coverage. This led to the identification of a median six CNVs per sample. Finally, we merged all the CNV calls across the samples, by considering any partial overlap. In total, we identified 134 regions containing CNVs across the core genome; 43 of these CNVs were present in at least five samples and were tested for association with the phenotype, as described above. For the 133 samples in the clinical study used for the survival analysis, we manually inspected each sample individually and scanned them for the presence of reads spanning the breakpoints.

Role of the funding source

The funders had no role in study design, data collection, analysis, interpretation, or report writing. The corresponding authors had full access to all data in the study and final responsibility for the decision to submit for publication.

Results

We analysed 486 P falciparum clinical isolates collected between July 9, 2010, and Dec 31, 2013, in three Cambodian provinces where artemisinin resistance is entrenched (Pursat), emerging (Preah Vihear), and uncommon (Ratanakiri; appendix).

, 

, 

To monitor the evolution of piperaquine resistance, we measured piperaquine IC50s for 297 P falciparum clinical isolates obtained directly from patients (ex vivo, 275 [93%] of 297) or following short-term culture (in vitro, 22 [7%] of 297) in 2011, 2012, and 2013. In all three provinces, IC50s increased significantly over time, especially when comparing province-stratified data from 2011 and 2013 (Kruskal-Wallis test, all p≤10−3; appendix). Despite parasites having comparable IC50s at all sites in 2011 (medians 20·0 nmol/L [IQR 13·7–29·0], 19·3 nmol/L [15·1–26·2], and 19·6 nmol/L [11·9–33·9] for Ratanakiri, Preah Vihear, and Pursat, respectively), they were remarkably differentiated at all sites by 2013, showing about two, three, and four times increases in IC50s (medians 39·2 nmol/L [IQR 32·8–48·1], 66·2 nmol/L [49·9–83·0], and 81·1 nmol/L [61·3–113·1], respectively). These regional differences are consistent with the relative prevalences of artemisinin-resistant parasites and recrudescent parasitaemias following dihydroartemisinin–piperaquine therapy.

, 

To investigate the genetic basis of piperaquine resistance in vitro, we did a GWAS analysis of the 297 samples (appendix) for which we had both piperaquine IC50s and whole-genome sequences (see Methods). We tested 11 630 SNPs that were well covered in this set of samples, and where one of the two alleles was present in at least ten samples. IC50s were treated as a continuous dependent variable in a linear mixed-model algorithm, and we treated the province of sample origin, presence of mdr1 amplification, and presence of kelch23 mutations as covariates; to correct for the confounding effect of population structure, we treated genetic similarity across samples as a random effect (appendix).

GWAS analysis identified one major locus on chromosome 13, containing two non-synonymous SNPs, that were significantly associated with raised IC50s (p ≤ 8·6 × 10−7; table 1, figure 1A, appendix). The strongest signal (2·55 times increase in IC50 after adjustment, p = 2·7 × 10−9) is from a non-synonymous SNP (referred to as exo-E415G) producing a Glu415Gly substitution in a putative exonuclease. The second strongest signal (p=2·3 × 10−7) from a non-synonymous variant is an SNP (referred to as mcp-N252D) within the same locus, producing an Asn252Asp substitution in a putative mitochondrial carrier protein 1. A locus on chromosome 12 also reaches genome-wide significance, but only with a single synonymous SNP. To exclude the possibility of a batch effect, we repeated the GWAS analysis by including the year of sample collection as a covariate in the linear mixed model (appendix). This analysis confirmed the locus on chromosome 13 as the only one significantly associated with the phenotype. When also correcting for year, mcp-N252D is the only SNP to reach genome-wide significance (p=7·3 × 10−7), although exo-E415G is just below the threshold (p=1·3 × 10−6). These data indicate that the chromosome 13 locus (referred to as L13) shows multiple significant GWAS hits that are robust to multiple corrections.

Table 1Genome-wide single-nucleotide polymorphisms (SNPs) most strongly associated with piperaquine 50% inhibitory concentrations (IC50s)

Significant SNPs (Bonferroni-corrected p≤8.6 × 10−7) and suggestive SNPs (p≤10−4) associated with piperaquine IC50s, according to increasing p value, are shown. For each SNP, we list the locus name, chromosome number, nucleotide position, gene ID, gene description, whether the SNP is non-synonymous (N) or synonymous (S), encoded aminoacid alteration, and association p value.

• Open table in a new tab

Figure 1Manhattan plot showing the statistical significance of single-nucleotide polymorphism (SNP) associations in the genome-wide association study (A), piperaquine 50% inhibitory concentrations (B), and exo-E415G frequency distribution (C)

Show full caption

(A) Each point represents one of the 11 630 SNPs with minor allele frequency greater than 0·033 in the set of 297 Plasmodium falciparum clinical isolates, alternately coloured red and blue according to chromosome. Genomic location is shown on the x-axis. The p value for each SNP association, calculated using a linear mixed model, is shown on the y-axis; point size is proportional to significance level. Province of sample origin, status of kelch23 (mutant vs wild type), presence or absence of mdr1 amplification, and a genetic relatedness matrix were added as fixed effects to the analysis. Four SNPs reached the Bonferroni-corrected, genome-wide significance level (p≤8·6 × 10−7; above horizontal blue line). Loci containing these significant SNPs, plus those containing suggestive SNPs reaching the significance level (p≤10−4; above horizontal green line), are listed in table 1. Vertical red lines mark known drug resistance loci: dhfr, mdr1, crt, dhps, and kelch23 on chromosomes 4, 5, 7, 8, and 13, respectively. Dashed grey vertical lines are plotted every 500 kb from the beginning of each chromosome. Solid grey vertical lines mark telomeric, sub-telomeric, and internal hypervariable regions. (B) Each point represents the piperaquine 50% inhibitory concentration (IC50) for a P falciparum clinical isolate carrying the wild-type A allele (WT) or mutant G allele (Mut) allele of exo-E415G (Pf3D7_13_v3:2504560). Bold and thin horizontal lines indicate the median and IQR of each distribution, respectively. Samples are divided into three coloured groups depending on their geographical origin. Coloured lines represent the least-squares linear regression of the phenotype on the two genotypes, calculated in each group separately. Shaded areas represent 95% CIs of the regression. (C) Coloured bars indicate the mutant allele frequencies in each of the three provinces over time (no samples were available from Preah Vihear in 2010). Error bars indicate 95% CIs of the estimation. Geographical coordinates are shown on the axes.

• View Large Image
• Figure Viewer
• Download Hi-res image
• Download (PPT)

To prioritise the multiple GWAS hits, we analysed their geographical spread, on the assumption that piperaquine resistance is emerging in Cambodia and neighbouring countries where dihydroartemisinin–piperaquine therapy is used (appendix). Taken together, the results of the GWAS and global allele frequency survey identify exo-E415G as the best candidate SNP marker of piperaquine resistance.

To gain a better picture of its frequency distribution in Cambodia, we used Sanger sequencing of PCR-amplified DNA fragments to genotype this SNP (ie, exo-E415G) in 168 additional samples for which whole-genome sequence data were not available, from the same initial cohort of 241 patients.

As expected, given the low p value of this SNP in the GWAS, the phenotype distribution differs significantly between the wild-type and the mutant alleles. When stratified by province and year, the association seems to be particularly strong in Pursat and Preah Vihear (figure 1B, appendix), and from 2012 onwards (appendix). Consistent with these results, the frequency of exo-E415G appears to increase in Pursat and Preah Vihear over time (figure 1C, table 2). In summary, exo-E415G is the best predictor of raised IC50s in the current dataset, although this finding alone is insufficient to infer a causal role for it in mediating piperaquine resistance.

Table 2Joint distribution of the two piperaquine resistance markers in Cambodia

The number of samples that are wild-type (WT) or mutant (Mut) for the exonuclease E415G mutation (exo-E415G) and have one (CN1) or multiple (CN2+) copies of plasmepsin 2–3, according to province and year of collection, are shown. Samples where data for one of the two markers were not available or the marker was present in a mixed infection were excluded.

• Open table in a new tab

Since gene CNVs have also been associated with drug resistance in P falciparum,

we investigated the association of piperaquine IC50s with 43 genomic regions exhibiting CNVs in at least five samples. We did a GWAS using the same method we used in the SNP GWAS, except that genotypes were the presence or absence of the CNV in each sample. We found that two CNVs strongly associated with IC50s (p≤2·3 × 10−4; figure 2A, appendix). As expected,

mdr1 amplification was associated (p=1·4 × 10−6) with low IC50s. However, the most significant association (p=7·6 × 10−7) was found for a novel amplification on chromosome 14, encompassing two of the ten plasmepsin genes in the P falciparum genome—plasmepsin 2 and plasmepsin 3—that encode aspartic proteases involved in the haemoglobin-to-haemozoin conversion pathway targeted by quinolines.

Figure 2Manhattan plot showing the statistical significance of copy number variation associations in the genome-wide association study (A) and piperaquine 50% inhibitory concentrations according to mdr1 and plasmepsin 2–3 copy number (B, C)

Show full caption

(A) Each point represents one of the 43 copy number variations (CNVs) present in at least five samples, alternately coloured red and blue according to chromosome. Genomic location is shown on the x-axis. The p value for each CNV’s association, calculated using a linear mixed model, is shown on the y-axis; point size is proportional to significance level. The province of sample origin, status of kelch23 (mutant vs wild type), and a genetic relatedness matrix were added as fixed effects to the analysis. Two CNVs reached the Bonferroni-corrected, genome-wide significance level (p≤2·3×10−4; above the horizontal blue line), one including plasmepsin 2 and plasmepsin 3, and one including mdr1. All 43 CNVs are marked by black lines at the top and are listed in the appendix. Dashed grey vertical lines are plotted every 500 kb from the beginning of each chromosome. Solid grey vertical lines mark telomeric, sub-telomeric, and internal hypervariable regions. (B, C) Each point represents the piperaquine 50% inhibitory concentrations (IC50s) for a Plasmodium falciparum clinical isolate carrying wild-type (WT) or amplified (Amplification) mdr1 (B) or plasmepsin 2–3 (C) genes. Bold and thin horizontal lines indicate the median and IQR of each distribution, respectively. Filled circles identify samples also carrying exo-E415G.

• View Large Image
• Figure Viewer
• Download Hi-res image
• Download (PPT)

A detailed analysis of sequencing reads aligned within the plasmepsin 2–3 locus revealed that the amplification boundaries were identical in all samples. The putative breakpoints lie near the 3′ end of both plasmepsin 1 and plasmepsin 3, so that each amplification creates an intact extra copy of plasmepsin 2 together with a new chimeric version of plasmepsin 3, with its 3′ end replaced by the 3′ end of plasmepsin 1 (appendix). Due to the degree of homology between plasmepsin 1 and plasmepsin 3, the aminoacid sequence of the chimeric plasmepsin 3 protein is identical to that of plasmepsin 3, except that an asparagine residue is replaced by two consecutive lysines.

The prevalence of plasmepsin 2–3 amplification in our cohort shows a steady increase over time in Pursat and Preah Vihear compared with Ratanakiri (appendix), and reflects our observations of increasing piperaquine IC50s and rising prevalence of exo-E415G in these two provinces (table 2, appendix). Despite being on different chromosomes, exo-E415G and plasmepsin 2–3 amplification are in significant linkage disequilibrium (r2=0·56, empirical p=1·6 × 10−5; appendix), and have a very similar allele frequency distribution in these Cambodian data. In particular, of the 462 samples where both markers were reliably typed, 72% (n=334) have neither marker, 19% (n=86) have both markers, and only 9% (n=42) have one of the two markers. Although the co-occurrence of the two markers in the population is interesting and surprising, the few samples where the two markers are found separately makes any conclusion regarding their relative effect on the phenotype difficult to support statistically (appendix).

To further investigate whether the exo-E415G and plasmepsin 2–3 amplification markers segregate with a newly described piperaquine resistance phenotype, we subjected a subset of 12 clinical isolates in triplicate to the in-vitro piperaquine survival assay.

26

• Duru V
• Khim N
• Leang R
• et al.

Plasmodium falciparum dihydroartemisinin-piperaquine failures in Cambodia are associated with mutant K13 parasites presenting high survival rates in novel piperaquine in vitro assays: retrospective and prospective investigations.

• Crossref
• PubMed
• Scopus (104)
• Google Scholar

In this assay, the survival of early ring-stage parasites following exposure to piperaquine for 48 h is assessed relative to non-exposed parasites tested in parallel. For these experiments, we selected two groups of parasites (appendix). In the first group, parasites did not recrudesce after dihydroartemisinin–piperaquine treatment, had low piperaquine IC50s ex vivo and after cultivation in vitro, were wild type for kelch23 and exonuclease, and carried single copies of mdr1 and plasmepsins 2–3. In the second group, parasites recrudesced after dihydroartemisinin–piperaquine treatment, had high piperaquine IC50s ex vivo and after cultivation in vitro, carried single copies of mdr1, and harboured the following mutations: kelch23-C580Y, exo-E415G, and plasmepsin 2–3 amplification. This second group of parasites showed significantly higher piperaquine survival rates (median 61·6 nmol/L [IQR 56·8–67·0]; n=6) than the first (2·4 nmol/L [1·6–2·9]; n=6; p=0·002). The relationship between raised piperaquine survival rate and the presence of either exo-E415G or plasmepsin 2–3 amplification were fully concordant.

Given that exo-E415G and plasmepsin 2–3 amplification are associated with raised piperaquine IC50s and survival rates, which, in turn, have been associated with parasite recrudescence,

, 

26

• Duru V
• Khim N
• Leang R
• et al.

Plasmodium falciparum dihydroartemisinin-piperaquine failures in Cambodia are associated with mutant K13 parasites presenting high survival rates in novel piperaquine in vitro assays: retrospective and prospective investigations.

• Crossref
• PubMed
• Scopus (104)
• Google Scholar

we directly tested their association with dihydroartemisinin–piperaquine failures. Of the 241 samples with clinical outcome data, we analysed a subset of 133 samples for which we had complete genetic and clinical information (appendix). Both exo-E415G and plasmepsin 2–3 amplification mutants showed a highly significant enrichment in recrudescent samples (p=1·6 × 10−8 and 1·8 × 10−11, respectively, Fisher’s exact test; appendix), with aHRs of recrudescence of 13·4 (95% CI 5·3–33·5; p=3·0 × 10−8) and 16·7 (5·8–48·1; p=1·7 × 10−7), respectively. Furthermore, the aHR for plasmepsin 2–3 amplification is still significant (5·2 [95% CI 1·5–17·7]; p=8·6 × 10−3) when only kelch23 mutants were considered, suggesting that whereas artemisinin resistance is certainly a risk factor, plasmepsin 2–3 amplification potentially has an additional, independent effect on piperaquine resistance and, ultimately, dihydroartemisinin–piperaquine treatment failure. To explicitly clarify the effect of the two markers on treatment success, we analysed the survival rate in samples with or without the two markers. We found that the adjusted survival rates at 63 days in the presence or absence of exo-E415G were 0·38 (95% CI 0·25–0·51) and 0·93 (0·85–0·97; figure 3A), respectively. Samples with or without plasmepsin 2–3 amplification had survival rates of 0·41 (95% CI 0·28–0·53) and 0·95 (0·87–0·98; figure 3B), respectively. Considering that non-recrudescent infections are not necessarily indicators of drug sensitivity, these values are likely underestimated. Together, these data identify exo-E415G and plasmepsin 2–3 amplification as strongly predictive markers of dihydroartemisinin–piperaquine treatment failures in Cambodia.

Figure 3Treatment success rate stratified by exo-E415G (A) and plasmepsin 2–3 copy number (B)

Show full caption

Adjusted survival curves showing the proportion of recrudescence cases following dihydroartemisinin–piperaquine treatment, according to the presence of the exo-E415G single-nucleotide polymorphism (A) and plasmepsin 2–3 amplification (B).

• View Large Image
• Figure Viewer
• Download Hi-res image
• Download (PPT)

Discussion

In GWAS analyses of piperaquine IC50 phenotypes, we identified two genetic markers of piperaquine resistance in vitro and of dihydroartemisinin–piperaquine treatment failures in patients: a non-synonymous SNP encoding a Glu415Gly mutation in a putative exonuclease (exo-E415G), and amplification of the plasmepsin 2 and plasmepsin 3 genes (plasmepsin 2–3 amplification). The prevalence of these two markers has increased substantially in recent years in Pursat and Preah Vihear, where artemisinin resistance is prevalent and where dihydroartemisinin–piperaquine has been the front-line artemisinin combination therapy since 2008 and 2010, respectively. In a global dataset of SNP allele frequencies, the exo-E415G allele is observed only in eastern mainland southeast Asia, where dihydroartemisinin–piperaquine is used, and completely absent where this artemisinin combination therapy has not been used.

Of the associated variations found in this study, plasmepsin 2–3 amplification is the strongest candidate for mediating piperaquine resistance because of the role of plasmepsins in haemozoin synthesis pathways in the parasite food vacuole. Since piperaquine, like chloroquine, is believed to inhibit the conversion of toxic haem moieties to non-toxic haemozoin crystals during haemoglobin digestion, it is possible that piperaquine targets plasmepsin 2, plasmepsin 3, or both, and that overproducton of these plasmepsins counteracts the drug’s action. The SNP markers, exo-E415G and mcp-N252D, do not lend themselves to an equally simple explanation and, until functional studies are done, it will not be possible to determine whether either is directly involved in mediating piperaquine resistance, whether their role is compensatory for lost fitness in piperaquine-resistant mutants, or whether they are associated with piperaquine resistance simply because of their strong linkage to plasmepsin 2–3 amplification or some other functional mutation in the Cambodian parasite populations.

Assigning a role to the single copy variant of mdr1, which is associated with lower sensitivity to piperaquine, is equally problematic. Although it is possible that mdr1 plays a functional part in piperaquine resistance (eg, by restricting the amount of drug that enters the parasite’s food vacuole), recent changes in Cambodia’s antimalarial drug policy have caused a decline in mefloquine pressure, which might have promoted the loss of mdr1 amplifications in the parasite populations. In view of the pronounced population structure following the emergence of artemisinin resistance in Cambodia,

it is possible that the association with a single copy of mdr1 is the result of piperaquine resistance having emerged in specific artemisinin-resistant, mefloquine-sensitive populations.

Taken together, these findings identify plasmepsin 2–3 amplification as the most likely candidate to be a causal variant of piperaquine resistance. In the samples studied here, it is strongly linked to exo-E415G, which is on a different chromosome and might represent some other functional component of the resistance phenotype. In practice, this means that the exo-E415G SNP can be presently used as a predictive marker of piperaquine resistance in Cambodia; however, this might not be the case elsewhere. Therefore, plasmepsin 2–3 amplification should be monitored in areas where dihydroartemisinin–piperaquine is used, although it is somewhat more laborious to genotype. Despite unresolved questions about causal mechanism, we now have tools to identify areas where treatment failures with dihydroartemisinin–piperaquine are likely to occur, and thereby to empower national malaria control programmes to make informed decisions about whether to switch to alternative artemisinin combination therapies for first-line treatment of P falciparum malaria.

Contributors

RA, PL, OM, DPK, and RMF designed the study. RA, PL, CA, DD, ATN, SSr, SSu, and ED collected data. RA, PL, OM, RDP, JA-G, DJ, JS, DPK, and RMF analysed and interpreted data. RA, PL, DPK, and RMF wrote the manuscript.

Declaration of interests

We declare no competing interests.

Supplementary Material

• View Large Image
• Download Hi-res image
• Download .PPT

Why is the heterozygous sickle cell genotype common only in malarial regions?

What is another name for a reproductive population?

Which of the following conditions contributes to evolution?

Which of the evolutionary forces is most likely to decrease variation between populations?