Artemisinin resistance without pfkelch13 mutations in Plasmodium falciparum isolates from Cambodia

Culture-adaptation and maintenance of TRAC parasites

All parasite samples were collected under protocols approved by ethical review boards in Cambodia, at Oxford University and at the Harvard T. H. Chan School of Public Health. Culture-adaptation of parasites was accomplished by thawing cryopreserved material containing infected red blood cells (iRBCs) that had been mixed with glycerolyte. Parasites were maintained in fresh human blood (O+) and Hepes buffered RPMI media containing 12.5% AB+ human serum (heat inactivated and pooled). Cultures were placed in modular incubators and gassed with 1% O₂/5% CO₂/balance N₂ gas and incubated with rotation (50 rpm) in a 37 °C incubator (Additional file 1).

Sample extraction

Genetic material were extracted from filter papers (Whatman) and culture-adapted material using Promega DNA IQ Casework Pro Kit for Maxwell 16 (Promega Corp., Madison, WI, USA) and Qiagen (QIAmp DNA Blood Mini Kit) commercial kits, respectively, according to manufacturer instructions.

Genotyping

Development of a high-resolution melt (HRM) assay to screen populations for mutations around amino acid position 580 in the pfkelch13 locus. The forward primer (5′-GGCACCTTTGAATACCC-3′), reverse primer (5′-CATTAGTTCCACCAATGACA-3′), and unlabeled, blocked probe (5′-AGCTATGTGTATTGCTTTTGAT-block-3′) were amplified asymmetrically at 0.5, 0.1, and 0.4 μM, respectively with 1 ng template DNA. After an initial 2-min hold at 95 °C, 5 or 10 μL reaction mixtures with 2.5 × HRM master mix (BioFire Defense, Salt Lake City, UT. USA) were PCR amplified for 55 cycles: 95 °C for 30 s, 66 °C for 30 s, and 74 °C for 30 s, followed by a pre-melt step of 95 and 28 °C for 30 s each. Products were melted from 45 to 90 °C on a BioFire Defense LightScanner-384 or -32 and analysed using the manufacturer software. Two plasmid controls containing the wild-type and mutant alleles were included as standards for every HRM run. Molecular barcoding was performed as described [24, 25] to identify monogenomic samples with unique parasite genotypes.

Whole genome sequences

This publication uses sequencing data generated by the Pf3k project [26]. The variant call files generated from this project were used to identify SNPs in each of the samples used in this study.

Pfkelch13 PCR sequencing strategy

The propeller domain of pfkelch13 was PCR amplified using Phusion HF DNA Polymerase kit and primers 3F and 1R′ in all 68 culture adapted parasites. DNA from KH001_024 was also amplified with primers 4F and 3R′. An aliquot of PCR product was resolved by gel electrophoresis to check for specificity and yield and the remaining product was purified using DNA Clean & Concentrator, ZymoResearch and sequenced using the same primers used to amplify the product by Genewiz. The full ORF of pfkelch13 was PCR amplified using primers 1F and 1R from 3D7 and individual culture-adapted parasites. The resulting PCR product of ~2.2 kb was purified by gel extraction (QIAquick Gel Extraction Kit, Qiagen) and sequenced at Genewiz using primers 1F, 2F, 3F, 1R, 2R and 3R. Primer sequences are as follows: 1F: 5′-ATGGAAGGAGAAAAAGTAAAAACAAAAGCAAATAG-3′; 2F: 5′-GGTAGGTGATTTAAGAATTACATTTATTAATTGGT-3′; 3F: 5′-CATTCCCATTAGTATTTTGTATAGGTG-3′; 4F: 5′-GTAGAGGTGGCACCTTTGAATACCCCTAGATCATC-3′ 1R: 5′-TTATATATTTGCTATTAAAACGGAGTGACCAAATCTG-3′; 1R′: 5′-TTA TAT ATT TGC TAT TAA AAC GGA GTG-3′; 2R: 5′-AGCCTTATAATCATAGTTATTACCACCAAAAACG-3′; 3R: 5′-TGTTGGTATTCATAATTGATGGAGAATTC-3′; 3R′: 5′-ATAAAATGTGCATGAAAATAAATATTAAAG-3′.

