Amplicon Deep Sequencing Reveals Multiple Genetic Events Lead to Treatment Failure with Atovaquone-Proguanil in Plasmodium falciparum

ABSTRACT Atovaquone-proguanil (AP) is used as treatment for uncomplicated malaria, and as a chemoprophylactic agent against Plasmodium falciparum. Imported malaria remains one of the top causes of fever in Canadian returning travelers. Twelve sequential whole-blood samples before and after AP treatment failure were obtained from a patient diagnosed with P. falciparum malaria upon their return from Uganda and Sudan. Ultradeep sequencing was performed on the cytb, dhfr, and dhps markers of treatment resistance before and during the episode of recrudescence. Haplotyping profiles were generated using three different approaches: msp2-3D7 agarose and capillary electrophoresis, and cpmp using amplicon deep sequencing (ADS). A complexity of infection (COI) analysis was conducted. De novo cytb Y268C mutants strains were observed during an episode of recrudescence 17 days and 16 h after the initial malaria diagnosis and AP treatment initiation. No Y268C mutant reads were observed in any of the samples prior to the recrudescence. SNPs in the dhfr and dhps genes were observed upon initial presentation. The haplotyping profiles suggest multiple clones mutating under AP selection pressure (COI > 3). Significant differences in COI were observed by capillary electrophoresis and ADS compared to the agarose gel results. ADS using cpmp revealed the lowest haplotype variation across the longitudinal analysis. Our findings highlight the value of ultra-deep sequencing methods in the understanding of P. falciparum haplotype infection dynamics. Longitudinal samples should be analyzed in genotyping studies to increase the analytical sensitivity.

