Blood infected with Plasmodium chabaudi AS parasites was kept in liquid nitrogen and thawed on ice upon use. Infected blood was inoculated into six-to-eight week old CD1 male mice and parasite growth was observed during time by inspection of Giemsa-stained blood smears by optical microscopy. The parasite clones used in this study were ( Figure 1 ): AS-ATN, which displays low level resistance to artesunate (Afonso et al., 2006); AS-ATNMF1, derived from AS-ATN through in vivo selection with ATN + MF administered simultaneously (Rodrigues et al., 2010); AS-ATN27P, a parasite that was obtained after twenty-seven consecutive sub-inoculations of AS-ATN into untreated animals (Rodrigues et al., 2010); clones AS-SENS (drug-sensitive) and AS-50SP (resistant to sulphadoxine and pyrimethamine).

All animal work was conducted according to relevant national and international guidelines: in Portugal, after approval by the Ethics Committee of the Instituto de Higiene e Medicina Tropical of Lisbon, Portugal, under PARECER 2/2006 from August 1^st, 2006; in the United Kingdom, in compliance with the United Kingdom Animals (Scientific Procedures) Act 1986.

Parasite DNA samples were extracted from blood by two different methods: (a) Fresh blood samples – Infected mice were exsanguinated by brachial artery puncture. Blood was processed, and DNA extracted as previously described (Grech et al., 2002). Host white blood cells were removed by CF11 cellulose (Whatman) and Plasmodipur filters (Eurodiagnostica); and (b) Blood samples preserved in filter paper – Blood samples were collected by a small puncture on the mouse tail vein. A drop was preserved in Whatman N° 4 filter paper and was extracted by boiling in Chelex-100 (Plowe et al., 1995).