In vitro 72-h drug susceptibility by SYBR green staining

Drug susceptibility was measured using the SYBR Green I method as previously described [27]. Briefly, tightly synchronized 0–6 h rings were grown for 72 h in the presence of different concentrations of drugs in 384-well plates at 1% haematocrit and 1% starting parasitaemia; and, growth at 72 h was measured by SYBR Green staining of parasite DNA. Except for PPQ and KH001_053 where a 24-point dilution series was used, a 12-point dilution series of each drug was carried out in triplicate and repeated with at least three biological replicates. DMSO stocks of drugs were dispensed by a HP D300 Digital Dispenser (Hewlett Packard Palo Alto, CA, USA) except for the CQ and PPQ stocks that were prepared in water and dispensed with a Velocity 11 Robot (Bravo). Relative fluorescence units (RFU) was measured at an excitation of 494 nm and emission of 530 nm on a SpectraMax M5 (Molecular Devices Sunnyvale, CA, USA) and analysed using GraphPad Prism version 5 (GraphPad Software La Jolla, CA. USA). EC₅₀ values were determined with the curve-fitting algorithm log(inhibitor) vs response—variable slope, except for PPQ and KH001_053. Due to the bimodal dose response of KH001_053 to PPQ, curve fitting didn’t give an accurate EC₅₀ value. The reported PPQ EC₅₀s for KH001_053 are estimates using biphasic curve fitting. Spearman correlation analysis was performed to assess the relationship between the anti-malarial EC₅₀ values and in vivo clearance half-life, ring survival assay value or pfmdr1 copy number. p values <0.05 were considered significant.

Copy number variation assays

To determine copy numbers for pfmdr1, plasmepsin II and the 63 kb amplicon genes (PF3D7_0520100, PF3D7_0520500, PF3D7_0520600, PF3D7_0520900 and PF3D7_0521000), real time quantitative PCR was performed on genomic DNA (extracted with QIAmp Blood Mini Kit, Qiagen) as previously described [28] with the following modifications: Amplification reactions were done in MicroAmp 384-well plates in 10 μL SYBR Green master mix (Applied Biosystems), 150 nM of each forward and reverse primer and 0.4 ng template. Forty cycles were performed in the Applied Biosystems ViiA™ 7 Real-time PCR system (Life Technologies). Forward and reverse primers used were as previously described to amplify the following loci: pfmdr1 (PF3D7_0523000) [29], the 63 kb region on chromosome 5 (PF3D7_0520100, PF3D7_0520500, PF3D7_0520600, PF3D7_0520900 and PF3D7_0521000 [20]) and plasmepsin II (PF3D7_140800) [22]. For the endogenous controls, β-tubulin forward and reverse primers [28] were used for pfmdr1, PF3D7_0520100 and PF3D7_0520900 while pfldh forward and reverse primers [30] were used for PF3D7_0520500, PF3D7_0520600 and PF3D7_0521000. Target primers used were validated to have the same PCR efficiencies as their endogenous control primers; and, average copy number values were calculated for each gene using data from three independent experiments.

Sequencing the pfcrt locus

The entire pfcrt locus was sequenced as previously described [20] with some modifications. Briefly, total RNA was extracted using RNeasy kit (Qiagen) and used to generate cDNA using Superscript III (Invitrogen). The resulting cDNA was then used as template for PCR amplification of pfcrt [20], followed by Sanger sequencing (GENEWIZ) [20]. Sequence data analysis was performed using MacVector.

Ring survival assay (RSA_0–3h)

The RSA_0–3h was performed as described previously [6]. Essentially, parasites were sorbitol synchronized twice at 40-h intervals, synchronous 40–44 h segmented schizonts were incubated for 15 min at 37 °C in serum-free media supplemented with heparin to disrupt agglutinated erythrocytes and late stages were purified with 35/65% discontinuous Percoll gradient. The segmented schizonts were washed and cultured with fresh RBCs for 3 h, after which late stages were removed by sorbitol treatment. Cultures with 0–3 h rings were adjusted to 2% haematocrit and 1% parasitaemia and seeded into a 24-well plate with 1 ml complete media per well. To these wells, either DHA at 700 nM or 0.1% DMSO were added immediately and incubated for 6 h at 37 °C, washed and incubated in drug free media. At 72 h from seeding, thin blood smears were made from control and treated wells and survival rates were measured microscopically by counting the proportion of next generation viable rings with normal morphology. Survival rates were expressed as ratios of viable parasitaemias in DHA-exposed and DMSO-treated controls. Parasites were counted from 10,000 RBC, and two separate individuals served as independent slide readers.