no parasites were observed. Eighteen days after her initial admission, and without subsequent travel, the patient was feeling unwell and readmitted to the hospital. She was febrile (40.2°C), developed headache, nausea, and vomiting. She was diagnosed with P. falciparum malaria after a blood smear revealed 2.6% parasitemia. Shortly after she was initiated on doxycycline (100 mg) and quinine sulfate (600 mg) for a 7-day course. She was discharged 2 days after her last admission with no further symptoms.
Genomic markers of resistance to atovaquone-proguanil. A total of 622,605 reads were sequenced for the three resistance markers. A read depth average 6 standard error of the mean (SE) per gene across all samples was computed, with a mean and SE depth of 1,719 6 319.3 for cytb, 2,190 6 865.7 for dhfr, and 1,864 6 999.2 for dhps. SNP calling was performed, with SNPs called only if 100 reads mapped to a given position and at least 1% of reads contained a SNP. A total of 7 SNPs and 5 deletions were observed across the three genes of interest over the 12 clinical samples. Of these, four SNPs previously correlated with treatment resistance were observed, including the A803G (coding for Y268C) in cytb conferring resistance to atovaquone, c1790g (coding for A437G) in dhps conferring resistance to sulfadoxine, and the dhfr a343t (N51I) and g514t (S108N), known to result in resistance to cycloguanil (5) ( Table 1). The dhps and dhfr SNPs associated with treatment resistance were observed across all 12 clinical samples (Fig. S1). The cytb mutant, on the other hand, was undetected across the seven samples prior to the treatment failure but became the dominant genotype during recrudescence. We thus found no evidence, within the limits of detection of this FIG 1 Haplotyping and SNP calling workflow. Whole blood samples (n = 12) were collected across the patient care. P. falciparum diagnostic and parasitemia was assessed using a Giemsa-stained thin blood smear. DNA was extracted from whole-blood samples followed by PCR amplification and Illumina sequencing of three genes associated with treatment failure in P. falciparum malaria (cytb, dhfr, and dhps). Haplotyping was performed using three different approaches by amplifying the following gene markers: msp2-3D7, msp2-FC27, and cpmp. SNP calling was performed using GATK. Traditional haplotyping approaches (2% agarose gel and capillary electrophoresis) were performed on msp2-3D7 and msp2-FC27. Amplicon deep sequencing haplotyping was performed with cpmp using DADA2. This figure was created with BioRender.com. approach, that the variant was preselected and already present in the first malaria episode (Fig. 2). During the episode of recrudescence, three out of five samples had a few reads mapped to the cytb WT gene. A total of 3 WT reads (0.20%) were mapped to the 19d 5h sample, 6 WT reads were mapped to the 20d 4h sample, and 13 WT reads were mapped to the 20d 14h sample. Overall, the cytb sequencing profile indicates a fixation of the cytb Y268C at the time of treatment failure.
Interestingly, the dhps c2879g UTR SNP was found at 100% prevalence in four of the seven samples before recrudescence (0 days 0 h, 1 day 6 h, 1 day 14 h, and 2 day 4 h), and two of the five samples during recrudescence (18 day 10 h, 19 day 5 h, 19 day 17 h) (Table S1). This pattern of alternating mutant prevalence was observed before and during recrudescence. No sequence variants in cytb, dhfr, or dhps were detected in the P. falciparum 3D7 positive control.
Complexity of infection and haplotyping analysis. To determine whether the cytb Y268C mutant observed during the treatment arose due to a de novo mutation or was simply a minor clone undetected during the first days of infection, we performed a haplotyping analysis using gel electrophoresis, capillary electrophoresis, and ADS ( Fig. 3 and Table 2). This approach also allowed us to determine if the change in frequency of Y268C resulted in a soft or hard selection sweep.
The first method for haplotyping consisted of estimating the number and size of msp2-3D7 and msp2-FC27 haplotypes via agarose electrophoresis. All the samples amplified the msp2 gene (n = 12). A range of 3 to 4 msp2-3D7 haplotypes was observed prior to the recrudescence, and 3 haplotypes were observed after treatment failure. All haplotypes were observed before and during the episode of recrudescence. (Table S2 and Fig. S2). The msp2-3D7 fragment sizes ranged from 370 bp to 668 bp. The 668 bp was the most prevalent haplotype across all samples. No msp2-FC27 family-specific alleles were observed. None of the msp2-3D7 haplotypes observed by this method were consistently found across all 12 samples. Similar to the agarose electrophoresis, we estimated the number and size of msp2-3D7 and msp2-FC27 haplotypes using capillary electrophoresis. A range of 6 to 8 msp2-3D7 haplotypes was observed across the samples prior to the treatment failure, and a range of 4 to 7 haplotypes was observed during the episode of recrudescence. The msp2-3D7 fragment sizes ranged from 346 bp to 638 bp (Table S3). As seen with gel electrophoresis, no msp2-FC27 family-specific alleles were observed. All haplotypes were observed before and during the recrudescence with the exception of one unique haplotype in an individual sample after the treatment failure (19 days 17 h). Similar to the agarose electrophoresis, none of the haplotypes were consistently observed across all samples. To determine the ability of this assay to detect minor haplotypes, we generated contrived samples of two haplotypes mixed at varying frequencies. We were able to detect the minor haplotype using .10 copies/ml of DNA for samples mixed 50%-50%, 75%-25%, and 90%-10%, and using .100 copies/ml of DNA for those mixed 95%-5% and 99%-1%.  The last approach consisted of the haplotype analysis of the cpmp gene using ADS (mean of 652,235 6 29,417) ( Fig. S4). Here, we observed a range of 4 to 6 cpmp haplotypes across the samples prior to the treatment failure, and a range of 5 to 6 haplotypes was observed during the recrudescence. All cpmp haplotypes were observed before and after recrudescence (Table S4). Four of the six haplotypes (66.66%) were found across all 12 samples. Overall, two major haplotypes were dominant across all samples, without any significant changes in their proportion, while the remaining four haplotypes were found at a lower prevalence (Fig. 4a). A range of 3 to 5 polymorphisms and a gap were observed between the haplotypes reported and the P. falciparum 3D7 cpmp reference gene (PlasmoDB: PF3D7_0104100) (Fig. 4b). Using single and mixed haplotype contrived controls, we observed a limit of detection of 0.1 copies/ml. No artifactual variants were observed in the control samples using the specified thresholds.
The three haplotyping profiles generated did not identify a minor clone that could have been preselected prior to the treatment failure. All haplotypes generated were identified before and during the episode of recrudescence (apart from one haplotype on 19d 17h by capillary electrophoresis), suggesting maintenance of genetic diversity during the episode of recrudescence, leading to a selective soft sweep.
None of the haplotypes observed in the clinical samples by agarose gel or capillary electrophoresis, nor the haplotypes estimated by the ADS approach were detected in the positive controls. This indicates msp2 and cpmp are good markers for haplotype discrimination.
To investigate the performance of haplotyping approaches, we compared the COI by methods before and during the episode of recrudescence (Fig. S3). When comparing the COI values before the episode of recrudescence, we observed significantly different results between the agarose electrophoresis versus the capillary electrophoresis (Dunn's adjusted P value = 0.01), the agarose electrophoresis versus the cpmp ADS (Dunn's adjusted P value = 0.002), and no significant differences between the capillary electrophoresis versus the cpmp ADS (Dunn's adjusted P value . 0.99). The COI values were also compared during the episode of recrudescence, with significantly different values between the agarose electrophoresis versus the capillary electrophoresis (Dunn's adjusted P value = 0.04), agarose electrophoresis versus the cpmp ADS (Dunn's adjusted P value = 0.01), and no significant differences between the capillary electrophoresis versus the cpmp ADS (Dunn's adjusted P value . 0.99). Thus, the same COI trends across the methods were observed before and during recrudescence.

