Decoding efficacy and resistance space at a drug binding site

Domain: neuroscience
Published: 2026-02-04
Editorial outcome: accept

Abstract
Assessing and understanding the impacts of all possible mutations at a drug binding site remain challenging. Here we use multiplex oligo targeting for mutational profiling, and computational modelling, to decode efficacy and resistance space at the otherwise native binding site for an anti-trypanosomal proteasome inhibitor. We saturation-edited twenty codons in the Trypanosoma brucei proteasome and subject the resulting libraries to stepwise drug selection and codon variant scoring, yielding dose-response profiles for >100 resistance-conferring mutants. Codon variant scores are predictive of relative resistance observed using a bespoke set of mutants, while fitness profiling reveals otherwise extensive constraints on mutational fitness and resistance space. The resistance profile is predictive of routes to spontaneous drug resistance observed within ‘accessible’, single nucleotide mutational space, while in silico predictions are closely aligned with impacts on drug resistance observed in cellulo. Thus, multiplex oligo targeting facilitates assessment of all possible mutations at a drug binding site.

Introduction
Interactions between drugs and their targets impact efficacy and, when altered by mutation, can result in resistance1–3, undermining the durability of chemotherapeutics used against infectious diseases or in oncology. An improved understanding of the impacts of all possible mutations at drug binding sites would therefore facilitate drug discovery and resistance surveillance, but assessing these impacts remains challenging4–8. Insights into drug – target interactions typically emerge from a combination of empirical experimental observation in cellulo, structural studies, and computational modelling and prediction in silico. Currently, although powerful approaches for saturation mutagenesis, combined with multiplexed assays, are available, in cellulo data describing comprehensive sets of mutations at native drug bindingsites are often lacking9, and confidence in in silico predictions can suffer from unquantified limitations as a result.

Several CRISPR-based editing approaches have been used to assess drug resistance associated mutations in eukaryotic cells. However, these approaches employ Cas9 nucleases that can yield undesirable and imprecise recombinants, base editors that typically introduce <50% of all possible edits6,10, or prime editors that currently display limited efficiency4,5,7. CRISPR-based approaches are also typically constrained by the availability of suitable protospacer-adjacent motifs in the editing window and necessitate the introduction of further changes at the target site to prevent multiple rounds of editing. These approaches also present multiplexing challenges in diploid cells since different edits will likely be introduced in each allele, complicating genotype to phenotype assessments; although this latter challenge can be mitigated by using haploid cells11. Deep mutational scans, widely used for multiplex analysis of variant effects12, have also been used to assess drug resistance mechanisms. However, these approaches typically employ ectopic expression systems that may alter target protein abundance and/or expression patterns13.

Alternative approaches for the delivery of codon edits to native genes are available. Indeed, editing using multiplex genome engineering14, or oligo targeting15, simply require the delivery of short single-stranded oligonucleotides. These approaches are DNA cleavage-free, and CRISPR-independent, and can deliver the full range of possible base edits, albeit at relatively low frequency. These approaches have been described in Saccharomyces cerevisiae14, in embryonic stem cells16, and in the trypanosomatids, Trypanosoma brucei, T. cruzi and Leishmania15; parasitic protozoa that are transmitted by insects, causing neglected tropical diseases in humans, and a range of veterinary diseases. Indeed, we recently reported the development of oligo targeting for the assessment of drug resistance mechanisms in the trypanosomatids15.

In this work, we scale the oligo targeting approach for saturation-editing of the codons encoding a drug binding site in the T. brucei proteasome β5 subunit. The proteasome is a promising drug target in the trypanosomatids17–21, in malaria parasites22, and against multiple myeloma23,24. Cryo-electron microscopy structures of the Leishmania proteasome reveal a common binding site for both GSK349424521 and LXE40825, and anti-trypanosomal proteasome inhibitors are currently in clinical development against both visceral leishmaniasis, and T. cruzi, which causes Chagas disease. Multiplex oligo targeting, amplicon sequencing, and codon variant scoring, combined with computational modelling to predict impacts on ligand affinity and proteasome function, provide unprecedented insights into the relationships between a drug and a mutated, but otherwise native, drug target. The readouts reveal how each binding site residue, and all possible binding site mutations, impact drug efficacy and drug resistance.

Results
Mutagenesis at a drug binding site
The proteasome is a high-priority drug target and was selected for the development and testing of multiplexed oligo targeting in trypanosomes. The proteasome comprises four stacked heteroheptameric rings with the β subunits contributing six proteolytic active sites that reside in the central chamber. These sites are autocatalytically activated in the assembled structure via cleavage of N-terminal propeptides; at the glycine preceding the catalytic threonine in the T. brucei β5 subunit (Tb927.10.6080)26. Anti-trypanosomal proteasome inhibitors selectively target the chymotrypsin-like site in this subunit17,21,25.

Oligo targeting in T. brucei involves delivery of single-stranded oligodeoxynucleotides (ssODNs) of ~50 b in length, with higher efficiency editing achieved when using ‘reverse-strand’ ssODNs15; a single allele is typically edited in T. brucei, which is diploid. To target the chymotrypsin-like site in the β5 subunit, we considered a pooled ssODNs delivery protocol but noted the potential for introduction of multiple edits in the same cell14, which would allow for enrichment of bystander edits, complicating codon-based genotype to phenotype assessments. We, therefore, tested co-editing efficiency by delivering pairs of ssODNs to wild-type T. brucei cells, one designed to introduce a resistance-conferring edit (G98S-GGA-AGC), equivalent to G98S in T. cruzi15, and a second designed to introduce a synonymous edit either ~100 (C63C, TGC-TGT) or 200 bp (V31V, GTT-GTA) upstream. Although an allele replacement frequency of only 0.003 % was reported previously in wild-type T. brucei15, not readily detectable using Sanger sequencing, we readily detected synonymous edits in drug-resistant populations generated using this approach (Supplementary Fig. 1), indicating efficient co-editing at adjacent sites, likely at a common transcription or replication fork. Accordingly, and to avoid co-editing, we subsequently delivered each degenerate ssODN individually.

We selected compound DDD247 (compound 7a in ref. 21), a low nanomolar potency inhibitor, for detailed analysis, and used cryo-EM structural data21 to identify amino acids located within 5 Å of the proteasome β5 subunit binding site (Fig. 1a). The cognate codons for these twenty amino acids were distributed over a region of ~500 bp in the β5 gene. To assemble a DDD247 binding site mutant library, we designed twenty reverse-strand, 53-b, ssODNs, each with a centrally located degenerate ‘NNN’ (N=A, C, T, G) codon (Fig. 1b; Source Data, sheet 1). Mismatch-repair, previously reported to supress editing efficiency ~50-fold15, was transiently knocked down using MSH2 RNA interference for 24 h and, to avoid co-editing at adjacent sites in individual cells (see above), each degenerate ssODN was delivered individually. We then pooled the resulting populations and split the pool to generate a pair of Multiplex Oligo Targeting (MOT) libraries. Given a calculated allele replacement frequency of >0.04% using this approach (Supplementary Fig. 2), we estimated an average yield of >1000 codon edits per ssODN.

To assess editing at each of the targeted codons, we extracted genomic DNA before, and six hours after library assembly, PCR-amplified the edited region in the proteasome β5 gene, and deep-sequenced the amplicons (Fig. 1b). A scan to quantify variant codons across the edited region revealed highly specific editing at all twenty targeted sites (Fig. 1c; Source Data, sheet 2). We concluded that all 1280 possible variants around the β5 subunit drug binding site were likely represented in our pooled MOT libraries. These results served to validate the edited MOT libraries and the codon variant scoring approach.

