Deep Mutational Scans as a Guide to Engineering High Affinity T Cell Receptor Interactions with Peptide-bound Major Histocompatibility Complex*

Proteins are often engineered to have higher affinity for their ligands to achieve therapeutic benefit. For example, many studies have used phage or yeast display libraries of mutants within complementarity-determining regions to affinity mature antibodies and T cell receptors (TCRs). However, these approaches do not allow rapid assessment or evolution across the entire interface. By combining directed evolution with deep sequencing, it is now possible to generate sequence fitness landscapes that survey the impact of every amino acid substitution across the entire protein-protein interface. Here we used the results of deep mutational scans of a TCR-peptide-MHC interaction to guide mutational strategies. The approach yielded stable TCRs with affinity increases of >200-fold. The substitutions with the greatest enrichments based on the deep sequencing were validated to have higher affinity and could be combined to yield additional improvements. We also conducted in silico binding analyses for every substitution to compare them with the fitness landscape. Computational modeling did not effectively predict the impacts of mutations distal to the interface and did not account for yeast display results that depended on combinations of affinity and protein stability. However, computation accurately predicted affinity changes for mutations within or near the interface, highlighting the complementary strengths of computational modeling and yeast surface display coupled with deep mutational scanning for engineering high affinity TCRs.

The process of increasing the affinity of a protein occurs naturally with antibodies, where somatic mutation within the variable region genes is followed by antigen-driven selection of B cells that express membrane-bound antibodies. In contrast, T cell receptors (TCRs) 3 do not undergo somatic mutations and bind to their antigen, a peptide-MHC (pepMHC), with low (micromolar) affinities. However, improvements in TCR affinity to the same levels of antibodies can be achieved by in vitro approaches involving the generation of mutant TCR libraries followed by antigen selection (1)(2)(3).
For therapeutic purposes, the affinity of a variety of proteinprotein interactions, and especially antibody-antigen interactions, has been enhanced using in vitro directed evolution approaches, including phage, yeast, ribosomal, and mammalian display (e.g. see . These methods rely on the generation of large libraries of mutants at residues within the proteinprotein interface, followed by several rounds of selection for desired parameters (such as affinity, stability, and expression levels) (8,9).
Although directed evolution using larger degenerate libraries has been successful, the most recent techniques involving deep sequencing of single-codon libraries have the potential both to provide mechanistic structural information about a binding site and at the same time to provide leads for affinity improvements. Sequence fitness landscapes have successfully been utilized to map protein-DNA interactions (10), protein-peptide interactions (11), and protein-protein interactions (12). Furthermore, using a PDZ protein domain as a model system, McLaughlin and colleagues were able to manipulate ligand-binding specificity through key mutations identified using sequence fitness landscapes (13). Additionally, a sequence fitness landscape of an influenza-binding protein inhibitor was used to enhance affinity and redirect specificity toward a single H1 hemagglutinin subtype (14).
A major goal in protein engineering is to be able to accurately identify mutations that yield improvements in stability or affinity. In addition to directed evolution and sequence fitness landscapes, when the protein structures are known, structure-based computational design has also been used to achieve these goals. Although there have been inspiring successes (reviewed in Ref. 15), advances in computational approaches require a thorough understanding of the rela-tionships between protein structural and physical properties as well as increases in the ability to rapidly and accurately sample different conformational and configurational states (16).
In the present study, we focused on TCRs because they have evolved to bind to a diverse repertoire of clinically relevant targets and thus represent a class of molecules with significant therapeutic potential. In addition, because of their naturally low affinities, they represent protein engineering targets for both stability and affinity. Previously, we reported successful affinity engineering of TCRs by directed evolution using yeast display (2,(17)(18)(19)(20)(21) and mammalian cell display (7,22). We also described structure-guided design strategies that estimated the binding energies of both favorable and unfavorable mutations and led to the design of additional high affinity TCRs (23,24).
More recently, we reported the use of single-codon libraries with two different TCRs to generate sequence fitness landscapes that allowed analysis of the impact of each residue on binding to their cognate peptide⅐HLA-A2 complexes (25). Sequence fitness landscapes offer a powerful perspective on protein-protein interactions not available from structural data alone by experimentally determining, on a residue-by-residue basis, which amino acids contribute to binding as well as the optimal amino acids at each position (11,26). Accordingly, the two higher affinity TCRs, A6-c134 and RD1-MART1 HIGH , that are specific for Tax⅐HLA-A2 and MART1⅐HLA-A2, respectively, were examined both structurally and by deep scanning mutagenesis to determine the basis of specificity and binding (25).
Here we further studied the mutations that were highly enriched in the sequence fitness landscape of the cancer antigen-specific TCR, RD1-MART1 HIGH , interacting with its target peptide⅐MHC by yeast surface display. We demonstrate that the mutations that exhibited the highest levels of enrichment acted both individually and in synergy to significantly increase the affinity and yeast surface levels of the RD1-MART1 HIGH TCR. We also compared strategies that involved individual combinations of mutations versus construction and selection of multicodon libraries, both based on the sequence fitness landscape. These approaches were further compared with computationally predicted affinity-enhancing mutations using a previously described in silico approach (24,27). Some of the most highly enriched residues identified in the sequence fitness landscape were successfully identified by the in silico approach, but these were not distinguished from many other substitutions that did not yield substantial enrichment. This is probably due in part to the sensitivity of yeast surface display to not only changes in binding but also changes in protein stability. The in silico approach nonetheless performed well when focused on mutations experimentally shown to improve binding and proved advantageous for providing structural interpretations. Overall, we show that deep sequencing methods combined with yeast display provides for significant opportunities for enhancing TCR affinity while also controlling for impacts on stability and demonstrate the utility of computational modeling in providing structural interpretations for affinity gains.

Approaches to Engineering a TCR-MHC Interaction for
Higher Affinity-The A6 TCR was used previously as a template for directed evolution of a novel peptide specificity, whereby the specificity was switched from Tax⅐HLA-A2 to the completely distinct cancer peptide, MART1⅐HLA-A2 (21). The resulting TCR, called RD1-MART1, was further affinity-matured to yield the RD1-MART1 HIGH TCR, with a K D value of about 100 nM. Recently, we generated a sequence fitness landscape of RD1-MART1 HIGH by deep sequencing analysis of antigen-selected, single-codon libraries (25). For the sequence fitness landscape and the work conducted in the present study, RD1-MART1 HIGH and various mutants were displayed on the surface of yeast as single-chain TCRs containing the V␤ and V␣ domains linked by a flexible linker (Fig. 1A). HA and c-Myc tags were attached at the N terminus and C terminus, respectively.
To understand the most expeditious strategy to generate higher affinity binding leads from a protein, here we explored three approaches. The first two were based on the substitutions from the single-codon libraries that yielded increases in the number of antigen-selected TCR, characterized as enrichment values (Fig. 1B). In the third, the results of these experimental approaches were compared with an in silico computational approach using the structure of the RD1-MART1 HIGH ⅐ MART1⅐HLA-A2 complex to predict the structural consequences and the energetic impacts of every individual substitution in the CDRs.
The individual enrichment values for the top substitutions of the RD1-MART1 HIGH TCR are plotted in Fig. 1C. To explore the robustness of these results from the sequence fitness landscape, we tested eight individual mutants (arrows in Fig. 1C) within this range of enrichment values, from 10-to 51-fold (T91Q␣, A99Y␤, V50D␤, G28E␣, S51G␣, Y50W␣, N28G␤, and K97Q␣). The top binders among this group were then tested in combination to determine whether there was binding additivity (Fig. 1B, left).
It is possible that the individual mutation approach, followed by combination of improved mutations, will not identify those collections of mutations that act additively or cooperatively to yield the highest affinity TCRs. Furthermore, positions where there were multiple enriched substitutions (e.g. positions 91␣ and 99␤) would not be fully sampled by analyzing only eight of the individual mutations. Accordingly, we took an alternative approach in which a combinatorial library at four residues with the highest enrichment values was generated. This combinatorial library was sorted with MART1⅐HLA-A2 by FACS (Fig. 1B, right) to yield clones that bound with higher affinity.
Site-specific Mutants Act Individually and in Combination to Enhance Affinity and Expression-Based on the sequence fitness landscape of RD1-MART1 HIGH , 8.3% of the mutations showed enrichment above the wild-type TCR (25). As indicated, eight mutations with enrichment values ranging from 10-to 51-fold (Fig. 1C, arrows) were expressed individually on the surface of yeast and tested for binding to the MART1⅐HLA-A2 ligand. Interestingly, the four residues with the highest enrichment values were each in different CDR loops ( Fig. 2A), raising the possibility that the combined mutations could yield additive effects. Mutants on yeast were stained with various concentrations of MART1⅐HLA-A2 monomer to determine affinities ( Fig. 3A and Table 1). Previous studies with TCR domains or scFv fragments have shown that K D values determined by such flow cytometry-based titrations are similar to K D values determined for the soluble proteins by surface plasmon resonance, at least in the nanomolar K D range (e.g. see Refs. 18, 28, and 29).
In some cases, we also stained yeast with an excess of MART1⅐HLA-A2 tetramer to compare relative surface levels of the TCR mutants (Table 1) because surface levels of properly folded TCRs in the yeast display system are known to be influenced by stability of the proteins (18, 30 -36), and the multimeric form would be able to assess total surface levels even with lower affinity mutants, due to avidity (37). We also used high levels of the multimeric ligand for this purpose, rather than antibodies to the C-terminal c-Myc tag, because the yeast quality control machinery may allow expression of full-length fusions even if misfolded or partially folded domains are pres-ent (38). This would allow us to address the possibility that the enrichment of TCRs with mutations in the single-codon libraries was due to enhanced affinity, higher surface levels of the folded protein, or both.
Of the eight mutations, five showed significant improvements in affinity compared with the RD1-MART1 HIGH TCR ( Fig. 3A and Table 1). The five mutants exhibited 3-fold (Y50W␣) to 11-fold (A99Y␤) increases in affinity above the wild-type protein. Although the other three mutants did not exhibit affinity increases, they did appear to be expressed at higher surface levels, which probably accounts for the enhancement values observed in the deep mutational scan.
To test whether the mutations could be combined to yield further improvements, four double mutants and a triple mutant were generated (Figs. 3 and 4 and Table 1). The double mutants yielded increases of 10 -20-fold, whereas the triple mutant yielded an increase of almost 100-fold, to an estimated K D value of 2 nM. Furthermore, combination of these mutations also led to increased TCR surface levels ( Table 1), indicating that the FIGURE 1. Methods for identifying affinity-enhancing mutations of RD1-MART1 HIGH . A, schematic of the single-chain TCR (scTCR) of RD1-MART1 HIGH with HA and c-Myc tags at the N terminus and C terminus, respectively. B, diagram of steps used in different methods for identifying affinity-enhancing mutations of RD1-MART1 HIGH . C, rank-ordered mutations from the RD1-MART1 HIGH sequence fitness landscape with enrichment scores of Ն10 are shown with their respective enrichment value. Mutations tested for binding by yeast display are indicated with a red arrow. mutations also may operate by stabilizing the TCR domains. A plot of affinities (K D values) versus enrichment value (Fig. 4C) supports the view that the highest enhancement values in a deep mutational scan can be used to reliably generate higher affinity variants of a protein. Indeed, when the eight mutants were viewed as changes in binding free energy relative to "wildtype" RD1-MART1 HIGH , there was a linear relationship between enrichment and ⌬⌬G 0 , with a correlation coefficient (R) of Ϫ0.82 (Fig. 4C, inset). This relationship is also apparent based on the finding that the three top mutations could be combined to yield an additive effect of 100-fold higher affinity.

Selection of a Combinatorial Library Yields Further Improvements and Correlates with the Sequence Fitness Landscape-To
compare with individual site-specific mutagenesis, we generated a combinatorial library at the four residues (28␣, 91␣, 50␤, and 99␤) with the most highly enriched mutations in the deep mutational scan (Figs. 1C and 2A). Selection of this library could yield improved combinations that were not identified in the one-by-one mutational approach. The combinatorial library was sorted three times with MART1⅐HLA-A2 monomer at concentrations ranging from 1 M to 1 nM, and the top 1% fluorescent population was collected (Fig. 5A). After the third sort, eight clones were sequenced, and there were three unique sequences identified (represented by clones S3-2, S3-7, and S3-8) ( Table 2). Titrations of these three clones were performed and compared with the triple mutant described above (Fig. 5, B and C). Mutants S3-7 and S3-8 were similar in affinity to the triple mutant, but clone S3-2 consistently appeared to have a slightly higher affinity and a higher surface level (Table 1). Based on this result, the combinatorial library/selection approach can yield modest but significant improvements above the single mutation approach (e.g. when there are positions where multiple substitutions yield high enrichment values).
The sequences of the mutants were interesting when considered with regard to the sequence fitness landscape and the triple mutant. For example, the T91Q␣ mutation had one of the highest affinities and expression levels among individual mutants (Table 1), and this mutation was present in all three clones from the combinatorial library. Position 50␤ had either an aspartic or glutamic acid in all three clones isolated, consistent with the nearly 4-fold improvement of V50D␤ (Table 1) and the deep mutational scan results (Fig. 1C). Although all three clones identified in the combinatorial library did not contain the A99Y␤ mutation that showed the highest affinity individually (Table 1), each had a different highly enriched mutation (A99D␤, A99E␤, and A99V␤) (Fig. 1C). Interestingly, in clones S3-2 and S3-8, the Gly-28␣ position was mutated to arginine despite this residue having a strong preference for glutamic acid in the deep mutational scan (Fig. 1C) and the improvement in affinity of the G28E␣ mutant. In fact, the G28R␣ mutation was not enriched compared with wild type in the deep mutational scan, suggesting that this mutation may have operated only in the context of one or more of the other three CDR mutations.  NOVEMBER 18, 2016 • VOLUME 291 • NUMBER 47

Engineering Protein-Protein Interactions
In conclusion, the combinatorial approach yielded the greatest improvement, as shown by the S3-2 mutant, which exhibited a K D value of 0.7 nM, an increase in affinity of about 250fold compared with RD1-MART HIGH . The improvement in affinity of the S3-2 mutant was also validated by a separate approach using soluble single-chain TCRs refolded from Escherichia coli inclusion bodies. Previous surface plasmon resonance measurements of soluble RD1-MART HIGH TCR yielded a K D value of 250 nM (21). Here we examined soluble S3-2 TCR by surface plasmon resonance, which yielded a K D value of 2.2 nM (Fig. 5D).
Evaluation of the Sequence Fitness Landscape and Individual Mutations via an in Silico Approach-We next evaluated the same residues from the sequence fitness landscape computationally, examining all 20 amino acids at the 39 CDR loop positions explored in the RD1-MART1 HIGH TCR. We used our previously described in silico structure-guided design approach and scoring function developed and trained using TCR-pep-MHC binding data (24,39). (Note that values of Ͻ1 are consid-ered candidates for improved affinity). There was a general, qualitative agreement between prediction and experiment, because the majority (92%) of the mutations predicted to weaken affinity were also selected against (Fig. 6A). However, 25% of the negatively enriched mutations were also predicted to improve binding.
The majority (78%) of the negatively enriched mutations predicted to improve binding were mutations that reduced polarity, with 62% involving introduction of a hydrophobic amino acid. This suggests that some of the discrepancy between computational prediction and sequence enrichment can be attributed to changes in protein stability. As noted above, this is a well characterized selection parameter in yeast display (18, 30 -36) but is not an output of a computational design approach focused solely on improving affinity. A good example is found in the mutation of Ala-99␤. Mutation of Ala-99␤ to Tyr improved affinity, as shown in Table 1. Mutation of Ala-99␤ to Trp was predicted to have an even more substantial impact on affinity because the modeling indicated the Trp side chain would pack more tightly against the peptide. However, substitution of Ala-99␤ with Trp is predicted to result in the added exposure of ϳ50 Å 2 of hydrophobic surface area when the TCR is displayed on the yeast surface (compared with substitution with Tyr). Using well recognized estimates of 20 -30 kcal/mol in free energy per Å 2 of exposed hydrophobic surface area (40,41), the resultant destabilization could be as high as ϳ1.5 kcal/mol.
A wider comparison of enrichment versus predicted changes in hydrophobic solvent-accessible surface further suggested that sequence enrichment in the fitness landscapes was at least in part associated with reduced exposed hydrophobic surface area, and vice versa (Fig. 6B). Of the residues with log 2 enrichment ratios Ն1, 68% are predicted to have reduced solvent-

TABLE 1 Affinity and MFI MAX ratios of RD1-MART1 HIGH mutations identified from sequence fitness landscape
The affinities of RD1-MART1 HIGH and its mutants are indicated, based on the average values determined from two or more titrations (n) for each TCR, with monomeric MART1⅐HLA-A2. MFI MAX ratio is the ratio of MFI (in mean fluorescence units) from staining mutants with either a high concentration of MART1⅐HLA-A2 tetramer or the highest concentration of monomer, both relative to MFI from RD1-MART1 HIGH . ND, not determined. Engineering Protein-Protein Interactions exposed hydrophobic surface. In contrast, for those residues with enrichment ratios Յ1, 66% are predicted to have increases in exposed hydrophobic surface. We attempted to compute changes in protein stability and incorporate them into binding predictions, as performed previously (14). We used multiple approaches to compute stability changes, including quantifying energetic costs of changes in hydrophobic surface exposure and computing stability using both our TCR-specific and more general scoring procedures. These approaches, however, did not provide additional insight beyond that obtained from qualitatively assessing hydrophobicity, as noted above and in Fig. 6C. This result may be related to recognized difficulties in computing protein stability (42), especially in loop structures, as well as an intrinsically low sta- . Binding constants were calculated using non-linear regression of the binding curve (one site-specific binding in GraphPad) based on the highest concentrations of pepMHC as an estimate for maximum MFI. For the N28G␤ and K97Q␣ mutants, a maximum MFI could not be determined strictly from pepMHC monomer staining; therefore, maximum MFI was estimated based on the ratio of tetramer binding to A99Y␤, for which a maximum MFI was reached with monomer. Inset, a plot of the ⌬⌬G values calculated for each mutant relative to the wild-type RD1-MART1 HIGH interaction, versus enrichment value. The dashed line shows a linear regression to the data. The correlation coefficient (r) is Ϫ0.82. NOVEMBER 18, 2016 • VOLUME 291 • NUMBER 47

JOURNAL OF BIOLOGICAL CHEMISTRY 24571
bility of the single-chain RD1-MART1 HIGH protein coupled with kinetic, rather than thermodynamic, stability effects that may lead to low yeast surface expression levels (43).
To test whether reduced yeast surface levels of the properly folded mutants were in part responsible for lower enhancement values, we examined four individual mutants. These mutants were predicted to have higher affinities, yet they were observed to have negative enrichment values (Fig. 6A, green). The four mutants (A99W␤, S94M␣, G28W␣, and G28F␣) were compared with the wild-type RD-MART1 HIGH TCR for binding to MART1⅐HLA-A2 monomers (see Fig. 7A for mutants and Fig.  3A for wild type) and tetramers (Table 1) using flow cytometry. Consistent with the notion that negative enrichment values could be due to lower surface levels of the properly folded proteins, all four mutants, and especially the S94M␣, G28W␣, and G28F␣ mutants, were expressed at significantly lower surface levels than wild-type RD-MART1 HIGH TCR. Although the A99W␤ mutant appeared to have a similar affinity to the wildtype RD-MART1 HIGH TCR, based on the monomer titrations (Fig. 7B), it was not possible to estimate affinities for the other three mutants due to their low surface levels.
Computational Analysis of Mutations That Yielded Higher Affinities-If we restrict the focus solely to mutations with binding measurements, after computing binding ⌬⌬G 0 relative to RD-MART1 HIGH using the values in Table 1, moderate agreement between computation and experiment is observed, with a correlation coefficient of 0.68 (Fig. 8A). This analysis, however, includes the N28G␤ and K97Q␣ mutations. Neither of these positions interacts directly with pepMHC either in the wild-type complex or in the models of the mutants through

Engineering Protein-Protein Interactions
either short or long range interactions (the closest predicted approaches for the N28G␤ and K97Q␣ mutants are 16 and 11 Å, respectively). Excluding these mutants from the comparison substantially improved the correlation between prediction and experiment (r ϭ 0.91). Prediction scores of the double and triple mutants also corroborated well with experiments, with a correlation coefficient of 0.87 when these were included in the comparison (Fig. 8B).
The accurate recapitulation of the relative binding energies of the RD1-MART1 HIGH mutants in Table 1 suggested that the structural models would be of use in rationalizing how the mutations acted to impact binding. The V50D␤ mutation, which improved binding nearly 4-fold, may function by introducing an electrostatic interaction with the side chain of Arg-65 of HLA-A2 (Fig. 2B). The T91Q␣ mutation occurs in the V␣/V␤ interface and may stabilize the CDR3␣ loop through an intraloop hydrogen bond and electrostatic interactions within the neighboring CDR3␤ loop (Fig. 2C). A99Y␤ is the most favored mutation and is predicted to improve packing with the peptide and the HLA-A2 ␣1 helix (Fig. 2D). The new interactions for A99Y␤ are largely non-directional contacts, perhaps accounting for why multiple substitutions here are both highly enriched and favorably predicted. G28E␣ is predicted to form new van der Waals and electrostatic interactions with Trp-167 of HLA-A2 (Fig. 2E).
It is notable that, with the exception of A99Y␤, the distance of these positions from the interface would tend to place them as a lower priority for investigation in typical structure-guided design studies, with the V50D␤ mutation serving as the best example. This highlights the strength of deep sequencing and yeast surface display in identifying important affinity-enhancing "second shell" mutations, which for TCRs could be important in preserving or even improving antigen specificity (25).

Discussion
Proteins often require affinity and stability enhancements for therapeutic purposes. To achieve the desired characteristics, various protein-engineering techniques that use in vitro directed evolution have been employed. These include ribosomal, phage, bacterial, yeast, and mammalian display of degenerate mutational libraries, followed by various selection

Amino acid residues at selected CDR positions of RD1-MART1 HIGH TCR mutants
The wild-type RD1-MART1 HIGH sequence at the 28␣, 91␣, 50␤, and 99␤ CDR residues is indicated. Mutants created based on the sequence fitness landscape and mutations of clones isolated from the combinatorial library are shown. The horizontal dashed line indicates a binding score of 0.92, above which binding is predicted to be weaker compared with wild type and below which binding is predicted to be stronger. The majority of the negatively enriched mutants are predicted to bind weaker, but there are many mutants predicted to bind stronger that are also negatively enriched. The four most highly enriched mutants are indicated in red and noted on the plot. Four mutants that were examined further in Fig. 7 are highlighted in green. B, of the amino acids with enrichment ratios Ն1, 68% are predicted to have reduced solvent-exposed hydrophobic surface. In contrast, for those residues with enrichment ratios Յ1, 66% are predicted to have increases in exposed hydrophobic surface. C, prediction score versus enrichment score with each data point colored by the predicted change in hydrophobic solvent-accessible surface area. NOVEMBER 18, 2016 • VOLUME 291 • NUMBER 47

JOURNAL OF BIOLOGICAL CHEMISTRY 24573
schemes. More recently, it has been possible to conduct deep mutational scans of every residue in a protein's binding site to construct a sequence fitness landscape. Here we used information from the sequence fitness landscape of a TCR (RD1-MART1 HIGH ) specific for the cancer antigen MART1⅐HLA-A2 as a guide to compare various approaches to engineer improvements in the affinity and/or stability. We were able to identify specific mutants that enhanced the affinity Ͼ200-fold and increased the yeast surface level expression by 6-fold. Furthermore, we compared these site-specific mutants with structurebased predictions using an in silico approach that modeled the effects of the mutations and predicted the impacts on binding.
To determine the robustness of identifying higher affinity mutations strictly from the results of the sequence fitness landscape, we generated single site-specific mutations with enrichment values ranging from 10-to 51-fold. These mutations were tested as individual and combined mutants. Site-specific mutagenesis was able to increase the TCR affinity by 100-fold when three of the highest affinity mutants were combined. A previous study that aimed to improve the affinity of an influenza inhibitor reported a maximum affinity improvement of 28-fold using a combination of eight mutations identified from a sequence fitness landscape (14). In our study, when enrichment values were Ͼ12, we observed a direct correlation with enrichment value and affinity enhancement. In addition to affinity increases, the yeast display approach also selects for TCRs that are expressed at higher levels due to the increased stability of the mutant (18, 30 -32). Accordingly, the present study identified mutations that increased surface levels up to about 6-fold. Similar observations have been made about the role of residues in the thermal stability of antibody scFv fragments (32,44,45).
Although the individual mutations with the highest enhancement values in the sequence fitness landscape yielded higher affinity, and their combinations provided even higher affinity (100-fold), we were interested in determining whether a directed evolution approach based on the same data could yield combined mutations with even higher affinity or stability. This combinatorial approach would include variants with multiple mutations that have additivity or positive cooperativity in binding. To test this, we generated a combinatorial library at the four most highly enriched residues (Gly-28␣, Thr-91␣, Val-50␤, and Ala-99␤). Three successive selections of this library yielded multiple clones that bound with higher affinity than RD1-MART1 HIGH to the MART1⅐HLA-A2 ligand. One of these clones (S3-2) had the highest affinity (700 pM) and greatest surface level (6-fold above RD1-MART1 HIGH ) among all of the mutants examined in the present study. The sequences of this and the other two unique clones revealed that three of the four positions in the library had mutations that were highly enriched in the sequence fitness landscape. Unexpectedly, the Gly-28␣  position had two clones, including S3-2, with an arginine substitution, whereas the G28R␣ substitution was not enriched in the sequence fitness landscape. Because G28R␣ is predicted to interact strongly with Glu-1 of the peptide as well as Glu-58 of HLA-A2, this mutation alone may destabilize the TCR, a consequence offset by the other mutations. Accordingly, the use of sequence fitness landscapes to generate higher affinity mutants can benefit by subsequent selections of combinatorial libraries in the codons of enriched amino acid positions.
Our recently described structure-guided computational design approach for engineering TCRs (24,39) provided an opportunity to evaluate the sequence fitness landscapes computationally and compare enrichment with predicted impacts on binding. Although there was qualitative agreement between those mutations predicted to weaken binding and sequence depletion, the in silico approach did not distinguish well between mutations predicted to improve binding and sequence enrichment, because 25% of mutations selected against were predicted to enhance binding. Much of this may reflect the sensitivity of yeast surface display to changes in protein stability, resulting in some substitutions that are either reduced or enhanced in surface levels. This is supported by our analysis of the modeled structures of each mutant for changes in exposed hydrophobic surface area; sequences enriched showed a much stronger trend to reduce hydrophobic solvent-accessible surface area, and vice versa. There can be a tendency for structureguided design to select for hydrophobic mutations (12,24), and our findings reiterate that improvements in in silico approaches may be found from continued attention to electrostatic features, such as polar solvation, hydrogen bonding, and interactions with formal charges (16, 46 -49).
For those mutants studied directly, our in silico analysis showed very good agreement between predicted and impacted effects on binding, permitting the use of the structural models in helping to assess how the mutations acted to improve binding. A key observation was that multiple "second shell" mutations appeared to be enriched, and these most likely act via indirect or long range effects (exemplified by the V50D␤ mutation). Such sites are not always considered in structure-guided computational design, highlighting a strength of sequence fitness landscapes and yeast surface display as well as further suggesting how in silico methods might be extended.
The single-codon approach used with sequence fitness landscapes also improves upon the many early studies using directed evolution and error-prone PCR techniques. Although the latter can sample the entire protein interface, it is limited to single-site mutations that are generated by a single base substitution. This is because the probability of having two substitutions in the same codon or in two codons with beneficial substitutions is quite low and often beyond the size of the libraries (or detrimental mutations at a higher error rate would obscure these mutations). As an example, the A99Y␣ mutant was one of the most highly enriched mutations of RD1-MART1 HIGH , but it would require two nucleotide substitutions to mutate from alanine to tyrosine.
Finally, it is worth considering whether information about binding affinities of mutants gained in the yeast display system could be extrapolated to these TCRs expressed in their normal context, T cells. The TCRs are in single-chain form (V␤-linker-V␣) as Aga2 fusions for yeast display, whereas in T cells, the V␣ and V␤ each contain transmembrane-spanning constant regions. Despite these differences, we have invariably observed with other TCRs that the mutations yielding higher affinity in the yeast display system also yielded higher affinity when used as full-length TCRs transferred into T cells (18,37,50,51).
Thus, although not formally tested, we have reason to believe that the mutations of RD1-MART HIGH TCRs identified in the present work would exhibit binding with enhanced affinities when transferred to the full-length TCR context.
In conclusion, deep mutational scans of the entire proteinprotein interface provide physical insights into binding, and they will allow improvements in structure-guided in silico analyses. Sequence fitness landscapes also serve as a robust guide for conducting site-specific mutagenesis to enhance affinity. In the future, one might be able to simply combine all of the top enriched mutations into a single mutant and thus bypass the need for any selections. This might be especially true where the residues are located at different locations in the interface. However, the combinatorial library/selection approach, based on the results of the sequence fitness landscape, yielded an improved mutant compared with the single-site mutant strategy. Thus, there may still be advantages to using a directed evolution, selection-based strategy for optimizing improvements.

Experimental Procedures
Reagents and Flow Cytometry-The procedure for induction of yeast to express single-chain TCRs has been described previously (21). HLA-A2 heavy chain and ␤ 2 -microglobulin were expressed as inclusion bodies in E. coli and refolded with a UVcleavable peptide (KILGFVFJV) (52). At the C terminus of HLA-A2 heavy chain, a biotinylation site was included for the addition of biotin with BirA enzyme (Avidity). To exchange peptide, refolded HLA-A2 with UV-cleavable peptide was incubated with a 100-fold excess of MART1 10-mer peptide (ELAGIGILTV) during 45 min of UV exposure.
Yeast stained with HLA-A2 monomer were incubated with various concentrations of MART1⅐HLA-A2 monomers (at the lowest concentration, a minimum of 2-fold excess monomer relative to the total number of TCRs, estimated to be 50,000 receptors/cell, was used). Cells were washed and stained with streptavidin-PE (1:300) and subsequently analyzed by flow cytometry. To determine binding affinities, the mean fluorescence intensity (MFI) of secondary reagent control was subtracted from the MFI of the yeast cells stained with monomer. Binding constants were calculated using non-linear regression of the binding curve (one site-specific binding in GraphPad). Maximum MFI values were the MFI achieved at the highest concentration of monomer (for titrations that had a plateau), or, in cases where saturation was not reached (e.g. N28G␤ and K97Q␣ mutants), maximum MFI values were estimated from the ratio of tetramer binding to a TCR run at the same time that did exhibit a plateau. Titrations were performed a minimum of two times, in completely separate experiments.
Single-site Mutants and Combinatorial Library Design-Single-site mutants of the RD1-MART1 HIGH TCR in the pCT302 yeast display vector were created using the QuikChange II sitedirected mutagenesis kit (Agilent Technologies). The identity of mutants was confirmed by sequencing. EBY100 yeast were subsequently transformed using a lithium acetate heat shock protocol (53). A four-codon combinatorial library of RD1-MART1 HIGH spanning four CDR domains was created by splice overlap extension (SOE) PCR. Pre-SOE products were created using the following primers: Pre-SOE products were combined by SOE PCR using the T7 and Splice 4L primers. Libraries were created using homologous recombination in EBY100 yeast as described previously (54). The theoretical diversity of the four-codon library is 32 4 , or 1 ϫ 10 6 . Based on the number of independent transformants, the estimated size of the combinatorial library was 2 ϫ 10 7 , which should provide coverage of all amino acid combinations of the library. The library was sorted three times with MART1⅐HLA-A2 monomer at 50, 0.5, and 0.1 nM, respectively, collecting the top 1% of fluorescent cells at each sort. After the third sort, eight high affinity clones were sequenced, of which three unique sequences were identified.
Characterization of Binding by Surface Plasmon Resonance-Recombinant single chain S3-2 TCR and MART1⅐HLA-A2 were generated from bacterially expressed inclusion bodies, refolded, and purified as described (25). Binding analysis was performed at 25°C using surface plasmon resonance on a Biacore T200 instrument. 150 resonance units of S3-2 TCR was coupled to a CM5 sensor surface using standard amine coupling. MART1⅐HLA-A2 analyte was injected in a series of concentrations at a flow rate of 45 l/min. Data were corrected for injections over a blank surface, and the responses versus time were fit globally to a 1:1 kinetic binding model using Biacore T200 evaluation software version 3.0.
Computational Prediction of Mutations-Structural modeling and energetic scoring followed closely the procedures described by Riley et al. (24). For modeling, Rosetta with the Talaris2013 score function was used (55)(56)(57), using the PyRosetta interface (58). The crystallographically absent linker region of the RD1-MART1 HIGH structure was constructed and modeled using kinematic loop closure (59). Following linker modeling, the entire structure was relaxed through multiple cycles of backbone minimization and rotamer optimization (60). For comparison with the experimental sequence landscapes, 19 amino acids were introduced separately into each position, and the mutants plus the wild-type structures (760 in total) were modeled from three independent Monte Carlobased simulated annealing trajectories. Modeling utilized Rosetta's Loop_Mover_Refine_CCD mover with three outer cycles and 10 inner cycles, with an initial metropolis acceptance criterion of 2.2 that decreased linearly to 0.6 (61). For double and triple mutants, mutations were introduced simultaneously, followed by a minimum of six independent minimization trajectories to account for additional structural impacts. For all models, the large number of packing operations introduced some minor variability when scoring. Therefore, the unweighted score terms for the three trajectories were averaged and stored for calculations (62).
Each model was scored in the Rosetta interface with the TR3 score function trained on TCR-specific data. 4 TR3 utilizes dynamic information from the wild-type structure as positional modifiers to account for protein flexibility. Accordingly, molecular dynamics simulations were calculated using the AMBER molecular dynamics suite (63). Starting coordinates for the complex used the relaxed structure described above, and the coordinates for the free TCR were generated by stripping away the pepMHC. After adding solvent and counterions, systems were energy-minimized and heated to 300 K with solute restraints. Solute restraints were gradually relaxed and followed with 2 ns of simulation with no solute restraints for equilibration, after which 100-ns production trajectories were calculated. Trajectory analysis for determining root mean square fluctuation (RMSF) values used ccptraj from the AMBER suite. The RMSF values were used along with the additional terms of the TR3 function to score each of the 760 modeled structures. The scores of the 114 wild-type control structures determined the baseline score of 0.92 energy units, where lower scores indicated stronger affinity and higher scores indicated weaker affinity. For double and triple mutants, RMSF descriptors were averaged between the positions for scoring purposes (i.e. for double mutant XY, the RMSF of position X is averaged with the position Y RMSF to give an RMSF descriptor for double XY). Solventaccessible surface areas were calculated using a 1.4-Å probe.
Calculations of protein stability included using structural models to calculate changes in exposed hydrophobic surface area and computing changes in stability using both the Talaris2013 and TR3 score functions. Computing binding probabilities as in Whitehead et al. (14) or by accounting for percentage of bound ligand using both folding and binding partition functions together did not yield more quantitative insight than simply considering changes in exposed hydrophobic surface area.