Characteristics of original and culture-adapted parasite isolates

A total of 157 cryopreserved TRAC study samples collected in 2011 from Pursat or Pailin, Cambodia where ART resistance is observed were obtained. These isolates were selected before the discovery of the pfkelch13 marker and were designed to comprise western Cambodian samples, with two-thirds of the samples coming from the upper segment of the clearance half-life distribution and one-third being isolates with the shortest clearance times. 68 parasites were culture-adapted for further investigation, without knowledge of their in vivo clearance data. A summary of the overall population and the sub-set selected for culture-adaptation and further characterization is provided in Additional files 2 and 3. Collectively, 63% of both the original and the adapted parasites were from Pursat and the remaining 37% were from Pailin (Additional file 3). Similarly, 67% of both parasite sets (original and adapted) had a delayed clearance of ≥5 h and 33% had a clearance half-life of <5 h; for one isolate there were insufficient data to determine an in vivo half-life. All culture-adapted samples were confirmed to harbour monogenomic infections by molecular barcode analysis [25], only those parasites with unique barcode signatures [25] were initially pursued to maximize genetic diversity and exclude highly similar parasites among the analysed isolates.

Pfkelch13 propeller mutations are associated with parasite clearance half-life ≥5 h

Next, parasites were genotyped for the pfkelch13 locus previously associated with delayed parasite clearance (defined as >5 h) and increased ring-stage survival (defined as RSA_0–3h ≥1%) [7]. An HRM genotyping assay [24, 25] for the common non-synonymous mutation conferring the C580Y change in PfKelch13 was developed and validated, among the 157 original samples with an unambiguous genotyping call, 71% (108/152) had the C580Y mutation (Additional file 2). These data, confirmed a significant positive association between the C580Y non-synonymous mutation and in vivo parasite clearance half-life values for these parasites [26]. This association was more pronounced among Pursat parasites (unpaired Student’s t test with Welch’s correction, p < 0.0001) as compared to parasites from Pailin (p = 0.0322, Additional file 4) for this sample set.

A combination of available whole genome sequence (WGS) data from the original samples [26, 31, 32], and pfkelch13-specific PCR-based sequencing identified both pfkelch13 mutations other than C580Y in this population and confirmed persistence of all pfkelch13 mutations in the culture adapted parasites further analysed. Overall there was good concordance between the molecular barcode derived from filter paper and cultured parasite samples, and between the pfkelch13 mutations identified in WGS or PCR re-sequencing data. Six different mutations were present in the pfkelch13 propeller domain: Y493H, R539T, I543T, C580Y, D584V, and H719N, and two pfkelch13 mutations outside the propeller region: H136N, E270K. The two parasites with the H136N change also contained C580Y, but E270K was seen by itself, in a single isolate (Additional file 2). All propeller domain mutations were positively associated with in vivo clearance values (Fig. 1a), including parasites with the double mutation H136N/C580Y genotype (Additional file 2). Additionally, two asparagine (Asn) insertions at codon 142 were identified in all isolates tested except those with the Y493H change and one parasite with the D584V genotype (Additional file 2). The location of this Asn insertion was distinct from that previously noted among Asian parasites [33]. No association between this Asn insertion and either parasite clearance half-life or RSA_0–3h phenotype was detected.

Pfkelch13 propeller mutations are positively associated with increased ring-stage survival by in vitro RSA_0–3h

There was a positive association between increased RSA_0–3h survival phenotypes and pfkelch13 mutations among 36 culture adapted parasites (Fig. 1b) and parasites with the D584V change previously associated with increased in vivo parasite clearance time [5] demonstrated an increased ring-stage survival by RSA_0–3h.

