E.coli CSH501 strain carries a mutation in the gyrA gene which abolishes inhibition and poisoning by CcdB (Bajaj et al., 2008). The EBY100 strain of Saccharomyces cerevisiae has the aga1 gene under the Gal1 promoter for inducible expression and a TRP1 auxotrophic mutation. The strain lacks the aga2 gene, so only Aga2p fused protein expressed from the plasmid, will form a complex with the Aga1p for yeast cell surface display (Boder and Wittrup, 2000). The ccdB gene was cloned in the pBAD24 plasmid for controllable expression in E. coli. ccdB mutants were cloned in the pPNLS shuttle vector for yeast cell surface expression (Najar et al., 2017).

ccdB mutants in pBAD24 were generated using three fragment Gibson assembly. Briefly, ccdB was amplified in two fragments using two sets of oligos. For each fragment one of the oligos binds to the vector and the other binds to the gene. The primer of both fragments which bind to the gene were completely overlapping and contained the desired mutation. The fragments were gel extracted and Gibson assembled with NdeI and HindIII digested pBAD24 vector. The Gibson assembled product was electroporated in E. coli CSH501 strain and positive transformants were selected on LB agar media containing ampicillin (100 μg/ml). The sequence was confirmed by Sanger sequencing. Sequence confirmed WT or mutant ccdB in pBAD24 vector was used as a template for PCR to amplify the ccdB gene by Vent DNA polymerase. The PCR amplified product was co-transformed with SfiI digested pPNLS vector in the EBY100 strain of Saccharomyces cerevisiae using LiAc/SS carrier DNA/PEG method for in vivo recombination (Gietz and Schiestl, 2007). Positive transformants were selected on SDCAA Tryptophan dropout media plates and the sequence was confirmed by Sanger sequencing.

WT and mutant CcdB was purified as described previously (Chattopadhyay and Varadarajan, 2019). Briefly, an overnight culture was diluted 100-fold in LB media containing ampicillin (100 μg/ml) and induced with L-arabinose (0.2% w/v) at an OD₆₀₀ of ∼0.5. Following induction for 3 h, cells were harvested and lysed by sonication. The soluble fraction was separated using centrifugation and incubated with CcdA peptide (residues 45–72^nd) coupled to Affigel-15 at 4°C. The unbound fraction was removed and the column was washed with bicarbonate buffer (50 mM NaHCO₃, 500 mM NaCl, pH 8.5). The bound protein was eluted with 200 mM glycine (pH 2.5) and collected in an equal volume of 400 mM HEPES buffer (pH 8) to neutralize the acidity of glycine.

GyrA14 was purified as described previously (Dao-Thi et al., 2004). Briefly, an overnight culture was diluted 100-fold in LB media containing ampicillin (100 μg/ml) and induced with IPTG (1 mM) at an OD₆₀₀ of ∼0.5. Following induction for 3 h, cells were harvested and resuspended in TES buffer (0.2 M Tris, pH 7.5, 0.5 mM EDTA, 0.5 M sucrose and 1 mM PMSF). Cells were lysed and the soluble fraction was separated using centrifugation. The soluble fraction was incubated with pre-equilibrated Ni-NTA beads for 2 h at 4°C. The unbound fraction was removed, and the column was washed with 100 column volumes of wash buffer (50 mM imidazole in 0.05 M Tris, pH 8, 0.5 M NaCl). The protein was eluted with 500 mM imidazole in 0.05 M Tris, pH 8, 0.5 M NaCl and dialysed against 1x PBS.

E.coli CSH501 strain, transformed with pBAD24 plasmid containing WT or mutant ccdB, was grown in media containing ampicillin for 16 h at 37°C and 180 RPM. A secondary culture was grown by diluting overnight grown culture 100-fold. Upon reaching an OD₆₀₀ of 0.4–0.5, CcdB variants were induced with Arabinose at a final concentration of 0.2% (w/v) for 3 h. The cells were harvested from 1.5 ml culture and lysed in 500 µl 1X PBS, using sonication. Supernatant and pellet fractions were separated by centrifugation at 13,000 RPM at 4°C. The pellet fraction was resuspended in 500 µl 1X PBS and equal volumes of pellet and supernatant fractions were loaded on Tricine-SDS-PAGE to measure the relative amounts of protein in each fraction.

