SARS-CoV-2 beta variant substitutions alter spike glycoprotein receptor binding domain structure and stability

The emergence of severe acute respiratory syndrome-related coronavirus 2 (SARS-CoV-2) and the subsequent COVID-19 pandemic have visited a terrible cost on the world in the forms of disease, death, and economic turmoil. The rapid development and deployment of extremely effective vaccines against SARS-CoV-2 have seemingly brought within reach the end of the pandemic. However, the virus has acquired mutations. and emerging variants of concern are more infectious and reduce the efficacy of existing vaccines. Although promising efforts to combat these variants are underway, the evolutionary pressures leading to these variants are poorly understood. To that end, here we have studied the effects on the structure and function of the SARS-CoV-2 spike glycoprotein receptor-binding domain of three amino-acid substitutions found in several variants of concern, including alpha (B.1.1.7), beta (B.1.351), and gamma (P.1). We found that these substitutions alter the receptor-binding domain structure, stability, and ability to bind to angiotensin converting enzyme 2, in such a way as to possibly have opposing and compensatory effects. These findings provide new insights into how these variants of concern may have been selected for infectivity while maintaining the structure and stability of the receptor binding domain.

The emergence of severe acute respiratory syndrome-related coronavirus 2 (SARS-CoV-2) in late 2019 and its subsequent spread around the world have caused the deadliest airborne pandemic in the United States, recently surpassing the 1918 influenza pandemic nearly a century ago (1). The international scientific community has risen to the challenge of combating SARS-CoV-2 and COVID- 19. The year 2020 ended with the fastest development of vaccine candidates, starting with the genetic sequence of the virus being reported (2) to human trials of novel mRNA-based vaccines within 3 months. Now, there are three SARS-CoV-2 vaccines approved for use within the United States and many more next-generation and pan-coronavirus vaccines currently in development. These advances have made substantial contributions to the control of the COVID-19 pandemic within the United States. Despite multiple manufacturers receiving emergency use authorization and an unprecedented vaccination campaign, significant challenges remain including uncertainty regarding durability, vaccination hesitancy, limited access to healthcare among disadvantaged persons, as well as the continued emergence of variants of concern (VOC). Our ultimate success in quelling this pandemic may lie in our ability, not only to characterize new variants, but also to be able to predict the emergence of new variants. Such advances will require an increased understanding of evolutionary pressures and constraints on viral variation.
Three SARS-CoV-2 lineages, the alpha variant lineage B.1.1.7 (or 501Y.V1) first identified within the United Kingdom, the beta variant lineage B.1.351 (or 501Y.V2) identified in South Africa, and the gamma variant lineage P.1 (or 501Y.V3) identified in Brazil, have been demonstrated to possess increased infectivity (3) and in the case, beta and gamma exhibit reduced neutralization by antibodies reacting with the cognate regions of the spike protein within the original Wuhan strain of SARS-CoV-2 (4-6). The alpha variant possesses the N501Y substitution within the spike glycoprotein receptor-binding domain (RBD) which has been shown to enhance binding to angiotensin converting enzyme 2 (ACE2), the entry receptor for SARS-CoV-2 (7)(8)(9). The beta and gamma variants possess N501Y as well as substitutions at two other sites within the RBD, E484K, and K417N in beta and K417T in gamma (10). These RBD substitutions present in the spike protein of the B.1.351 and P.1 variants have been shown to reduce the binding and neutralization of mRNA vaccine-induced antibodies as well as potent human monoclonal antibodies (11).
The consequences of the K417N, E484K, and N501Y substitutions on RBD-ACE2 interactions have also been examined, with the increased infectivity of the alpha variant resulting from the enhanced binding to ACE2 when the RBD N501Y substitution is present (9). The E484K substitution has been shown to enhance ACE2 binding (12) and reduce the efficacy of neutralizing antibodies (13). A recent study examined the effects of the K417N substitution on ACE2 binding and antibody interactions using molecular dynamics and found that K417N disrupts RBD-ACE2 interactions, as well as interactions with a monoclonal antibody (14). However, the effects of these substitutions on the structure of the RBD itself have not been examined. Based on the nature of these substitutions, including residue changes in charge or polar to nonpolar substitutions, we hypothesized that the K417N, E484K, and N501Y substitutions alter the RBD structure and stability as well as ACE2 binding interactions. We studied those changes in single-substitution RBD variants as well as in the RBD containing all three substitutions using molecular dynamics and biophysical approaches. Our data suggest that these VOC substitutions significantly alter RBD structure and stability, with consequences for ACE2 binding and proteolytic susceptibility, having potentially opposing consequences for the fitness of new variants. These findings have implications for viral evolution and the design of subunit vaccine candidates.

