Emergence and spread of Plasmodium falciparum PX1 polymorphisms associated with decreased susceptibility to antimalarials in Uganda

Abstract Artemisinin-based combination therapies are the cornerstone of malaria treatment and control. In Africa, artemether–lumefantrine is the most widely used first-line artemisinin-based combination therapy, but its efficacy in Uganda is increasingly threatened by the emergence of artemisinin partial resistance and reduced lumefantrine susceptibility. To identify loci contributing to this decreased susceptibility, here we assessed signatures of selection in 157 whole-genome sequences of Plasmodium falciparum from Uganda. Although extended haplotypes were observed around Kelch13 C469Y and A675V mutations, the strongest signal of recent selection was centered on a segment of chr. 7 encoding the phosphoinositide-binding protein (PX1, PF3D7_0720700). A haplotype, represented by three PX1 mutations (L1222P, M1701I and D1705N) and two deletions (designated PIN), was first seen in 2008 and rapidly increased, reaching a prevalence >50% in northern Uganda by 2016 and eastern Uganda by 2023. PIN-carrying parasites showed significantly decreased ex vivo susceptibilities to lumefantrine, mefloquine and dihydroartemisinin, an active metabolite of artemether. A parasite strain in which px1 was disrupted in vitro showed increased susceptibility to the three drugs. Thus, PX1 polymorphisms appear to impact on the susceptibilities of African malaria parasites to key drugs. Subjects Main Artemisinin-based combination therapies (ACTs) are the primary drugs used to treat uncomplicated Plasmodium falciparum malaria1, making them a cornerstone of malaria control2,3,4. Unfortunately, the efficacy of ACTs is under threat due to the emergence of artemisinin partial resistance (ART-R), defined clinically as delayed parasite clearance following treatment with an artemisinin or in vitro as increased parasite survival after artemisinin exposure in the ring-stage survival assay (RSA). This resistance is associated with treatment failures in parts of Southeast Asia5,6,7,8. ART-R first emerged in Southeast Asia and appears to be primarily mediated through a number of different mutations in Kelch13 (K13), resulting in decreased parasite uptake of hemoglobin and reduced artemisinin activation9,10,11,12,13. Recently, the emergence and rapid spread of K13 mutations validated to mediate ART-R have been reported in Rwanda, Uganda and the Horn of Africa14,15,16,17,18,19,20. The most widely used ACT in Africa is artemether–lumefantrine (AL), which was rolled out as the first-line antimalarial for uncomplicated malaria in Uganda in 2006. The mechanism of action of lumefantrine is unknown2,3. Although clinical resistance has not been identified, the ex vivo susceptibility of P. falciparum to lumefantrine in Uganda has decreased in recent years21,22. Decreased susceptibility to lumefantrine has been associated with polymorphisms in P. falciparum multidrug resistance protein 1 (MDR1), including N86, 184F, 500N, 1042N and D1246 and gene duplication, and the chloroquine resistance transporter (CRT) K76 allele22,23,24,25,26,27,28,29. However, these variants mediate only modest changes in drug susceptibility, and validated genetic determinants of lumefantrine resistance have yet to be clearly established. Uganda is an epicenter for emerging ART-R. Five validated or candidate K13 mutations have been detected at concerning prevalences, with A675V and C469Y being the most common, particularly in northern and eastern Uganda. These two mutations were first reported in northern Uganda in 2016 and have now spread across much of the country14. In addition, parasitological surveillance has demonstrated decreasing susceptibility to lumefantrine, first in northern Uganda, and more recently in eastern Uganda21,22,25,30. Decreased susceptibility to lumefantrine has been accompanied by reversion to wild-type MDR1 N86 and CRT K76 alleles, although these genotypes are associated with only small decreases in susceptibility. The clinical consequences of ART-R and decreased lumefantrine susceptibility in Uganda remain uncertain. Recent therapeutic efficacy studies have shown corrected treatment efficacies for AL < 90% at some sites31,32, but decreased AL efficacy was not seen at the site with the highest prevalence of ART-R K13 mutations32, and interpretation of the results is confounded by the difficulty of distinguishing recrudescence and new infections after therapy in high transmission sites. Studies of K13 flanking microsatellite haplotypes have identified independent emergences of Ugandan ART-R parasites14, but few ART-R African isolates have been characterized by whole-genome sequencing (WGS). To directly address this gap, we performed WGS on carefully selected Ugandan P. falciparum isolates and systematically interrogated the genome to identify signals of recent directional selection potentially underlying both ART-R and declining lumefantrine susceptibility17,20,33. Results Whole-genome sequencing of Ugandan isolates Plasmodium falciparum samples collected in northern and eastern Uganda from people with uncomplicated malaria as part of ongoing molecular and parasitological surveillance14,22,34 were subjected to selective whole-genome amplification followed by WGS. A total of 190 low complexity of infection (COI) samples were selected to include temporally matched parasites having (1) the common Ugandan K13 C469Y and A675V mutations, (2) relatively low lumefantrine and dihydroartemisinin (DHA) ex vivo susceptibility or (3) K13 wild-type sequences or relatively high lumefantrine and/or DHA ex vivo susceptibility. To evaluate sample quality before high-throughput sequencing, an Illumina iSeq run was performed. Of the 190 samples, 158 (83%) with fewer than tenfold excess human reads were subsequently sequenced across two Illumina NovaSeq runs (Supplementary Table 1). The majority of these samples (88%) were sequenced to a mean genome-wide coverage of ≥50-fold. One sample that did not reach acceptable coverage was excluded from downstream analyses. Supplementary Table 2 summarizes the distribution of sequenced samples by site, region, k13 genotype and ex vivo drug susceptibility category. Variant calling using the optimized GATK4 pipeline35, followed by recalibration and filtering, yielded 55,100 high-quality single nucleotide polymorphisms (SNPs) with minor allele frequencies (MAF) ≥2% (ref. 36) that were obtained for downstream analysis. Of the 158 sequenced infections, 118 (75%) were estimated to be monogenomic. Evidence of positive selection around the K13 C469Y and A675V mutations To better understand the evolutionary history of the K13 C469Y and A675V mutations, the predominant K13 mutations in northern and eastern Uganda, we compared the extended haplotype homozygosity (EHH) of each mutation to that of the wild-type allele for monogenomic samples and for dominant strains in polygenomic infections (n = 157). The K13 C469Y and A675V mutations each had a strong EHH signal, indicative of a sweep due to positive selection (Extended Data Fig. 1a). To determine the number of unique haplotypes associated with each mutation, we visualized and clustered the flanking variation around the k13 gene using monogenomic C469Y (n = 15) and A675V (n = 11) samples. We identified a single haplotype for C469Y as well as one major and two minor haplotypes for A675V (Extended Data Fig. 1b). Principal component analysis (PCA) showed no clustering by K13 mutation status, indicating substantial outcrossing and arguing against clonal expansion of mutant parasites (Supplementary Fig. 1a,b). Furthermore, neither PCA (Supplementary Fig. 1a) nor the pairwise identity-by-descent (IBD) network (Supplementary Fig. 2) revealed clustering by geographic region, supporting the absence of population structure or lineage-specific background effects. A haplotype centered around px1 gene shows strong signals of positive selection To identify genomic variations potentially responding to AL pressure, we performed genome-wide scans of two complementary measures of positive selection: (1) the isoRelate statistic (iR) based on allele-level pairwise IBD fractions, and (2) the integrated haplotype homozygosity score (iHS). iR analysis detected seven significant peaks: one each on chr. 5, chr. 8, chr. 12 and chr. 13 (the k13 region), and three on chr. 7 (Fig. 1a). Analyses of samples stratified by geographical origin (north, n = 116; east, n = 42) (Extended Data Fig. 2) or by K13 genotype (675V mutant, n = 31; 469Y, mutant n = 35 and wild-type, n = 91) (Extended Data Fig. 3) showed variation in the consistency of these signals by region and by genetic background. iHS analysis identified many of the same regions as iR, including peaks on chr. 7, chr. 8 and chr. 12 (Fig. 1b and Extended Data Fig. 4). Additional peaks not detected by iR were seen on chr. 1 and chr. 10. To focus on mutations potentially driving recent selection, we excluded SNPs that were (1) already common (MAF ≥ 5%) in global P. falciparum samples based on the MalariaGEN Pf6 dataset (Pf6), or (2) under balancing selection (Tajima’s D > 1). With the 15,137 SNPs remaining after filtering, iHS analysis detected 11 nonsynonymous mutations with significant selection signals (adjusted P < 10−5, false discovery rate (FDR)) (Fig. 1b). Overall, the strongest iR signal of selection (Figs. 1a and 2a) involving previously unidentified, rapidly increasing mutations (Fig. 2b) was found on chr. 7 between positions 728081 and 988719 and contained 69 genes (Supplementary Table 3). This candidate sweep was detected consistently in both geographical regions (although the signal was stronger in the north) (Extended Data Fig. 2) and in A675V, C469Y and wild-type K13 genetic backgrounds (Extended Data Fig. 3). Notably, of the 11 filtered SNPs with significant iHS P values, 5 fell within this chr. 7 region (Fig. 2b). The maximal iHS signal co-located with the px1 gene (PF3D7_0720700)37. Among all nonsynonymous SNPs identified across the candidate sweep, PX1 L1222P and D1705N had the highest pairwise IBD fraction (0.016) and the largest delta change in allele frequency from 2017 to 2022 (14.5% and 14.2%, respectively) (Fig. 2b). The iHS signal for PX1 D1705N (P = 1.3 × 10−5, FDR) was the highest in this genomic region and was among the strongest across the entire genome. Genes immediately proximal and distal to the peak signal in px1 lacked significant candidate SNPs based on iHS, apart from PF3D7_0716700 and PF3D7_0721000, genes with unknown functions carrying one (R705I) and two (N2024S and E2044K) SNPs with significant iHS scores (Fig. 2b) and EHH plots (Fig. 2c), respectively38,39,40. These SNPs had lower delta changes in allele frequency compared to px1 variants (Fig. 2b) (2.5% for R705I, 5.5% for N2024S and 7.6% for E2044K). To further assess the five SNPs detected in the candidate sweep, we generated recombination maps for this region using all monogenomic samples from our study (n = 118) and 96 monogenomic controls from the Tanga region in Tanzania available through the Pf7 dataset. In the Tanzanian controls, all three genes fell within regions of high chromosomal crossover activity (recombination rate ρ > 10) (Extended Data Fig. 5). By contrast, among Ugandan samples, px1 showed markedly reduced recombination (ρ = 0.78–1.97), representing the lowest rate in the candidate sweep region. The other genes in the sweep exhibited substantially higher recombination rates, with ρ values up to 38.8 for PF3D7_0716700 and 21.6 for PF3D7_0721000 (Extended Data Fig. 5). A screen of monogenomic samples for structural variations in the region using local reference-free assembly did not detect any large deletions or duplications41 (Supplementary Fig. 3). We also assessed genome-wide linkage disequilibrium (LD) to evaluate potential associations between the sweep region and other genomic loci. Both LD profiles (Extended Data Fig. 6a) and decay (Extended Data Fig. 6b) analyses showed no evidence of long-range associations, indicating that the sweep is restricted to the region surrounding the px1 locus. Considering regions identified by iR, the chr. 12 peak corresponded to the gene encoding the cell-traversal protein for ookinetes and sporozoites (celtos, PF3D7_1216600), a major vaccine candidate42. The two non-PX1 iR peaks on chr. 7, corresponding to erythrocyte binding antigen 175 (eba-175, PF3D7_0731500) and exported protein family 1 gene (epf1, PF3D7_0713200), were significant only in eastern Uganda (Extended Data Fig. 2) and were not significant when stratified by K13 mutation (Extended Data Fig. 3). The chr. 5 signal was no longer significant upon stratification by region or k13 genotype. The chr. 8 peak corresponded to pfa55-14 (PF3D7_0809200), which encodes the asparagine-rich antigen that was previously reported to be under directional selection in sub-Saharan Africa and Southeast Asia40. This signal was present in both regions of Uganda (Extended Data Fig. 2) but appeared to be limited to parasites carrying the wild-type k13 allele (Extended Data Fig. 3). The signal on chr. 12 was limited to eastern Uganda and neared significance only in parasites carrying the wild-type k13 allele (Extended Data Figs. 2 and 3). Genetic regions only identified by iHS similarly corresponded to loci previously reported to show signs of positive selection in African parasite populations, consisting mainly of genes encoding known antigens, such as exported protein family 1 (EPF1) and PF3D7_0710200 (refs. 40,43,44,45). Thus, proteins encoded by genes other than px1 with iHS and iR selection signals were less compelling as candidate mediators of drug resistance, because they had characteristics consistent with evolving antigenic variation or were not consistently significant across populations, and most of them have previously been shown to be under selection in multiple countries in whole-genome analyses40,43,44,45. A haplotype block centered around the px1 gene Visualization and clustering of variations in the candidate sweep region (~260 kb, Pf3D7_07_v3:728081–988719) were performed using sequences from monogenomic samples with no coverage gaps (L1222P and D1705N mutant, n = 46; wild-type, n = 21) and indicated that all implicated SNPs belonged to a single shared haplotype (Supplementary Fig. 4). Although the length of this haplotype varied by isolate, the most conserved segment across samples (Pf3D7_07_v3:873736–918738) coincided with the peak of IBD fractions (Figs. 2a,b and 3a and Supplementary Fig. 4) including 12 genes with px1 at the center. We named this shared haplotype PIN for the three nonsynonymous mutations, L1222P, M1701I and D1705N, with the highest frequency changes and measures of EHH. The PIN haplotype was found in 58.2% of isolates sequenced and the wild-type haplotype, LMD (L1222, M1701 and D1705), in 28.5%. All other detected haplotypes (LID, LIN, PID, PMD and PMN) were represented in <10% of isolates. Two other px1 SNPs (D384A and S1673N) were weakly associated with the PIN haplotype: the D384A mutation was present in 82.8% of PIN samples (77 of 93) but showed a nonsignificant iHS score (P > 10−5, FDR), and the S1673N mutation was present in all PIN samples and showed a nonsignificant iHS score (P = 0.18, FDR). Of note, the S1673N mutation was also present in 66.7% of non-PIN haplotype samples and appeared at high frequency in the Pf6 dataset. Overall, sampling locations and K13 mutations were not associated with clustering in the PIN genomes (Supplementary Fig. 4). The shared haplotype appeared to shorten over time, particularly in northern Uganda (Supplementary Fig. 4). The px1 gene contains four exons and is ~9 kb long, with conserved C-terminal and N-terminal regions (Fig. 2b). Each of the SNPs discussed above is located in exon 3, which measures 5.78 kb (nucleotides 892083–897457) and also encodes a repetitive region with evidence for deletions (Supplementary Fig. 4). Oxford Nanopore Technologies (ONT) long-read sequencing tiling the px1 gene resolved 5 in-frame deletions in the coding sequence of 24 representative monogenomic samples (Supplementary Fig. 5). Among the deletions in the PIN haplotype, two (delV1680–N1685 and delN811–Y822, deletions of 18 and 36 nucleotides at positions 897247 and 894638, respectively) were uncommon (MAF < 5%) in the Pf6 dataset and showed high EHH signal in our Uganda WGS data (Extended Data Fig. 7). Prevalence of the px1 PIN haplotype has been increasing rapidly in Uganda To assess the spatiotemporal prevalence of the px1 PIN haplotype, we genotyped 1,598 available samples collected from eastern Uganda in 2004, 2008 and 2012, and from eastern and northern Uganda from 2016 to 2024. Using one of the primer pairs designed to tile the px1 gene, a 2.6-kb region (894851–897457) spanning all the key variants was amplified and ONT sequenced (Supplementary Table 4 and Supplementary Fig. 5). A total of 1,436 samples had sequencing coverage ≥25× and underwent haplotype calling. In 2004, before ACTs were recommended for treatment of malaria in Uganda, the most common haplotypes seen were LMD (76.9%) and LID (12.1%); the PIN and PMN haplotypes were not found (Fig. 4). The PIN haplotype was first seen in 2008, significantly increasing prevalence over time in both regions, reaching 55% in eastern Uganda (P = 0.009, Mann–Kendall test) and 84% in northern Uganda (P = 0.009, Mann–Kendall test) in 2024 (Fig. 4). The same trend was observed after analyzing the prevalence over time by sampling site using a Bayesian model (Supplementary Fig. 6). We evaluated the co-occurrence of K13 and putative ACT partner drug resistance mutations, as well as PIN haplotypes, in samples from northern Uganda, where K13 mutations first emerged and were more prevalent than in other regions14. Compared to wild-type K13, the C469Y and A675V mutations were consistently more prevalent in parasites carrying the PIN haplotype over time, although statistical assessment was limited by small sample sizes upon stratification (Supplementary Fig. 7). The proportions of MDR1, CRT, dihydrofolate reductase and dihydropteroate synthase mutations associated with resistance to other antimalarials were similar in PIN and LMD parasites, suggesting no co-occurrence of these markers with the PIN haplotype (Supplementary Fig. 8). The px1 PIN haplotype was associated with reduced ex vivo susceptibility to DHA and lumefantrine We evaluated associations between the px1 haplotypes and ex vivo drug susceptibility22, considering samples with ≥25× sequencing depth in which the called allele represented ≥90% of reads to reduce ambiguity introduced by mixed genotype infections (n = 465). Overall, the PIN haplotype was associated with decreased DHA, lumefantrine and mefloquine susceptibilities (increased half-maximal inhibitory concentration (IC50)) compared to the wild-type (LMD) haplotype (P < 0.001 for all comparisons, Wilcoxon test) (Supplementary Fig. 9). When considering samples with wild-type K13, the PIN haplotype was associated with decreased susceptibility compared to the LMD haplotype for lumefantrine (median IC50 [interquartile range (IQR)]: 14.3 nM [9.8–24.0], n = 110 versus 6.2 nM [4.4–10.5 nM], n = 240), mefloquine (17.1 nM [10.4–26.7], n = 109 versus 11.0 nM [7.7–16.0], n = 236), and DHA (3.7 nM [2.4–5.5], n = 110 versus 1.8 nM; [1.2–2.6], n = 240) (Fig. 5a). Similar trends were seen with samples carrying the K13 C469Y and A675V mutations, although the small number of isolates with the PIN haplotype and K13 mutations (n = 26 for C469Y and n = 11 for A675V) limited assessment of statistical significance (Fig. 5a and Supplementary Fig. 10). Notably, no significant differences in ex vivo RSA survival between haplotypes were detected (n = 169), consistent with our recent observation that ex vivo RSA results from 2019–2024 did not correlate with other markers of drug susceptibility in Uganda22. Similar results were found when stratifying samples based on year or geographical origin (Supplementary Fig. 11a,b). Mixed-effects modeling controlling for site, parasitemia, COI, sampling year, and k13 genotype revealed that the PIN haplotype was associated with reduced drug sensitivity for both DHA and lumefantrine (Extended Data Table 1). For DHA IC50, the PIN haplotype was the most influential biological factor, accounting for 12.3% of the total phenotypic variance (ΔR2 = 0.123) in K13 wild-type samples. Parasites carrying the PIN haplotype exhibited a significantly higher estimated marginal mean DHA IC50 (3.1 nM) compared to the LMD haplotype (2.3 nM; P = 3.0 × 10−8), with a moderate Cohen’s effect size (0.44). Notably, the substantial gap between marginal (R2marginal = 0.131) and conditional (R2conditional = 0.264) R2 values indicates that geographic and temporal factors (site and year) contributed an additional 13.3% to the variance in DHA susceptibility. In lumefantrine and mefloquine IC50 assays, the PIN haplotype also contributed significantly to the model, explaining 4.9% (ΔR2 = 0.049) and 6.5% (ΔR2 = 0.065) of the variance, respectively. PIN-carrying parasites showed a markedly higher marginal mean lumefantrine IC50 (14.2 nM) than those with the LMD haplotype (8.8 nM; P = 2.5 × 10−3). For lumefantrine, the total model accounted for 18.7% of the variance (R2conditional = 0.187), suggesting that although the genetic effect of the PIN haplotype is clear, susceptibility is also likely shaped by broader contextual factors. By contrast, the PIN haplotype was not associated with RSA survival (R2conditional = 0.093, P = 0.13). To determine whether the influence of px1 variant on drug susceptibility was dependent on the presence of K13 mutations (C469Y and A675V), we added an interaction term between the two genes to our linear mixed-effects model. For lumefantrine and mefloquine IC50 values, and RSA, no significant interaction was observed (P > 0.65), suggesting the effects of these variants are independent (Extended Data Table 2). However, for DHA IC50, the interaction term approached statistical significance (P = 0.061), hinting at a potential synergy or modification of the phenotype when both variants are present. In vitro loss of px1 leads to decreased lumefantrine, mefloquine and DHA IC50 values Building on the previously published genetic and cellular characterization of PX1 using selection-linked integration-based targeted gene disruption37, we next sought to functionally validate the contribution of PX1 to antimalarial susceptibility. We tested two independently derived 3D7 px1-disrupted (knockout) clones (PX1-KO-1 and PX1-KO-2) in standardized 72-hour in vitro susceptibility assays against chloroquine, cycloguanil, mefloquine, lumefantrine and DHA. Chloroquine IC50 values for both knockouts were comparable (P > 0.05, one-way analysis of variance (ANOVA) with Tukey’s multicomparison) (Fig. 5b) to the unmodified 3D7 parent and the PX1–GFP integration control, indicating no effect of px1 loss on chloroquine sensitivity. As expected, both PX1–GFP and the knockout lines showed equivalent increased resistance to cycloguanil relative to 3D7 presumably because of the human dihydrofolate reductase cassette rather than px1 disruption (P > 0.05, one-way ANOVA with Tukey’s multicomparison) (Fig. 5b). By contrast, both knockout clones displayed marked hyper-susceptibility to lumefantrine, DHA and mefloquine, with significantly reduced IC50 values relative to 3D7 (all comparisons by one-way ANOVA with Tukey’s multicomparison) (Fig. 5b). Relative to 3D7, PX1-KO-1 and PX1-KO-2 showed 1.87-fold (P = 0.0237) and 1.75-fold (P = 0.0357) reductions in DHA IC50, respectively. For lumefantrine, hyper-susceptibility was even more pronounced, with PX1-KO-1 showing a 3.16-fold (P = 0.0226) and PX1-KO-2 a 3.89-fold (P = 0.0145) decrease in IC50 relative to 3D7. Similarly, PX1-KO-1 and PX1-KO-2 exhibited 3.56-fold (P = 0.0030) and 3.23-fold (P = 0.0038) reductions in mefloquine IC50, respectively. Importantly, PX1–GFP showed no significant differences from 3D7 for these drugs, confirming that the hyper-susceptibility phenotype was attributable to loss of PX1 function rather than plasmid integration effects. The px1 PIN haplotype was rare in global parasite populations in the Pf6 dataset The Pf6 WGS dataset, comprising 6,388 samples collected between 2001 and 201546, was interrogated to assess historical global distribution of px1 alleles35. px1 allele frequencies varied geographically, and the PIN haplotype was rare, occurring in only 5 of 3,570 samples, all from Africa (Fig. 6 and Supplementary Figs. 12 and 13). These samples were collected in the Democratic Republic of Congo (DRC, two in 2012 and one in 2014, n = 364) and Kenya (two in 2014, n = 110), countries neighboring Uganda. The PIN haplotype was present as a mixed genotype in all but one case. Notably, the PIN haplotype was not detected in the few Ugandan samples collected in 2010 (n = 13) in this dataset (Fig. 6). The PMN haplotype was found in one polygenomic sample from Kenya in 2014. The LID haplotype was at high prevalence world-wide (Supplementary Fig. 13). The LIN haplotype was at a maximum prevalence of 9.3% (34 of 364) in the DRC. In two high-quality samples from the DRC, we confirmed the two associated deletions were present (Supplementary Fig. 14). Discussion The decreased lumefantrine susceptibility and emergence of ART-R reported in Uganda14,21,22,25,30,47,48,49 are a major concern for malaria control in Africa, prompting a comprehensive search for genomic regions of directional selection that might underlie these changes. We leveraged WGS of Ugandan P. falciparum isolates with varied K13 mutations and ex vivo drug susceptibilities to identify a locus associated with decreased susceptibility to DHA and lumefantrine, the most widely used antimalarials in Africa. Based on iR and iHS analyses, the strongest signal of recent selection was on chr. 7, centered on the px1 gene. The px1 PIN haplotype carrying L1222P, M1701I, D1705N mutations and two deletions showed the strongest signal of selection, with a long-shared flanking haplotype. This haplotype was not detected in samples collected in 2004, before rollout of AL in Uganda, but increased in prevalence thereafter, and appears to have emerged before K13 mutations. Importantly, isolates with the PIN haplotype had significantly decreased ex vivo susceptibilities to lumefantrine and DHA compared to parasites with the wild-type px1 haplotype. A role for PX1 in modulating susceptibilities to these drugs is further supported by the increased in vitro susceptibilities demonstrated in px1 knockouts. Thus, remarkably, PX1 mutations are associated with and may mediate decreased susceptibility to both components of AL. Our genomic findings suggest that px1 has been the target of strong directional selection in Ugandan malaria parasites, leading to a steady increase in frequency of the PIN haplotype. The candidate sweep signal encompassed a large region (~260 kb) containing 69 genes, among which a core haplotype (Pf3D7_07_v3:873736–918738) of 12 genes, including px1, showed limited recombination. Our identification of px1 as the selection target was based on an array of evidence. First, px1 was in the center of the iR peak, a P. falciparum-optimized measure of selection50. Second, this peak had the highest IBD fraction across the genome. Third, the PX1 D1705N mutation had the strongest iHS and EHH signals among all SNPs located in the candidate sweep region. Fourth, L1222P and D1705N mutation frequencies had the highest delta changes over time. Fifth, there was an abundance of synonymous and noncoding SNPs surrounding the px1 gene, a strong indicator of genetic hitchhiking resulting from rapid selection. Sixth, px1 had the lowest recombination rates among all the genes located in this candidate sweep. Emergence of the PIN haplotype in Uganda has important implications for our understanding of the spread of ART-R and decreased lumefantrine susceptibility in the country. Our findings suggest that PIN emerged in Uganda before the emergence of the K13 C469Y and A675V mutations, because it was found in samples collected in 2008, whereas K13 mutations linked to ART-R were first seen in samples collected in 2016. Based on the high proportion of PIN-carrying parasites from northern Uganda with C469Y and A675V, the haplotype may have been part of the background on which K13 mutations emerged, and its role in facilitating the emergence of K13-mediated ART-R warrants further investigation. Although the PIN haplotype was often seen in K13 mutant parasites, decreased lumefantrine and DHA susceptibilities were linked to this haplotype independent of K13 mutations associated with ART-R and of MDR1 mutations linked to decreased lumefantrine susceptibility. The PIN haplotype is now present in the vast majority of parasites in northern Uganda and it appears to be quickly rising in frequency in eastern Uganda. In addition, emergence of the PIN haplotype was temporally and geographically associated with decreasing susceptibility to lumefantrine, first seen in northern and later in eastern Uganda21,22,25,30. Overall, these findings suggest that heavy exposure to AL, beginning in about 2006, led to selection of the px1 PIN haplotype and then the K13 mutations that mediate ART-R. Whether the selection of PIN is representative of what has happened elsewhere in Uganda will require additional sequencing efforts, because current surveillance panels do not cover the locus. In addition, broad surveillance is needed in neighboring countries because, given the observed selection and high-levels of human movement across borders, it is likely that the PIN haplotype has spread to other countries. However, these data are currently lacking and it is possible that differences in transmission intensity, treatment coverage and drug use patterns or parasite population structure may have limited the geographical spread of the haplotype1,51,52. The PIN haplotype was associated with decreased ex vivo susceptibility to lumefantrine, mefloquine and DHA. In vitro testing of px1-disrupted (knockout) parasites showed hypersensitivity to these drugs. Together, these findings support a role for PX1 in modulating susceptibility to lumefantrine, mefloquine and DHA. Prior studies demonstrated impaired growth and reduced hemoglobin transport to the parasite food vacuole in px1 knockout parasites37. The role for PX1 in facilitating efficient hemoglobin trafficking to the digestive vacuole37 appears analogous to that of K13 (ref. 12), and is consistent with PX1 loss mediating artemisinin resistance by decreasing hemoglobin processing. However, despite their functional parallels, PX1- and K13-mediated artemisinin resistance are likely to differ mechanistically. Notably, px1 knockout parasites display opposing phenotypes in stage-specific assays: increased survival in 0–3-hour RSA indicative of ART-R37, but hyper-susceptibility in 72-hour growth-inhibition assays, suggesting that the stage specificity of PX1 involvement with artemisinin action differs from that of K13. The ex vivo version of the RSA detected a modest but nonsignificant difference in drug response between PIN and LMD parasites, potentially reflecting limited statistical power (n = 169; Extended Data Table 1) and the use of unsynchronized field isolates. Further work is required to define the mechanistic basis of PX1-mediated drug responses and the functional linkage between susceptibility to DHA, lumefantrine and mefloquine. Although px1 knockout data support a role in drug responses, targeted studies are needed to determine how the PIN haplotype modulates these effects. Most importantly, clinical investigations are essential to determine the therapeutic consequences of the px1 PIN haplotype. In vitro drug susceptibility assays have yet to be shown to directly predict treatment outcomes. Of note, our ex vivo assays were performed with Albumax serum substitute, without plasma lipoproteins that bind lumefantrine, resulting in increased free drug concentrations and decreased IC50 values compared to assays supplemented with serum. Whether our measured shift from a lumefantrine IC50 of ~6 nM for parasites with the wild-type LMD haplotype to ~14 nM for parasites with the PIN haplotype is sufficient for changes in treatment responses is unknown. However, reports of treatment failures in nonimmune travelers and >10% polymerase chain reaction (PCR)-corrected failures in multiple clinical trials testing the efficacy of AL suggest that marginal shifts in susceptibility may be clinically relevant32,47,49. Collectively, our findings suggest that px1 polymorphisms contribute to decreased susceptibility to both components of AL, the most widely used ACT in sub-Saharan Africa. The rapid increase in prevalence of the PIN haplotype and its association with reduced ex vivo drug susceptibility underscore the need for its inclusion in molecular surveillance panels and for clinical studies to determine its impact on treatment outcomes. Methods Study design and sample collection For WGS activities, we leveraged clinical samples collected as part of ongoing health facility-based molecular14,34,53 and parasitological21,22 surveillance activities. Briefly, molecular surveillance samples were collected from individuals aged >6 months diagnosed with malaria by rapid diagnostic test or microscopy at up to 16 health facilities across Uganda between 2016 and 2024. Following written informed consents, dried blood spots were collected by finger prick. For parasitological surveillance, individuals aged >6 months diagnosed with high parasitemia malaria by microscopy at health facilities near parasitology laboratories in Tororo, Tororo District in eastern Uganda and Kalongo, Agago District in northern Uganda between 2016 and 2024 were consented and up to 5 ml of blood was collected into heparin tubes by venipuncture. From these samples, we used pre-existing genotyping and ex vivo drug susceptibility data14,22 to select low COI samples with K13 mutations, low lumefantrine and/or DHA susceptibility or high RSA survival. Each of these samples was then matched by collection year and site with a low COI sample encoding a wild-type K13 allele or having unremarkable lumefantrine and DHA susceptibility profiles. To estimate prevalences of px1 genotypes over time, we performed long-read ONT sequencing on a random subset of 50 samples that had undergone ex vivo drug susceptibility assessment from each site for each year of surveillance activity (50 each year from 2016 to 2024 for the eastern region and 50 each year from 2021 to 2024 for the northern region). For each year when ex vivo samples were not available from the north, we sequenced a random subset of 50 molecular surveillance samples collected from the Patongo health facility. Finally, to provide an understanding of changes in px1 diversity in the early stages of AL utilization, we evaluated 91 pretreatment samples collected as part of a 2004 therapeutic efficacy study54 and 92 samples collected per year in 2008 and 2012 from children (aged <5 years) enrolled in a cohort study55,56. Both studies were conducted in Tororo. Library preparation, whole-genome sequencing and variant calling Genomic DNA extracted from dried blood spots underwent two rounds of specific whole-genome amplification, as previously described57. The amplified products were combined, and WGS libraries were prepared using the Watchmaker DNA Library Kit with Fragmentation (Watchmaker Genomics). The resulting libraries were pooled and sequenced using Illumina 2 × 150 bp chemistry at Psomagen on an Illumina X Plus (Psomagen). After sequencing, Trimmomatic was used to trim off adapters and select properly paired reads before mapping. Reads were competitively mapped onto a hybrid reference genome obtained from the concatenation of P. falciparum 3D7 (version 3.1) and human genome assembly (version GRCh38) using BWA-MEM (version 0.7.17-r1188). We used Samtools (version 1.19.2) and GATK (version 4.3.0.0) to select and clean reads that specifically mapped to the P. falciparum genome. Samples with a human-to-parasite read ratio <10 were retained and those among these with low sequencing depth (first quartile of read depth <35×) were repooled and rebalanced for another NovaSeq X Plus run. Cleaned binary alignment map files from different sequencing runs were merged before variant calling using a P. falciparum-optimized GATK4 pipeline (https://github.com/Karaniare/Optimized_GATK4_pipeline/tree/main), as previously described35. An accurate in silico positive training dataset built in the pipeline was used for machine learning variant recalibration accounting for multiple mapping parameters including read depth, mapping quality and strand bias. Variants that failed this filtering were removed, as were samples and variants with genotype missingness >10% and >20%, respectively. Subtelomeric and internal hypervariable regions that are hard to map were excluded from the variant call format (VCF) file to focus the downstream analysis on the core genome as previously defined58. The fraction of reads supporting the alternate allele was added in the format field to enable detection of major alleles in mixed infection samples. Estimation of complexity of infection We selected high-quality SNPs with MAF ≥ 2% and <10% genotype missingness to estimate COI using The REAL McCOIL package59 as implemented in the MIPTools pipeline. The total number of Markov chain Monte Carlo (MCMC) was set to 2,000 with 500 burn-in iterations. Selection analysis The rehh R package (version 3.2.2)60 was used to estimate the EHH around specific makers and to scan the genome for allele-specific iHS signals from filtered VCFs. An initial analysis was performed with all the SNPs at MAF ≥ 2% in monogenomic samples. Tajima’s D analysis was also performed in this SNP set using VCF-kit (version 0.2.9)61 to identify balancing selection signals, likely reflecting immune-mediated pressure. For the iHS scan, SNPs with Tajima’s D > 1 or present at MAF ≥ 5% in the Pf6 dataset from samples collected up to 2015 were excluded to enrich for signals of recent directional selection. Raw iHS values were standardized within derived-allele frequency bins (bin size = 50) to correct for allele-frequency dependence. Statistical significance was assessed assuming a standard normal distribution of the standardized iHS values, and two-sided P values were calculated as P = −log10(2Φ(−|iHS | )), where Φ is the cumulative distribution function of the standard normal distribution, corresponding to a two-sided Z-test with an empirically derived null distribution60. A minimum of four haplotypes was required for evaluation at each locus. Multiple testing correction was performed using the Benjamini–Hochberg procedure as implemented in the R function p.adjust to control the FDR. isoRelate (version 0.1.0) iR statistics was also used to scan the genome for recent positive selection signals based on IBD50. IBD segments were inferred allowing a genotyping error of 0.001, using all SNPs with MAF ≥ 2% following recalibration filtering. Only shared IBD segments ≥50 kb in length and supported by at least 20 SNPs between sample pairs were retained. The isoRelate function getIBDiR was used to compute pairwise iR statistics per SNP across the full sample set and after stratification by K13 mutation status or sampling regions. Statistical significance of iR was evaluated using empirical two-sided P values derived from the genome-wide distribution of iR statistics, defined as the proportion of loci with absolute iR values greater than or equal to the observed value. Multiple testing corrections for iR were conducted as for iHS. To further characterize the selection sweeps detected across the genome, SnpEff annotations were used to identify whether SNPs are nonsynonymous or synonymous or from noncoding regions. The delta changes of allele frequencies for these SNPs over time were also calculated. For more robust analysis of spatiotemporal change in PX1 mutations, prevalences of detected haplotypes were calculated from 2008 to 2024 in the east and from 2016 to 2024 in the north. A linear regression model was used to measure the increase of haplotype frequencies over time in each region. Haplotype visualization A subset of the quality-filtered WGS VCF containing SNPs with MAF ≥ 2% from the px1 flanking region was extracted. Polygenomic and low coverage (read depth <50) samples were removed. The VCF subset was converted into a genotype table with two values (0 and 2) containing samples in the rows and SNP positions in the columns. The genotype table was visualized using ComplexHeatmap R package (version 2.21.1) (https://jokergoo.github.io/ComplexHeatmap-reference/book/a-single-heatmap.html). The complete-linkage clustering method was applied to cluster samples based on their relatedness. Metadata variables were added as bar plots to the heatmap and included PX1 mutations, K13 mutations, region and year of sample collection. Recombination analysis The LDhat package (version 2.2a)62 was used to estimate per locus recombination rates using SNPs with MAF ≥ 1% and monogenomic samples. A likelihood lookup table of 192 sequences was generated to compute sample recombination rate profiles using composite likelihood and piecewise constant model62,63. A total of 10,000,000 MCMC iterations and a background block penalty of 5 were used, as well as 100,000 burn-in iterations and 2,000 MCMC iterations between samples. px1 genotyping using Oxford Nanopore Technologies To resolve px1 haplotypes, we designed primers (Supplementary Table 5) tiling across the gene using the Multiply2 package (https://github.com/JasonAHendry/multiply) developed for multiplex PCR panel design for ONT sequencing. Minimum and maximum amplicon sizes were 1,605 and 2,319 bp (Supplementary Table 5), respectively. Two rounds of long-range PCR were performed with GoTaq master mixes (Promega). The first PCR (30 cycles) amplified targeted regions from the genomic DNA template using each primer set in separate simplex reactions (Supplementary Table 5). Forward (5′GACTCGCCAAGCTGAAGNNNN3′) and reverse (5′ACGTGTGCTCTTCCGATCTNNNN3′) primers were attached to linkers, oligo-sequences containing binding sites for Illumina barcoding primers. The second PCR (ten cycles) was performed to barcode the products of the first PCR using Illumina primers attached to indexes. The primer set P2_Px1Block3_v2_F/P2_Px1Frag3_v8_R (894851–897457) spanning the PIN haplotype variants was used for large scale genotyping. After barcoding, all samples were pooled and bead cleaned for ONT library preparation. A ligation sequencing kit (SQK-LSK114) was used without ONT barcoding, because samples were already barcoded using Illumina indexes. The library pool was sequenced using a PromethION 2 Solo device. Duplex super-accurate base calling was performed from pod5 files using dorado (https://github.com/nanoporetech/dorado). The sam file obtained from the base calling step was converted into a fastq file before demultiplexing using the elucidator package (https://github.com/nickjhathaway/elucidator). The reads were mapped onto Plasmodium falciparum 3D7 reference genome (version 3.0) using minimap2.0 (https://github.com/lh3/minimap2) and the bcftools package (version 1.13) was used for variant calling with ONT-optimized settings (https://samtools.github.io/bcftools/howtos/variant-calling.html). Samples with read depth <25× (10 of 1,608) were removed from the analysis. Samples carrying mutations that are supported by fewer than 50% of reads (88 of 1,608) were also excluded from the analysis to prevent any ambiguity due to contamination during the PCR steps. In vitro 72-hour drug susceptibility assay of genetically manipulated parasites px1 genetically manipulated parasites37 were grown in complete RPMI medium in the presence of 2.5 nM WR99210 (Jacobus Pharmaceuticals) and washed human red blood cells (Interstate Blood Bank) at 5% crit at 37 °C under a gas mixture of 5% CO2, 5% O2 and 90% N2. Incomplete medium consisted of RPMI 1640 with l-glutamine (Gibco), supplemented with 50 mg l−1 hypoxanthine (Calbiochem) and 25 mM HEPES (Corning). Complete medium was prepared by adding 0.5% Albumax II (Gibco), 10 mg l−1 gentamicin (Gibco) and 0.225% NaHCO3 to incomplete culture medium. Cultures were enriched for schizonts by Percoll gradient centrifugation and left to invade RBCs for 6 h. Parasites at 2% hematocrit and 0.2% parasitemia were grown for 72 h in the presence of different concentrations of drugs in 96-well plates. Growth at 72 h was measured by SYBR Green (Invitrogen) staining of parasite DNA on a plate reader (Flexstation 3). A dilution series of the drugs was carried out in three biological replicates. Relative fluorescence units were measured at an excitation of 490 nm and emission of 525 nm on a plate reader and analyzed using GraphPad Prism version 10 (GraphPad Software). IC50 values were determined with the curve-fitting algorithm log(inhibitor) versus response-variable slope. Geospatial mapping of px1 haplotype frequencies The VCF obtained after reanalyzing the Pf6 dataset with the optimized GATK4 variant calling pipeline was used to detect px1 genotypes in multiple parasite populations. The prevalences of px1 haplotypes were calculated for each country and plotted on the map (https://www.naturalearthdata.com/) using igraph (version 2.3.0), sf (version 1.0-15), ggrepel (version 0.9.8), ggspatial (version 1.1.10) and rnaturalearth (version 1.2.0) R packages. Population structure analysis For allele-based population structure analysis, we selected quality-filtered SNPs with MAF ≥ 2% and LD < 0.20. PLINK (version 2.0.0-a.6.9) was used to calculate the pairwise variance-standardized genetic relationship matrix on which we performed PCA using the FactoMineR (version 2.14) package in R, as previously implemented64. We used factoextra to visualize the results of the PCA. Identity-by-descent network analysis Genetic relatedness among samples was assessed by estimating pairwise IBD using filtered genome-wide SNPs. Analyses were restricted to major alleles and SNPs with a MAF ≥ 2%. IBD inference was performed using the hmmIBD package (version 2.0.0)65, with a maximum of 5 fitting iterations, a requirement of at least 200 informative sites per sample pair and an assumed genotyping error rate of 0.1%. IBD networks were constructed using the igraph package (version 2.3.0). Copy number variation analysis Copy number variation was analyzed using PathWeaver (version 1.0), a de novo assembler optimized for P. falciparum41 (https://github.com/nickjhathaway/PathWeaver). Accessible regions of chr. 7 were identified by excluding short repetitive sequences using a tandem repeat finder function. PathWeaver uses a de Bruijn graph-based assembly strategy with iterative recruitment of unmapped reads to improve assembly accuracy and provides region-specific summary metrics, including read depth used for each assembly. Per-base coverage across assembled regions was used to assess structural variation along chr. 7. For each sample, per-base coverage values were first normalized by the mean coverage across all assembled regions. These sample-normalized values were subsequently re-normalized by the median normalized coverage for each region across all samples. Linkage disequilibrium analysis Filtered genome-wide SNPs with MAF ≥ 2% were used to estimate LD among SNPs in the candidate sweep and between these SNPs and the remainder of the genome using PLINK66. Pairwise LD was quantified as the squared correlation coefficient (r2), computed using the default maximum pairwise distance of 10,000 kb. Pairwise r2 values between candidate sweep SNPs and SNPs across the genome were visualized to assess LD profiles along each chromosome. To characterize LD decay along chr. 7, mean r2 values were calculated in bins of increasing pairwise genomic distance, with bin sizes incremented in 100-bp intervals. Locally estimated scatterplot smoothing was applied to fit a smooth curve to the mean r2 values as a function of genomic distance. Ex vivo drug assays For the published ex vivo data leveraged, drug susceptibilities were assessed using a 72-hour growth-inhibition assay with SYBR Green detection, as previously described22,25. Briefly, threefold serial dilutions of 10 mM (50 mM for pyrimethamine) stocks of chloroquine, monodesethylamodiaquine (the active metabolite of amodiaquine), piperaquine, pyronaridine, mefloquine, lumefantrine, DHA, quinine and pyrimethamine, in Albumax-supplemented complete medium were placed in 96-well microplates (50 μl per well), including drug-free and parasite-free controls. Parasites were diluted with uninfected erythrocytes and added to diluted drugs to a final culture volume of 200 μl at 0.2% parasitemia and 2% hematocrit. Plates were maintained at 5% CO2, 5% O2 and 90% N2 for 72 h at 37 °C in a humidified modular incubator. After 72 h, cells were lysed and stained with 100 μl of SYBR Green lysis buffer, incubated for 1 h in the dark at room temperature and fluorescence (485 nm excitation and 530 nm emission) was measured. IC50 values were derived from plots of fluorescence intensity versus log(drug concentration) and fit to nonlinear curves using a four-parameter Hill equation in Prism. The ex vivo RSA was performed as previously described25. Briefly, clinical samples were diluted with uninfected erythrocytes to no more than 1% parasitemia and were incubated with 700 nM DHA or, for controls, 0.1% dimethyl sulfoxide for 6 h. Cells were then washed to remove the drug and cultures were maintained for an additional 66 h. Giemsa-stained thin smears were then prepared for exposed and control parasite cultures, parasitemias were counted and RSA survival was calculated as the proportion of viable parasites in the DHA-treated cultures relative to controls. Genotyping of key antimalarial drug resistance markers All dried blood spot samples (n = 1,035) were subject to DNA extraction using a Chelex-Tween protocol as previously described67,68. Comprehensive genotyping of existing antimalarial drug resistance markers was performed using both molecular inversion probe and Sanger sequencing. For molecular inversion probe, the DR2 panel was used as previously described67,68. Briefly, drug resistance targets were captured by oligo-probes into DNA circles which are enriched by digesting linear DNA using exonucleases. A final PCR amplification step was performed to open the circular DNA and add indexes for Illumina sequencing using NextSeq 550. After sequencing, read processing and Freebayes variant calling steps were done using the MIPTools package (https://github.com/bailey-lab/MIPTools). Drug resistance mutation data was filtered and analyzed in R using the miplicorn package (https://github.com/bailey-lab/miplicorn). Dideoxy sequencing was performed on a subset of samples to supplement k13 genotypes as previously described53. Sequences were evaluated using CodonCode Aligner version 9.0.1 (CodonCode Corporation). Public whole-genome sequencing dataset Raw WGS reads from MalariaGEN’s Pf6 release46 were previously downloaded from the Sequence Read Archive and analyzed with the optimized GATK4 pipeline by our laboratory35. The VCF obtained from this variant calling was leveraged to detect common alleles across the genome in global parasite populations. The same data was utilized to estimate and map (Supplementary Methods) frequencies of PX1 mutations in older samples from malaria-endemic countries. Statistical analysis Microsoft Excel (versions 97-2003 and 2016) was used for data entry. Data analysis was performed using R (version 4.3.1). Data visualization was performed using R and Adobe Illustrator CC (version 17.0.0). We used the Mann–Kendall test for trend analysis of haplotype prevalence over time in northern and eastern Uganda. The trend by sampling site was analyzed using the Bayesian model with default priors and number of MCMC draws as previously described69. To evaluate genotype–phenotype associations, we utilized all samples with lumefantrine, DHA or mefloquine IC50 values or RSA data that underwent px1 and k13 sequencing. Wilcoxon rank-sum test (independent sample sets, two-sided) was used to compare phenotype scores (IC50 values or RSA survival rates) between px1 haplotypes (PIN and LMD) with stratifications by different variables including K13 wild-type, 675V and C469Y alleles, region and year. To assess a potential contribution of confounders in drug susceptibility–px1 haplotype associations, a linear mixed-effects regression model was fitted with site, parasitemia, COI, year of sample collection and K13 mutation status as covariates. For each drug, we calculated the estimated marginal means of IC50 and its confidence interval, P values and Cohen’s d effect size. Values of P < 0.05 indicated statistically significant differences. We estimated both marginal and conditional R2 values as previously described70. Marginal R2 was used to estimate the proportion of variance explained solely by the px1 PIN haplotype. By contrast, conditional R2 was computed to represent the total variance explained by the entire model. We also evaluated the potential synergy between the px1 haplotype and K13 mutations using the mixed-effects model. Each drug assay was modeled with an interaction term between px1 and k13 treated as fixed effects, allowing us to determine whether the phenotypic effect of the PIN haplotype was modified by the K13 C469Y and A675V backgrounds. Complete-linkage hierarchical clustering method was used to cluster haplotypes flanking the px1 gene based on the genotype matrix in monogenomic samples. Inclusion and ethics statement This study arose from a long-standing scientific collaboration among the University of California, San Francisco (UCSF), Brown University, The University of North Carolina at Chapel Hill (UNC) and the Infectious Diseases Research Collaboration (IDRC) in Uganda. The work leveraged biological samples collected through ongoing, health facility-based molecular and parasitological malaria surveillance activities in Uganda, together with archived samples from prior studies. Surveillance activities were led by Ugandan investigators and conducted at up to 16 public health facilities across eastern and northern Uganda, selected based on established clinical and laboratory infrastructure. Ugandan study team members played central roles in conducting the study and led key components of the work in collaboration with UCSF investigators, including participant enrollment, sample collection, ex vivo drug susceptibility assays and targeted sequencing using molecular inversion probes developed by the Brown University team and Sanger methods. WGS was performed at UNC, and downstream genomic and statistical analyses were conducted by the Brown University team. Targeted ONT long-read sequencing of historical samples was jointly planned and led by Brown University and UCSF investigators. All collaborating institutions jointly agreed on data ownership, intellectual property and authorship in advance of the research. All studies were conducted in accordance with local and international ethical standards. Samples were collected under approved protocols with informed consent, including consent for future use of biological specimens. Genomic and phenotypic data were analyzed in de-identified form, with access restricted to authorized study personnel. The use of archived and prospectively collected samples was designed to maximize scientific value while minimizing additional risk to participants, and no new human participant recruitment was undertaken specifically for the genomic analyses reported here. Findings from previous studies conducted in the same regions were used to inform sample selection for WGS and to contextualize the impact of the px1 haplotype identified in this work. Relevant prior studies are cited accordingly. Ethics For all the molecular and parasitological studies, consent for future use of biological samples was given for all samples and ethical approval was obtained from the Makerere University Research and Ethics Committee, the Uganda National Council for Science and Technology, and the University of California, San Francisco, Human Research Protection Program (approval numbers: 24-41284, SBS-2023-497, HS4481ES, 16-19084 and 10-03144). Reporting summary Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article. Data availability The WGS data for 158 P. falciparum isolates generated in this study have been deposited in the NCBI Sequence Read Archive under BioProject accession PRJNA1298911. Targeted long-read amplicon sequencing data used for the spatiotemporal analysis of the px1 PIN haplotype have been deposited in the same repository under BioProject accession PRJNA1450510. PCR primers used in this work are available in Supplementary Table 5. The ex vivo data has been previously published as source data in ref. 22. The publicly available dataset used in this paper was previously published by the MalariaGEN Consortium46. Source data are provided with this paper. Code availability All the codes used for the genomic and statistical analyses are available via Zenodo at https://doi.org/10.5281/zenodo.21283869 (ref. 71). All codes are present in this repository to replicate the results. References Global Malaria Programme. World malaria report 2025 (World Health Organization, 2025); https://www.who.int/teams/global-malaria-programme/reports/world-malaria-report-2025 Shibeshi, W., Alemkere, G., Mulu, A. & Engidawork, E. Efficacy and safety of artemisinin-based combination therapies for the treatment of uncomplicated malaria in pediatrics: a systematic review and meta-analysis. BMC Infect. Dis. 21, 326 (2021). Makanga, M. & Krudsood, S. The clinical efficacy of artemether/lumefantrine (Coartem®). Malar. J. 8 (Suppl. 1), S5 (2009). Balint, G. A. Artemisinin and its derivatives: an important new class of antimalarial agents. Pharmacol. Ther. 90, 261–265 (2001). Dondorp, A. M. et al. Artemisinin resistance in Plasmodium falciparum malaria. N. Engl. J. Med. 361, 455–467 (2009). Witkowski, B. et al. Novel phenotypic assays for the detection of artemisinin-resistant Plasmodium falciparum malaria in Cambodia: In-vitro and ex-vivo drug-response studies. Lancet Infect. Dis. 13, 1043–1049 (2013). Amaratunga, C. et al. Dihydroartemisinin–piperaquine resistance in Plasmodium falciparum malaria in Cambodia: a multisite prospective cohort study. Lancet Infect. Dis. 16, 357–365 (2016). Phyo, A. P. et al. Declining efficacy of artemisinin combination therapy against P. falciparum malaria on the Thai–Myanmar border (2003–2013): the role of parasite genetic factors. Clin. Infect. Dis. 63, 784–791 (2016). Mok, S. et al. Population transcriptomics of human malaria parasites reveals the mechanism of artemisinin resistance. Science 347, 431–435 (2015). Rocamora, F. et al. Oxidative stress and protein damage responses mediate artemisinin resistance in malaria parasites. PLoS Pathog. 14, e1006930 (2018). Rahman, A. et al. Artemisinin-resistant Plasmodium falciparum Kelch13 mutant proteins display reduced heme-binding affinity and decreased artemisinin activation. Commun. Biol. 7, 1499 (2024). Birnbaum, J. et al. A Kelch13-defined endocytosis pathway mediates artemisinin resistance in malaria parasites. Science 367, 51–59 (2020). Rosenthal, P. J., Asua, V. & Conrad, M. D. Emergence, transmission dynamics and mechanisms of artemisinin partial resistance in malaria parasites in Africa. Nat. Rev. Microbiol. 22, 373–384 (2024). Conrad, M. D. et al. Evolution of partial resistance to artemisinins in malaria parasites in Uganda. N. Engl. J. Med. 389, 722–732 (2023). Uwimana, A. et al. Emergence and clonal expansion of in vitro artemisinin-resistant Plasmodium falciparum kelch13 R561H mutant parasites in Rwanda. Nat. Med. 26, 1602–1608 (2020). Wernsman Young, N. et al. High frequency of artemisinin partial resistance mutations in the great lakes region revealed through rapid pooled deep sequencing. J. Infect. Dis. 231, 269–280 (2025). Juliano, J. J. et al. Prevalence of mutations associated with artemisinin partial resistance and sulfadoxine-pyrimethamine resistance in 13 regions in Tanzania in 2021: a cross-sectional survey. Lancet Microbe 5, 100920 (2024). Mihreteab, S. et al. Increasing prevalence of artemisinin-resistant HRP2-negative malaria in Eritrea. N. Engl. J. Med. 389, 1191–1202 (2023). Fola, A. A. et al. Plasmodium falciparum resistant to artemisinin and diagnostics have emerged in Ethiopia. Nat. Microbiol. 8, 1911–1919 (2023). Ishengoma, D. S. et al. Evidence of artemisinin partial resistance in northwestern Tanzania: clinical and molecular markers of resistance. Lancet Infect. Dis. 24, 1225–1233 (2024). Tumwebaze, P. K. et al. Decreased susceptibility of Plasmodium falciparum to both dihydroartemisinin and lumefantrine in northern Uganda. Nat. Commun. 13, 6353 (2022). Okitwi, M. et al. Changes in susceptibility of Plasmodium falciparum to antimalarial drugs in Uganda over time: 2019–2024. Nat. Commun. 16, 7353 (2025). Venkatesan, M. et al. Polymorphisms in Plasmodium falciparum chloroquine resistance transporter and multidrug resistance 1 genes: parasite risk factors that affect treatment outcomes for P. falciparum malaria after artemether–lumefantrine and artesunate–amodiaquine. Am. J. Trop. Med. Hyg. 91, 833–843 (2014). Fançony, C. et al. Artemether–lumefantrine treatment selects Plasmodium falciparum multidrug resistance 1 (pfmdr1) increased copy number among African malaria infections. J. Infect. Dis. 231, e1119–e1128 (2025). Tumwebaze, P. K. et al. Drug susceptibility of Plasmodium falciparum in eastern Uganda: a longitudinal phenotypic and genotypic study. Lancet Microbe 2, e441–e449 (2021). Mungthin, M. et al. Association between the pfmdr1 gene and in vitro artemether and lumefantrine sensitivity in Thai isolates of Plasmodium falciparum. Am. J. Trop. Med. Hyg. 83, 1005–1009 (2010). Kubota, R., Ishino, T., Iwanaga, S. & Shinzawa, N. Evaluation of the effect of gene duplication by genome editing on drug resistance in Plasmodium falciparum. Front. Cell. Infect. Microbiol. 12, 915656 (2022). Laird, V. R. et al. Plasmodium falciparum multidrug resistance 1 gene polymorphisms associated with outcomes after anti-malarial treatment. Malar. J. 24, 186 (2025). Price, R. N. et al. Molecular and pharmacological determinants of the therapeutic response to artemether–lumefantrine in multidrug-resistant Plasmodium falciparum malaria. Clin. Infect. Dis. 42, 1570–1577 (2006). Rasmussen, S. A. et al. Changing antimalarial drug sensitivities in Uganda. Antimicrob. Agents Chemother. 61, e01516-17 (2017). Ebong, C. et al. Efficacy and safety of artemether–lumefantrine and dihydroartemisinin–piperaquine for the treatment of uncomplicated Plasmodium falciparum malaria and prevalence of molecular markers associated with artemisinin and partner drug resistance in Uganda. Malar. J. 20, 484 (2021). Kamya, M. R. et al. Efficacies of artemether–lumefantrine, artesunate–amodiaquine, dihydroartemisininpiperaquine and artesunate–pyronaridine for the treatment of uncomplicated Plasmodium falciparum malaria in children in Uganda: a randomized, open-label phase IV clinical trial. Social Science Research Network https://doi.org/10.2139/ssrn.5226564 (2025). Awor, P. et al. Indigenous emergence and spread of kelch13 C469Y artemisinin-resistant Plasmodium falciparum in Uganda. Antimicrob. Agents Chemother. 68, e0165923 (2024). Asua, V. et al. Changing prevalence of potential mediators of aminoquinoline, antifolate, and artemisinin resistance across Uganda. J. Infect. Dis. 223, 985–994 (2021). Niaré, K., Greenhouse, B. & Bailey, J. A. An optimized GATK4 pipeline for Plasmodium falciparum whole genome sequencing variant calling and analysis. Malar. J. 22, 207 (2023). Amambua-Ngwa, A. et al. Chloroquine resistance evolution in Plasmodium falciparum is mediated by the putative amino acid transporter AAT1. Nat. Microbiol. 8, 1213–1226 (2023). Mukherjee, A. et al. A phosphoinositide-binding protein acts in the trafficking pathway of hemoglobin in the malaria parasite Plasmodium falciparum. mBio 13, e0323921 (2022). Zhu, L. et al. The origins of malaria artemisinin resistance defined by a genetic and transcriptomic background. Nat. Commun. 9, 5158 (2018). Ravenhall, M. et al. Characterizing the impact of sustained sulfadoxine/pyrimethamine use upon the Plasmodium falciparum population in Malawi. Malar. J. 15, 575 (2016). Ocholla, H. et al. Whole-genome scans provide evidence of adaptive evolution in Malawian Plasmodium falciparum isolates. J. Infect. Dis. 210, 1991–2000 (2014). Hathaway, N. J. et al. Interchromosomal segmental duplication drives translocation and loss of P. falciparum histidine-rich protein 3. eLife 13, RP95334. Tang, W. K. et al. Multistage protective anti-CelTOS monoclonal antibodies with cross-species sterile protection against malaria. Nat. Commun. 15, 7487 (2024). Lefebvre, M. J. M. et al. Population genomic evidence of adaptive response during the invasion history of Plasmodium falciparum in the Americas. Mol. Biol. Evol. 40, msad082 (2023). Shah, Z. et al. Whole-genome analysis of Malawian Plasmodium falciparum isolates identifies possible targets of allele-specific immunity to clinical malaria. PLoS Genet. 17, e1009576 (2021). Mobegi, V. A. et al. Genome-wide analysis of selection on the malaria parasite Plasmodium falciparum in West African populations of differing infection endemicity. Mol. Biol. Evol. 31, 1490–1499 (2014). MalariaGEN et al. An open dataset of Plasmodium falciparum genome variation in 7,000 worldwide samples. Wellcome Open Res. 6, 42 (2021). van Schalkwyk, D. A. et al. Treatment failure in a UK malaria patient harboring genetically variant Plasmodium falciparum from Uganda with reduced in vitro susceptibility to artemisinin and lumefantrine. Clin. Infect. Dis. 78, 445–452 (2024). Sutherland, C. J. et al. Pfk13-independent treatment failure in four imported cases of Plasmodium falciparum malaria treated with artemether–lumefantrine in the United Kingdom. Antimicrob. Agents Chemother. 61, e02382 (2017). Oliveira, R. et al. Artemether–lumefantrine resistant falciparum malaria imported from Uganda. J. Travel Med. 32, taaf013 (2025). Henden, L., Lee, S., Mueller, I., Barry, A. & Bahlo, M. Identity-by-descent analyses for measuring population dynamics and selection in recombining pathogens. PLoS Genet. 14, e1007279 (2018). Espié, E. et al. Efficacy of fixed-dose combination artesunate–amodiaquine versus artemether–lumefantrine for uncomplicated childhood Plasmodium falciparum malaria in Democratic Republic of Congo: a randomized non-inferiority trial. Malar. J. 11, 174 (2012). Ferrari, G. et al. An operational comparative study of quinine and artesunate for the treatment of severe malaria in hospitals and health centres in the Democratic Republic of Congo: the MATIAS study. Malar. J. 14, 226 (2015). Asua, V. et al. Changing molecular markers of antimalarial drug sensitivity across Uganda. Antimicrob. Agents Chemother. 63, e01818-18 (2019). Yeka, A. et al. Artemisinin versus nonartemisinin combination therapy for uncomplicated malaria: randomized clinical trials from four sites in Uganda. PLoS Med. 2, e190 (2005). Arinaitwe, E. et al. Artemether–lumefantrine versus dihydroartemisinin–piperaquine for falciparum malaria: a longitudinal, randomized trial in young Ugandan children. Clin. Infect. Dis. 49, 1629–1637 (2009). Wanzira, H. et al. Longitudinal outcomes in a cohort of Ugandan children randomized to artemether–lumefantrine versus dihydroartemisinin–piperaquine for the treatment of malaria. Clin. Infect. Dis. 59, 509–516 (2014). Oyola, S. O. et al. Whole genome sequencing of Plasmodium falciparum from dried blood spots using selective whole genome amplification. Malar. J. 15, 597 (2016). Miles, A. et al. Indels, structural variation, and recombination drive genomic diversity in Plasmodium falciparum. Genome Res. 26, 1288–1299 (2016). Chang, H.-H. et al. THE REAL McCOIL: a method for the concurrent estimation of the complexity of infection and SNP allele frequency for malaria parasites. PLoS Comput. Biol. 13, e1005348 (2017). Gautier, M. & Vitalis, R. rehh: an R package to detect footprints of selection in genome-wide SNP data from haplotype structure. Bioinformatics 28, 1176–1177 (2012). Cook, D. E. & Andersen, E. C. VCF-kit: assorted utilities for the variant call format. Bioinformatics 33, 1581–1582 (2017). Auton, A. & McVean, G. Recombination rate estimation in the presence of hotspots. Genome Res. 17, 1219–1227 (2007). McVean, G. A. T. et al. The fine-scale structure of recombination rate variation in the human genome. Science 304, 581–584 (2004). Niaré, K. et al. Highly multiplexed molecular inversion probe panel in Plasmodium falciparum targeting common SNPs approximates whole-genome sequencing assessments for selection and relatedness. Front. Genet. 16, 1526049 (2025). Schaffner, S. F., Taylor, A. R., Wong, W., Wirth, D. F. & Neafsey, D. E. hmmIBD: software to infer pairwise identity by descent between haploid genotypes. Malar. J. 17, 196 (2018). Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007). Verity, R. et al. The impact of antimalarial resistance on the genetic structure of Plasmodium falciparum in the DRC. Nat. Commun. 11, 2107 (2020). Aydemir, O. et al. Drug-resistance and population structure of Plasmodium falciparum across the Democratic Republic of Congo using high-throughput molecular inversion probes. J. Infect. Dis. 218, 946–955 (2018). Meier-Scherling, C. P. G. et al. Selection of Plasmodium falciparum kelch13 mutations in Uganda in comparison with southeast Asia: a modelling study. Lancet Microbe 6, 101027 (2025). Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods Ecol. Evol. 4, 133–142 (2013). Niare, K. Karaniare/PfWGS-DRscanner: genome-wide scans for locus association with antimalarial drug responses in Plasmodium falciparum (version v1.0.0) [computer software]. Zenodo https://doi.org/10.5281/zenodo.21283869 (2026). Acknowledgements We are grateful to the study participants, caregivers and research teams at the Infectious Diseases Research Collaboration (IDRC) in Uganda. The funders had no role in study design, data collection, analysis, interpretation or the decision to publish. We thank S. Garg for performing molecular assays; P. K. Tumwebaze, O. Byaruhanga, E. Muhanguzi, S. Opio, I. Tibagambirwa, P. Angutoko, J. Asiimwe, Y. Tarema and T. Katairo for performing ex vivo assays; G. Dorsey for sharing metadata for historical samples; and C. Meier-Scherling for sharing codes for the Bayesian analysis. We thank L. Checkley, D. Shoue and D. Gagnon for helping with the in vitro IC50 experiments. Funding This work was funded by the National Institutes of Health (NIH)/National Institute of Allergy and Infectious Diseases (NIAID) (R01AI173557 to M.D.C. and K24AI134990 to J.J.J.). Sample collection and phenotyping was funded by NIAID (R01AI075045, U19AI089674, R01AI117001 and R01AI139179); the Medicines for Malaria Venture (RD/15/0001); and the Gates Foundation (INV-035751). In vitro drug evaluation was funded by NIAID (R01AI189911-01). Author information Authors and Affiliations Contributions K.N., M.D.C., J.J.J., J.A.B. and P.J.R. designed and conducted the study. K.N. analyzed data and interpreted results. M.D.C., J.J.J. and J.A.B. supervised the analysis and interpretation of the data. M.D.C. obtained primary funding and led the project. J.J.J. and J.M.S. performed the library preparation and rebalancing for WGS. K.N. conducted WGS quality control and coordinated library rebalancing. K.N. and B.T. performed Nanopore sequencing. M.T., O.K., V.A. and J.L. performed genotyping assays. J.M. performed the scan for large structural variations in chr. 7. M.O. and S.O. performed ex vivo assays. S.L.N., V.A. and A.Y. led lab activities and sample collections in Uganda. A.M. led the in vitro studies evaluating drug resistance phenotypes associated with px1 disruption, wrote the relevant section and analyzed the results in M.T.F.’s laboratory. D.R. supervised the px1 disruption experiments. K.N. wrote the primary draft of the paper. K.N., M.D.C., J.J.J., J.A.B. and P.J.R. critically reviewed and revised the paper. All authors revised the paper. Corresponding authors Ethics declarations Competing interests The authors declare no competing interests. Peer review Peer review information Nature Medicine thanks Olivo Miotto and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available. Primary Handling Editor: Lia Parkin, in collaboration with the Nature Medicine team. Additional information Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Extended data Extended Data Fig. 1 Signal of selection around K13 mutations C469Y and A675V. A) Decay of the extended haplotype homozygosity around each of the two markers. Mutant (red) represents 469Y or 675 V. Analysis included dominant alleles with minor allele frequency ≥2%. Numbers of unique P. falciparum isolates analyzed are indicated (n = 106 and n = 102) B) Visualization and hierarchical clustering of flanking haplotypes based on genotypes. Region size was based on the maximal extended haplotype heterozygosity region found in (A) and measured ~200 kb. Positions of key mutations are highlighted by red and blue vertical lines. Numbers of unique P. falciparum isolates analyzed are indicated (n = 11, n = 15). Extended Data Fig. 2 Selection signals by region. Two-side P values from IsoRelate’s chi-square-distributed test statistics for IBD (iR) were corrected for multiple testing using FDR (Benjamini–Hochberg method). Significance threshold of -Log10(FDR-corrected P value) of 5 was indicated by the red dotted horizontal line. Red dashed rectangles demarcate candidate sweep and k13 peak. The analysis included SNPs with minor allele frequency ≥2% and all samples (n = 157 unique P. falciparum isolates). Genes corresponding to peaks are indicated. Sample size by region is indicated in the figures (n = 116 and n = 41). Extended Data Fig. 3 Selection signal by K13 mutation status. Two-sided P from IsoRelate’s chi-square-distributed test statistics for IBD (iR) values were corrected for multiple testing using FDR (Benjamini–Hochberg method). Significance threshold of -Log10(FDR-corrected P value) of 5 was indicated by the red dotted horizontal line. The analysis included SNPs with minor allele frequency ≥2% and all samples (n = 157 unique P. falciparum isolates). Genes corresponding to peaks are indicated. Red dashed rectangles demarcate candidate sweep and k13 peak. Sample size by K13 mutation status (n = 66 and n = 91)is indicated in the figures. Extended Data Fig. 4 Unfiltered iHS signals. Two-sided P values from rehh’s iHS test were corrected for multiple testing using FDR (Benjamini–Hochberg method). Significance threshold of -Log10(FDR-corrected P value) of 5 was indicated by the red dotted horizontal line. Red dashed rectangles demarcate candidate sweep and k13 peak. The analysis included SNPs with minor allele frequency ≥2% and mono-genomic samples (n = 118 unique P. falciparum isolates). epf1: exported protein 1; celtos: cell-traversal protein for ookinetes and sporozoites; pfa55-14: asparagine-rich antigen. Extended Data Fig. 5 Genetic recombination map of the candidate sweep region. This analysis included SNPs with minor allele frequency ≥2% and mono-genomic samples (unique P. falciparum isolates) in our study (n = 118) and MalariaGEN Pf7 (n = 96). The px1 gene had lowest ρ values, ranging between 0.78-1.97. PF3D7_0721000 had high ρ values, ranging between 0.99-21.60. Shaded areas correspond to key genes in the candidate sweep with significant iHS. Extended Data Fig. 6 Linkage desequilibrium scan for the PIN haplotype. A) Profile of the Linkage Desequilibrium (LD) between the PIN haplotype and SNPs across the genome (n = 157 individual samples). LD R2 between the PIN haplotype and SNPs outside the candidate sweep are low and below 0.2 (indicated by the red dashed line). B) LD decay of the SNPs located in the candidate selective sweep region versus the rest of chromosome 7. Extended Data Fig. 7 Extended haplotype homozygosity profile for PX1 top variants (SNPs and indels). WT: wild-type. The analysis included SNPs and indels with minor allele frequency ≥2% and mono-genomic samples (n = 118). The two deletions (delV16Z80-1685 and delN1289-1294) were rare in global populations and belonged in the PIN haplotype (1222 P, 1701I, 1705 N and deletions delV1680-N1685 and delN811-Y822). Supplementary information Supplementary Information (download PDF ) Supplementary Tables 1–5 and Figs. 1–14. Supplementary Data 1 (download XLSX ) Statistical source data for Supplementary Fig. 8. Source data Source Data Figs. 4 and 5 and Extended Data Tables 1 and 2 (download XLSX ) Statistical source data for Figs. 4 and 5 and Extended Data Tables 1 and 2. Rights and permissions Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/. About this article Cite this article Niaré, K., Tafesse, B., Treat, M. et al. Emergence and spread of Plasmodium falciparum PX1 polymorphisms associated with decreased susceptibility to antimalarials in Uganda. Nat Med (2026). https://doi.org/10.1038/s41591-026-04590-5 Received: Accepted: Published: Version of record: DOI: https://doi.org/10.1038/s41591-026-04590-5