The thermal shift assay was conducted in an iCycle iQ5 Real Time Detection System (Bio-Rad, Hercules, CA). A solution of total volume 20 μl containing 10 μM of the purified CcdB protein and 2.5X Sypro orange dye in suitable buffer (200 mM HEPES, 100 mM glycine), pH 7.5 was added to a well of a 96-well iCycler iQ PCR plate. The plate was heated from 15°C to 90°C with a 0.5°C increment every 30 s. The normalized fluorescence data was plotted against temperature and T_m measured as described (Niesen et al., 2007; Tripathi et al., 2016).

Saccharomyces cerevisiae EBY100 cells containing WT ccdB or mutant in pPNLS plasmids were grown in 3 ml SDCAA media (glucose 20 g/L, yeast nitrogen base 6.7 g/L, casamino acid 5 g/L, citrate 4.3 g/L, sodium citrate dihydrate 14.3 g/L) for 16 hours. Grown cells were diluted to an OD₆₀₀ of 0.2 in 3 ml SDCAA media and grown till the OD₆₀₀ reached two. Thirty million cells were harvested using centrifugation and resuspended in 3 ml SGCAA induction media (galactose 20 g/L, yeast nitrogen base 6.7 g/L, casamino acid 5 g/L, citrate 4.3 g/L, sodium citrate dihydrate 14.3 g/L) for 16 hours at 30°C, 250 RPM (Chao et al., 2006). One million cells were used for flow cytometric analysis. The amount of total protein expressed on the yeast cell surface was estimated by incubating the induced cells in 20 μl FACS buffer (1X PBS and 0.5% BSA), containing chicken anti-HA antibodies from Bethyl labs (1:600 dilution) for 30 min at 4°C. This was followed by washing the cells twice with 100 μl FACS buffer at 4°C. Washed cells were incubated with 20 µL FACS buffer containing goat anti-chicken antibodies conjugated to Alexa Fluor 488 (1:300 dilution), for 20 min at 4°C. Fluorescence of yeast cells was measured by flow-cytometric analysis. The total amount of active protein on the yeast cell surface was estimated by incubating the induced cells in 20 μl FACS buffer containing 100 nM GyrA14 for 45 min at 4°C. Cells were washed and incubated with 20 µl mouse anti-FLAG antibodies (1:300). This was followed by washing the cells twice with FACS buffer, followed by incubating with 20 µl rabbit anti-mouse antibodies conjugated to Alexa Fluor 633 (1:1,600 dilution). The flow-cytometric analysis was carried out on BD Accuri or BD Aria III instruments.

Previously, an SSM library of ccdB was generated in the pBAD24 vector (Adkar et al., 2012; Tripathi et al., 2016). The library was PCR amplified using primers having homology to the pPNLS vector. The PCR amplified library was gel extracted and cloned in pPNLS vector using yeast in vivo recombination.

A similar protocol was used for sample preparation of the library for FACS as described above for the single mutants with slight modifications. Briefly, ten million cells were taken for FACS sample preparation and the reagents were used in 10X higher volumes compared to the earlier flowcytometric analysis. Two different concentrations of GyrA14 (100 nM, 5 nM) were used for sorting CcdB mutants based on the binding in the 1D histogram. The cells were sorted in 11 and 10 different populations (bins) in case of binding with GyrA14 at concentrations of 100 and 5 nM respectively. Additionally, 11 different populations (bins) were sorted from the expression histogram. The experiment was repeated in a biological replicate. The sorting of CcdB libraries was performed using a BD Aria III cell sorter.

Sorted populations were grown on SDCAA agar plates for 48 h. Colonies were scraped and plasmids were extracted from the cells. The ccdB gene was PCR amplified using primers which bind upstream and downstream of the ccdB sequence and had multiplex identifier (MID) sequence to segregate the reads from different sorted bins. The DNA was amplified for 15 cycles using PCR and the amplified product was gel extracted and purified. Equal amounts of DNA from each sorted population were pooled, and the library was generated using the TruSeq™ DNA PCR-Free kit from Illumina. The sequencing was done on an Illumina HiSeq 2,500 250 PE platform at Macrogen, South Korea after incorporating 20% ϕX174 DNA in the library.

