Systematic in vitro evolution in Plasmodium falciparum reveals key determinants of drug resistance

Abstract

Surveillance of drug resistance and the discovery of novel targets-key objectives in the fight against malaria-rely on identifying resistance-conferring mutations in Plasmodium parasites. Current approaches, while successful, require laborious experimentation or large sample sizes. To elucidate shared determinants of antimalarial resistance that can empower in silico inference, we examined the genomes of 724 Plasmodium falciparum clones, each selected in vitro for resistance to one of 118 compounds. We identified 1448 variants in 128 recurrently mutated genes, including drivers of antimalarial multidrug resistance. In contrast to naturally occurring variants, those selected in vitro are more likely to be missense or frameshift, involve bulky substitutions, and occur in conserved, ordered protein domains. Collectively, our dataset reveals mutation features that predict drug resistance in eukaryotic pathogens.

Fig. 1.. Compound-selected single nucleotide variants (SNVs) and indels have different characteristics compared to naturally occurring variants.

(A-D) Comparison between the set of 1,448 compound-selected variants (SNVs and indels) in the core genome from this study (left axis) and a set of 17,613 control variants obtained by aligning Dd2 to the 3D7 reference genome. Axes have been normalized to display differences in proportions of different categories for compound-selected (left) and control (right) variants. (B) Reference base and alternate base for 1,141 unique core SNVs (coding and noncoding) compared to 15,793 control variants. (C) Analysis of amino acid changes for 640 core missense and 6,600 core missense control SNVs with arrows showing relative increases or decreases. (D) 905 evolved and 13,763 control protein coding variants were analyzed for location within a defined InterPro domain obtained from PlasmoDB v61. In all **** indicates P < 0.0001; *** < 0.001 and ** < 0.01 calculated using a Chi-square test.

Fig. 2.. Copy number variants (CNVs) are a frequent driver of antimalarial compound resistance.

(A) Manhattan plot showing the number of times each gene was amplified as part of an independent CNV in our dataset. Genes that are likely drivers of the selective advantage conferred by CNVs are annotated, including compound targets and multidrug resistance genes. (B) Heatmap visualization of amplifications containing pfmdr1 for evolved clones (pfmdr1 CNVs that independently involved in the same compound selection were omitted for brevity). Denoised log₂ gene copy ratios are plotted for each clone; a contiguous segment of high log₂ copy ratio (yellow) suggests an amplification including those genes. For CNVs identified as tandem duplications, more precise CNV boundaries are indicated by red markers. (C) Comparison of genomic classification distribution between 284 independent tandem duplication CNV breakpoints and all sites in the core genome. (D) AT content around tandem duplication breakpoints, averaged over the 284 independent tandem duplications. AT content for each relative position was computed as the proportion of bases that are A or T over all sequences beginning at that position with a sliding window of five bp. The plot is shown as a running average with a window of 100 bp.

Fig. 3.. Enriched genes mutated in compound resistance selections.

(A) Classification of evidence types supporting assignment of target or resistance mediator for 118 compounds (Data S5). “Other evidence” indicates the target/resistance gene was identified but its representation did not reach statistical significance (hypergeometric test); other experimental evidence was required for confirmation. “Ambiguous CNVs” were cases in which a CNV was found in a majority of selected clones but a clear target or resistance gene was not identified. (B) Enrichment P values (hypergeometric mean function) for genes that were recurrently mutated in selections for a given compound (Table S1). Structures of compounds that gave rise to SNV/indel mutations in select genes are shown in Fig. S7, including the Tanimoto chemical similarity score. Not all compounds are shown for pfatp4 due to undisclosed structure. (C) Condensed version of a network linking genes for which the same pair of mutations arose in distinct compounds, filtered to show gene pairs with at least ten shared variant pair-compound pair occurrences (full network shown in Fig. S6). In cases of more than one variant pair-compound pair existing between two genes, this multiplicity was encoded as edge weight. All seven digit gene IDs have PF3D7_ removed. Circles represent known proteins; diamonds, conserved proteins; squares, putative proteins. Node size is proportional to node degree and edge color maps to occurrence of disruptive mutations from low (light blue) to high (dark blue) after adjustment for non-missense mutations. (D) dN for 449 genes that had more than one non-singleton SNV in the compound-selected dataset plotted against dN among Pf6 samples (top) and dS vs. dN of 4,938 core genes in the Pf6 dataset of worldwide variation (bottom).

Fig. 4.. Compound susceptibility assays comparing compound-selected mutations and field variants in pfcarl, pfmdr1, and pfatp4.

(A) Heatmap showing mean IC₅₀ fold changes relative to 3D7-A10 from dose-response experiments for four laboratory lines and 15 clinical isolates from Uganda and Senegal with pfcarl variants against three compounds that select for PfCARL mutations. Pfcarl and pfat1 variants were confirmed by WGS; the left heatmap shows allele frequencies for the tested lines. For reference, fold change based on previously reported IC₅₀ values are shown for resistant lines with PfCARL substitutions, highlighted in orange. (B) Mean IC₅₀ fold changes relative to isogenic parent for five edited lines with pfmdr1 mutations. Parasite lines highlighted in orange have mutations identified from compound selections, while the rest occur in the field. (C) Mean IC₅₀ fold changes relative to Dd2–2D4 for eight clinical isolates from Senegal and a KAE609-pressured resistant mutant, 3D7-ATP4^T416N, highlighted in orange. Pfatp4 allele frequencies confirmed by WGS are displayed in the left heatmap. (D-F) Ligand-filled models of PfCARL and PfATP4 were obtained from AlphaFill; PfMDR1 homology model was constructed with SWISS-MODEL (43) using 7a69 (see Fig. S8 for details). Mutated residues in our compound-selected dataset (with the exception of PfCARL Ile-1139, which was found in (21)) are colored magenta; mutated residues in the field with global major allele frequency (gMAF) > 0.002 were obtained from (39) and are colored green, while those found in both sets are colored blue. Labels are colored orange and green, respectively, for compound-selected vs. field mutations phenotyped in (A-C). Other compound-selected mutations are also labeled in black.