More than 100 resistance-conferring base edits
To select drug-resistant mutants, the MOT libraries were grown with DDD247 at 8 nM; twice the EC50 (Effective Concentration of drug to inhibit growth by 50%). The selected libraries were then split three ways and grown with DDD247 at either 40, 200 or 1000 nM (Fig. 2a). We extracted genomic DNA following selection at each drug concentration, PCR-amplified the edited region in the proteasome β5 gene, deep-sequenced the products (Fig. 2a), and quantified variant codons. Principle Component Analysis of variant codon read counts revealed similar profiles for the replicate MOT libraries (Fig. 2b). An assessment of edits generating stop/nonsense codons (n=60) or synonymous codons (n=61), neither of which were expected to yield resistant cells, revealed negligible read counts across the dose-response profile (Fig. 2c). In contrast, a number of edits generating mis-sense (non-synonymous) codons yielded elevated read counts across the dose-response profile, at the Y113 codon for example (Fig. 2c).

The heatmap shown in Fig. 2d (upper panel) shows relative representation of all 1280 codon variants in the unselected library and across the dose-response profile (Source Data, sheet 3). We observed >100 resistance-conferring edits that were enriched following drug selection (Fig. 2d). By averaging read counts for sets of synonymous codons, we derived a second heatmap showing the behaviour of all 400 possible amino acid variants across the dose-response profile, revealing forty-six resistance-conferring edits that were enriched following drug selection (Fig. 2d, lower panel).

Those edits associated with drug resistance show that oligo targeting can be used to introduce all twelve possible base edits, with no apparent bias in terms of single-nucleotide edits (SNEs), double-nucleotide edits (DNEs), or triple-nucleotide edits (TNEs). Eight of the twenty sites targeted yielded edits associated with drug resistance, with G98, Y113, V128 and S132 emerging as resistance ‘hotspots’ (Fig. 2d). Conversion to asparagine yielded resistance in six cases and glutamine yielded resistance in four, while only conversion to proline or arginine failed to yield resistance.

We next selected a set of nine drug-resistant mutants to interrogate the predictive power of the MOT-library profile. We used the sequencing read count data to reconstruct virtual dose-response curves, and to derive virtual EC50 values for these mutants (Fig. 3a; Source Data, sheet 4; see Methods). As can also be observed in Fig. 2d, codons for the same amino acid edit yielded highly consistent results. A bespoke set of ssODNs (Source Data, sheet 1) was used to introduce each selected codon-edit in otherwise wild-type cells, and a panel of mutants was assembled, a pair of independent clones for each edit. The mutants, comprising six SNEs (G98A, GGA-GCA; Y113H, TAC-CAC; D115N, GAC-AAC; D116E, GAC-GAA; V128M, GTG-ATG; and Y136F, TAT-TTT), a DNE (V128Q GTG-CAG) and two TNEs (Y113G, TAC-GGA, S132M TCC-ATG), were validated by sequencing, confirming the expected edits in every case (Fig. 3b; Supplementary Fig. 3); these results also indicated that heterozygous edits were sufficient to yield resistance in every case. The mutant panel was then used to derive conventional dose-response curves and EC50 values (Fig. 3c; Supplementary Fig. 4). A comparison with the virtual EC50 data revealed an excellent correlation (R2=0.98), indicating that codon variant scores from MOT libraries were indeed predictive of relative resistance observed using the bespoke set of mutants, and across a wide range of EC50 values (Fig. 3d). Thus, we identified forty-six distinct drug-resistant mutants, which displayed up to 70-fold increased resistance relative to wild-type cells (Source Data, sheet 4), and which included all four previously know drug resistance mutations in trypanosomatid proteasome β5 subunits; Y113F in T. brucei25, D115N in T. cruzi27, and G98C and G98S in Leishmania21. We concluded that MOT-library profiling can be used to identify large numbers of distinct drug resistance associated edits in trypanosomes, and also that sequence read counts for variant codons across a dose-response profile are excellent predictors of relative resistance.

As a further test of MOT-library profiling, we screened the libraries used above with bortezomib, which binds a site immediately adjacent to the DDD247 binding site (Supplementary Fig. 5a). Bortezomib is a modified Phe-Leu dipeptide with a boronic acid warhead25; the compound is approved for cancer chemotherapy, but resistance is a significant problem23. Selection of the MOT-library with 4 nM bortezomib, approximately twice the EC50 for wild-type T. brucei, yielded resistant cells, while selection with higher drug concentrations failed to do so. Amplicon sequencing and codon variant scoring revealed only two conservative amino acid edits associated with bortezomib resistance, M97I and M97V (Supplementary Fig. 5b; Source Data); edits that notably failed to yield DDD247 resistance. These results highlight the ability of MOT-library screens to discriminate between inhibitors that engage immediately adjacent sites within the same target.

Fitness constraints on resistance space
More than 90% of edits fail to yield drug resistance in our screens, either because they fail to impact drug binding affinity, or simply because they yield defective proteasomes. Thus, mutational resistance space will be constrained by mutational fitness space at any given drug binding site. In order to assess mutational fitness space at the proteasome β5 subunit drug binding site, we deleted one copy of the β5 gene, assembled MOT libraries in the diploid and haploid strains, and derived codon variant scores after six hours and four days, without any drug selection (Fig. 4a; Source Data). To validate the approach, we first assessed edits producing stop codons, which were severely depleted in the haploid strain, but could be retained in the diploid strain, as expected (Fig. 4b). This analysis also revealed dominant-negative impacts of more distal stop codons in the diploid strain, likely due to the expression of longer, albeit truncated, β5 peptides.

We next turned our attention to edits generating mis-sense (non-synonymous) codons, comparing relative codon-variant scores from the diploid and haploid strains, and for fifteen sites for which we obtained sufficient data. Mutational tolerance was reduced at all fifteen of these sites in the haploid strain relatively to the diploid strain, indicating widespread loss-of-fitness associated with mutant proteasomes (Fig. 4c). Beyond this feature, three notable categories of edited codons emerged. First were codons for which most edits were poorly tolerated in both strains, consistent with dominant-negative defects (T1, G48, S127); second were codons for which many edits were tolerated only in the diploid strain, consistent with loss of proteasome function (A46, G47, M97, D116, V128, S130, Y136, A169); and third were codons for which many edits were tolerated in both strains, consistent with retained proteasome function (G23, S96, Y113, S132) (Fig. 4c). Dominant-negative edits were not expected to yield drug resistance, and indeed we saw no evidence for resistance-associated edits at these three sites, which included the catalytic threonine at T1 (Fig. 2d). The second loss-of-function category edits were also not expected to yield robust drug resistance, and indeed resistance-associated edits were limited at these eight sites, typically with failure to tolerate >40 nM drug selection (Fig. 2d); V128 did emerge as a resistance hotspot above but the mean virtual EC50 for edits at this site was only 23 nM (n=8). In contrast, edits in mutational fitness space, those associated with retained proteasome function, have the potential to yield robust resistance, and indeed we identified both Y113 and S132 as resistance hotspots above (Fig. 2d). Thus, although only a minority of edits around the proteasome β5 subunit drug binding sites sustain fitness, a substantial proportion of these edits yield robust drug resistance. We concluded that MOT-library profiling can be employed to elucidate both drug resistance space at a given drug binding site, and constraints imposed on resistance space by fitness space.

Predicting resistance in SNP-accessible space
The codon-degenerate ssODNs used above were able to introduce the full range of SNEs, DNEs, and TNEs, and to access all amino acid variants. The vast majority of naturally occurring mutations are limited to single nucleotides, however, and can access only nine alternative codons, or five to seven alternative amino acids (see Fig. 5a, lower left segment). Using the MOT-library profiling data to predict naturally occurring mutations that confer compound DDD247 resistance, we found that only fifteen of forty-six amino acid edits linked to drug resistance above were accessible via an SNE (Fig. 5a, upper right segment, green squares); six are at the Y113 resistance ‘hotspot’, while only DNEs and TNEs were linked to resistance at the S132 ‘hotspot’ (Fig. 5b).