Deep sequencing data for the ccdB mutants obtained from the Hiseq 2,500 platform was processed using a pipeline developed by adopting certain aspects from an already existing in-house protocol (https://github.com/skshrutikhare/cys_library_analysis). The latter method involved the alignment with wild type sequence followed by merging of the paired-end reads, while in the modified protocol, the reads are first merged and then aligned with the wild-type sequence. The present methodology consists of the following steps: assembling the paired end reads, quality filtering, binning, alignment and mutant identification. All these steps were incorporated in a pipeline and made executable from a single command using a parameter file unique to a given data-set. In the first step, paired end reads were assembled using the PEAR v0.9.6 (Paired-End Read Merger) tool (Zhang et al., 2014). The “quality filtering” step involved deletion of terminal “NNN” residues in the reads, and removal of reads, not containing the relevant MID and/or primers, along with the reads having mismatched MID’s. Finally, only those reads having bases with Phred score ≥20 are retained. A binning step involved further filtering, which eliminated all those reads having incorrectly placed primers, truncated MIDs/primers (due to quality filtering) and shorter/longer sequences than the length of the wild type sequences. The remaining reads were binned according to the respective MIDs. In the alignment step, reads were aligned with the wild type ccdB sequence using the Water v6.4.0.0 program (Smith and Waterman, 1981) and reformatted. The default values of all parameters, except the gap opening penalty, which was changed to 20, was used. In the final step of “substitution”, reads were classified based on insertions, deletions and substitutions (single, double etc mutants).

Reads of each mutant were normalized across different bins individually (Equation 1), and the fraction of each mutant (Xi) distributed amongst the different bins was calculated (Equation 2). The reconstructed MFI for an individual mutant was calculated by the summation of the product, obtained upon multiplying the fraction (Xi) of the mutant in a particular bin 1) with the MFI of the corresponding bin obtained from the FACS experiment (Fi), across the various bins populated by the respective mutant (Equation 3). Normalized read of mutant in   bin     i   ( Ni )   =   No .   of reads of mutant   i   in   bin   i   ∑ reads   in   bin     i Equation 1 Fraction of mutant in each gate ( Xi ) = Ni ∑ 1 n Ni Equation 2 Reconstructed MFI   =   ∑ 1 n F i ∗ X i Equation 3

The MFI_seq of the biological replicates were different so the MFI_seq of one of the replicates was adjusted using “m” and “c” obtained from the correlation between the replicates and then averaged. Average MF I seq   =   MFIseq   ( replicate   1 ) + ( m ∗ MFIseq  ( replicate   2 ) + C ) 2

Reads of each mutant were normalized within and across the bins. The fraction of each mutant (Xi), distributed amongst the different bins, was calculated as explained in the above section. The fraction (Xi) was multiplied with a scaling factor to convert the data into integers as this is required by the package below. The mlMFI was calculated using a maximum likelihood method using the fitdistrplus R package as explained earlier (Starr et al., 2020). The “fitdistcens” function in the fitdistplus R package helps in the estimation of fluorescence values for such observations using a maximum likelihood approach, where the values are transformed into a data frame of two columns left and right, describing each observed value as an interval and assuming a normal distribution of values. The left column contains the left bound of the interval and the right column contains the right bound of the interval for interval-censored observations, based on the fluorescence boundaries of each bin. The maximum likelihood approach was used to estimate the MFI of binding and expression for each mutant, based on its distribution of reads across the sorted bins, and the fluorescence boundaries of each sorted bin.

The bins were merged following which mlMFI amd MFI_seq were calculated for GyrA14 binding (100 nM) for replicate 1. The fraction of each mutant in each bin was calculated as explained in the earlier sections. To merge bins for a given mutant, fractions present in each of the bins to be merged were added arithmetically. For mlMFA calculation, the minimum and maximum fluorescent boundary of the merged bin was set at the lowest and highest value of the fluorescent boundary for that set of bins. The mlMFI of CcdB mutants was calculated as explained above. In the case of MFI_seq, the mean fluorescent intensity of merged bins was determined by making a new bin spanning the set of merged bins. The MFI_seq of CcdB mutants was then calculated as explained above.

