SARS-CoV-2 Omicron subvariant spike N405 unlikely to rapidly deamidate

The RGD motif on the SARS-CoV-2 spike protein has been suggested to interact with RGD-binding integrins αVβ3 and α5β1 to enhance viral cell entry and alter downstream signaling cascades. The D405N mutation on the Omicron subvariant spike proteins, resulting in an RGN motif, has recently been shown to inhibit binding to integrin αVβ3. Deamidation of asparagines in protein ligand RGN motifs has been demonstrated to generate RGD and RGisoD motifs that permit binding to RGD-binding integrins. Two asparagines, N481 and N501, on the Wild-type spike receptor-binding domain have been previously shown to have deamidation half-lives of 16.5 and 123 days, respectively, which may occur during the viral life cycle. Deamidation of Omicron subvariant N405 may recover the ability to interact with RGD-binding integrins. Thus, herein, all-atom molecular dynamics simulations of the Wild-type and Omicron subvariant spike protein receptor-binding domains were conducted to investigate the potential for asparagines, the Omicron subvariant N405 in particular, to assume the optimized geometry for deamidation to occur. In summary, the Omicron subvariant N405 was primarily found to be stabilized in a state unfavourable for deamidation after hydrogen bonding with downstream E406. Nevertheless, a small number of RGD or RGisoD motifs on the Omicron subvariant spike proteins may restore the ability to interact with RGD-binding integrins. The simulations also provided structural clarification regarding the deamidation rates of Wild-type N481 and N501 and highlighted the utility of tertiary structure dynamics information in predicting asparagine deamidation. Further work is needed to characterize the effects of deamidation on spike-integrin interactions.

Asparagine residues in RGN/NGR motifs have been discovered to undergo spontaneous deamidation through formation of a succinimide intermediate and subsequent hydrolysis to generate isoDGR motifs (Fig. 1D) that may also interact with RGD-binding integrins, such as avb3 [11,12]. Two solvent-exposed asparagine residues on the Wild-type spike protein receptor-binding domain (RBD), N481 and N501, have been experimentally determined to have deamidation half-life rates of 16.5 and 123 days, respectively [13]. Thus, deamidation of the D405N mutation may restore the ability of the spike protein to interact with RGD-binding integrins. Although molecular dynamics simulations have been previously utilized to detect changes in structural conformations that permit deamidation of asparagine residues, no study has yet investigated such tertiary structural effects on coronavirus spike proteins [14].
Thus, herein, the conformations of asparagine residues, particularly D405N, in the WT and O 0 RBDs were analyzed to predict the structural potential for succinimide formation. Hydrogen bonding with downstream E406 was found to stabilize the O 0 N405 residue away from the neighboring backbone nitrogen atom, suggesting that the N405 is unlikely to rapidly undergo deamidation. Nevertheless, even low half-life rates may result in a small number of deamidated sites per virion, thus potentially restoring the ability to interact with RGD-binding integrins. The experimentallydetermined deamidation rates of asparagines on the WT RBD were further explained and supported by the data, and the impact of tertiary structure dynamics on additional asparagine residues on the WT and O' spike RBDs was clarified. Additional structural validation is necessary to better understand deamidation rates on coronavirus spike proteins and potential changes in interactions with integrins.