Since generation of spontaneous proteasome inhibitor-resistant trypanosomatid mutants typically takes several months18,21,27, we took advantage of transient knockdown of mismatch-repair15 to rapidly generate a panel of resistant T. brucei mutants. We sub-cloned the inducible mismatch-repair knockdown strain, induced knockdown for 24 h, subjected each clone to either 40 nM or 100 nM drug pressure for 6-9 days, sub-cloned each resistant population, and sequenced the β5 gene in a single clone from each population. Sequence analysis revealed that twenty of twenty-three of these independent clones had mutations in the β5 gene, all of which were single nucleotide mutations (Fig. 5c), as expected. Thirteen were V128M (GTG-ATG) mutants, which may have arisen frequently due to the presence of a native editing template; we identified a potential A-mismatched template sequence (TGAAATTTTTAGTAT) on four distinct chromosomes. Notably, this mutant also registered a particularly strong resistance signal in Fig. 2d. Although there are twice as many possible transversions (A or G-C or T) as transitions (A-G or C-T; see Fig. 5a, b), transitions, like the V128M mutation, are more frequent. Indeed, assessment of 355,892 polymorphisms28 suggested that transitions are 1.8 times more frequent than transversion in T. brucei. Accordingly, the other seven β5 subunit mutants comprised six transitions (Y113C, TAC-TGC; Y113H, TAC-CAC; D115N, GAC-AAC*4) and a transversion (V128G, GTG-GGG) (Fig. 5c, Supplementary Fig. 6a), while the remaining three resistant clones had single nucleotide mutations in the proteasome β4 subunit (Tb927.10.4710), all at F24L (Supplementary Fig. 6b), which also interacts with compound DDD24721.

The resistance-associated mutation profile observed in T. brucei using MOT-library profiling (Fig. 2d) was also predictive of spontaneous drug resistance mutations reported previously in trypanosomatid proteasome β5 subunits; Y113F in T. brucei, which conferred resistance to LXE40825, D115N in T. cruzi, which conferred cross-resistance to DDD248, GNF6702, and TCMDC-14319427, and G98C and G98S in Leishmania, which conferred resistance to DDD24821. All of these mutations are accessible via single nucleotide changes; with G98C (GGC-TGC) and G98S (GGC-AGC) only accessible via single nucleotide mutations in Leishmania due to distinct codon-usage at this site, which is GGA in T. brucei. In summary, all spontaneous DDD247-resistant T. brucei clones had a single nucleotide mutation in a proteasome subunit, with 87% in the β5 subunit; all these β5 subunit mutations were predicted by MOT-library profiling (Figs. 2d and 5a); and no β5 subunit mutations were identified in residues more than 5 Å from the drug binding site (Fig. 1a). These observations, taken together, indicate that MOT-library profiling in T. brucei can facilitate accurate prospective prediction of spontaneous drug resistance mutations in multiple trypanosomatids, and for multiple anti-trypanosomal proteasome inhibitors.

Modelling ligand affinity
For computational modelling of ligand-proteasome binding affinity, we built a homology model of the T. brucei proteasome (Fig. 6a), revealing a high degree of similarity between the ligand binding sites in Leishmania21 and T. brucei (Fig. 6b). Docking calculations focused on the cavity formed by residues within 5 Å of the DDD247 ligand (Supplementary Fig. 7a, b) and yielded excellent correspondence with the published Leishmania - GSK3494245 structure (RMSD average: ~0.35 Å)21, while docking of compound DDD247 yielded a similar pose, with very good alignment to the aromatic and aliphatic rings (RMSD average: ~0.35 Å) (Supplementary Fig. 7c). Molecular dynamics simulations of the complex with compound DDD247 revealed stabilisation of ligand-protein interactions by hydrogen bonding with β5 residues, Y113, G129, and S130, in particular, and with other hydrophobic interactions also present along the simulation, with β5-Y113 and β4-F24; while the morpholine group displays evidence of being rather flexible and solvent exposed (Supplementary Fig. 8a, b).

We next considered all possible mutations in those twenty β5 subunit residues within 5 Å of the DDD247 ligand. After clustering snapshots from the molecular dynamics simulation trajectory (Supplementary Fig. 9), selected clusters were used for ΔΔG Molecular Mechanics General Born Surface Area (MM/GBSA) calculations29, revealing G98, Y113, V128, and S132 as hotspots where mutations were predicted to decrease ligand affinity (Fig. 6c). Remarkably, these were the same four hotspots identified by in cellulo MOT library profiling above (Fig. 2d). Among these hotspots, we observe ligand displacement at G98 and S132 (Fig. 6d). Indeed, the resistance-associated mutations we observe in cellulo all increase the side-chain mass at these sites; by 58+/−31 Da at G98 (n=7) and by 35+/−8 Da at S132 (n=7). Exemplar images illustrate how G98 and Y136 mutations have substantial or minimal impacts on ligand binding, respectively (Fig. 6e).

Hydrogen-bonding and lipophilicity are also predicted to have substantial impacts of protein-ligand affinity (Fig. 6f). Both features are predicted to contribute to the impact of Y113 mutations, since this residue can form a hydrogen bond with the ligand and also interact with in a lipophilic environment. Mutations at G98, although they do not form a hydrogen bond with the ligand, all disrupt hydrogen bonding between the adjacent Y113 residue and the ligand, while impacts on lipophilicity appear to be the major drivers of drug resistance for V128 mutations. In the case of S132, molecular dynamics clusters often placed this residue distal from the DDD247 ligand binding pose. Cluster 8, however, indicated a direct hydrogen bond between S132 and the ligand, and this bond is lost in mutants associated with drug resistance. Accordingly, we suggest that this hydrogen bond does occur but is not observed in the cryo-EM structure. By comparing these ΔΔG (MM/GBSA) calculations with the full set of forty-six virtual EC50 values derived using MOT library profiling above, we obtained an R² value of 0.42 (Fig. 6g). As an example of excellent correlation between observation in cellulo and prediction in silico, Y113 mutations incorporating highly charged residues (Y113K, Y113D, Y113Q, Y113S) are associated with a higher degree of resistance relative to mutations that conserve hydrophobicity (Y113F, Y113L, Y113H, Y113M).

Modelling drug resistance
We elucidate above the mutational fitness space that restricts mutational resistance space at the DDD247 binding site; whereby non-functional mutant proteasomes fail to yield drug-resistant cells (Fig. 4). We, therefore, sought to combine in silico predictions of loss of protein function with our ligand affinity predictions. We used Evo 2, a DNA language model trained on >9 trillion DNA base pairs from all domains of life30, to predict loss of proteasome function due to codon changes, considering the full set of 1280 mutations surveyed by MOT-library profiling above. Evo2 Δ_scores (Supplementary Fig. 10) were used to assign penalties relative to loss-of-function for each mutated β5 subunit residue. Both the affinity scores (Fig. 7a, left-hand panel) and functional scores (Fig. 7a, right-hand panel) were then combined to generate a heatmap showing resistance predictions (Fig. 7a, lower panel); for those mutations that both disrupt ligand affinity and retain fitness. The resistance hotspots described above, G98, Y113, V128 and S132, were retained, and the computationally predicted impact of each individual mutation can now be seen in the combined heatmap. Although the analysis predicted reduced ligand affinity due to mutations at A46 and G47, loss of proteasome function is also predicted due to these mutations, meaning that resistance is predicted to be minimal (Fig. 7a). Indeed, mutations at these sites registered a substantial loss-of-fitness in cellulo (Fig. 4c) and failed to register as resistance mutations (Fig. 2d) in our in cellulo screen. Further highlighting the potential for added value when combining affinity and function parameters to predict drug resistance, Pearson’s r for the full dataset was increased from 0.46 to 0.51 after applying penalties for loss-of-fitness; Pearson’s r for loss-of-fitness alone was 0.01. A circular plot shows the full in cellulo resistance profile and in silico predicted resistance profile side-by-side for comparison, and a correlation plot shows how effectively our in silico analyses predicted the in cellulo derived resistance hotspots (Fig. 7b, Pearson’s r=0.9).