DISCUSSION
To our knowledge, this is the first described observation of an early clinical failure of AP treatment caused by multiple cytb Y268C de novo mutations with distinct genetic backgrounds (COI . 3). In this study, three haplotyping approaches and a multilocus ADS approach were used to determine the genetic longitudinal profiles of 12 samples before and during an episode of recrudescence. Data presented in this work support prior reports of de novo Y268 cytb mutations in returning travelers and in malariaendemic populations while being treated with AP (25)(26)(27)(28)(29).
The three methodologies agree on a similar haplotyping pattern across all samples, suggesting no clearance of any of the clones during the first 3 days of treatment. Additionally, no Y268C reads were observed during the seven samples prior to recrudescence, eliminating the possibility of an existing preselected mutant. The maintenance of genetic diversity by the haplotyping results suggests a soft selective sweep occurred within a window of 17 days and 16 h of at least four independent clones. Current studies highlighting AP treatment failure have previously identified one to two haplotypes independently acquiring the same de novo mutation (4,14,(25)(26)(27)30). No studies have shown evidence of multiclonal (COI . 3) de novo mutations at the 268th position of the cytb gene.
Despite having a consistent haplotype profile before and during the episode of recrudescence, the capillary and gel electrophoresis failed to identify a consistent haplotype across all 12 samples. ADS, on the other hand, was able to capture 4 out of 6 haplotypes consistently. Despite recording nonsignificant differences in COI by capillary electrophoresis and cpmp ADS, we observed a higher variability in haplotype calls from capillary electrophoresis. In addition, both the capillary electrophoresis and the ADS methodologies estimated a higher COI, suggesting a higher sensitivity as previously discussed (3,7). Similar to our findings, patterns of unique haplotypes have been previously recorded in longitudinal studies using msp1 and msp2 as haplotyping markers (28). This suggests the COI in the parasite population could be underestimated when using individual samples. This is of paramount importance for therapeutic efficacy trials and genomic surveillance.
Current reports highlighting AP treatment failure are limited by using Sanger sequencing, agarose gel haplotyping of length-polymorphic markers, or a small sample size (4,9,14,25,(29)(30)(31)(32). This could underestimate the haplotype profile and/or the type of selective sweep in de novo mutations. Indeed, previous studies have highlighted the low mutant detectability of Sanger sequencing in contrived samples, suggesting a limit of detection as low as 70% (17). Additionally, length-polymorphic markers estimate haplotypes by looking at the size of alleles between treatment and follow-up samples. These methods have several limitations, such as the preferential amplification of shorter alleles, higher limits of detection (LOD), decreased specificity, and formation of spurious bands (21,22,33,34). Despite this, recent studies still use length polymorphic markers (msp1, msp2, glurp) and microsatellites for haplotyping (35,36). This could lead to a bias in their findings, especially if no posterior Bayesian correction for recrudescence or reinfection was applied or if a unique sample was taken at a given time (28,37,38).
Classical models of treatment failure consist of a beneficial allele being part of a population and expanding under directional drug selection (30). Remarkably, no cytb Y268C allele was detected in any of the samples prior to the recrudescence within the limit of detection of ADS. Because no parasitic recombination occurs during the erythrocytic stage, it is likely that the cytb Y268C arose independently in multiple clones with distinct genetic backgrounds while maintaining their genetic diversity after the treatment failure. An alternative explanation is that Y268C selection occurred due to a long period of subtherapeutic AP exposure that could have been competing with other genotypes under low AP pressure. It is possible but less likely that low levels of the mutant were not detected in ADS due to a skew in the read distribution by the 30 cycles of amplification during the library preparation. Furthermore, all the clinical samples from this study carried two dhfr mutations (S108N, N51I) that render resistance to cycloguanil (5), which have been reported in southwest Uganda (39). Amplicon deep sequencing protocols should be used to determine the nature of the selective sweep when evaluating the genetic emergence of resistance in Plasmodium, in order to understand the implications behind the propensity of AP resistance development.
Surprisingly, the time to recrudescence was shorter than the estimated intervals for AP therapeutic failure previously reported by Sutherland (7,40). Previous studies have observed early AP treatment failures but were not associated with a point mutation in the cytb gene (25,41). Here, we hypothesize that subtherapeutic levels of atovaquone in the patient's plasma from the chemoprophylactic agent prior to the first hospital admission led to increased selective pressure. This could have decreased the time period to for the symptoms to reappear. Indeed, suboptimal drug levels have been correlated with a atovaquone resistance in animal models (15).
Canadian guidelines for the treatment of uncomplicated P. falciparum consist of the administration of AP, or quinine and a second drug (e.g., doxycycline) (1). AP is a less favorable option to fixed-dose artemisinin combination treatments (ACTs) due to the risk of de novo resistance as observed in this study. Given the results reported in this study and other instances of treatment failure with AP in the literature, there is mounting evidence that revision of the treatment guidelines in Canada may be in order. However, access to alternative options need to be in place. Overall, this is the first study to observe an AP early treatment failure consisting of multiple P. falciparum haplotypes present before and during the episode of recrudescence. The data suggests a soft selective sweep of the cytb Y268C mutant that occurred 17 days after the initiation of treatment. No resistant mutants were observed in any of the samples prior to the episode of recrudescence. Amplicon deep sequencing is a useful tool in fully exploring the genetic emergence of resistance in P. falciparum.

MATERIALS AND METHODS
Included samples. Twelve whole-blood samples were drawn and preserved at 280°C during patient care from a patient diagnosed with P. falciparum malaria returning from a 9-week trip to Uganda and South Sudan. The samples were collected before, during, and after treatment failure at day 0, 1, 2, 18, 19, and 20, and were subsequently analyzed for SNP calling and haplotyping (Fig. 1). Ethical approval was obtained from the Conjoint Health Research Ethics Board of the University of Calgary (REB15-1160). Patient written informed consent was obtained.
DNA extraction, PCR, and sample purification. DNA extraction was performed on frozen whole blood using the QIAamp DNA Blood minikit (Qiagen, Germany) following the manufacturer's protocol. PCR was performed on a total of five genes: merozoite surface protein 2 (msp2), cpmp, dhfr, dhps, and cytb. The full-length cytb, dhfr, and dhps genes were amplified by adapting previous protocols (20,42) (Tables S5 to 10). The cpmp PCR protocols were adapted from previous publications (20) (Tables S11 to 14).
The msp2 gene was used for haplotyping using agarose gel electrophoresis and capillary electrophoresis, whereas the cpmp gene was used for ADS haplotyping. A primary PCR was performed on the msp2 gene followed by a nested PCR on the 3D7 and FC27 alleles by adapting previous protocols (43) (Tables S15 to 20). The thermocycling conditions of the nested reaction were taken from the original protocol (43). The primer sequences for the primary and nested reactions were taken from the original protocol as well (43). The nested reaction was performed twice on each sample for the haplotyping analysis to prepare amplicons for agarose and capillary electrophoresis. For the latter, the 6-FAM fluorophore was added to the 59 end of the reverse 3D7 and FC27 primers.
The cpmp, cytb, dhfr, and dhps PCR products (referred to further as "amplicons") were run on a 1% agarose gel for 1 h at 110 V. Each amplicon was visualized in the gel using GelRed Nucleic Acid Gel Stain and exposed to blue LED light in a Blue Light Transilluminator (Maestrogen). Amplicons were then and purified using the QIAquick Gel Extraction kit (Qiagen, Germany). Amplicons with a length below 1,000 bp were compared against a 100 bp ladder (Quick-Load Purple 100 bp DNA Ladder, NEB number N0551), while the remaining were compared against a 1 kb ladder (Quick-Load Purple 1 kb Plus DNA Ladder, NEB number N0550).
Library preparation and sequencing. The cytb, dhfr, and dhps amplicons of each sample were pooled in an equimolar fashion prior to the library preparation step. Two library preparations were generated. The first library preparation was performed across the cytb, dhf, and dhps amplicons by tagmenting using the Illumina DNA prep kit (Illumina, San Diego, CA, catalogue number 20060060) and barcoded with the Nextera XT v2 indexes (Illumina, catalogue number: FC-131-1024). The second library preparation was performed following the amplification of the cpmp gene by following Illumina 16S profiling protocol without tagmentation and barcoded with the Nextera XT v2 indexes (Illumina, catalogue number: FC-131-1024). Sequencing was carried out in two runs for each library using an Illumina MiSeq v3, 600 cycles, a 2% PhiX spike-in, and a maximum output of 25 million reads. A set of contrived samples was included in order to tune the software parameters for the cpmp ADS haplotype calling. This set consisted of six mixtures of the MRA-914 and MRA-1001 strains (100%:0%, 99%:1%, 95%:5%, 90%:10%, 75%:25%, 50%:50%) across six 10-fold dilutions (0.001 copies/mL to 100 copies/mL).
Resistance markers SNP calling. The sequencer fastq files were mapped to the corresponding resistance gene marker references from the PlasmoDB database (44) (cytb: PF3D7_0417200; dhfr: PF3D7_0417200; dhps: PF3D7_0810800) using the Burrows-Wheeler aligner (BWA) (45), followed by duplicate read flagging with Picard (46). The SNP calling parameters and file consolidation were adapted from previously validated methods (47). Briefly, SNP calling was performed with the Genome Analysis Toolkit (GATK) software (48) with SNPs called only if 100 reads mapped to a given position and at least 1% of reads contained a SNP with a quality Phred score of 30. File validation, manipulation, and indexing were done with Picard (46) and Samtools (49).
Gel and capillary electrophoresis haplotyping. Three techniques were employed to determine the haplotype profile before and after treatment failure: (i) agarose gel haplotyping on the family-specific alleles of msp2, (ii) capillary electrophoresis haplotyping of the family-specific alleles of msp2, and (iii) amplicon deep sequencing of the highly heterozygous cpmp gene marker.
The length-polymorphic msp2-3D7 and msp2-FC27 alleles were haplotyped using a 2% agarose gel and capillary electrophoresis. Both methods were adapted from previously described protocols (43). Briefly, the allele-specific amplicons were visually distinguished from each other based on fragment size assessed using a Quick-Load Purple 100 bp DNA Ladder (NEB number N0551) after being separated by a 2% agarose gel electrophoresis for 2 h at 110 V. The gels were visualized using a ChemiiDoc Touch Imaging System (Bio-Rad, version 2.3.0.07). The size and number of bands was estimated using Image Lab software (Bio-Rad, version 6.0.0). Alleles in each family were considered to be identical if the fragment sizes were within 20 bp of one another, as previously described (50,51).
Capillary electrophoresis was performed using the GeneScan 1200 LIZ standard (Applied BioSystems) on a 3730 DNA Analyzer (Applied BioSystems) according to the manufacturer's manual at the University of Guelph. Two replicates of each sample were run along with a negative control. In addition, a set of six mixtures of the MRA-914 and MRA-1001 strains (100%:0%, 99%:1%, 95%:5%, 90%:10%, 75%:25%, 50%,50%) across six 10-fold dilutions (0.001 copies/mL -100 copies/mL) were used to determine the parameters for peak discrimination. To reduce false-positive calls from stutter peaks, companion peaks, and primers, four thresholds were applied: (i) allele-specific fragments below 200 bp were excluded (52); (ii) allele peaks within 3.3 bp of one another were considered to be the same (52); (iii) the allele peaks had to be present in both replicates within 0.4 bp of one another to be considered real; (iv) peaks below 550 relative fluorescent units (RFU) or below 20% of the highest peak were excluded. The analysis of the electropherograms was performed with the Open Source Independent Review and Interpretation System (OSIRIS) software (NCBI, version 2.16) (53).
Amplicon deep sequencing haplotyping. The cpmp gene marker was used for amplicon deep sequencing haplotyping using DADA2 (version 1.22) (54) implemented in R (version 4.1.3). The raw fastq files were demultiplexed. A Trimmomatic (55) sliding window of four nucleotides with a Phred score of ,15 was used to remove low-quality reads. The primers of each read were trimmed with DADA2. The samples were filtered and truncated, followed by error rate estimation using the in-house DADA2 machine learning algorithm (54). Haplotype sample inference was performed on the filtered and trimmed data, followed by read merging and chimera removal with the RemoveBimeraDenovo function in DADA2. These putative haplotypes were then further filtered to prevent the call of false positives by excluding haplotypes (i) with a read depth below 1% (unless present in other independent samples at a minimum read depth of 0.1%), (ii) deviated from the expected nucleotide length of 430 bp, and (iii) a minimum of 100 reads. MAFFT (56) was used in L-INS-i mode to align the identified haplotypes. The complexity of infection was calculated and defined as the number of unique haplotypes within each independent sample.
Data availability. Sequencing data have been deposited at the NCBI BioProject under the accession PRJNA875487. The electropherograms and raw gel picture are available at https://github.com/dcastaneda5/ pillai_lab.git.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, DOCX file, 3.1 MB.