In vitro RSA_0–3h phenotypes were evaluated as a measure of ART resistance for 36 adapted parasites exhibiting a range of in vivo clearance phenotypes and harbouring a variety of pfkelch13 mutations. Consistent with previous reports [6], the distribution of RSA_0–3h phenotypes spanned two log₁₀ ranges for all parasites analysed. Overall, there was a significant difference between wild-type (RSA_0–3h = 0.1–6.0%, n = 7) parasites and those with pfkelch13 mutations (RSA_0–3h = 2.1–38.6%, n = 29) in terms of RSA_0–3h parasite survival rate (Student’s t test p = 0.0018). RSA_0–3h survival rates also differed between parasites with individual pfkelch13 mutations; whereby parasites with either a C580Y or R539T mutation exhibited the highest survival rates (Fig. 1b). One isolate with a PfKelch13 E270K mutation outside the propeller domain exhibited a clearance of >4 h, yet had an RSA_0–3h value of 0.4%. Among the parasites with RSA_0–3h values >1%, 19% (7/36) had in vivo clearance values between 4 and 5 h, which suggests that for this specific population a clearance value of 4 h or greater is more consistent with ART resistance based upon an RSA_0–3h phenotype (Additional file 5).

All pfkelch13 propeller polymorphisms found among these parasites were previously identified [7, 33]. However, a novel positive association was demonstrated between the D584V mutation and in vitro ring survival values (8.0 h clearance and 8.5% RSA survival for D584V), not previously tested for in vitro ART resistance (Fig. 1b). These data confirm that all pfkelch13 propeller mutations found in these TRAC isolates are associated with an RSA_0–3h phenotype ≥1%, thus considered in vitro ART resistant.

Identification of discordant parasites that exhibit increased ring-stage survival but lack Pfkelch13 mutations

ART resistance in vitro has been suggested to correlate with an RSA_0–3h survival value of ≥1% [6], and all isolates in this study with a mutation in pfkelch13 had an RSA_0–3h ≥1%. All wild-type pfkelch13 parasites had in vivo clearance half-lives below 5 h (Additional file 2); however, there was a range of in vitro RSA_0–3h survival phenotypes from 0.1 to 6%, including four isolates with RSA_0–3h of approximately ≥1% (0.8, 1.1, 3, and 6%) (Fig. 1b). Parasites with an RSA_0–3h of ≥0.8% that lack pfkelch13 mutations are considered discordant. These data suggest changes outside pfkelch13 could also confer an increased ring-stage survival phenotype indicative of ART resistance in this parasite population.

Using WGS data previously generated [26] for 33 of the RSA_0–3h phenotyped parasites, the question was explored whether mutations outside pfkelch13 were evident in the four parasites with RSA_0–3h ≥0.8% values that lack pfkelch13 ORF mutations. This small data set was underpowered to identify novel mutations using a genome-wide strategy, but a candidate gene approach was undertaken to ask whether any known drug resistant loci or other previously identified secondary mutations could potentially explain the high RSA_0–3h survival values (Additional file 6). Variant positions were identified across 23 genes previously associated with ART resistance among the set of 33 RSA_0–3h phenotyped parasites [7, 15, 26], which included both sensitive and resistant parasites. Examination of both sensitive and resistant parasites would allow detection of mutations that may account for the differences in RSA_0–3h survival values and identification of mutations that confer high RSA_0–3h survival values independent of pfkelch13.

This analysis identified 37 SNPs (30 non-synonymous, four synonymous, two intronic) spread across 19 unique genes (Additional file 7). Eleven of these positions were invariant among these four isolates that lack pfkelch13 ORF mutations but exhibited an RSA_0–3h survival value ≥0.8%. These 37 SNP positions are found in genes known to be associated with drug resistance including dihydrofolate reductase (dhfr, PF3D7_0417200) [34] and dihydropteroate synthase (dhps, PF3D7_0810800) [35]. Another locus, phosphatidylinositol-4-phosphate 5-kinase (PIP5K, PF3D7_0110600) [36], is involved in cellular signaling pathways, including synthesis of PIP2, a substrate for PI3K to produce PIP3, which activates the AKT family of serine/threonine kinases, a pathway implicated in an ART resistance mechanism [15]. The putative NLI Interacting factor-like phosphatase (NIF4, PF3D7_1012700) [36–39] is a protein phosphatase 2C protein. Several members of the Plasmodium phosphatome contain kelch domains and have been shown to be key regulators of parasite development and differentiation [37, 40], but NIF4 lacks an intact DxDx(T/V) motif, and thus may be catalytically inactive [39]. Pfmdr2 (PF3D7_1447900) [41, 42] is implicated as a secondary locus associated with delayed parasite clearance [26]. No other changes associated with ART resistance reported in the literature [43, 44] were detected among these parasites (Additional file 6).