Results
Sequence-based calculations using NGOME-Lite [15] suggest deamidation half-lives of the O 0 N405 and N481 residues at 23.7 and 29.3 days, respectively. These results, however, should be taken as rough estimates since the WT N481 and N501 residues were predicted to have half-lives of 27.0 and 20.2 days, respectively, and were experimentally verified at 16.5 and 123 days, respectively, at 37 C [13]. Thus, initial sequence-based predictions indicate that the O' N405 may be deamidated during a time relevant for infection.
To examine the structural potential for deamidation of the asparagine residues, molecular dynamics simulations were used to determine changes in residue geometry that may permit succinimide formation [16e18]. Dihedral phi and psi angles between 40-60 and 35e50 , respectively, have been suggested to be predictive of succinimide formation, while residue flexibility and solvency have been debated regarding their effects on predicting deamidation [19,20]. The WT spike RGD motif has been previously suggested to be solvent accessible [21]. The WT spike RBD (PDB: 6m0j) N481 and N501 residues were analyzed alongside the O' spike RBD (PDB: 7zxu) N405 and N481 residues as controls and to compare with published experimental data [13].
Three replicates of 1 ms-long all-atom molecular dynamics simulations were conducted for both the WT and O 0 RBDs. The root mean square deviation (RMSD) values were, first, calculated to establish the time of convergence of the simulations ( Fig. 2A and B). All runs, with the exception of the third WT run, converged after 250 ns. Subsequently, the root mean square flexibility (RMSF) for each residue in the RBDs was assessed. The runs reported similar RMSF values with the exception of the third simulation WT RBD residues 440e450 and residues near the C-and N-terminal residues. Since the RMSF calculated for all residues except 440e450 during the third WT simulation were consistent with the other WT two runs, the data from all six productions were included ( Fig. 2C and D). The WT and O 0 N481 residues were found to be in highly flexible regions with average RMSF values of 3.5 ± 0.54 and 3.80 ± 0.18 Å, respectively, while both WT N501 and O 0 N405 were discovered to be largely inflexible with RMSD values of 1.19 ± 0.25 and 0.88 ± 0.01 Å, respectively. These data suggest that the WT and O 0 N481 residues may more likely assume the optimal geometries for deamidation, while the WT N501 and O' N405 may only undergo rapid deamidation if naturally found in the correct conformation.
The dihedral phi and psi angles of the productions were calculated, although only the dihedral psi angles were found to be indicative of the geometries conducive to succinimide formation when compared with experimentally-determined deamidation. Thus, the dihedral psi angles were used as a measure for deamidation potential. The WT and O 0 N481 residues were found to adopt the preferred geometry for, on average, 9.4 ± 8.34 and 56.31 ± 35.95 ns, respectively ( Fig. 3B and C and 3K-L, respectively; Table 1). The WT N501 residue was not captured in the optimal state at any time during the simulation replicates ( Fig. 3E and F). Visualization of the simulations reveal that the WT and O' N481 residues are facing outwards in a flexible region with minor steric constraints on adoption of the preferred conformation ( Fig. 3A and J, respectively). Inspection of the N501 residue revealed that the side chain amine may form hydrogen bonds with the backbone carbonyl of Y505, thus stabilizing the side chain group away from the Nþ1 amine necessary for deamidation (Fig. 3D). The WT results reflect the conclusions by Lorenzo et al. that the WT N481 may undergo deamidation during the viral life cycle while N501 most likely does not and such deamidation has no particular advantageous role e which is further evident in the N501Y mutation [22]. These data also highlight the predictive capability of structural dynamics for better understanding asparagine deamidation potential.
Interestingly, the O 0 N405 residue was found to be in the state that undergoes succinimide formation during the first 50e250 ns of the simulations and, subsequently, assumed a state that primarily did not allow such conformations for the remainder of the simulations e only 0.027 ± 0.0058 ns ( Fig. 3H and I). Visualization of the simulations revealed that the carboxylic group of the downstream E406 eventually forms a hydrogen bond with the amine of the N405 residue, thus stabilizing the group outside of the permitted range for deamidation (Fig. 3G). The starting conformations are likely to be artefacts of the crystallization process [23]. A simulation of the E406A O 0 mutation was, thus, run to test for the effect of this putative hydrogen bonding. The resulting flexibility of O 0 N405 was slightly increased to 0.95 Å, which reflects a loss of stabilization due to hydrogen bonding. The mutation was found to increase the time for which N405 side group was at an the optimal psi angle to 0.88 ns. These data suggest that both hydrogen bonding with E406 and backbone constraints hinder N405 from structurally adopting the preferred deamidation geometries (Fig. 4AeC). The pKa of E406 is predicted to be 3.21 and 3.3 using the DelPhiPKa [24] and Hþþ 3.0 [25] web servers, respectively, which are below physiologic and lysosome pH. Thus, potential protonation of E406 is not likely to occur or affect hydrogen binding with N405. These results suggest that the O 0 N405 does not rapidly undergo deamidation, as in the case of WT and O' N481, and its half-life rate more likely resembles that of the WT N501Y residue.
The geometries of all other surface-exposed asparagines near the receptor-binding motif were measured to note additional evolutionary differences in deamidation between the WT and O 0 RBDs (Fig. 4D) [26]. Comparing WT and O 0 RBD asparagine content not discussed above, one mutation, N440K, removed and two mutations, K417N and S477N, added an asparagine to the spike RBD, leading to one additional asparagine in the O 0 spike RBD. The reported dihedral psi angles revealed that the WT N440 and WT/O 0 N439, N450, and N487 residues are likely to undergo deamidation while the WT/O 0 N437, N448, and N460 and O 0 N417 may not (Table 1). Most strikingly, although N450 was largely found to be able to adopt the preferred geometry in the WT RBD, it was diminished in the O 0 RBD. Inspection of the WT/O 0 N450 and surrounding residues revealed that the Q498R mutation may affect the conformation of the backbone of N450 (Fig. 4E). Hydrogen bonding between WT Q498 and Y449 pulls the backbone of the N450 residue further from potential steric clashes with the Y451 residue. The O 0 R498 side chain hydrogen bonds, instead, with the backbone carbonyl of V445 and permits the backbone of N450 to pack closer to the RBD, thus imposing clashes with Y451 with respect to the geometry required for succinimide formation. Furthermore, hydrogen bonding between the sidechain amide of WT N448 and the backbone carbonyl of D442 was found to stabilize the amide away from the neighboring backbone amine. This hydrogen bonding was also disrupted by the backbone of N450 backbone moving inwards as a result of the Q498R mutation, thus increasing the potential for the O' N448 to be deamidated. The discovery of such putative local and long-range interactions further emphasizes the effect of tertiary structure on deamidation potential.