RBD beta variant substitutions alter RBD structure in silico
We began our studies on the RBD-VOC substitutions with molecular dynamics to investigate whether single substitutions or all three substitutions found within the beta variant would alter the RBD structure. We used homology modeling with residues 319 to 541 of the trimeric spike glycoprotein as a template (PDB 6VXX) to generate structures of the RBD (15,16). The resulting model of the WT/Wuhan strain RBD was used as a template for further modeling of the RBD variants containing single substitutions and all three substitutions (K417N/E484K/N501Y) found within the beta variant. Molecular dynamics simulations were performed with the GRO-MACS 2020.5 package (17) using the CHARMM36 or ff14SB force fields (18,19). After solvation and neutralization, the systems were minimized and equilibrated before undergoing a 25-nanosecond production run within the NPT ensemble. From these trajectories, we observed no large differences in RMSD relative to the starting structure for any of the RBD variants ( Fig. 1, A and B). In the CHARMM36 force field, the RMSD for the N501Y and K417N/E484K/N501Y RBD variants began to diverge relative to the other RBD variants after about 15 ns. The RMSD for the K417N/E484K/N501Y variant quickly converged with the others by the end of the simulation run, whereas in the ff14SB force field, the RMSD was consistent across all the RBD variants. Changes in residue-specific fluctuations were observed for K417N, N501Y, and K417N/ E484K/N501Y RBD variants in the CHARMM36 force field (Fig. 1C) but were absent in the ff14DB force field except for K417N/E484K/N501Y exhibiting a slight increase in the overall residue fluctuations (Fig. 1D). We also observed a decrease in the residue RMSF values around the region of residues 468 to 488 in the ff14SB force field for the E484K, N501Y, and K417N/E484K/N501Y variants that were not seen for the WT and K417N RBD. Overall, in the CHARMM36 force field, we observed VOC substitutions causing more variation in the RMSD and RMSF than in the ff14SB force field.
Next, we examined hydrogen bonding and radius of gyration for all the RBD variants in both force fields. Hydrogen bond content was lower in the CHARMM36 force field ( Fig. 2A) than in the ff14SB force field (Fig. 2B) whereas in both the force fields, the K417N substitution increased average hydrogen bond content relative to all the other RBD variants except for E484K in the ff14SB force field. In the CHARMM36 force field, all the RBD variants exhibited an increase in the radius of gyration, a measure of the compactness of the molecule, relative to the WT RBD throughout the production run with the greatest increases observed for the N501Y and K417N/E484K/N501Y RBD variants (Fig. 2C). In the ff14SB force field, the radius of gyration was overall consistent for all the RBD variants throughout the production run with some slight increases observed for the K417N and K417N/E484K/ N501Y RBD variants within the initial 10 ns of the simulation (Fig. 2D). Taken together, these observations suggest that the CHARMM36 force field predicts that VOC substitutions will have a more disruptive effect on the RBD structure than the ff14SB force field.
To identify common structural changes and differences in the residue contacts, we extracted frames from the simulation trajectories and examined the RBD structure near the sites of the beta variant substitutions at the beginning (100 ps, green) and end (25,000 ps, cyan) of the production run in both force the fields. In the CHARMM36 force field for the WT RBD, we observe K417 and E484 interacting with the solvent whereas N501 forms a hydrogen bond with the polypeptide backbone as well as a hydrogen bond between R403 and E406 ( Fig. 3A). At the end of the simulation, we observed changes in the secondary structure and repositioning of the 468 to 488 loop. When K417 is substituted for asparagine, we observed N417 interacting with E406 at the beginning of the simulation but by the end, N417 has formed solvent interactions whereas a significant rearrangement of the 468 to 488 loop is observed with E484 forming a hydrogen bond with R403 (Fig. 3B). This change in conformation occurs within the first 5 ns of the simulation as shown by the rapid decrease in distance between residues 403 and 484 (Fig. 3G) and residues 417 and 484 (Fig. 3H). For the E484K variant, both the early and late structures are very similar with some changes in the 468 to 488 loop secondary structure (Fig. 3C). For the N501Y variant, we observed that the substitution at position 501 eliminated the backbone interactions observed for N501 in the WT structure and changes in the position of the 468 to 488 loop (Fig. 3D). We also observed changes in the positioning of K417, E406, and R403 whereas the hydrogen bond between R403 and E406 was maintained. For the K417N/E484K/N501Y RBD variant, we observed the same R403, E406 interactions we observed in the other RBD variants (Fig. 3E). However, N417 maintained an interaction with E406 that was not observed in the K417N trajectory. By the end of the simulation for K417N/E484K/ N501Y, we observed the 468 to 488 loop backbone carboxyl groups interacting with arginine 457 (Fig. 3F), a shift in conformation that occurs after about 14 ns in the simulation trajectory and can be seen by changes in the distances between residues 403 and 484 (Fig. 3G) as well as residues 417 and 484 (Fig. 3H). Overall, the observed changes in residue contacts for the K417N and K417N/E484K/N501Y RBD variants reduced the RMSD, relative to the starting structure of the 468 to 488 loop compared with the WT, E484K and N501Y variants (Fig. 3I). In the CHARMM36 force field, we observed stabilization of the 468 to 488 loop in both the K417N and K417N/ E484K/N501Y variants by hydrogen bonding between R403 for K417N and R457 for K417N/E484K/N501Y. We next compared structures extracted from trajectories in the ff14SB force field. For the WT RBD, we observed K417 forming a hydrogen bond with E406 at both time points (Fig. S1A). Here, we also observed a transition in the 468 to 488 loop from a more solvent-exposed position to one forming primarily backbone contacts rather than a hydrogen bond between E484 and R403, which remain relatively far apart throughout the simulation run (Fig. S1F). However, this conformation does not appear to be occupied for a long period of time as the RMSD for this region is elevated compared with K417N and K417N/E484K/N501Y (Fig. S1H). We do observe N417 participating in electrostatic interactions with R403 and E406 (Fig. S1A). For the K417N variant in ff14SB, we do not observe a similar conformational change to that observed in CHARMM36 but rather, we observed the formation of a hydrogen bond between E484 and R457, and this region exhibits a reduced RMSD relative to the other variants except K417N/E484K/N501Y (Fig. S1, B and H). Changes in the structure of the E484K variant in ff14SB are similar to those observed for the WT RBD with the exception of the 468 to 488 loop which does not form the same backbone contacts seen for the WT RBD (Fig. S1C). For the N501Y variant in ff14SB, we do not observe the same K417 -E406 interaction seen in the WT RBD (Fig. S1D). For the K417N/E484K/N501Y variant in ff14SB, we observed changes in the secondary structure within the 468 to 488 loop that may explain the reduced RMSD in this region relative to the other variants except K417N (Fig. S1E). We also observed electrostatic interactions between N417, R403, and E406. As in CHARMM36, in the ff14SB force field, we observed stabilization of the 468 to 488 loop for the K417N and K417N/E484K/N501Y RBD variants by different mechanisms. For K417N, a hydrogen bond between R457 and E484 whereas for K417N/E484K/N501Y, it appears that helical formation in the 468 to 488 loop is responsible for its relative stability compared with the WT RBD. Beta variant substitutions alter RBD stability Beta variant substitutions alter RBD stability RBD beta variant substitutions alter RBD stability and ACE2 binding affinity The molecular dynamics simulations of the RBD variants suggested that VOC substitutions alter the RBD structure and hydrogen bonding. Therefore, we hypothesized that these RBD substitutions alter the RBD stability and resistance to unfolding. To test this, we measured guanidine-induced unfolding of RBD variants using a fluorescence-based unfolding assay (20,21). 4 0 ,6-diamidino-2-phenylindole (bis-ANS) fluorescence was measured for the RBD variants in the presence of increasing guanidine-HCl concentration (Fig. 4A). Fluorescence data were normalized and the unfolding curves were fit by nonlinear regression to estimate the free energy of unfolding ( Fig. 4B) (22). We observed that the E484K substitution significantly destabilized the RBD whereas the N501Y and K417N substitutions significantly stabilized the RBD relative to WT (Fig. 4C). We observed no difference in unfolding energy when all three beta variant substitutions were present in the RBD compared with the WT protein. Changes in m-values for the RBD variants corresponded with changes in folding free energy. We observed a decrease in the m-value for the E484K variant suggesting a reduction in the difference between solvent accessible surface area between the folded and unfolded states compared with the WT RBD (Fig. 4D) (23). For both the K417N and N501Y variants, we observed increased m-values suggesting that the unfolded and folded states exhibit greater differences in exposed surface area upon unfolding relative to the WT RBD. A summary of the parameters fit by nonlinear regression are shown in Table 1. Changes in resistance to denaturant-induced unfolding did not correspond . RBD variant of concern substitutions alter RBD stability. A, the raw bis-ANS fluorescence data as a function of guanidine-HCl concentration for the indicated RBD variant tested. B, unfolding curves for the indicated RBD variant generated by monitoring bis-ANS fluorescence as a function of guanidine-HCl concentration. Nonlinear regression was used to analyze unfolding data. The data points indicated the mean measurement of three replicates whereas error bars indicate SD. C, the calculated ΔG 0 values from nonlinear regression analysis for each RBD variant were compared by one-way ANOVA and Tukey's test for multiple comparisons. p value less than 0.05 is indicated by *, p value less than 0.01 is indicated by **. The error bars indicate SD. D, nonlinear regression m-values for each RBD variant were compared by one-way ANOVA and Tukey's test for multiple comparisons. p value less than 0.05 is indicated by *, p value less than 0.01 is indicated by **, and p value less than 0.001 is indicated by ***. The error bars indicate SD. Guanidine unfolding data are representative of three technical replicates using the same preparation of RBD. E, temperature-induced unfolding of the RBD variants. Bis-ANS fluorescence data as a function of temperature is plotted for the indicated RBD variant. The peak of each curve was taken as the melting temperature. Melting curves were arbitrarily shifted on the y-axis to better visualize differences between the individual plotted curves. F, the melting temperatures for the indicated RBD variant were compared by one-way ANOVA and Tukey's test for multiple comparisons. p value less than 0.05 is indicated by *. The error bars indicate SD. Each RBD variant was measured in triplicate using a unique preparation of RBD. bis-ANS, 4,4 0 -dianilino-1,1 0 -binaphthyl-5,5 0 -disulfonic acid; RBD, receptor-binding domain.
with protein melting temperatures as measured by a fluorescence melting assay that we have applied to the study of ovalbumin variants (Fig. 4E) (20). We observed a reduction in the melting temperature of the beta variant RBD (Fig. 4F).
As discussed above, it has been reported that the substitutions to the spike protein observed in the alpha and beta variants exhibit a higher binding affinity for ACE2. We sought to replicate these observations and examine the effects of the individual K417N and E484K substitutions on ACE2 binding affinity. It has also been reported that K417N reduces the RBD binding affinity for ACE2 by both computational and experimental investigation (14,24). To test the effects of the RBD substitutions on ACE2 binding, we used a binding competition assay where the RBD variants compete with horseradish peroxidase-labeled WT RBD (HRP-RBD) for binding to ACE2 adsorbed to a microplate. As the RBD variants are diluted, HRP-RBD can outcompete for binding to ACE2 resulting in increased absorbance at 450 nm after plate development (Fig. 5A). These data were analyzed by nonlinear regression to calculate LogIC50 values. We found that the K417N substitution significantly reduced ACE2 binding as less dilution was needed for HRP-RBD to outcompete for ACE2 binding (Fig. 5B). E484K exhibited a slightly higher affinity for ACE2 than the WT RBD but this did not reach the level of statistical significance. The N501Y substitution resulted in the increased ACE2 affinity consistent with previous reports, and the presence of all three variant of concern substitutions exhibited the greatest increase in ACE2 affinity compared with the WT RBD (Fig. 5B). These results were consistent with previously reported observations of the effects of variant of concern substitutions on RBD-ACE2 binding (8,9,14,24).

