Comprehensive Deep Mutational Scanning Reveals the Immune-Escaping Hotspots of SARS-CoV-2 Receptor-Binding Domain Targeting Neutralizing Antibodies

The rapid spread of SARS-CoV-2 has caused the COVID-19 pandemic, resulting in the collapse of medical care systems and economic depression worldwide. To combat COVID-19, neutralizing antibodies have been investigated and developed. However, the evolutions (mutations) of the receptor-binding domain (RBD) of SARS-CoV-2 enable escape from neutralization by these antibodies, further impairing recognition by the human immune system. Thus, it is critical to investigate and predict the putative mutations of RBD that escape neutralizing immune responses. Here, we employed computational analyses to comprehensively investigate the mutational effects of RBD on binding to neutralizing antibodies and angiotensin-converting enzyme 2 (ACE2) and demonstrated that the RBD residues K417, L452, L455, F456, E484, G485, F486, F490, Q493, and S494 were consistent with clinically emerging variants or experimental observations of attenuated neutralizations. We also revealed common hotspots, Y449, L455, and Y489, that exerted comparable destabilizing effects on binding to both ACE2 and neutralizing antibodies. Our results provide valuable information on the putative effects of RBD variants on interactions with neutralizing antibodies. These findings provide insights into possible evolutionary hotspots that can escape recognition by these antibodies. In addition, our study results will benefit the development and design of vaccines and antibodies to combat the newly emerging variants of SARS-CoV-2.

INTRODUCTION SARS-CoV-2, which causes viral pneumonia in humans, is the cause of COVID-19 (Lai et al., 2020). Under an electron microscope, the virus exhibits crown-like morphology ("corona") and is thus named coronavirus (Gui et al., 2017). The World Health Organization declared COVID-19 as a pandemic. In April 2021, there were 142.5 million confirmed cases of COVID-19, including 3,043,707 deaths (daily online worldwide data about COVID-19 1 ). Common symptoms of SARS-CoV-2 infection include diarrhea, dry cough, fever, nasal congestion, respiratory problems, and sore throat (Baj et al., 2020). In severe cases, kidney failure, severe acute respiratory syndrome, and pneumonia may ensue, eventually leading to death (Lai et al., 2020).
SARS-CoV-2, a single-stranded positive-sense enveloped RNA virus, consists of an RNA sequence of approximately 30,000 bases (Naqvi et al., 2020). This viral genome has 10 open reading frames (ORF) (Tsai et al., 2020). Of these, ORF1ab encodes polyprotein lab (pp1ab), which is cleaved by the proteases 3CL pro and PL pro to yield multiple proteins associated with viral RNA replication and transcription (Graham et al., 2008;Moustaqil et al., 2021) as well as 16 non-structural proteins, creating the replicationtranscription complex of SARS-CoV-2 (Romano et al., 2020). In addition, ORFs 2-10 encode four structural proteins: spike (S), membrane (M), nucleocapsid (N), and envelope (E). The N protein is critical for packing the RNA genome, and the S, M, and E proteins are essential for viral coating. The S protein is a large oligomeric transmembrane protein responsible for the entry of the virus into the host cell (Lan et al., 2020). It comprises two functional domains: S1 and S2; the S1 domain comes into contact directly with the angiotensin-converting enzyme 2 (ACE2) receptor on the host cell (Wrapp et al., 2020), whereas the S2 domain mediates cell membrane fusion Wrapp et al., 2020).
SARS-CoV-2 enters the host cell through ACE2; thus, the S protein partly determines its transmissibility and infectivity (Hoffmann et al., 2020). The receptor-binding domain (RBD) of the S1 subunit directly interacts with ACE2 (Lan et al., 2020;Yang et al., 2020). Thus, some antiviral drugs targeting RBD were developed. Small molecules, such as chloroquine, hydroxychloroquine, ivermectin, and azithromycin, have been reported to target the S protein-ACE2 interface (Pandey et al., 2020;Batalha et al., 2021;Mirtaleb et al., 2021). Moreover, novel drug-like compounds DRI-C23041 (Rajgor et al., 2020) and DRI-C91005 (Lan et al., 2020) have been observed to inhibit the S protein-ACE2 interaction, with low micromolar activity. The S protein is immunogenic; hence, several approaches have targeted it for viral neutralization. Neutralizing antibodies targeting RBD have also been developed (Pinto et al., 2020;Rogers et al., 2020;Xiaojie et al., 2020;Liu et al., 2021;Lu et al., 2021). Some antibody-based antiviral therapeutics have demonstrated high specificity, potency, and modularity. However, RNA viruses continually change through mutations, leading to the emergence of new variants (Pachetti et al., 2020;Nagy et al., 2021;Wang et al., 2021). This antigenic evolution 1 https://www.worldometers.info/coronavirus/ leading to RBD mutations overcomes the established neutralizing antibody immunity (Eguia et al., 2021;Greaney et al., 2021a). It is therefore critical to systematically monitor the antigenic evolution and investigate viral mutations that can impair the immune response conferred by neutralizing antibodies.
Here, we comprehensively estimated the RBD mutations that destabilize the binding of five representative neutralizing antibodies: the H11-D4 and VH1-2-15 nanobodies, MR17 and SR4 sybodies, and P2B-2F6 Fab, which target RBD's receptorbinding motif. We employed complex structures retrieved from the Protein Data Bank to calculate binding stability through detailed mutational scanning, in which a single residue was replaced by all other 20 amino acids in RBD to systematically investigate the hotspots that affect binding. The resulting heatmaps demonstrated that mutations at R403, K417, G447, N448, Y449, N450, L452, Y453, L455, F456, E484, G485, F486, Y489, F490, P491, L492, Q493, S494, Y495, and G496 were unfavorable for binding with antibodies. Notably, the E484K and L452R mutants are also present in the P.1 viral lineage. Moreover, F456 variants have reduced binding to neutralizing antibodies, and L455, F486, and F490 have substantial antigenic effects. The N501Y mutant is present in emerging viral lineages, such as B1.1.7 and B.1.351. All the aforementioned clinical and experimental reports support our findings. Thus, the interactive residues of RBD (Y449, L452, L455, E484, Y489, F490, L492, Q493, and S494) identified in this study can be hotspots for further antibody engineering or vaccine developments to combat potential variants of SARS-CoV-2.