Methods
T. brucei growth and manipulation
Bloodstream form Lister 427 Trypanosoma brucei wild-type and MSH2 RNAi15 strains were cultured in HMI-11 (Gibco) supplemented with 10% fetal bovine serum (Sigma) at 37 °C and with 5% CO2 in a humidified incubator. Genetic manipulations was carried out by electroporation using a Nucleofector (Lonza), with a Human T-cell kit (Lonza), or with custom-made transfection buffer34, with the Nucleofector set to Z-001 (Amaxa).

Oligo targeting
Oligo targeting was carried out using single-stranded oligodeoxynucleotides (ssODNs, Thermo Fisher Scientific) essentially as described15 (Source Data, sheet 1). Briefly, for routine assays, we used 40 μg of each ssODN in 10 μl of 10 mM Tris-Cl, pH 8.5, mixed with 25 million cells in 100 μl transfection buffer. Selection with DDD247 or DDD248, compounds 7a and 7, respectively in Wyllie et al.21, was applied 6 h later, as appropriate, and cells were subjected to limiting dilution in 96-well plates as required. Genomic DNA was extracted using a Qiagen DNeasy Kit followed by PCR amplification of the β5 gene (25 cycles). PCR products were purified using a PCR Purification Kit (Qiagen). To identify individual edits, amplicons were generated using Phusion polymerase (NEB) and amplicons were subjected to Sanger sequencing (Genewiz, Azenta Life Sciences).

Assembly and screening of multiplexed oligo targeting libraries
For site saturation mutagenesis using a set of twenty degenerate ssODNs, we used 40 μg of each ssODNs in 10 μl of 10 mM Tris-Cl, pH 8.5, mixed with 6.5 million cells in 100 μl transfection buffer. MSH2 knockdown was induced 24 h prior to transfection by the addition of 1 μg/mL tetracycline (Sigma). Each ssODN was transfected individually and the cultures were then pooled in 150 ml of medium to generate the library, comprising ~65 million cells based on an electroporation survival rate of ~50%. After 6 h, cells from 20 ml of the culture were cryo-preserved, cells from another 20 ml were collected for (pre-selection) DNA extraction, and the remainder were split into two cultures and subjected to selection with 8 nM compound DDD247. Following 4–5 d, 50 ml of each culture was collected for DNA extraction, and the remainder was split into three cultures that were subjected to selection with either 40 nM, 200 nM or 1 μM compound DDD247. Following 2–3 d, each culture was collected for DNA extraction.

For selection with bortezomib (Sigma), a cryo-preserved library sample was thawed, allowed to recover for 24 h, and split into two parallel cultures. Bortezomib selection was applied at 4 nM (EC50 = 2.2 nM). Following 7–9 d, cultures were collected for DNA extraction and the remainder was split into three cultures that were subjected to selection with higher concentrations of bortezomib; none of these latter cultures yielded resistant cells.

For fitness profiling, we first generated a T. brucei MSH2 RNAi strain with a single β5 allele. A construct was synthesized comprising a BLA (blasticidin deaminase) selectable marker, flanked by T. brucei tubulin mRNA processing sequences and ~300 bp of β5 gene flanking sequence to promote homologous recombination. This cassette was excised from the plasmid using AscI and used to transfect MSH2 RNAi cells. BLA selection was applied at 10 μg/ml after 6 h and positive clones were assessed using PCR assays. Multiplex oligo targeting libraries were then assembled in both strains, either haploid or diploid for the β5 gene, as described above, and again using the same set of twenty degenerate ssODNs. In this case, cells were collected for DNA extraction after 6 h and 4 d without any drug selection. To identify multiplexed edits, PCR amplicons from the β5 gene were generated using KOD Hot Start DNA Polymerase (Merck). We used 0.5 μg of genomic DNA template to assess library coverage and 0.1 μg of genomic DNA template to assess drug-selected samples; 25 cycles in all cases. Amplicons were purified using a PCR Purification Kit (Qiagen).