RBD beta variant substitutions alter RBD proteolytic susceptibility
The molecular dynamics simulations and results from unfolding studies presented here lead us to the conclusion that the RBD-VOC substitutions alter the RBD structure and stability. The changes in protein structure and stability are associated with changes in proteolytic resistance (20,21,25,26). Based on these previous observations and those we have reported here thus far, we hypothesized that the RBD substitutions alter proteolytic susceptibility in accordance with changes in stability. Substitutions that stabilize the RBD increase proteolytic resistance, whereas destabilizing substitutions will decrease proteolytic resistance. To test this, we performed limited proteolysis experiments using the lysosomal protease cathepsin S at pH 5.6. The proteolysis reactions were sampled after 0, 15, 30, and 60 min of incubation and analyzed by SDS-PAGE and Coomassie staining (Fig. 6A). The gels were imaged and intensities of the RBD bands were measured with gel analysis software. Band intensities were normalized to the 0 min time point and compared by two-way ANOVA. We observed that the K417N substitution significantly increased proteolytic resistance at all time points compared with the WT RBD (Fig. 6B). The E484K substitution significantly decreased proteolytic resistance compared with the WT RBD except after 60 min of incubation. The N501Y substitution exhibited increased resistance to proteolysis by cathepsin S only after 60 min of incubation relative to the WT RBD. The presence of all three substitutions resulted in a large increase in cathepsin S proteolytic susceptibility after 15 and 30 min, but degradation was similar to that seen for the WT RBD after 60 min (Fig. 6B). Taken together, these results demonstrate that the RBD-VOC substitutions significantly alter RBD proteolytic susceptibility.