Discussion
In summary, structural investigations suggest that the N405 does not undergo rapid deamidation. These data also provide mechanistic details at atomic resolution into the experimentallydetermined deamidation rates of asparagines on the WT RBD. Notably, two asparagines, N481 and N501, are positioned directly upstream of a glycine, which is thought to be predictive of deamidation potential, although only N481 was found to rapidly deamidate [27]. The tertiary positioning of downstream charged residues, such as glutamate in the context of O 0 N405, has been previously found to impact deamidation potential [28]. For instance, a human interleukin 1a asparagine is directly upstream of an aspartate, but the aspartate is pointing inwards towards the core of the protein and away from the asparagine, thus having no effect [29]. Additionally, presence in an alpha helix, such as in the case of O' N405, has been found to impede deamidation, which is further  supported by these data [14]. These findings and those in this study further highlight the importance of inspecting tertiary structure and its dynamics to understand and predict asparagine deamidation potential. Furthermore, rapid deamidation of residues found in solvent-exposed epitopes, such as the N477, may provide an evolutionary benefit through antibody evasion [13]. Notably, other factors, such as pH and solvent catalysts, may have an impact on deamidation rates and, thus, should be investigated further. Additional work is necessary to characterize the extent to which coronavirus spike proteins undergo deamidation and its effect on protein-protein interactions. Nevertheless, it should be noted that even a small deamidation rate per receptor-binding domain, may not necessarily imply functional irrelevance. This is because there are approximately 40 spike proteins per virion, with three receptor-binding domains per spike protein, i.e., about 120 receptor-binding domains per virion [30]. Assuming that the conversion of the RGN to an RGD or RGisoD sequence follows a first order exponential equation, the fraction of receptor-binding domains that contains an RGN motif after time t (F (RGN)t ), is given by: F (RGN)t ¼ 2 -(t/t1/2) , where t 1/2 is the half-life for deamidation. Or, equivalently, the fraction of RGN motifs converted to RGD/RGisoD after time t (F (RGisoD)t ) is given by: F (RGisoD)t ¼ 1-2 -(t/ t1/2) . Assuming that the deamidation of RGN to RGD/RGisoD occurs independently and stochastically and that this individual probability per unit time is low (given that the molecular dynamics simulations suggest a long half-life), the conversion can be modelled by Poisson statistics. Hence, the fraction of viruses with at least one receptor-binding domain that undergo deamidation at this site after time t is: 1 -e (À120 * F(RGisoD)t) . At time ¼ t 1/2 , by definition, only half of the receptor-binding domains have deamidated at this site. But the fraction of viruses that has at least one receptorbinding domain with an RGD/RGisoD sequence is: 1 -e À60 , which is essentially 1. SARS-CoV-2 has been found to be stable on several surfaces, such as plastic, glass, and latex gloves, for at least seven days [31]. Calculating the longest reported half-life of 123 days, based on the findings by Lorenzo et al., for seven days, then: .0039. Thus, after seven days, the fraction of viruses that are predicted to have at least one RGD/RGisoD is: 1-e (À7 * 120)~0 .99 (Table 2). Even after one day, approximately half of a given sample population of viruses will have at least one deamidated N405.
Thus, even the potential slow deamidation of O' N405 may generate a mixture of RGN, RGisoD, and RGD motifs on the virion spike proteins that may restore the ability to bind with integrins. Further work is needed to characterize the capacity for deamidated SARS-CoV-2 subvariant spike proteins to interact with RGD-binding integrins. The investigation of potential non-RGD-based interactions between integrins and SARS-CoV-2 subvariant spike proteins remains warranted [21,32]. These findings re-emphasize  the relevance of integrin-based therapies for SARS-CoV-2 and other viral infections [9,33,34]. Monitoring of mutations and their potential impacts on protein-protein interactions may reveal novel insights into viral pathogenesis and tropism [35].

Molecular dynamics simulations and analyses
All-atom molecular dynamics simulations were prepared and run using GROMACS 2021.3 [45] using the University of Cambridge High Performance Computing resources. Periodic boundary conditions were established in a cubic box with the smallest distance between the protein and box edge set at 2 Å. All systems were solvated in 150 mM NaCl at a zero net charge using the TIP3P model [46] and were run at 300 K [47]. The Particle-mesh Ewald method was used to calculate long-range electrostatic interactions, and the cutoff for Coulomb interactions and van der Waals interactions were set to 10 Å [48]. Following a 100 ps steepest descent minimization, all systems were subjected to an NVT equilibration ensemble with temperature coupling using velocity rescaling [49] for 100 ps and an NPT equilibration ensemble with pressure coupling using the Parrinello-Rahman method [50] for 100 ps. All simulations were run using the OPLS-AA/L all-atom force field with 2 fs time steps [51]. Triplicate productions of 1 ms each were run for both the WT and O 0 RBDs, and one 1 ms production was run for the O' E406A mutant model. VMD was used to visualize snapshot structures from the productions [52]. The last 250 ns e after convergence of the simulations as evidenced by the RMSD values e of the runs were analyzed to account for the influence of crystallization on protein structure [23,53]. The RMSD, RMSF, and dihedral phi and psi angles were extracted from the simulation trajectories using the gmx rms, gmx rmsf, gmx rama GROMACS modules, respectively [54]. The data were plotted using R version 3.6.3.

Data availability
Data reported in this manuscript is available upon request.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.