Preparation of Protein Structures
The structures of the SARS-CoV-2 S protein in complex with nanobodies, sybodies, Fabs, and ACE2 were retrieved from the Protein Data Bank (PDB IDs: 6M0J, 6YZ5, 7BWJ, 7C8V, 7C8W, and 7L5B). The CHARMm Polar H forcefield was applied to all complex structures before computations.

Calculation of Mutational Binding Stability
The mutational binding stability of RBD with its targets was estimated by Discovery Studio (DS) 3.5 (Accelrys, San Diego, CA, United States), MutaBind2 , FoldX (Schymkowitz et al., 2005), and mCSM-PPI2 (Rodrigues et al., 2019). For the prediction in DS, the Calculate Mutation Energy (Binding) protocol was used to estimate the mutational binding stability. The complex structure was used as the "input typed molecule, " and the complexed partners of RBD were employed as the "ligand chain." Additionally, interactive residues of RBD, making direct contact with its targets (within a maximum distance of 5 Å from the targets' interface), were selected for a mutational study. Furthermore, a single mutation was used as "mutation sites, " and all 20 amino acids were chosen as the "mutations" parameter." Additionally, the dielectric constant [The dielectric constant of a molecular interior corresponds to the measure of electric potential energy. Therefore, the induced polarization is stored within a given volume of substance under the action of an electric field. It is expressed as the ratio of the dielectric permittivity of the material to that of a vacuum or dry air (Ilic, 2001)], solvent dielectric constant, maximum structures to save, and maximum number of mutants were set to 10, 80, 25, and 25, respectively. All other parameters were set to default. For the prediction in FoldX, the binding stability was determined according to a previous report (Teng et al., 2021). In addition, the calculations in MutaBind2 and mCSM-PPI2 were followed by the online prediction instructions. The mutational binding stability resulted from the calculations of the effect of mutations on the binding affinity. We performed combinatorial scanning mutagenesis in protein complexes, depending on the selected mutation sites. For each single mutant, the differences in the free energy of binding between the wild-type and mutated structures are calculated. Mutational energy is the total free energy difference between the wild-type and mutated structures. It is calculated as a weighted sum of the VDW, electrostatic, entropy, and non-polar terms. Accordingly, the negative and positive values of mutational energies represent the stabilized and destabilized binding stabilities.

Generation of Heatmaps of Mutational Binding Energy
The binding energy values upon single mutations were subjected to heatmap generation by using Excel. The residues for mutations were used as the y-axis, and the amino acid types of single mutations were used as the x-axis. The binding energy values of all mutations were selected for conditional formatting. The blue and red colors are used to indicate the values of binding energies, respectively, with the gradient to a range of cells for the color scales.

Ligplot Analyses
The molecular interactions of each complex structure of the S protein with its targets were analyzed using a ligplot (Wallace et al., 1995;Laskowski and Swindells, 2011), and its DIMPLOT module was employed for plotting protein-protein and domaindomain interactions.

Mutational Binding Stability of RBD Variants Interacting With the H11-D4 Nanobody
Before investigating the effects of RBD mutations on the binding of antibodies, we searched for the dissolved complex structure of RBD with antibodies in the Protein Data Bank. The H11-D4 nanobody has high affinity to the S protein and blocks its attachment to ACE2. The structure of the H11-D4 nanobody in complex with RBD is illustrated in Figure 1A. The epitope of RBD to H11-D4 was analyzed with a ligplot and is illustrated in Figure 2A. The RBD residues N450, E484, F490, Q493, and S494 predominantly form hydrogen bonds, and residues R346, Y449, L452, F456, V483, Y489, and L492 make hydrophobic contacts with the H11-D4 nanobody. These interactive residues were further mutated into all 20 amino acids to estimate their effects on the binding stability against the H11-D4 nanobody. Additionally, RBD residues that make direct contact with the H11-D4 nanobody (within a maximum distance of 5 Å from the nanobody interface) were selected for the mutational study. The resultant binding stability of the RBD variants was plotted into a heatmap, with the y-axis as the residues mutated and the x-axis as the amino acid types of single mutations. The calculated binding stability values are colored with the gradient of a range between blue (stabilized binding) and red (destabilized binding). The results revealed that mutations at Y449, N450, L452, E484, Y489, F490, P491, L492, Q493, and S494 were mostly not favorable for binding stability ( Figure 2B and Table 1).

Changes in the Binding Stability of RBD Variants Targeting P2B-2F6 Fab
The mutational effects on the interactions of SARS-CoV-2 RBD variants with the plasma antibody P2B-2F6 Fab were also examined. The complex structure of P2B-2F6 Fab and SARS-CoV-2 RBD is illustrated in Figure 1B. P2B-2F6 Fab primarily interacts with RBD through its heavy chain. As illustrated in the ligplot (Figure 6A), the epitope residues are in the RBM of RBD, including Y449, N450, and E484 and K444, G446, G447, N448, L452, V483, G485, F490, and S494. We systematically analyzed their mutational effects on the binding stability of RBD toward P2B-2F6 Fab. The results revealed that the amino acid replacements of the interactive residues G447, Y449, N450, L452, V483, E484, and F490 resulted in unstable binding between RBD and P2B-2F6 Fab ( Figure 6B and Table 1).

DISCUSSION
Epidemics have been affecting humans for decades (Khan et al., 2020). Millions of people have died from epidemics, including the 6th and 7th "EI Tor" cholera (Zhang et al., 2014), influenza A (H2N2) (Joseph et al., 2015), and HIV/AIDS FIGURE 2 | The mutational binding stabilities of RBD variants interacting to H11-D4 nanobody. (A) The molecular interactions of SARS-CoV-2 RBD with H11-D4 nanobody analyzed by ligplot. The chains E and F correspond to RBD and H11-D4 nanobody, respectively. (B) The heatmap of interactive residues of RBD derived from the calculated mutational binding stabilities by using Discovery Studio 3.5 (DS), Mutabind2, FoldX, and mCSM-PPI2. The boxes of each mutations were colored with the gradient of a range between blue (stabilized binding) and red (destabilized binding). In all panels, the solid and hollow circles denote significant and moderate decreases of the binding stabilities, respectively.
Frontiers in Microbiology | www.frontiersin.org . c Mutations at other structurally adjacent sites in the RBD's receptor-binding ridge can also have substantial antigenic effects.
The bold values are residue number in SARS-COV-2 RBD.
Frontiers in Microbiology | www.frontiersin.org The heatmap of interactive residues of RBD derived from the calculated mutational binding stabilities by using Discovery Studio 3.5 (DS), Mutabind2, FoldX, and mCSM-PPI2. The boxes of each mutations were colored with the gradient of a range between blue (stabilized binding) and red (destabilized binding). In all panels, the solid and hollow circles indicate significant and moderate decreases of the binding stabilities, respectively. (Gayle and Hill, 2001). At present, the COVID-19 pandemic presents a major global threat. Its causative agent, SARS-CoV-2, is an RNA virus and can thus mutate and evolve (Zheng, 2020).
RNA viruses have higher mutation rates than DNA viruses (Sanjuan and Domingo-Calap, 2016). SARS-CoV-2 has been evolving at a rate of one-two mutations every month in the The heatmap of interactive residues of RBD derived from the calculated mutational binding stabilities by using Discovery Studio 3.5 (DS), Mutabind2, FoldX, and mCSM-PPI2. The boxes of each mutations were colored with the gradient of a range between blue (stabilized binding) and red (destabilized binding). In all panels, the solid and hollow circles denote significant and moderate decreases of the binding stabilities, respectively.
Several experimental studies made effort to investigate the mutational escape from neutralizing and convalescent antibodies (Andreano et al., 2020;Li et al., 2020;Starr et al., 2020;Weisblum et al., 2020a;Garcia-Beltran et al., 2021a;Greaney et al., 2021a,b). Weisblum et al. (2020b) employed a recombinant chimeric VSV/SARS-CoV-2 reporter virus to investigate the mutations in the RBD, which confer resistance to monoclonal antibodies or convalescent plasma. They found E484K, F490L, and Q493K/R occured at high freequency during recombinant chimeric VSV/SARS-CoV-2 passage in the presence of neutralizing antibodies or plasma. In addition, the mutants E484K and Q493R caused apparently complete resistance to monoclonal antibodies (Weisblum et al., 2020b). Li et al. (2020) investigated 80 RBD variants for the infectivity and reactivity to a panel of neutralizing antibodies and sera from convalescent patients. They reported that most variants were less infectious, but L452R and F490L became resistant to some neutralizing antibodies (Li et al., 2020). Also, the immune escape of lineage B.1.351 of South Africa was examined and revealed that variants K417N/T, E484K, and N501Y were highly resistant to neutralization (Garcia-Beltran et al., 2021b). Greaney et al. (2021a) have mapped all the mutations to the RBD that escape binding by antibodies isolated from convalescent plasma. Their yeast-display deep mutational scanning revealed that antibodies were escaped by mutations to sites K417, N450, L452, L455, F456, E484, F486, F490, and Q493 (Greaney et al., 2021a). In our study, we analyzed that RBD residues Y449, L452, L455, E484, Y489, F490, L492, Q493, and S494 were hotspots with significantly destabilizing effects on binding to neutralizing antibodies. Especially, the identified hotspots, L452, L455, E484, F490, and Q493 are well consistent with the reported immune escape mutations from neutralizing and convalescent antibodies. Moreover, our findings that impaired binding stability of RBD and antibodies resulted from the mutational effects at L455, F456, G486, F486, and F490 (Table 1) are in good concordance with those of Greaney et al. (2021a) -mutations at sites near the structurally adjacent site of RBD's receptor-binding ridge (e.g., L455, F456, G485, F486, and F490) have substantial antigenic effects. All these consistencies indicate the precision and reliance of our computational study in revealing the hotspots of SARS-CoV-2 RBD-specific neutralizing antibodies. In addition to the antibody immune escape of RBD, its binding and interaction to ACE2 were also characterized by in vitro and in silico mutational studies. Starr et al. (2020) systematically changed every amino acid in the RBD of the SARS-CoV-2 spike protein and determined the effects of the substitutions on RBD expression, folding, and ACE2 binding. They found that there are handful of sites where ACE2 binding imposes strong constrain (e.g., Y489, G502, and Y505), and mutations at interface residues (Y449, L455, F486, and Y505) enhance RBD expression but destabilize the effect of surface-exposed hydrophobic patches required for ACE2 binding. Mutations that enhance ACE2 binding affinity of RBD are notable at sites Q493, Q948, and N501. Besides, Teng et al. (2021) have conducted a computational study to investigate the effects of mutations on SARS-CoV-2 RBD-ACE2 binding affinity. They reported that mutations on residues G476, V483, Q498, T500, G496, and G502 exert apparently destabilizing impacts in RBD-ACE2 complex. In our study, we also analyzed the mutational effects of RBD on interacting with ACE2 (Supplementary Figure 1), in which variants of Y449, L455, F456, F486, N501, G502, and Y505 conspicuously impaired the binding stability, corroborating with Starr's and Teng's findings.
To explore the possible mechanism of action of the identified hotspots, we further investigated the molecular interactions of mutational variations in the structural complex of RBDantibody. The residue K417 of RBD connected with D101 of MR17 by charge-charge interactions. Thus, RBD-MR17 complex could be destabilized when K417 was mutated to non-charged N417 and T417 (Figures 7A,B). Structurally, the residue L452 of RBD contacted with W99 and V102 of SR-4 and H11-D4, respectively, through hydrophobic interactions. These hydrophobic contacts are probably disrupted as L452 is substituted by an arginine residue, causing the decreased binding affinity of RBD with SR-4 and H11-D4 (Figures 7C,D). In addition, the charge-charge interactions can also be seen in residue E484, which interacts with R52, R59, and R12 of H11-D4, MR17, P2B-2F6, individually. As well, E484 formed hydrogen bond with Y33 of VH1-2-15. All these molecular interactions were abolished while E484 was replaced by a lysine residue. This could also detract from the binding affinity of RBD with neutralizing antibodies (Figures 7E-H). Moreover, the aromatic residue F490 of RBD makes contributions in binding to antibodies by exerting hydrophobic and cation-pi interactions (Figures 7I-M). These molecular interactions were disrupted when F490 was mutated to S490, therefore destabilizing the binding stability of RBD and antibodies. Notably, the mutation S494P could breakdown the hydrogen bond interactions of S494 (RBD) with N101 and V102 of H11-D4, thus weakening the binding affinity ( Figure 7N). The variant N501Y showed extra hydrophobic interactions with Y41 and K353, significantly increasing the binding stability of RBD to ACE2 (Figures 7O,P), explaining the high concurrency of the N501Y variant in several countries. Taken together, our data revealed all the possible 2 cov-lineages.org hotspots that may substantially impair the binding of SARS-CoV-2 RBD to both ACE2 and neutralizing antibodies. Our findings will benefit the development and engineering of new and potent antibodies and vaccines against SARS-CoV-2.

CONCLUSION
In this study, we explored the possible hotspots of SARS-CoV-2 RBD that can enable virus escape from recognition by neutralizing antibodies. Computational analyses demonstrated that specific variants of RBD significantly impair binding to neutralizing antibodies. Particularly, the RBD residues Y449, L452, L455, E484, Y489, F490, L492, Q493, and S494 were found to be hotspots because their variants could markedly destabilize the binding to neutralizing antibodies. Notably, the hotspots K417, L452, L455, E484, F490, and S494 were supported by evidence from the literature. The hotspots Y449, L455, and Y489 were commonly observed to disrupt the binding to ACE2 and neutralizing antibodies. Conclusively, our data provide insights into the putative impacts of the possible immuneescaping hotspots on interactions with neutralizing antibodies, which can help develop new therapeutic agents against potential variants of SARS-CoV-2.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
TS-T and KC-T designed the experiments and wrote the manuscript. TS-T, KC-T, and YC-L performed the experiments and analyzed the data. All authors contributed to the article and have approved the submitted version.

ACKNOWLEDGMENTS
We thank the National Center for High-Performance Computing (NCHC) for providing computational and storage resources.