Convergent evolution of artemisinin and chloroquine resistance in Ethiopian Plasmodium falciparum parasites

The emergence of antimalarial drug resistance has threatened the control and elimination of malaria in Africa. Ethiopia, once a success story in case reduction, is now facing resurgence. This study examined the key drug resistance genes (Pfmdr1, Pfcrt, Pfk13, Pfdhfr, and Pfdhps) and mitochondrial genomes of 605 Plasmodium falciparum isolates from 15 districts across Ethiopia, varying in transmission intensity and P. vivax co-endemicity. A dominant PfMDR1 NFSND haplotype, associated with reduced lumefantrine susceptibility, was identified alongside near-fixation of the PfCRT CVIET chloroquine-resistant haplotype in specific areas. Concerningly, PfK13 variants associated with partial artemisinin resistance – R622I (10%), A675V (1.7%), and P441L (1.1%) – were expanding. Drug resistance markers were primarily found in settings with low transmission of P. falciparum and high levels of P. vivax coendemicity. Along with a distinct parasite lineage and limited gene flow, these findings suggest that local evolutionary and therapeutic pressures shape resistance. This underscores the urgent need for targeted genomic surveillance to guide tailored interventions and contain the spread of multidrug-resistant P. falciparum.

Methods
Surveillance sites and sample collection

Fifteen sites across Ethiopia's administrative districts were selected for genomic surveillance of antimalarial resistance (Fig. 1A). Site selection prioritized settings that reported malaria outbreaks [31,68,69], markers of drug-resistant parasites [34],high pfhrp2/3 deletion rates [34], representation of the transmission strata of the country (ranging from low and unstable to high, Fig. 1B), varying P. vivax co-endemicity with P. falciparum (Fig. 1C), and areas prone to importation of parasites from neighboring countries.

A total of 605 blood samples were collected from symptomatic P. falciparum between October 2019 and January 2023. Finger-prick blood samples (0.5 mL) collected in EDTA microtainer tubes were used for malaria diagnosis by microscopy and to prepare dried blood spots (DBS) on Whatman 3MM filter paper.

Ethical Consideration

The study protocol was approved by the Armauer Hansen Research Institute and ALERT Hospital Ethics Review Committee (PO/46/20) and the University of Notre Dame Institutional Review Board (approval no. 20-12-6350). Written informed consent was obtained from all study participants and their parents or guardians before recruitment. The study adhered to the Declaration of Helsinki principles.

DNA extraction, P. falciparum species confirmation, parasite quantification

Genomic DNA was extracted from 6 mm diameter punches of DBS using a MagMAX magnetic bead-based DNA multi-sample kit on a Kingfisher Flex robotic extractor machine (Thermo Fisher Scientific) following the manufacturer’s protocol and eluted in 150 µL low-salt elution buffer. Quantitative polymerase chain reaction[P1]  (qPCR) targeting the P. falciparum 18S rRNA small subunit gene for P. falciparum was performed using primer and probe sequences as previously described [70], using TaqMan Fast Advanced Master Mix (Applied Biosystems). P. falciparum parasites were quantified using standard curves generated from a serial dilution of NF54 ring-stage parasites (10^6 to 10^3 parasites/mL). Samples with parasitemia of 50 parasite/μL or higher were selected for drug resistance marker genotyping. 

Amplification of Antimalarial Drug Resistance Genes, Sample and Library Preparation, and Targeted Deep-Amplicon Sequencing