Discussion
In this study, we tested the hypothesis that SARS-CoV-2 beta variant substitutions within the spike glycoprotein RBD alter the RBD structure, stability, and ACE2 binding affinity. We studied the RBD by molecular dynamics simulation and found that the K417N substitution alone as well as the presence of all three substitutions, K417N/E484K/N501Y within the beta variant RBD may alter the flexibility of a loop region spanning residues 468 to 488 that forms contacts with ACE2. A similar molecular dynamic simulation analysis of this region was reported for the WT RBD and the E484K, N501Y, and the K417N/E484K/N501Y substitutions, albeit with different observations to those reported here (27). The authors observed that the E484K and N501Y substitutions alone increased flexibility relative to WT whereas both E484K and N501Y together along with all three substitutions (K417N/E484K/ N501Y) did not alter flexibility compared with the WT RBD. Our results differ here, in the CHARMM36 force field, we observed minor changes in the RMSF values for the 468 to 488 region for all the RBD variants whereas in the ff14SB force field, the K417N substitution did not result in large RMSF changes whereas the E484K, N501Y and K417N/E484K/ N501Y substitutions reduced RMSF values relative to the WT RBD which, except for the K417N/E484K/N501Y RBD variant, were not consistent with the results from the previous study. Although the previous study performed their simulations in the ff14SB force field, differences in the equilibration parameters and production run time may account for these discrepancies. It is challenging to say which results are more plausible and comparison to experimental data is difficult as there are no available structures of the RBD in the unbound state.
The major differences observed here between CHARMM36 and ff14SB are the distances between residues 403 and 484 as well as 417 and 484. In CHARMM36 for K417N, and less so for K417N/E484K/N501Y, these residues quickly form close contacts that are not observed for the other variants whereas in   K417N/E484K/N501Y variant in the CHARMM36 force field.
Others have recently reported that the CHARMM36 force field favors more disordered protein conformations whereas the ff14SB force field favors the native state and a stabilization of alpha helical structures which may explain the helix formation within the 468 to 488 loop of the K417N/E484K/ N501Y in ff14SB (28). The biases reported for the CHARMM36 and ff14SB force fields may also explain our observations of the differences in RMSD and radius of gyration, we observed between these two force fields. CHARMM36 predicted greater changes in RMSD and radius of gyration resulting from VOC substitutions than ff14SB. In a binding competition assay, the N501Y and K417N/ E484K/N501Y variants all exhibited a higher binding affinity for ACE2 whereas the K417N substitution alone exhibited reduced binding affinity. The E484K substitution caused a modest but insignificant increase in affinity as measured by the competition assay employed here, and other studies that have reported the E484K substitution results in a negligible increase in affinity (29). Nonetheless, our observations are all in agreement with numerous previous reports on the effects of these substitutions on RBD-ACE2 binding (9,12,14,24,(30)(31)(32)(33). The reduced binding to ACE2 in K417N may result from a decrease in the flexibility of the 468 to 488 region that undergoes a conformational change when bound to ACE2 in the RBD variant containing only the K417N substitution. Analysis of the simulation trajectories in the CHARMM36 force field suggests that glutamate 484 forms a hydrogen bond interaction with arginine 403 when only the K417N substitution is present as these residues are much further apart in simulation data generated from the WT and E484K-RBD variant. The presence of the hydrogen bond between E484 and R403 results in a decrease in the RMSD of the loop spanning residues 468 to 488 in K417N when compared with the WT and E484K RBD variants. In the ff14SB force field, the 468 to 488 loop in the K417N appears to be stabilized by a hydrogen bond between R457 and E484. The 468 to 488 loop region has been reported to undergo a structural transition that is important for ACE2 binding (34). Therefore, a hydrogen bond between E484 and R403 or E484 and R457 that restricts the flexibility of the 468 to 488 loop may be responsible for the reduction in ACE2-binding affinity we observed here but more experiments would be necessary to determine this. More extensive simulation experiments in a previous study have shown that the K417N substitution may abolish a salt bridge of K417 with an aspartate residue at position 30 within ACE2 and therefore reduce the strength of binding interactions (14). If this alone is enough to explain the observed reduction in ACE2 binding by the K417N variant or if it simply contributes remains to be determined.
The utility of the K417N substitution to the virus remains unclear. It has not been reported alone in other variants and this may be due to its negative effects on ACE2 binding reported here and by others (14). K417N also appears to reduce the strength of antibody-RBD interactions and may have been selected for immune evasion despite its potential negative effects on transmissibility. Our data suggest that K417N alone alters the RBD structure in such a way that conformational changes necessary for ACE2 binding may be disfavored relative to the WT or other variants. We also observed that the E484K substitution alone significantly destabilizes the RBD by both denaturant-induced unfolding and limited proteolysis. Therefore, the stabilizing effects of the K417N substitution may be necessary to offset the negative structural effects of the E484K substitution. The N501Y substitution also stabilized the RBD in our studies and the iota variant (B.1.526) has been identified in New York that possesses both the N501Y and E484K substitutions (35). The K417N and E484K substitutions have not yet been reported to exist alone in a SARS-CoV-2 variant, and it is possible their beneficial effects for immune evasion in the case of both K417N and E484K or ACE2 binding in the case of E484K are outweighed by negative structural effects. However, the caveat here is that the RBD does not exist in isolation but rather as part of the much larger spike glycoprotein, but the effects of these substitutions on structure and stability may still be relevant for understanding the evolution of SARS-CoV-2.
The presence of K417N alone and N501Y alone significantly stabilized the RBD relative to WT and it is unclear from the simulation data how this was achieved for the N501Y substitution found in the alpha variant. Surprisingly, the E484K substitution alone significantly destabilized the RBD relative to WT. Our simulation data suggest that E484K alone reduces overall hydrogen bond content in the RBD. Reduced hydrogen bond content may explain the observed decrease in stability and suggests that the elimination of specific interactions between residue 484 and other RBD residues alone cannot explain the observed reduction in RBD stability. There are, however, limitations to the bis-ANS unfolding assay used in this study. This assay is based on the measurement of bis-ANS fluorescence as a function of denaturant-induced protein unfolding. The fluorescent dye bis-ANS is quenched by water and binds to hydrophobic regions of proteins increasing fluorescence. Owing to the binding nature of bis-ANS, the dye may unbind before the protein is completely unfolded as has been seen for the unfolding of ovalbumin (20) and was observed here. Fluorescence signal peaked at 1.25 to 1.5 M guanidine HCl and began to decrease at the RBD became more solvent-exposed. As denaturant concentration increases eventually, the hydrophobic regions of proteins will become completely exposed to water solvent and the dye will unbind and be quenched. However, as was seen for ovalbumin and was reported in this study, bis-ANS binding in the presence of denaturant can still be informative of protein stability. It is also possible that amino acid substitutions may alter the binding affinity of the dye rather than protein stability. We do not consider that to be a major factor for our results as the substitutions studied here were of polar residues or, in the case of N501Y, surface-exposed, and therefore we would not expect those substitutions to alter dye-binding sites. Furthermore, the results from the bis-ANS unfolding experiments were corroborated by limited proteolysis with cathepsin S. The stabilized K417N and N501Y variants were more resistant to cathepsin S proteolysis than the destabilized E484K variant.
In summary, we have reported here that the RBD VOC substitutions alter the RBD structure and stability with K417N and N501Y increasing stability whereas E484K reduces stability. All three substitutions together, as found in the beta variant, exhibit similar stability to the WT RBD albeit with a possibly more open conformation and significantly higher ACE2-binding affinity that is greater than the sum of its parts. Taken together, our findings support the notion that the evolution of the SARS-CoV-2 RBD has been guided by the pressure for increased ACE2 binding and immune evasion within the constraints of maintaining RBD structure for optimal interactions with the ACE2 receptor.

Molecular dynamics simulations
Molecular dynamic simulations were performed with the GROMACS 2020.5 package (17,36) with the CHARMM36 allatom force field (18) or the ff14SB force field (19), and the TIP3P water model (37). The WT SARS-CoV-2 spike glycoprotein RBD (residues 319-541) and the variants were modeled from the full-length spike glycoprotein Cryo-EM structure PDB entry 6VXX (16) using SWISS-MODEL (15). The RBD models were solvated in a dodecahedral box with a minimum protein to edge distance of 1.5 nm and 37,162 water molecules. The system charge was neutralized with approximately 105 Na + and 113 Cl − ions at a concentration of 0.15 M. The assembled systems were minimized for 10 ps and then equilibrated for 100 ps in the NVT ensemble followed by further equilibration for 1000 ps with harmonic position restraints on heavy protein atoms (1000 kJ mol −1 nm −2 ) and Berendsen coupling (38) to maintain temperature and pressure (P 1 bar and T 300 K). Production runs were then performed for each RBD variant for 25 ns in the NPT ensemble. Trajectory data were saved every 1 ps for CHARMM36 simulations and every 0.5 ps for ff14SB simulations. The analysis was performed with GROMACS, VMD 1.9.3 (39), PyMOL, and GraphPad Prism.

Protein expression and purification
The vector pCAGGS containing the receptor-binding domain (residues 319-541) of the spike glycoprotein from SARS-CoV-2, Wuhan-Hu-1 was obtained from the BEI resources. Individual K417N, E484K, and N501Y substitutions were introduced to the coding sequence by site-directed mutagenesis. The RBD variant containing three substitutions (K417N/E484K/N5101Y) was synthesized and cloned by Genscript into the pcDNA3.1 vector and subcloned into pCAGGS. The vectors were used to transform DH5α bacteria for transfection DNA preparation. Transfection grade DNA was purified using the PureLink plasmid midiprep kit (Invitrogen), filtered through a 0.22 μm filter and stored at −20 C. The SARS-CoV-2 spike glycoprotein RBD was expressed in Freestyle 293 cells (ThermoFisher) grown in Freestyle 293 media at 37 C with 8% CO 2 and shaking at 135 rpm in a humidified incubator. The cells were transfected with plasmid DNA and 25 kDa polyethyleneimine (Polysciences) at 1 μg/ml in Freestyle 293 media. Media was harvested four days after transfection and clarified by centrifugation at 4000g for 10 min. His-tagged RBD protein was purified from the cell culture media by nickelnitrilotriacetic acid chromatography and an AKTA pure system. Columns were washed with buffer A (20 mM Sodium Phosphate pH 7.2 and 500 mM NaCl) and the bound protein was eluted with buffer B (buffer A + 250 mM Imidazole). The protein was then desalted using a HiPrep 26/10 desalting column and PBS pH 7.4 and concentrated using a 10 kDa centrifugal concentrator (Milipore-Sigma). Protein concentration was determined by absorbance at 280 nm and an estimated extinction coefficient of 33,850 M −1 cm −1 . The concentrated protein was aliquoted, snap-frozen in liquid nitrogen, and stored at −80 C.

Unfolding experiments
Guanine unfolding experiments were performed by fluorescent dye assay using bis-ANS (Tocris), as described previously (21) with minor modifications. The protein (5 μM) was mixed with bis-ANS (10 μM) in PBS and guanidine HCl (prepared in PBS) ranging from 0 M to 2.5 M in a 96-well black plate in duplicate. The mixtures were incubated for 1 h at room temperature before fluorescence readings were taken using a GloMax Explorer (Promega) plate reader with an excitation filter at 405 nm and emission filter at 500 to 550 nm. Fluorescence data were used to calculate the fraction folded at a given guanidine concentration and analyzed by nonlinear regression (22) to estimate the free energy of unfolding and the best-fit values were compared by one-way ANOVA. Data presented are representative of three independent experiments using the same preparation of RBD. The temperature induced unfolding experiments were performed, as described previously for ovalbumin (20). Briefly, 5 μM RBD protein was mixed with 50 μM bis-ANS in 50 μl of 10 mM sodium phosphate buffer at pH 7.2 in a 200 μl PCR plate, which was then sealed with a plate sealer. The plate was heated from 10 C to 90 C at a rate of 0.016 C per second and bis-ANS fluorescence was recorded. Fluorescence data was plotted as a function of temperature, and the peak was taken as the melting point. Tm values were compared by one-way ANOVA.

RBD-ACE2 binding competition assay
RBD-ACE2 binding competition assay was developed using the SARS-CoV-2 surrogate virus neutralization test kit (Genscript). First, a 5-fold dilution series of RBD variant starting at 10 μM was prepared in sample dilution buffer in duplicate. Next, the serum/antibody incubation step was replaced by incubating RBD-HRP 1:1 with diluted RBD variant at 37 C for 30 min. Then, the kit procedure was followed, as described in the manufacturer's instructions. After development, the absorbance at 450 nm was measured with a GloMax Explorer plate reader. The absorbance values were plotted as a function of RBD-variant concentration and analyzed by nonlinear regression to estimate LogIC50 values. The best-fit values were compared by one-way ANOVA. The data presented are representative of three independent experiments.

Limited proteolysis
Limited proteolysis experiments were conducted with recombinant human cathepsin S (Milipore-Sigma) diluted in phosphate-citrate buffer at pH 5.6 and 37 C (21). Proteolysis reactions were prepared in 150 μl phosphate-citrate buffer pH 5.6, containing 0.5 μg/μl of RBD variant, 0.025 μg/μl Cathepsin S, and 2 mM DTT. The reactions were prepared in triplicate and sampled after 0, 15, 30, and 60 min of incubation. Proteolysis was halted by mixing with SDS-PAGE loading buffer containing 150 mM 2-mercaptoethanol and incubation at 95 C for 5 min. The proteolysis reactions were analyzed by SDS-PAGE and staining with Coomassie blue. The gels were imaged with a Chemidoc MP imaging system (Bio-Rad) and analyzed with ImageJ software. Band intensities of the intact RBD were extracted from images and normalized to the 0-min band intensity before analysis by two-way ANOVA.

Data availability
Molecular dynamics trajectory files are available upon request to the corresponding author Daniel L. Moss (dmoss2@ tulane.edu). All other data are contained within the article.
Supporting information-This article contains supporting information.