Detection of PPQ resistance phenotype

To further evaluate the relationship between RSA_0–3h phenotype and pfkelch13 mutations, a sub-set of eight isolates with a range of ring-stage survival phenotypes was selected. This sub-set included the two most sensitive (in green) and the two most resistant (discordant, in blue) parasites without pfkelch13 mutations, as well as four resistant parasites that contain pfkelch13 mutations (in red) (Fig. 1b). First standard in vitro drug tests [27] were used to evaluate anti-malarial drug responses among these parasites (Fig. 2). EC₅₀ values for ART and its derivatives did not consistently correlate with either the in vivo clearance or in vitro RSA_0–3h phenotypes, as previously reported [6, 7] (Additional file 8; Fig. 2). The parasite clearance half-life and RSA_0–3h survival phenotypes did not consistently associate with responses to other anti-malarial compounds tested (chloroquine (CQ), MFQ, lumefantrine (LUM), atovaquone (ATV), quinine (QN) and PPQ, Additional file 8). A significant positive correlation was observed between EC₅₀ values for ART (r = 0.881, p = 0.0072), artesunate (AS, r = 0.719, p = 0.0368), MFQ (r = 0.7619, p = 0.0368) and LUM (r = 0.857, p = 0.0107) and increased pfmdr1 copy number (Spearman correlation, Additional file 8 and Fig. 2), consistent with previous reports [30, 45, 46].

Interestingly, one isolate, KH001_053, with a slow clearance time (8.7 h) and high ring survival (30%) exhibited a bimodal response to PPQ that indicated parasite survival at high PPQ concentrations (Fig. 3a). Compared to other isolates, KH001_053 also demonstrated a high survival phenotype (50–70%) when grown at high concentrations of PPQ (125 and 250 nM) for 72 h [47] (Fig. 3b). This PPQ-resistant parasite also demonstrated hyper-susceptibility to MFQ and LUM (Fig. 2), explicable at least in part by a single copy of the pfmdr1 gene (Additional file 8).

KH001_053 was then tested for evidence of known molecular markers associated with PPQ resistance [20], including C101F mutation in pfcrt, an amplification of 63 kb region on chromosome 5, and amplification of plasmepsin II [21, 22]. All eight isolates were evaluated for SNPs in the pfcrt locus by PCR amplification and sequencing of cDNA across pfcrt. Eight non-synonymous mutations (M74I, N75E, K76T, A220S, Q271E, N326S, I356T, R371I) were found that were fixed in these eight isolates. No evidence of the previously reported C101F haplotype in the KH001_053 parasite was detected [20]. Since PPQ resistance has been associated with an increased copy number of a 63 kb amplicon on chromosome 5 (including PF3D7_0520100, PF3D7_0520500, PF3D7_0520600, PF3D7_0520900 and PF3D7_0521000, 825600–888300 bp) and a decreased number of pfmdr1 in vitro [20, 30] as well as an increased copy number of plasmepsin II (PF3D7_140800) locus in vivo [21, 22], copy number assays were performed on those loci. In concordance with to a previous study [30] no increased copy number on the 63 kb amplicon region on chromosome 5 was detected. However, KH001_053 did harbour two copies of plasmepsin II consistent with previous reports where the presence of an increased plasmepsin II and III copy number was predictive of recrudescence in DHA/PPQ-treated patients [21, 22]. These findings strengthen the hypothesis that increased plasmepsin II copy number may play a role in PPQ resistance.