DNA obtained from fresh blood samples was processed using standard methods according to the manufacturer’s recommendation (Illumina, 2010). AS-ATNMF1 genome sequencing was performed using the Illumina^® (Solexa) platform, with 50 bp reads, using a paired-end read approach at the Genepool facilities at the University of Edinburgh, UK. Sequences were aligned against an isogenic reference genome [AS-WTSI, curated by the Pathogen Sequencing Group, Welcome Trust Sanger Institute and available at Sanger Institute webpage (Sanger Institute, 2009). The updated reference genome is available at https://plasmodb.org/plasmo/app/record/dataset/DS_5dab6e4a9f, with the following accession number: GCA_900002335.2] by two software packages, Mapping and Assembly with Quality (MAQ) (Li et al., 2008) and Sequence Search and Alignment by Hashing Algorithm (SSAHA2), using default settings (Ning, 2001).

SNP detection was performed using in-built algorithms, with a read depth of 3 set as the minimum threshold for SNP detection. In order to filter out SNPs arising in the sensitive progenitor AS-SENS, the list of SNPs obtained for AS-ATNMF1 was compared against a list previously obtained for AS-SENS (Hunt et al., 2010). The two filtered lists (for MAQ and SSAHA2) were then combined. Small indels (< 3bp) were detected using the internal algorithm of SSAHA2 only. The list of small indels was filtered against a similarly obtained list for AS-SENS (Hunt et al., 2010).

Larger indels (≥ 3bp) and copy number variation (CNV) detection was performed using SSAHA2 as previously described (Hunt et al., 2010). Two approaches were used to indicate areas of significantly decreased or increased coverage (Hunt et al., 2010). Since AS-ATNMF1 was sequenced using paired-end reads but AS-SENS was originally sequenced using single-end reads (Hunt et al., 2010), this approach could not be adopted using this clone, due to the excessive variation introduced by the different sequencing strategies adopted. Instead, a clone selected for sulphadoxine resistance (AS-50SP) and re-sequenced using paired-end reads (Martinelli et al., 2011) was selected to act as a “filter”, with the caveat that false positives and negatives are likely to appear. Due to a complete sequencing and analysis of other clones in the AS lineage (Hunt et al., 2010), we can rely on the data previously obtained for resolving any uncertainties and discovering false negatives. The positions of all the mutations described here are in accordance to an older version of the P. chabaudi annotated genome available on Welcome Trust Sanger Institute website as indicated above.

All SNPs were treated as potentially genuine mutations pending verification by di-deoxy sequencing. In the case of small and large indels, the vast majority is still expected to be false positive calls (Hunt et al., 2010; Martinelli et al., 2011; Kinga Modrzynska et al., 2012) and limited di-deoxy sequencing was performed based on the potential biological role attributed to the proteins coded by the genes where mutations were present or, in the case of intergenic mutations, the role of genes immediately close (either upstream or downstream) to the identified mutations.

A sub-set of mutations was chosen for verification by di-deoxy sequencing. For that purpose, a region of about a 1.000 bp flanking each mutation was used for designing oligonucleotide primers ( Supplementary Table 1 ) and the fragments were amplified by PCR. DNA samples were obtained from dried blood spots and 1 μl of DNA was used as template in 50 μl reactions. The other reagents were added to 10x Green GoTaqR Flexi Buffer (Promega) to a final concentration of 1.5 mM MgCl₂, 0.2 mM deoxynucleoside triphosphate (each), 0.2 pmol/μl forward and reverse primers, and 1.25 U of Go Taq Flexi DNA Polymerase (Promega). Cycling conditions were 93°C for 3 minutes, followed by 35 repeats of denaturation at 93°C for 30 seconds, annealing for 45 seconds, and elongation at 72°C for 1 minute. A final elongation step was performed at 72°C for 10 minutes. PCR products thus obtained were sequenced using the commercial services of STABVIDA Laboratories (Oeiras, Portugal). Chromatograms were inspected using Chromas 2.33 software (Technelysium).

Predicted protein sequence and function was retrieved from PlasmoDB (Aurrecoechea et al., 2009) for all genes where a mutation was confirmed. Additionally, in cases where the predicted protein sequence did not contain a functional annotation in PlasmoDB, it was inspected for conserved domains and functional sites using InterProScan, a tool that combines different protein signature recognition methods into one resource (Zdobnov and Apweiler, 2001; Quevillon et al., 2005; Mulder and Apweiler, 2007).

The wild-type 3D structure of 26S proteasome subunit Rpn2 of P. falciparum (PF3D7_1466300) was built by threading approach through the I-Tasser server (Zhang, 2008; Yang et al., 2015). After modeling, the structure underwent refinement through the KobaMin server (Rodrigues et al., 2012), utilizing a knowledge-based potential of mean force for minimization (Summa and Levitt, 2007). Hydrogen atoms were added using the MolProbity server (Hintze et al., 2016). The mutant 3D structure E694K was derived from the wild-type structure by substituting the side chain of residue 694, using GROMACS v.5.1.2 (Van Der Spoel et al., 2005; Abraham et al., 2015). For the molecular dynamics (MD) simulations, the systems were solvated in a cubic box (solute box distance of 1.2 nm) and charge neutrality was achieved by introducing Na⁺ and Cl⁻ ions at a concentration of 0.15 mol/L. All simulations were conducted in GROMACS v.5.1.2, utilizing the AMBER FF99SB-ILDN force field (Cornell et al., 1995). Water molecules were described by the TIP3P model (Jorgensen et al., 1983).

The systems were prepared as mentioned above and then submitted to 1.000 steps of energy minimization using the steepest descent method, with harmonic restraints applied to all heavy atoms, to remove highly repulsive contacts or geometric strain. The minimizations were completed when the tolerance of 1000 kJ/mol was no longer exceeded. Subsequently, the following steps were performed to equilibrate the systems: (a) MD simulations using the NVT ensemble with harmonic restraints applied to all of the heavy atoms (1 ns); (b) MD simulations using the NPT (isothermal-isobaric) ensemble with harmonic restraints on all of the heavy atoms (1 ns); (c) MD simulations using the NPT ensemble without any restraints (1 ns). After these preparatory steps, MD trajectories of 100 ns for each system were generated as the production phase, in duplicates. The simulations were performed through periodic boundary conditions, using a cut-off of 1.0 nm, at a temperature of 300 K and 1 atm pressure. Electrostatics interactions were evaluated using the particle mesh Ewald algorithm (Leach, 2001). All RMSF and β-factor were calculated using the GROMACS gmx-toolbox (https://manual.gromacs.org/current/user-guide/cmdline.html). Visual Molecular Dynamics program (VMD) (Humphrey et al., 1996) was employed for visualizing MD trajectories, calculating residue-residue distances and rendering the molecular images. Additional molecular images were created using Python Molecular Viewer v.1.5.6 (Morris et al., 2009), while the graphs were generated using XMgrace software (http://plasma-gate.weizmann.ac.il/Grace/).

The genome of the resistant clone AS-ATNMF1 was sequenced and 23,732,208 reads were produced with 21,816,439 reads (approximately 91% of all reads) mapped to the reference genome by SSAHA2. The average coverage of the whole genome was ~50-fold and 43% of the genome had at least 40-fold coverage, 92% was covered by at least 10 reads, whereas less than 0.7% of genome was not covered by any read. A subset of mutations was chosen for verification by di-deoxy sequencing, as will be described in detail below. It is worth noting that all the mutations analyzed by di-deoxy sequencing in AS-ATNMF1 and AS-ATN were also investigated in the parasites sub-inoculated twenty-seven times in untreated animals, AS-ATN27P. Unsurprisingly, AS-ATN27P showed the same genotype as AS-ATN.

The 26S proteasome comprising a 20S proteolytic core and two 19S regulatory particles (19S-RP) plays a pivotal role in the proteolysis of ubiquitylated proteins in the cell. The 19S-RP subunits Rpn1 and Rpn2 are responsible for recognizing ubiquitylated proteins and priming them for proteolysis by the 20S core. The Rpn2 subunit is assembled from a rod-like N-terminal domain, a toroidal proteasome/cyclosome (PC) domain, and a globular C-terminal domain (He et al., 2012) ( Figure 2A ). The PC repeat domain is a concentric arrangement of α helices: two antiparallel α helices (axial helices), wrapped in a double layer of α helices, consisting of the inner and outer helices ( Figure 2B ). We built the wild-type structure of the P. falciparum Rpn2 subunit based on P. falciparum primary sequence (Uniprot ID Q8IKH3), through a fold recognition (threading) approach. Analysis of the model’s Ramachandran plot ( Supplementary Figure 1 ) revealed that 87.2% of the residues occupy the most favorable regions which is close to the ideal of 90% (Laskowski et al., 1993), suggesting that the built model is acceptable. Based on wild-type structure modeled, the E694K mutant was computationally built. This mutation on P. falciparum corresponds to the E738K mutation on the PCHAS_133430 gene of P. chabaudi ( Supplementary Figure 2 ). Considering the Rpn2 wild-type tertiary structure, we observe that E694 is surrounded by positively charged residues, such as H689, H693, K695 and R698. Hence, the substitution of glutamate with lysine (E694K) leads to a significant alteration in the residue’s standard interaction network.

To gain a deeper understanding of the dynamic behavior of the Rpn2 subunit in both the wild-type (WT) and E694K mutant (MT) models, we conducted two independent molecular dynamics simulations (MDs) (100 ns each, in duplicate) of the WT and MT structures. Figures 2C, D present snapshots of the most prevalent E694 network interactions (WT-MD) and K694 network interactions (MT-MD), respectively, as observed during the MD simulations. Hydrogen bond analysis of residues interaction, in WT and MT-MD simulations, are available as Supplementary Data ( Supplementary Figure 3 ). In WT-MD simulations, E694 network interactions were extended from inner to axial helices, through salt bridges ( Figure 2C ). In MT-MD simulations, K694 interacts with outer helices residues/water molecules, and there is a rearrangement of neighborhood residues (conformations and interactions) ( Figure 2D ). The residue R698 underwent a conformational change, deviating from its WT state by disrupting interactions with inner helix residues and establishing new associations with axial helix residues ( Figure 2D ). This shift in conformation affects the interplay between the axial-inner and inner-outer α helices within the proteasome/cyclosome (PC) domain, subsequently influencing the toroidal organization of the PC domain.