Amplicon sequencing targeting the entire length of six antimalarial drug resistance genes (Pfcrt, Pfmdr1, Pfk13, Pfdhfr, and Pfdhps) and the mitochondrial genome were employed as described before [71] with few modifications as detailed in Supplementary Information; Single gene PCR amplification assays were modified into two multiplex PCR reactions: Multiplex-I for Pfk13, Pfmdr1, and Pfcytb genes; and Multiplex-II for PfMit, Pfcrt, Pfdhps and Pfdhfr genes amplification. Normalized amplicons pooled from each patient sample were used for library preparation using the Nextera DNA Flex Library Preparation Kit (Illumina, USA). Seven PCR cycles were used for indexing followed by standard purification of the final libraries. Library quantification was performed using a Qubit 4.0 fluorometer with a DNA HS kit (Invitrogen, USA), and fragment size distribution was assessed using a Bioanalyzer 2100 with a DNA HS assay kit (Agilent, Germany). All libraries were normalized to 4 nM, pooled, and sequenced on a NextSeq 500/550 (Illumina, Singapore) using 2 × 149 bp paired-end reads. In total, 605 P. falciparumisolates were successfully amplified and pooled, and sequencing data were obtained for 604 isolates.

Sequence alignment, variant and haplotype identification

Raw Illumina paired-end reads were quality-controlled using the next-generation Sequence Analysis Toolkit  (https://github.com/yyr4/Nf-NeST/). Reads with a Phred score below 20 and length below 100 bp were removed. SNPs were identified using the same pipeline, which employs an ensemble of four variant callers (Samtools, freeBayes, HaplotypeCaller, and VarDict_v1.8.2) to increase the detection confidence. Known drug-resistance-associated SNPs were filtered for a minimum read coverage of 5´ at the corresponding position. Novel SNPs were reported only if they were detected by at least two variant callers and were present in at least ten samples. After quality control and removal of low-depth coverage samples, 586 isolates were retained for downstream analysis. Haplotypes were constructed for Pfcrt (codons 72-76), Pfmdr1 (codons 86, 184, 1034, 1042, and 1246), Pfdhps (codons 431, 436, 437, 540, 581, and 613), and Pfdhfr (codons 16, 51, 59, 108, and 164) using a variant allele frequency threshold of ≥95% for mutant alleles and ≤5% for wild-type alleles [73]. Samples that did not meet these criteria, which are potentially indicative of polyclonal infections, were excluded from the haplotype analysis.

Statistical and Population genetic analyses

Field data were initially collected using paper-based forms and subsequently double-entered into RedCap (mobile application version 5.20.11). Any discrepancies identified during the data entry process were resolved by cross-referencing the original paper forms and reaching a consensus among the data managers.

Associations between haplotype frequencies, the annual parasite index (API), and the fraction of P. vivax and total infections were analyzed using Spearman’s correlation. To adjust for small sample sizes and to achieve robust estimates, 200 bootstrap resampling iterations were used. The API was based on the total confirmed malaria cases and was calculated as follows: (total confirmed malaria cases per year ÷ population at risk) × 1,000. Co-occurrence patterns of key missense mutation pairs were evaluated using chi-square tests, highlighting statistically significant associations. Geographic variation in SNP and haplotype frequencies was assessed using chi-squared tests (statistical significance: p < 0.05), and the Wilson score interval method was used to calculate 95% confidence intervals. The median and interquartile range (IQR) were used to summarize the age distribution across sites.

For haplotype analysis of P. falciparum drug resistance genes (Pfcrt, Pfmdr1, Pfk13, Pfdhps, Pfdhfr), full-gene amplicon sequences were aligned using MAFFT software [74]. Haplotype networks were constructed using R package pegas [75]. SNP-based principal component analysis (PCA) was conducted using the scikit-allel Python package [76]. Publicly available P. falciparum genomic sequences were obtained from previous studies in East Africa and the Horns of Africa [77–79]. A neighbor-joining tree was created using pairwise genetic distance matrices and ape R package [75]. Identity-by-descent (IBD) was used to assess relatedness between isolates by estimating the pairwise fraction of shared ancestry between the sequenced segments. This was done using the hmmIBD software, which applies a hidden Markov model-based approach to estimate shared ancestry while taking recombination into account [81]. Network graphs were created using the igraph R package (v.1.3.5) [82]. The codes used for the bioinformatics and statistical analyses are available at https://github.com/leenvh/EMAGEN.