Amplicon sequencing and codon variant scoring
Amplicons were subjected to deep sequencing using DNBseq (BGI Genomics); the controls and drug-selected samples shown in Figs. 1 and 2 were sequenced in parallel. Filtering of sequencing reads was performed using the SOAPnuke software with parameters: “-n 0.001 -l 10 -q 0.4 --adaMR 0.25 --ada_trim --minReadLen 150” for drug selection experiments and “-n 0.01 -l 20 -q 0.3 --adaMR 0.25 --ada_trim --polyX 50 --minReadLen 150” for the haploid v diploid experiment. Codon variant scoring was performed with the OligoSeeker (0.0.5) Python package35 designed to process paired FASTQ files and count occurrences of specific codons. To visualise codon variant scores, we performed a normalization step by dividing each codon variant score by the total reads for that position and converting the fraction of reads to a percentage. This was followed by background correction whereby we subtracted the values for the control sample. Negative values were replaced with 0, and average values for the duplicate libraries were then calculated to give the codon variant scores. Codons registering elevated read-counts in the control sample, including all single nucleotide variants, were excluded from the analyses shown in Figs. 1c and 4c. To derive virtual EC50 values, we derived a weighted average of codon variant scores across the DDD247 compound dose range. Briefly, the drug concentration associated with the highest read-count was assigned a weighting of ‘1’ and relative weighting was assigned accordingly. For example, the virtual EC50 for G98A is the average of 0.48*8, 1*40, 0.34*200, and 0.02*1000 = 32 nM. To visualize relative codon frequency using the pyCirclize python package (https://github.com/moshi4/pyCirclize, v1.9.1), raw codon counts were normalized by the maximum values and replicates were then averaged.

Assembly of a panel of spontaneous drug-resistant mutants
We exploited the MSH2 RNAi strain to rapidly generate a panel of spontaneous and independent DDD247-resistant mutants. To ensure that each mutant was independent, we first used limiting dilution to isolate twenty clones from the MSH2 RNAi cell line. MSH2 knockdown was induced in each clone for 24 h, tetracycline was removed by washing in HMI-11 medium, and compound DDD247 was added at 40 nM or 100 nM. We observed cell death followed by recovery and obtained twenty-three independent DDD247-resistant cultures using this approach. We sub-cloned each culture by limiting dilution, and isolated one clone from each for DNA extraction, PCR amplification of the β5 gene and Sanger sequencing (Genewiz, Azenta Life Sciences). We also subsequently PCR amplified and sequenced the proteasome β4 gene in three clones for which no mutations were identified in the β5 gene.

Dose-response assays
To determine the Effective Concentration of drug to inhibit growth by 50% (EC50), cells were plated in 96-well plates at 1×103 cells/ml in a 2-fold serial dilution of selective drug. Plates were incubated at 37 °C for 72 h, 20 μl resazurin sodium salt (Alamar Blue, Sigma) at 0.49 mM in PBS was then added to each well, and plates were incubated for a further 6 h. Fluorescence was determined using an Infinite 200 pro plate reader (Tecan) at an excitation wavelength of 540 nm and an emission wavelength of 590 nm. EC50 values were derived using Prism (GraphPad).

Homology modelling, docking and molecular dynamics
A T. brucei proteasome homology model was generated using a cryo-EM structure of the Leishmania proteasome in complex with the GSK3494245 inhibitor (PDB: 6qm7)21, and the homology modelling tool implemented in Maestro (Schrodinger inc.)36. Docking was performed with Glide XP37, using GSK3494245 for validation, followed by compound DDD247, with RMSD comparisons calculated using an MCS protocol. The complex was placed in a Cubic box with waters and ions (0.15 M) using Maestro and the system was subjected to 80 ns of molecular dynamics with Desmond38. The trajectory was clustered by RMSD (resulting in 9 clusters), and representative frames were selected for further analysis. Those clusters were classified with interaction fingerprints39.

Free energy estimations
The ensemble of conformations was generated from the clusters produced by the simulation to estimate the binding free energy. The residue scanning tool implemented in Maestro was used to compute approximate binding energy scores between the ligand and wild-type/mutant proteins38. This tool implements a high-throughput workflow for testing mutations with the ΔΔG associated (calculated with the MM/GBSA method)29 using wild-type as reference (baseline). MM/GBSA is a free energy method that gives an approximation of relative binding free energy between states (such as ligand bound vs unbound). The method decomposes the calculation into terms; Van der Waals, Coulombic, hydrogen bonding, solvation and hydrophobic, for contributions to the final ΔG. All residues tested experimentally were mutated to build single mutants, with minimisation around 7 Å as the final step of preparation. Results were plotted with Seaborn and Matplotlib. For virtual EC50 values derived from MOT-library profiling data, a ΔpvEC50 was calculated for each mutant, and plotted against the ΔΔG values.

Protein function prediction
Machine learning approaches can effectively predict the impact of mutations on protein function40 and we used an Evo 2 model to predict loss of proteasome function here (evo2_7b_base)30. Taking native nucleotide sequence for the β5 subunit gene, we considered all possible mutations in those codons surveyed by MOT-library profiling. After the inference, the native sequence was used as a reference, and sequences with mutations were scored and subtracted from the reference score to obtain Evo2 Δ_scores. These scores were then used to apply penalties to the affinity predictions. Values were normalised to allow comparison of in cellulo and in silico predictions and the results were visualised in a circular plot with Seaborn and Matplotlib.

Data availability
The high-throughput sequencing data generated for this study have been deposited at the Sequence Read Archive under accession code PRJNA1295514 (https://www.ncbi.nlm.nih.gov/bioproject/1295514). Source data for Figs. 1c, 2b–d, 3a, 4b, c and 5a, b are available as a Source Data file and for Figs. 6c, d, f, g, 7a, b, Supplementary Figs. 7c and 10a, b in the Github repository (https://github.com/velocirraptor23/Decoding-efficacy-and-resistance-space-at-a-drug-binding-site-/tree/main) and Zenodo https://zenodo.org/records/18195459. The homology model used in Figs. 6 and 7 is available in Model Archive (https://modelarchive.org/doi/10.5452/ma-cxis8).

Code availability
The code used to develop the homology model, perform analyses and generate results in this study is publicly available and has been deposited in Github at https://github.com/velocirraptor23/Decoding-efficacy-and-resistance-space-at-a-drug-binding-site-/tree/main, under MIT license. The specific version of the code associated with this publication is archived in Zenodo and is accessible via https://zenodo.org/records/1819545941.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this work, the authors present a systematic mapping of the proteasome β5 subunit mutations conferring resistance to the inhibitor DDD247 using Multiplex Oligo Targeting (MOT). The key contributions lie in 1. the accurate design of the pooled MOT libraries with massive codon edits across residues in β5 subunit binding pockets. 2. High-resolution experimental profiling of resistance, including virtual EC50 estimation directly from pooled deep mutational scans. 3. Identification of four major resistance hotspots (e.g., Y113, G98, V128, S132), with codon-level granularity of tolerated vs. non-tolerated substitutions. 4. Integration of Evo2 predictions (sequence-based functional tolerance) with docking-based ligand affinity modelling to propose a combined predictive framework for resistance.\nFrom the computational side, the integration of Evo2, along with other computational approaches, is novel but the results are quite qualitative and are not significant and solid enough to support a unique novelty. Adding the following details would significantly strengthen the article:\n1. Virtual EC50 derivation:\n* Please provide more detail on how virtual EC50 values were computed from pooled sequencing data.\n* How was the Hill slope (n) handled, given that only four drug concentrations were tested? Was it fixed at 1? How sensitive the correlation is to this hyperparameter.\n2.The bespoke mutant validation set appears to be enriched for strong resistance alleles. This partly explains the near-perfect correlation. Please discuss this validation bias and, if possible, include weaker or borderline variants to test the robustness of virtual EC50 inference.\n3. Evo2 evaluation: the contribution of Evo2 is argued qualitatively.\nI encourage adding quantitative benchmarks comparing: a. Docking-only predictions, b. Evo2-only predictions, and c. Docking+Evo2 combined. Suggested metrics include correlation to (virtual) EC50, AUROC/AUPRC against a binary resistance definition, and hotspot overlap with experimental data.\n4. Figure 7 lacks quantitative metrics. It is difficult to judge whether predictions and experiments align. For example, a scatter plot of average predicted vs experimental resistance per residue would be helpful. Figure 6h is on the right track. However, the analysis is limited to selected variants and only reports R2. I encourage the authors to extend this comparison to the full dataset and include additional metrics such as Spearman correlation, AUROC/AUPRC, and site-level hotspot overlap (i.e., Experimentally observed hotspots (Y113, G98, etc) vs predicted hotspots from docking+Evo2. )\nAlthough the computational analyses could be strengthened, the genetic and experimental work of the paper is strong and novel.\n(Remarks on code availability)\nI have the following comments about the repository:\nPlease provide detailed descriptions for the MD_simulation, data, and MD_clusters subfolders, including:\n* Purpose and contents of each directory\nThe main README file should include:\n* A comprehensive introduction to each subfolder's role in the analysis pipeline\n* A step-by-step execution guide that walks through the entire workflow", "reviewer_2": "Altmann et al have conducted a comprehensive and impressive study employing oligo targeting to perform nearly complete saturation mutagenesis at the native drug-binding site of the Trypanosoma brucei proteasome beta-5 subunit, successfully generating detailed resistance and fitness landscapes. By combining high-throughput mutational screening, resistance profiling, and computational modelling, they demonstrate both the accuracy of their approach and its potential to predict spontaneous resistance mutations, providing valuable insights into the constraints shaping resistance evolution. The work significantly advances the methodology for mutagenesis in protozoan parasites and offers novel biological and pharmaceutical insights, especially regarding the identification of resistance hotspots and the interplay between fitness and resistance. Overall, this represents a major step forward for drug target characterization in neglected parasitic diseases.\nKey findings include the successful application of multiplex oligo targeting to saturate 20 residues within the drug-binding pocket, the identification of over 100 resistance-conferring mutations with associated fitness effects, and the validation that these mutations mirror spontaneous resistance mutations observed in vitro and across species. Further, the integration of structural modelling and AI-driven predictions supports the empirical data, highlighting the potential for this approach, though the correlation between computed and experimental data remains moderate (R² = 0.42). But this limitation is acknowledged.\nI have only minor comments:\n1. Lines 50-70: The authors effectively discuss the limitations of CRISPR-based editing and the advantages of oligo targeting. However, it would be helpful to include additional information about the scalability of this approach. Specifically, how large a library can realistically be generated, and how many sites can be targeted simultaneously, given current system limitations? A brief comparison of this scalability to what is achievable with CRISPR-based or base-editing screens would also strengthen this section.\n2. Figure 2a should include the untreated control condition for clarity and ease of interpretation.\n3. It would also be useful to briefly restate the rationale for using the MSH2 RNAi cell line, even if this was covered in Altmann et al. (2022), to make the manuscript more self-contained.\n4. A concise summary paragraph highlighting the key new findings and advances compared with previous literature would help readers clearly understand the novel contributions of this study. This is particularly important given the extensive existing literature on the T. brucei proteasomes.\n5. Lines 510-517: While the authors provide clear evidence that fitness profiling is feasible using their approach, it would be helpful for others if more methodological detail were included on how this was achieved without additional selection and despite the relatively low transfection efficiencies previously reported for oligo targeting. For example, could the authors clarify what the effective representation of the transfected library was (i.e., the number of cells per construct) and how many total cells were transfected to achieve this level of representation? I am also uncertain whether the library used for the fitness screen was actually subjected to 8 nM DDD247 selection. In the Methods section, the authors describe a general library transfection followed by compound selection, yet in the main text (line 197) they state that fitness profiling was performed “without any drug selection”. It appears that the fitness profiling did not involve drug treatment, but if so, sequencing depth must have been very high to achieve sufficient resolution of fitness profiles. Providing these details would substantially improve clarity.\n6. Finally, the library design does not seem to include a “non-targeting” control, which is typically incorporated into CRISPR or similar screening libraries to demonstrate assay specificity. While the authors do validate specificity by profiling the library with bortezomib and present other supporting data that increase confidence in the results, it would be useful if they could briefly comment on why such a non-targeting control was not included or whether they would include such controls in the future.\n(Remarks on code availability)\nWhile I have not tested the code, the GitHub repository appears well organized and clearly documented. It should be straightforward to reproduce the presented data and to follow the provided instructions for applying the workflow to similar datasets.", "reviewer_3": "Altmann et al describe the utilization of multiplex oligo targeting, deep sequencing and computational modeling to identify mutations that confer resistance and loss-of-function phenotype by carrying out saturation mutagenesis of selected interactions of the molecule (inhibitor) binding to target. They demonstrate the use of this workflow by utilizing known proteasome inhibitor and elegantly identify key mutations that are leading to robust resistance hotspots. This work integrates empirical resistance/fitness landscapes with in-silico docking, molecular dynamics clustering, and Evo2 functional predictions, demonstrating a combined framework to predict relative resistance and constrain mutational space by fitness.\nBelow are some of the potential points to consider and discuss\n1. Generalizability of findings: The work is centered on a single compound (DDD247) and one parasite species T. brucei; resistance landscapes may differ across related inhibitors or species. It would have been interesting to show, if these mutations can be validated by another proteasome inhibitor with a different chemical structure. Although, there is high likelihood that these mutations do confer resistance to other structurally different inhibitor.\n2. Co-editing artifacts: Authors have pointed out initial pooled ssODN approach can introduce multiple edits per cell, potentially confounding genotype–phenotype links. Although, the switch to individual delivery helps, but residual co-editing risk and selection bottlenecks may need to be carefully controlled.\n3. Heterozygous versus homozygous: The technology used is very helpful for targets which depend heterozygous mutations are dominant thereby showing resistance. If the resistance needs homozygosity (mutations in both alleles), this method might not be that useful. Authors could highlight this limitation in their discussion.\nAllele-specific effects: T. brucei is diploid and the current technology, typically edits a single allele; heterozygosity may mask fitness costs or resistance phenotypes compared to full homozygous states.\n4. Selection biases: Stepwise drug selection can enrich fitter mutants while missing rare but relevant variants; I am hoping this issue might not be of big concern given the extensive efforts made by authors using varying drug concentration and deep sequencing in order to isolate all the possible mutations.\n5. Modeling assumptions: Combining Evo2 functional scores with docking/free-energy predictions involves thresholds and normalization choices that may influence classification of loss-of-function vs resistance; Very heartening to see significant correlation between virtual EC50s and measured EC50s.\n6. Structural inference: Homology modeling and docking rely on a related proteasome structure; uncertainties in local side-chain positioning could affect predicted interaction fingerprints.\n7. The current methodology may not be able to capture structural changes that may affect compound binding, if the mutations far (>5Å), thereby leading to compound resistance and loss-of-fitness.\nOverall, the work appears methodologically solid and well-documented, with clear resistance hotspots identified. The main caveats are the scope (single target/compound/species) and potential biases from editing and selection workflows, which are common challenges in mutational scanning studies.\nThe MOT and computational modeling has potential not only to enable drug design but also to monitor resistance surveillance.\nMinor comments\nLine 23-24: the sentence ending “however”. May be you can start with “However, ….”\nLine 171-172: Authors refer to “forty-six distinct drug-resistant mutants, which displayed up to 70- fold increased resistance”. I was unable to find this tables. Please refer to Table.\n(Remarks on code availability)"}, "author_response": "Altmann et al.,\nDecoding efficacy and resistance space at a drug binding site\nNature Communications\nResponses to Reviewer Comments\nReviewer #1:\nIn this work, the authors present a systematic mapping of the proteasome β5 subunit\nmutations conferring resistance to the inhibitor DDD247 using Multiplex Oligo Targeting\n(MOT). The key contributions lie in 1. the accurate design of the pooled MOT libraries with\nmassive codon edits across residues in β5 subunit binding pockets. 2. High-resolution\nexperimental profiling of resistance, including virtual EC50 estimation directly from pooled\ndeep mutational scans. 3. Identification of four major resistance hotspots (e.g., Y113, G98,\nV128, S132), with codon-level granularity of tolerated vs. non-tolerated substitutions. 4.\nIntegration of Evo2 predictions (sequence-based functional tolerance) with docking-based\nligand affinity modelling to propose a combined predictive framework for resistance.\nFrom the computational side, the integration of Evo2, along with other computational\napproaches, is novel but the results are quite qualitative and are not significant and solid\nenough to support a unique novelty. Adding the following details would significantly strengthen\nthe article:\n1.1: Virtual EC50 derivation:\n* Please provide more detail on how virtual EC50 values were computed from pooled\nsequencing data.\n* How was the Hill slope (n) handled, given that only four drug concentrations were tested?\nWas it fixed at 1? How sensitive the correlation is to this hyperparameter.\nR1.1: We’ve added further detail to the Methods section, which now states “Briefly, the drug\nconcentration associated with the highest read-count was assigned a weighting of ‘1’ and\nrelative weighting was assigned accordingly. For example, the virtual EC for G98A is the\n50\naverage of 0.48*8, 1*40, 0.34*200, and 0.02*1000 = 32 nM”. We’ve also added “see Methods”\nat the relevant point in the Results section.\n1.2: The bespoke mutant validation set appears to be enriched for strong resistance alleles.\nThis partly explains the near-perfect correlation. Please discuss this validation bias and, if\npossible, include weaker or borderline variants to test the robustness of virtual EC50 inference.\nR1.2: The mean virtual EC is 94 nM (n = 46) and the mean conventional EC is 91 nM (n =\n50 50\n9); see Supplementary Data File (sheet 4). Thus, the mutant validation set is not significantly\nenriched for strong resistance alleles and there’s no significant validation bias (p = 0.94).\n1.3: Evo2 evaluation: the contribution of Evo2 is argued qualitatively.\nI encourage adding quantitative benchmarks comparing: a. Docking-only predictions, b. Evo2-\nonly predictions, and c. Docking+Evo2 combined. Suggested metrics include correlation to\n(virtual) EC50, AUROC/AUPRC against a binary resistance definition, and hotspot overlap\nwith experimental data.\nR1.3: We thank the reviewer for this suggestion and have now used quantitative data for the\ncontribution of Evo2 data to the predictions shown in Fig. 7, with further details included in\nSupplementary Fig. 10. We did not expect to see a strong correlation between Evo2 only\npredictions and our in cellulo resistance profile, partly because we believe that many “edits fail\nto yield drug resistance in our screens … because they fail to impact drug binding affinity”.\nIndeed Pearson’s r is 0.01 in this case. We were reassured, however, to see that predicted\nloss-of-fitness was significantly lower for the set of 46 mutations we linked to drug resistance,\nrelative to all 334 other possible mis-sense mutations (p = 0.01; Student’s t-test). We’ve now\nadded “Pearson’s r for the full dataset was increased from 0.46 to 0.51 after applying penalties\nfor predicted loss-of-fitness” to the Results text.\n1\n1.4: Figure 7 lacks quantitative metrics. It is difficult to judge whether predictions and\nexperiments align. For example, a scatter plot of average predicted vs experimental resistance\nper residue would be helpful. Figure 6h is on the right track. However, the analysis is limited\nto selected variants and only reports R2. I encourage the authors to extend this comparison\nto the full dataset and include additional metrics such as Spearman correlation,\nAUROC/AUPRC, and site-level hotspot overlap (i.e., Experimentally observed hotspots\n(Y113, G98, etc) vs predicted hotspots from docking+Evo2. )\nR1.4: We’ve now included a scatter plot for average predicted vs experimental resistance per\nresidue, and added “a correlation plot shows how effectively our in silico analyses predicted\nthe in cellulo derived resistance hotspots (Fig. 7b, Pearson’s r = 0.9)” to the Results text. The\nFig. 7 legend also now states “The correlation plot shows the same data averaged for all mis-\nsense mutations at each site, with resistance hotspots identified”.\nAlthough the computational analyses could be strengthened, the genetic and experimental\nwork of the paper is strong and novel.\nRemarks on code availability:\n1.5: I have the following comments about the repository:\nPlease provide detailed descriptions for the MD_simulation, data, and MD_clusters\nsubfolders, including:\n* Purpose and contents of each directory\nThe main README file should include:\n* A comprehensive introduction to each subfolder's role in the analysis pipeline\n* A step-by-step execution guide that walks through the entire workflow\nR1.5: Updated as requested.\nReviewer #2:\nAltmann et al have conducted a comprehensive and impressive study employing oligo\ntargeting to perform nearly complete saturation mutagenesis at the native drug-binding site of\nthe Trypanosoma brucei proteasome beta-5 subunit, successfully generating detailed\nresistance and fitness landscapes. By combining high-throughput mutational screening,\nresistance profiling, and computational modelling, they demonstrate both the accuracy of their\napproach and its potential to predict spontaneous resistance mutations, providing valuable\ninsights into the constraints shaping resistance evolution. The work significantly advances the\nmethodology for mutagenesis in protozoan parasites and offers novel biological and\npharmaceutical insights, especially regarding the identification of resistance hotspots and the\ninterplay between fitness and resistance. Overall, this represents a major step forward for drug\ntarget characterization in neglected parasitic diseases.\nKey findings include the successful application of multiplex oligo targeting to saturate 20\nresidues within the drug-binding pocket, the identification of over 100 resistance-conferring\nmutations with associated fitness effects, and the validation that these mutations mirror\nspontaneous resistance mutations observed in vitro and across species. Further, the\nintegration of structural modelling and AI-driven predictions supports the empirical data,\nhighlighting the potential for this approach, though the correlation between computed and\nexperimental data remains moderate (R² = 0.42). But this limitation is acknowledged.\nI have only minor comments:\n2.1: Lines 50-70: The authors effectively discuss the limitations of CRISPR-based editing and\nthe advantages of oligo targeting. However, it would be helpful to include additional information\nabout the scalability of this approach. Specifically, how large a library can realistically be\ngenerated, and how many sites can be targeted simultaneously, given current system\n2\nlimitations? A brief comparison of this scalability to what is achievable with CRISPR-based or\nbase-editing screens would also strengthen this section.\nR2.1: We felt that it was more appropriate to address this point in the Discussion. Several\nrecent studies report CRISPR-based editing of the EGFR gene, but only Kim et al (Nature\nBiotechnology, 2025) provide precise figures for library size, and we therefore use this report\nas a comparator. We now include the following paragraph in our Discussion:\n“The current T. brucei libraries incorporated 1,280 alternative codons encoding 400\nproteasome b5 subunit variants. In this case, our goal was to introduce every possible\nalternative codon at each of twenty targeted sites. Based on an allele replacement frequency\nof approx. 0.04 % following transient suppression of mismatch repair, we estimated 26,000\nedited cells among 65 million cells in a 150 ml culture volume, which equates to an average\ntwenty-fold redundancy for each alternative codon. For comparison, a library generated\nrecently using prime editors in human cells incorporated 2,476 single nucleotide edits\nencoding 1,726 epidermal growth factor receptor variants (Kim et al., 2025). In terms of\nscaling, our results suggest that larger libraries could readily be generated and screened using\nthe oligo targeting approach by simply using larger numbers of degenerate ssODNs.\nSubstantial further scaling would present challenges in terms of assessing every alternative\ncodon at each targeted site, but use of larger cell numbers and culture volumes would reduce\nbottlenecks and maintain library complexity”.\n2.2: Figure 2a should include the untreated control condition for clarity and ease of\ninterpretation.\nR2.2: The untreated control library is described in Fig. 1. To clarify, we added “described\nabove” to the Fig. 1 legend. We also now note in the Methods section that “the controls and\ndrug-selected samples shown in Fig. 1-2 were sequenced in parallel”.\n2.3: It would also be useful to briefly restate the rationale for using the MSH2 RNAi cell line,\neven if this was covered in Altmann et al. (2022), to make the manuscript more self-contained.\nR2.3: We previously restated the rationale on lines 114-115, but have now added ‘MSH2’ to\nthis text to clarify: “Mismatch-repair, previously reported to supress editing efficiency approx.\n50-fold 14, was transiently knocked down using MSH2 RNA interference for 24 h”.\n2.4: A concise summary paragraph highlighting the key new findings and advances compared\nwith previous literature would help readers clearly understand the novel contributions of this\nstudy. This is particularly important given the extensive existing literature on the T. brucei\nproteasomes.\nR2.4: We thank the reviewer for this suggestion and have added the following text to the\nsummary paragraph at the end of the Discussion:\n“We developed multiplex oligo targeting in T. brucei and used the approach to saturation edit\na native drug binding site in the proteasome b5 subunit. We profiled both fitness and\nresistance space and identified resistance hotspots and >100 resistance-conferring mutants,\nwhich allowed us to accurately predict spontaneous resistance mutations in multiple\ntrypanosomatids. Assessment of the empirical data using computational modelling provided\nunprecedented insights into anti-trypanosomal proteasome inhibitor structure-activity\nrelationships”.\n2.5: Lines 510-517: While the authors provide clear evidence that fitness profiling is feasible\nusing their approach, it would be helpful for others if more methodological detail were included\non how this was achieved without additional selection and despite the relatively low\ntransfection efficiencies previously reported for oligo targeting. For example, could the authors\nclarify what the effective representation of the transfected library was (i.e., the number of cells\nper construct) and how many total cells were transfected to achieve this level of\nrepresentation? I am also uncertain whether the library used for the fitness screen was actually\nsubjected to 8 nM DDD247 selection. In the Methods section, the authors describe a general\nlibrary transfection followed by compound selection, yet in the main text (line 197) they state\n3\nthat fitness profiling was performed “without any drug selection”. It appears that the fitness\nprofiling did not involve drug treatment, but if so, sequencing depth must have been very high\nto achieve sufficient resolution of fitness profiles. Providing these details would substantially\nimprove clarity.\nR2.5: We’ve now added further detail in the Methods section; “the library, comprising approx.\n65 million cells based on an electroporation survival rate of approx. 50 %”, and also more\ndetail in the Discussion, including “we estimated an average twenty-fold redundancy for each\nalternative codon” (see R2.1 above). For the fitness profiling section, we added more detail to\nclarify in the Methods section; “In this case, cells were collected for DNA extraction after 6 h\nand 4 d without any drug selection”. We’ve also added detail in the Figure legends for\nsequencing depth, which was indeed higher for the controls shown in Fig. 1c (“average of >33\nM reads mapped per site”) and for the fitness profiling shown in Fig. 4c (“average of >27 M\nreads mapped per site”), than for resistance profiling shown in Fig. 2d (“average of >6.5 M\nreads mapped per site”).\n2.6: Finally, the library design does not seem to include a “non-targeting” control, which is\ntypically incorporated into CRISPR or similar screening libraries to demonstrate assay\nspecificity. While the authors do validate specificity by profiling the library with bortezomib and\npresent other supporting data that increase confidence in the results, it would be useful if they\ncould briefly comment on why such a non-targeting control was not included or whether they\nwould include such controls in the future.\nR2.6: Non-targeting control sgRNAs are typically included in CRISPR experiments to establish\na baseline for off-target editing by Cas9. We feel that a non-targeting control oligo would be of\nlimited value for oligo targeting, since a ‘zero’ baseline is already set by several targeting\noligos, those that fail to yield any drug resistance associated edits.\nRemarks on code availability:\nWhile I have not tested the code, the GitHub repository appears well organized and clearly\ndocumented. It should be straightforward to reproduce the presented data and to follow the\nprovided instructions for applying the workflow to similar datasets.\nReviewer #3:\nAltmann et al describe the utilization of multiplex oligo targeting, deep sequencing and\ncomputational modeling to identify mutations that confer resistance and loss-of-function\nphenotype by carrying out saturation mutagenesis of selected interactions of the molecule\n(inhibitor) binding to target. They demonstrate the use of this workflow by utilizing known\nproteasome inhibitor and elegantly identify key mutations that are leading to robust resistance\nhotspots. This work integrates empirical resistance/fitness landscapes with in-silico docking,\nmolecular dynamics clustering, and Evo2 functional predictions, demonstrating a combined\nframework to predict relative resistance and constrain mutational space by fitness.\nBelow are some of the potential points to consider and discuss\n3.1: Generalizability of findings: The work is centered on a single compound (DDD247) and\none parasite species T. brucei; resistance landscapes may differ across related inhibitors or\nspecies. It would have been interesting to show, if these mutations can be validated by another\nproteasome inhibitor with a different chemical structure. Although, there is high likelihood that\nthese mutations do confer resistance to other structurally different inhibitor.\nR3.1: We had previously noted that our data can predict spontaneous drug resistance\nmutations observed in other trypanosomatids (see lines 257-262) and have now added more\ndetail to the current text regarding other inhibitors; “Y113F in T. brucei, which confers resistance\nto LXE408 23, D115N in T. cruzi, which confers cross-resistance to DDD248, GNF6702, and\nTCMDC-143194 25, and G98C and G98S in Leishmania, the latter of which confers resistance\nto DDD248 7”; also now noting that our data can predict “spontaneous drug resistance\n4\nmutations in multiple trypanosomatids, and for multiple anti-trypanosomal proteasome\ninhibitors”.\n3.2: Co-editing artifacts: Authors have pointed out initial pooled ssODN approach can\nintroduce multiple edits per cell, potentially confounding genotype–phenotype links. Although,\nthe switch to individual delivery helps, but residual co-editing risk and selection bottlenecks\nmay need to be carefully controlled.\nR3.2: A co-editing artefact would be expected to produce a background signal across tolerated\ncodons, including synonymous codons. We see no evidence for this in Fig. 2c-d. We have\nnow added more detail on library complexity and mitigation against bottlenecks, in both the\nMethods and Discussion sections (see R2.1 and R2.5 above).\n3.3: Heterozygous versus homozygous: The technology used is very helpful for targets which\ndepend heterozygous mutations are dominant thereby showing resistance. If the resistance\nneeds homozygosity (mutations in both alleles), this method might not be that useful. Authors\ncould highlight this limitation in their discussion.\nAllele-specific effects: T. brucei is diploid and the current technology, typically edits a single\nallele; heterozygosity may mask fitness costs or resistance phenotypes compared to full\nhomozygous states.\nR3.3: We’ve added the following text to the Discussion: “We do not observe homozygous edits\nusing oligo targeting here, and this may be a limitation for investigating some drug resistance\nmechanisms”.\nHeterozygosity does not appear to mask resistance phenotypes in our study, since\n“heterozygous edits were sufficient to yield resistance in every case” in the bespoke set of\nnine pairs of mutants. Fitness costs, on the other hand, are ‘masked’ by heterozygosity (see\nFig. 4b-c) in the sense that defective alleles are tolerated in a diploid background, but such\ndefective alleles would not be expected to yield resistance since “mutational resistance space\nwill be constrained by mutational fitness space”.\n3.4: Selection biases: Stepwise drug selection can enrich fitter mutants while missing rare but\nrelevant variants; I am hoping this issue might not be of big concern given the extensive efforts\nmade by authors using varying drug concentration and deep sequencing in order to isolate all\nthe possible mutations.\nR3.4: Mutations associated with loss of proteasome function and substantial loss-of-fitness\nunder drug pressure will inevitably have been lost from the population during screening. On\nthe other hand, we’re confident that our approach does not miss relevant variants since\n“codons for the same amino acid edit yielded highly consistent results”.\n3.5: Modeling assumptions: Combining Evo2 functional scores with docking/free-energy\npredictions involves thresholds and normalization choices that may influence classification of\nloss-of-function vs resistance; Very heartening to see significant correlation between virtual\nEC50s and measured EC50s.\nR3.5: We agree and have now used quantitative data from Evo2 predictions and added more\ndetail on this analysis (see R1.3-1.4).\n3.6: Structural inference: Homology modeling and docking rely on a related proteasome\nstructure; uncertainties in local side-chain positioning could affect predicted interaction\nfingerprints.\nR3.6: We do comment on ‘limitations’ associated with computational modelling in the\nDiscussion but also find it reassuring that our models “yield effective in silico predictions of\ndrug target affinity and relative resistance”, meaning that in silico analysis can be fruitful even\nwhen structural data are limited.\n3.7: The current methodology may not be able to capture structural changes that may affect\ncompound binding, if the mutations far (>5Å), thereby leading to compound resistance and\n5\nloss-of-fitness.\nR3.7: The current MOT libraries do indeed incorporate only mutations affecting residues <5Å\nfrom the drug binding site. Although results may differ depending upon the target, our analysis\nof spontaneous mutations in the proteasome b5 subunit indicated that “all [spontaneous] β5\nsubunit mutations were predicted by MOT-library profiling; and no β5 subunit mutations were\nidentified in residues more than 5 Å from the drug binding site”. In terms of loss-of-fitness, and\nas noted above, mutations associated with loss of proteasome function and substantial loss\nof fitness under drug pressure will inevitably have been lost from the population during\nscreening (see R3.4).\n3.8: Overall, the work appears methodologically solid and well-documented, with clear\nresistance hotspots identified. The main caveats are the scope (single\ntarget/compound/species) and potential biases from editing and selection workflows, which\nare common challenges in mutational scanning studies. The MOT and computational\nmodeling has potential not only to enable drug design but also to monitor resistance\nsurveillance.\nR3.8: We thank the reviewer for these comments and certainly intend to expand the scope in\nthe future.\nMinor comments\n3.9: Line 23-24: the sentence ending “however”. May be you can start with “However, ….”\nR3.9: We considered both alternatives, but prefer the original version.\n3.10: Line 171-172: Authors refer to “forty-six distinct drug-resistant mutants, which displayed\nup to 70- fold increased resistance”. I was unable to find this tables. Please refer to Table.\nR3.10: We now refer to “(Supplementary Data sheet 4)”.\nReviewer 1:\nThe authors have addressed my main concerns regarding the computational validation and\nthe methodology. The quantitative analysis is now sufficiently clear. Specifically:\n• The clarification of the virtual EC_{50} calculation as a weighted average is reasonable and\naddresses my concern about curve fitting parameters2.\n• I reviewed the supplementary data provided in response to the selection bias concern. I am\nsatisfied that the high correlation in Figure 3 is robust and not an artifact.\n• I appreciate the authors' transparency regarding Evo2's low standalone predictive power (r\n= 0.01). Please ensure this specific limitation is explicitly stated in the main text so the\nmodel's contribution is not overstated.\nWe’ve added “Pearson’s r for loss-of-fitness alone was 0.01” to the last section in the\nResults text.\nRemarks on code availability:\nThe documentation for the code is improved this time. However, for b) Free Energy\nMMGBSA/Calculations, the current instructions provided are minimal ('run externally using\nthe Schrodinger suite'). It is recommended to include more detailed information.\nWe’ve now added the following detail to the repository: \"Free energy MMGBSA/Calculations\nwere performed for each of the clusters with mutations. For these calculations, we provide\nPDB files from each cluster; however, the calculation needs to be run externally using the\nSchrodinger suite (Residue Scanning Tool, 2023-2). We provide the final csv file with results\n('all_combined_final.csv'). This tool implements a high-throughput workflow for testing\nmutations with the free energy associated (calculated with the MM/GBSA method) using\nwild-type as reference (baseline). For more details, please, refer to the paper associated\nwith this repository.\"\nReviewer 2:\nAll my previous comments have been addressed. I recommend the manuscript for\npublication.\nReviewer 3:\nAuthors have addressed all the comments raised by me. I am satisfied with the changes\nmade.\nRemarks on code availability:\nI am not expert in reading the codes. Hence, I am unable to comment on the code provided\n1"}]

Decision letter