Depth was calculated using the server DEPTH (Chakravarty and Varadarajan, 1999; Tan et al., 2011). Accessibility was calculated using the program NACCESS (Hubbard SJ, 1993). In both cases, the input co-ordinates were homodimeric CcdB (PDB ID 3VUB). RankScore and MS_seq are measures of mutational sensitivity in E. coli. Values were obtained from Adkar et al. (Adkar et al., 2012). Buried residues were those with <10% accessibility in 3VUB. Active-site residues were those with ΔASA>0. ΔASA difference between the solvent accessible surface area of CcdB residues in the free (3VUB) and GyrA14-bound forms (1X75) respectively (Aghera et al., 2020).

The deep mutational scanning data was taken from a recent report (Starr et al., 2020) in which two independent libraries of RBD were generated and sorted in four different bins based on expression or binding to ACE-2. In the MFI of binding and expression for individual mutants was reconstructed in that study using a maximum likelihood method using fitdistrplus R package. The expression MFI [Sortseq (expr)] data was shared by the authors in a repository (https://github.com/jbloomlab/SARS-CoV-2-RBD_DMS). We reconstructed the binding MFI [Sortseq (bind)] at an ACE-2 concentration of 100 pM (TiteSeq_09). For Sortseq (bind) estimation we used the script provided by the authors (https://github.com/jbloomlab/SARS-CoV-2-RBD_DMS/blob/master/results/summary/compute_expression_meanF.md). The authors used data from both single and multiple mutants, together with a model to account for epistatic effects to infer the MFI values for individual mutants. We modified the script to change the input data required to calculate Sortseq (bind). For both Sortseq (bind) and Sortseq (expr), we analyzed only single mutant data to avoid any artifacts that might arise from the epistatic model and took the average of delta Sortseq MFI {log (Sortseq (WT))—log [Sortseq (mutant)]} of mutants which had multiple barcodes. The Sortseq MFI values of mutants were averaged between the two libraries and the antilog was calculated for delta Sortseq MFI to analyse the ratio of Sortseq (bind) or Sortseq (expr) of mutants with respect to WT.

The correlations and p values for its significance were calculated using the GraphPad Prism software 9.0.0 (* indicates p < 0.05, ** indicates p < 0.01, **** indicates p < 0.0001). The weighted correlations were calculated using the weights function of R. For the computation of weighted correlation, a weight of 1/(σ/µ) was used on the mean values of replicates.

Yeast surface display (YSD) has become an increasingly popular tool for protein engineering and library screening applications (Pepper et al., 2008). Aga2p mating adhesion receptor of Saccharomyces cerevisiae is used as a fusion protein for yeast surface display. For surface expression, we used a vector in which CcdB is fused at the C-terminus of Aga2 (Sahoo et al., 2015). We generated (Supplementary Figure S1) and individually characterized 30 CcdB variants on the yeast cell surface. Most CcdB mutants had similar levels of expression to the WT protein (Figure 1A). However, the mutants showed different amounts of active protein as assayed by binding to the FLAG tagged GyrA14 compared to the WT protein (Figure 1B). Previously, we have characterized the in vitro thermal stability and in vivo solubility of several CcdB mutants (Tripathi et al., 2016). The amounts of total and active protein were estimated using antibodies against the HA-tag at the N-terminal of the yeast surface displayed CcdB and the C-terminal FLAG tag of GyrA14 respectively. The correlation coefficient (r) between amount of total protein on the yeast cell surface with in vivo solubility or T_m of the corresponding purified protein were 0.31 and 0.70 respectively (Figures 2A,B). It is unclear why mutants which have very low solubility in E. coli are highly expressed on the yeast cell surface. It was previously hypothesized that the protein folding quality control system in yeast is not as effective as in mammalian systems, therefore partially folded/molten globule/aggregated protein may exist on the surface of yeast (Park et al., 2006). A correlation of r = 0.80 was found between the amount of active protein on the yeast cell surface with its in vivo solubility determined in E. coli (Figure 2C). We also found a better correlation (r = 0.90) between amount of active CcdB protein on the yeast cell surface and its in vitro thermal stability (Figure 2D), compared to that between total CcdB protein on the yeast cell surface and thermal stability.

To extend these results, an SSM library of ccdB was expressed on the yeast cell surface. Different populations based on extent of binding to gyrase or cell surface expression were sorted. A total of 32 different populations were sorted at two different concentrations of GyrA14 (100 nM, 5 nM) as a function of either surface expression level or the extent of binding to GyrA14 (Supplementary Figure S2). The lower concentration of GyrA14 was chosen to be around the K_D of CcdB-GyrA binding (Supplementary Figure S3), the higher concentration was one where WT CcdB approaches saturation in binding with GyrA14 on the yeast cell surface. We hypothesized that at lower concentrations of GyrA14, the binding on the yeast cell surface will be a function of both stability as well as binding affinity. However, at saturating concentration of GyrA14, the binding on the yeast cell surface will largely be a function of amount of correctly folded protein that in turn might be a function of protein stability, rather than the K_d of the mutant(s). MFI was calculated for each mutant as explained in the Methods section. The MFI was calculated at different stringencies (where the stringency refers to the sum of reads for a given mutant over each gate of the histogram), namely 25, 50, 100, 150, and 200 reads. All mutants with a total read number less than the stringency value were removed from the analysis. As the stringency increased, the pairwise correlation between the biological replicates increased (Supplementary Figure S4, Supplementary Table S1). The data was analysed with a stringency of 50 reads, since at higher stringencies, correlation did not improve significantly, but the number of mutants reduced. Reconstructed Binding and Expression MFI from deep sequencing data are hereafter referred to as MFI_seq (bind) and MFI_seq (expr) respectively.

A few published studies have described estimation of MFI values using deep sequencing of sorted populations and are therefore similar to our experimental strategy. However, the procedure for MFI reconstruction in these reports was relatively complicated compared to that used here (Sharon et al., 2012; Noderer et al., 2014; Peterman and Levine, 2016; Cambray et al., 2018). In those studies, the fractions of reads were calculated in each bin for all the mutants and MFI (mlMFI) of mutants were calculated by fitting the data to a maximum likelihood distribution of the histogram. We found that if mutants are present in only one bin (highly destabilized or nonsense mutants) then this method is unable to perform the MFI calculation (Starr et al., 2020). For the remaining mutants we found a good correlation between MFI_seq and mlMFI for binding at 5 and 100 nM GyrA14, and for expression (Supplementary Figure S5). For mutants with over 50 reads, we could calculate the MFI of 11,153 mutants using the maximum likelihood method and 11,436 mutants using our method. We also found that progressively reducing the number of bins from eleven to six, does not significantly affect the estimated MFI values, however a further reduction to four bins results in a noticeable change in the estimated values using either method (Supplementary Figure S6). A good correlation was also found between the MFI of individually analysed mutants and their corresponding MFI_seq values, validating our approach of MFI reconstruction (Supplementary Figure S7A, 7B). Individually analysed mutants showed a good correlation between the amount of active protein on the cell surface and in vitro measured thermal stability of the purified protein. Similarly, we also found a good correlation between MFI_seq (bind) of mutants inferred from deep sequencing, and thermal stability as well as in vivo solubility for the selected mutants (Supplementary Figure S7C, 7D).

For the exposed residues (>10% accessibility) (Supplementary Figure S8), mutations did not affect the degree of surface expression and binding to GyrA14 (Figures 3A,B). Expression was also unaffected by mutations in the active-site residues (identified from PDB ID:1X75) (Figure 3C, Supplementary Figure S8). However, many buried site mutants showed very low expression, possibly because of aggregation and degradation inside cells or during export (Figure 3C). In the case of binding for buried and active-site residues, a very high mutational sensitivity was found (Figure 3D) similar to the previous report of CcdB mutants in E. coli (Tripathi et al., 2016). We also found a very high mutational sensitivity of binding for a few non-interacting residues in the loop connecting beta strands S2 and S3 at both 5 and 100 nM GyrA14 concentration (Supplementary Figure S9). The residues I24, I25 and D26 in this loop are directly involved in interacting with Gyrase and mutation at non-interacting residues (22, 23 and 27) in the loop might restrict or alter the conformation of the loop, thus reducing the affinity of CcdB mutants to GyrA14. However, there was no effect on the expression of the mutants in this loop, indicating that the mutant proteins are not destabilized (Supplementary Figure S9). We did not find a high correlation between MFI_seq (bind) and either accessibility or depth, because many mutations at both buried and active-site residues have high mutational sensitivity (Supplementary Table S2). The previously described parameter RankScore, is a measure of mutant activity in E. coli (Adkar et al., 2012) with high RankScore denoting lower activity. We found a poor correlation between the MFI_seq (bind) values of CcdB mutants at both exposed non active-site as well as active-site residues, and RankScore. In E. coli, most of the exposed non active-site residues do not show any mutational sensitivity, i.e., they have the same RankScore values as WT. However, in the present case many such CcdB mutants show lower binding to GyrA14 compared to WT. The loss of binding could be attributed to the decrease in the affinity between CcdB and Gyrase, or destabilization due to mutation. We defined a new parameter MrMFI (mean residue MFI) which is the mean of the MFI values of all the mutants at a certain position. MrMFI (expr) and MrMFI (bind) at 100 nM GyrA14, show a good correlation with RankScore (Supplementary Table S2). MrMFI (expr) also showed good correlation with Depth which is a structural measure of residue burial (Chakravarty and Varadarajan, 1999). However, in the case of binding at 5 nM, a weaker correlation of MrMFI (bind) with the aforementioned parameters was observed (Supplementary Table S2). In previous studies, identification of the active-site residues solely from the deep sequencing data was not very efficient (Adkar et al., 2012; Bhasin and Varadarajan, 2021), this is presumably because in vivo activity is often governed by threshold effects, and because mutations at buried residues also affect activity. The current methodology removes such drawbacks. We could distinguish between buried and active-site residues by comparing the MFI_seq (bind) and MFI_seq (expr). Most buried site residues showed low values of both MFI_seq (bind) and MFI_seq (expr) compared to WT. However, the active-site residues showed low MFI_seq (bind) but similar MFI_seq (expr) compared to WT. We found that the average MFI_seq values of charged residues are a good predictor to discriminate between buried and active-site residues. For calculating MrMFI_charged of charged WT residues, we only consider mutants with opposite charge. For some mutants at buried positions, we found a very low MrMFI_charged (expr) but the mutants were absent in MrMFI_charged (bind). We found that such mutants had very high reads, suggesting that the values of MrMFI_charged (expr) are correct. We anticipated that such mutants lack binding and are therefore present only in the bin which had a background level of binding signal, the presence of mutant in only that gate led to the removal of such mutants due to the stringency set for the analysis. Hence, such mutants were assigned a MrMFI_charged (bind) similar to other buried positions. MrMFI_charged had a bimodal distribution (Supplementary Figure S10), so k-means clustering was performed to identify the mean (µ) and standard deviation (σ) of each distribution. The distributions were named D1 (higher MrMFI_charged) and D2 (lower MrMFI _charged). Buried site residues were assigned to be those which have MrMFI_charged (bind) and MFI_seq (expr) less than the set threshold (µ+0.5*σ) for distribution D2. Active-site residues were assigned as those which had MrMFI_charged (bind) less than (µ+σ) of the D2 distribution and MFI_seq (expr) higher than (µ−2*σ) of distribution D1 (Figure 4). The accuracy, specificity and sensitivity of prediction of exposed non active-site, buried and exposed active-site residues are mentioned in Supplementary Table S3. We also compared our prediction results derived from saturation mutagenesis phenotypes with those of an in silico predictor, PROF (Rost and Sander, 1994). For a residue to be classified as buried by PROF, the relative solvent accessibility cut-off used is < 12. We observed a slightly lower specificity and accuracy for CcdB, and lower sensitivity in the case of RBD when predictions were made using PROF (Supplementary Table S4), relative to our predictions. We also examined the performance of PROF with other proteins and found that the specificity of the predictions was higher than 0.8 in all the cases except for CcdB. However, the sensitivity of the predictions was lower than 0.8 in all the cases except for CcdB, Gal4 and Ubiquitin. The accuracy for the PROF prediction was 0.77 and 0.78 for CcdB and RBD respectively, comparable but slightly lower than the corresponding values of 0.92 and 0.8 for CcdB and RBD respectively, from the saturation mutagenesis predictions in this work.

In the previous section, we discussed the correlation between protein biophysical properties such as thermal stability and in vivo solubility with either the amount of active protein or the ratio of active protein to total protein on the yeast cell surface for a few (30) mutants. However, most of these mutants were destabilized with respect to the WT protein. To confirm whether this correlation also holds for mutants that have stability similar or greater than WT, we selected a few CcdB mutants based on either the MFI_seq (bind) or MFI_seq (ratio) [MFI_seq (bind)/MFI_seq (expr)] for in vitro characterization of thermal stability. We examined the average and standard deviation of expression for all mutants and selected only those mutants based on MFI_seq (ratio) which cross a minimum cut-off (µ+0.5*σ) for MFI_seq (expr) to remove the bias created by mutants which have very low expression. No threshold for expression was set for selection of mutants based on their MFI_seq (bind). No selection of the mutants was performed based solely on the MFI_seq (expr).

Six mutants were characterized using the criteria MFI_seq (bind) at 5 nM GyrA14, none of them showed a higher T_m than WT (Figure 5A); whereas two of the mutants selected on the basis of MFI_seq (ratio) showed a significantly higher T_m than WT (Figure 5B). A subset of seven mutants was selected based on MFI_seq (bind) at 100 nM GyrA14, none of the mutants showed higher stability than WT CcdB (Figure 5C). Ten mutants were selected based on MFI_seq (ratio) and characterized, four showed higher stability, two mutants were similar to WT and the remaining four were less stable than WT CcdB (Figure 5D). We therefore hypothesize that if the stability of a mutant crosses a threshold then its expression will not increase further. To confirm this hypothesis, we measured the amount of active protein on the yeast cell surface for seven individual mutants which had T_m’s ranging from 60°C to 70°C, and found that the expression and binding for these mutants are similar to each other and to WT (Supplementary Figure S11).

For destabilized mutants we observed a good correlation between MFI_seq (bind) and T_m of individual mutants (Supplementary Figure S7D). Using this correlation, we next predicted the T_m of each mutant for an additional set of (n = 28) previously described CcdB mutants (Tripathi et al., 2016) based on their MFI_seq (bind). We found a good correlation (r = 0.82) between predicted and in vitro measured T_m for this set of CcdB mutants as well (Supplementary Figure S12A). This now allows us to identify putative destabilized mutants and accurately predict the extent of destabilization for all such mutants in the CcdB YSD library. We also predicted the thermal stability of CcdB mutants using the in silico predictor, HoTMuSiCv1.0 (Pucci et al., 2020), however, we did not find a good correlation between measured and predicted T_m (Supplementary Figure S12B). It has been shown that in vitro protein thermal stability and free energy of unfolding are correlated (Chen et al., 2000; Prajapati et al., 2007; Tripathi et al., 2016). We therefore predicted the free energy of unfolding for CcdB mutants using SDM (Pandurangan et al., 2017), mCSM (Pires et al., 2014b), PoPMuSiC (Dehouck et al., 2011), DynaMut (Rodrigues et al., 2018), DUET (Pires et al., 2014a), MAESTROweb (Laimer et al., 2016), DeepDDG (Cao et al., 2019), CUPSAT (Parthiban et al., 2006), PremPS (Chen et al., 2020) and INPS-MD (Savojardo et al., 2016). We found moderate correlations, with DeepDDG performing the best (r = 0.59), but still poorer compared to our prediction from YSD data (r = 0.82). For a more detailed comparison we analysed the predictions of stability by DeepDDG, since this showed the highest correlation with measured stability of individual mutants at non active-site residues. We excluded residues 21, 22, 23 and 27 as these positions behaved like active-site residues. We found that trends for ΔΔG predicted by DeepDDG for exposed non active-site residues are similar to those obtained from MFI_seq (bind) (Figures 6A,B). However, we observed some mutant specific differences at residues 8, 16, 50, 53 and 96. Mutations at residues 50 and 96 have highly deleterious effects which reduced GyrA14 binding to yeast surface displayed protein, these are only partially predicted by DeepDDG. In the case of charged and polar mutations at residue 8, 16 and 53 we did not observe a reduction in binding, but the software predicted them to be destabilizing. In the case of buried positions, we found mutation specific effects at 35, 52 and 94 where DeepDDG predicted changes were significantly smaller than the experimentally observed ones. We also found that most of the phenylalanine, tryptophan and arginine mutations were highly destabilizing and the mutants did not bind to GyrA14, however the software gave a lower stability penalty for these substitutions (Figures 6C,D). Our MFI based measurements suggested greater destabilization for several mutants relative to DeepDDG prediction. While the overall trends were similar, as discussed above, there are several differences between MFI based and DeepDDG based stability predictions.

To examine the generality of our approach, we also analyzed recently reported deep mutational scanning data of the SARS-CoV-2 receptor binding domain (Starr et al., 2020). In this study two separate libraries were generated and individually sorted based on expression and binding to ACE-2. The binding [Sortseq (bind)] or expression [Sortseq (expr]) MFIs relative to WT for barcoded mutants were calculated from the deposited NGS data as explained in the Methods section. Additionally, we analyzed binding at only one concentration of ACE-2 (100 pM, TiteSeq_09) at which the binding started to saturate. Buried residues were those with <10% side chain accessibility in chain C of PDB ID 7KMH (Jones B. E. et al., 2020). ACE-2 binding (active-site) residues were assigned as those contacting ACE-2 (Malladi et al., 2021). To identify the active-site and buried residues from Sortseq data, we calculated the MrMFI_charged for each position. Similar to CcdB, we observed a bimodal distribution for both MrMFI_charged (bind) and MrMFI_charged (expr) (Supplementary Figure S13) and k-means and standard deviation were calculated for both the distribution D1 (higher MrMFI_charged) and D2 (lower MrMFI_charged). As described above for CcdB, buried residues were identified as those which had MrMFI_charged (bind) and MrMFI_charged (expr) less than the set threshold (µ+0.5*σ) for distribution D2. The active-site positions were identified as those which had MrMFI_charged (bind) lower than the set threshold (µ+σ) for population D2 and MrMFI_charged (expr) values higher then (µ-2*σ) for population D1. We accurately identified most of the buried residues, however there were some false positive and false negative predictions relative to the crystal structure information (Figure 7). We found 21 positions to be false negative buried positions. We categorized these false negatives into two categories, namely, glycine and the side chains which are pointing towards the surface. The accessibility calculated by DEPTH server for glycine was zero and we therefore expected glycine to fall into the false negative buried category. Thirteen positions out of twenty-one false negative were glycine. Another six positions, 336, 348, 361, 443 and 480 had their side chains pointing towards the protein surface. We also found similar false negative buried residues in CcdB where the side chain hydrophilic group was pointing towards the protein surface. Position 363 and 365 in RBD had accessibility <10% and were pointing towards the core of the protein in the PDB (7KMH) used to calculate accessibility. However, we found that these positions have high accessibility (>30%) in another structure (PDB ID 7D2Z). All the available RBD structures are in complex with other molecules, this might be responsible for variation in the accessibility of residues in different RBD structures. We found 17 false positive buried residue predictions, seven of them were aromatic, seven are charged or polar, two are prolines and one is an aliphatic residue. These positions have both reduced expression and binding for charged residue substitutions (Supplementary Figure S14A, 14D) similar to the buried residues (Supplementary Figure S14B, 14E). The specificity, sensitivity and accuracy of prediction is mentioned in Supplementary Table S3. Active site residues were identified with very high accuracy (Supplementary Table S3), though there were a few false negative and false positive predictions. Additionally, we found several positions which had Sortseq (expr) like WT, however, they had very low Sortseq (bind) (Supplementary Figure S14A, 14D). We hypothesize that these positions are also assisting in the maintenance of proper RBM conformation and enabling its binding to ACE-2. Residues 447, 448, 473 and 476 which gave false positive results, 447 and 476 are part of the receptor binding motif (RBM) and contain glycine in a conformation which is available only for glycine. Hence mutation to a non-Gly residue will likely disrupt the conformation of the RBM thus decreasing binding to ACE-2. Mutations at positions 446, 453, 493 and 498 gave false negative results. Of these false negative positions, 446 is again glycine. We found that the Arg mutants at N493 and N498 positions have very little effect on expression and binding (Supplementary Figure S14C, 14F). We hypothesized that these positions may not have the most optimal WT residue, or they may show no mutational penalty for binding to ACE-2. A recent report showed that the affinity of Q498R to ACE-2 is higher than WT RBD (Xue et al., 2020) and was enriched as double mutant Q498R/N501Y when selection was performed for RBD mutants having high affinity towards ACE-2 (Zahradník et al., 2021). It has also been reported that when chimeric virus evolved in the presence of neutralizing antibodies C121 and C141, this enriched for the Q493R mutation. The mutant virus grows to high PFU titers similar to WT, and infectivity is also inhibited by a chimeric ACE-2 analog, similar to WT (Weisblum et al., 2020). The specificity, sensitivity and accuracy of prediction is mentioned in Supplementary Table S3.