key: cord-1001053-8mcylboq
authors: Chen, Ping; Wang, Jingfang; Xu, Xintian; Li, Yuping; Zhu, Yan; Li, Xuan; Li, Ming; Hao, Pei
title: Molecular dynamic simulation analysis of SARS-CoV-2 spike mutations and evaluation of ACE2 from pets and wild animals for infection risk
date: 2021-12-01
journal: Comput Biol Chem
DOI: 10.1016/j.compbiolchem.2021.107613
sha: 85bacb858ed67afb105c32979ccd4cb780e897d2
doc_id: 1001053
cord_uid: 8mcylboq

Coronavirus Disease 2019 (COVID-19) is an ongoing global health emergency that has caused tremendous stress and loss of life worldwide. The viral spike glycoprotein is a critical molecule mediating transmission of SARS-CoV-2 by interacting with human ACE2. However, through the course of the pandemics, there has not been a thorough analysis of the spike protein mutations, and on how these mutants influence the transmission of SARS-CoV-2. Besides, cases of SARS-CoV-2 infection among pets and wild animals have been reported, so the susceptibility of these animals requires great attention to investigate, as they may also link to the renewed question of a possible intermediate host for SARS-CoV-2 before it was transmitted to humans. With over 226,000 SARS-CoV-2 sequences obtained, we found 1573 missense mutations in the spike gene, and 226 of them were within the receptor-binding domain (RBD) region that directly interacts with human ACE2. Modeling the interactions between SARS-CoV-2 spike mutants and ACE2 molecules showed that most of the 74 missense mutations in the RBD region of the interaction interface had little impact on spike binding to ACE2, whereas several within the spike RBD increased the binding affinity toward human ACE2 thus making the virus likely more contagious. On the other hand, modeling the interactions between animal ACE2 molecules and SARS-CoV-2 spike revealed that many pets and wild animals’ ACE2 had a variable binding ability. Particularly, ACE2 of bamboo rat had stronger binding to SARS-CoV-2 spike protein, whereas that of mole, vole, Mus pahari, palm civet, and pangolin had a weaker binding compared to human ACE2. Our results provide structural insights into the impact on interactions of the SARS-CoV-2 spike mutants to human ACE2, and shed light on SARS-CoV-2 transmission in pets and wild animals, and possible clues to the intermediate host(s) for SARS-CoV-2.

resonance (SPR) results, the SARS-CoV-2 spike protein had a nearly 15 nM affinity to the human ACE2, which was approximately 20-fold higher than that of the SARS-CoV spike protein binding to the human ACE2 [12] .

First, as the pivotal viral molecule for human infection, any changes with the SARS-CoV-2 spike protein sequence, especially in its RBD region, could alter the interaction between SARS-CoV-2 and the human ACE2. Thus it is critical for us to investigate the mutational changes with the SARS-CoV-2 spike protein, and understand how the mutations potentially affect the viral infection in humans.

Throughout COVID-19 global pandemics, more than 96,000,000 human infections were recorded. However, only a small fraction of the SARS-CoV-2 viruses from human patients have been sequenced. While mutations have been found in all SARS-CoV-2 coding genes, there has not been a thorough analysis on the spike protein mutations, and on how these mutations affect the interaction between spike protein and human ACE2 molecule, hence as a consequence changing SARS-CoV-2 infect risk.

Second, while many cases of animal infections, including pets and wild animals, by SARS-CoV-2, have been reported [13] [14] [15] [16] , the susceptibility of these animals requires great attention to investigate, as they may also link to the renewed question of what is the possible intermediate host(s) for SARS-CoV-2 before it transmitted in human populations. Almost all the coronaviruses causing human illness have animal origins. Early studies indicated that SARS-CoV and MERS-CoV were transmitted from palm civet and dromedary camel to human, respectively [17] [18] [19] . HCoV-OC43 and HKU1 are considered to have originated from rodents [20, 21] To answer these important questions, in the current study we used more than 220,000 SARS-CoV-2 virus sequences collected from patients worldwide and identified over 2700 mutations in the spike protein gene. Using molecule dynamic simulation, we have analyzed the non-synonymous mutations and assessed their effect on interactions with the human ACE2. Then to study the interactions of ACE2 from pets and wild animals to SARS-CoV-2 spike protein, we constructed structural models of ACE2 from twenty-five different animals and used molecular dynamic simulation to optimize their binding complex models. To assess the binding affinity of these complexes, we also employed the molecular mechanics Poisson-Boltzmann surface area (MM-PB/SA) approach to calculate the binding free energies of the SARS-CoV-2 spike to the animals' ACE2 molecules. Our results provide structural insights into the interactions of the animal ACE2 to SARS-CoV-2 spike mutants and shed light on the transmission of SARS-CoV-2 in pets and wild animals, and possible Table S1 ).

The spike gene sequences from complete SARS-CoV-2 sequences were extracted and mapped to the SARS-CoV-2 spike gene reference, for which 2,782 mutations were found within the spike gene (Additional file 2: Table S2 ). These mutations were further analyzed, and 1573 of them were found to be missense mutations (Additional file 3: Table S3 ). Among the missense mutations, 226 are in the critical RBD region of spike protein (viral genome position 22517-23185) ( Figure 1 and Additional file 3: Table S3 ). The RBD domain is located at the C-terminal of the S1 subunit that directly interacts with host receptors hence facilitating virus infection [29] .

To detect the structural influence of the aforementioned 226 RBD missense mutations, we modeled these missense mutations based on the crystal structure of SARS-CoV-2 spike RBD [10]. Among the RBD missense mutations, 74 are found to be located in the region of interaction interface between RBD and ACE2 (Additional file 3: Table   S3 ). Modeling and MM-PB/SA free energy calculations indicated that most of the missense mutations in the interface had little impact on the spike binding to ACE2.

Interestingly, two mutations, Y449F and Y453F, caused the loss of hydrogen bond interactions between RBD and ACE2 ( Figure 2A) . However, the binding free energies J o u r n a l P r e -p r o o f 7 / 28 for Y449F (-59.55 kcal/mol) and Y453F (-59.17 kcal/mol) did not show a significant reduction, compared to that of the wild-type spike protein (-60.64 kcal/mol), which helps explain why they maintained infection in people. These results indicate a compensation mechanism may exist for the viral spike protein interacting with human ACE2 molecule to maintain the stability of the complex. Interestingly, three mutations, E471Q (binding free energy: -62.34 kcal/mol), S477R (-62.62 kcal/mol), and N501Y (-63.48 kcal/mol), were found to have a decrease in binding free energies, which made them bind more tightly to human ACE2. Notably, N501Y is one of the critical mutations recently found in the dominant SARS-CoV2 variants circulating in the UK (B.1.1.7) and in South Africa (501Y.V2) that became more infectious among people.

However, the E484K mutation found in the South Africa variant 501Y.V2 did not significantly alter binding free energy (-60.06 kcal/mol) by itself. So over the course of the pandemics, several missense mutations within the spike RBD region increased the binding affinity toward human ACE2 thus making the virus likely more contagious, while most of the mutations had little impact on the spike ACE2 interaction.

Both experimental and theoretical evidence indicated that SARS-CoV-2 shared similar receptor recognition and binding mechanism with SARS-CoV due to the similarity in the spike protein sequence [6, 12] . To investigate the interaction of ACE2 from pets and wild animals with SARS-CoV-2 spike protein, we first analyzed the crystal structure of the spike-ACE2 complex (PDB entry 6m0j) [10] for critical residues for ligand recognition. In the binding interface between the SARS-CoV-2 J o u r n a l P r e -p r o o f 8 / 28 spike and human ACE2, three groups (left, middle, and right) of critical residues for ligand recognition with the hydrogen bonding interactions were present ( Figure 2A ).

The hydrogen-bonding interactions between the ACE2 residues Asp38 and Tyr41, and the spike residues Tyr449 and Gln498 defined the recognition group on the left side of the binding interface. These hydrogen bonds are also known to exist in the binding interface of the SARS-CoV spike and human ACE2. The middle recognition group was composed of two hydrogen-bonding interactions, i.e. Tyr453-His34 and Gln493-Glu35. These two hydrogen bonds that also existed in the SARS-CoV spike-ACE2 complex, were essential for the transmission of SARS-CoV from animal to human [4, 30] . Mutations in these residues, for example, N479L and N479R in the SARS-CoV spike protein, significantly reduced the binding affinities to the human ACE2 [31, 32]. The recognition group on the right side contained two hydrogen bond interactions between the ACE2 residues Ser19 and Tyr83, and the spike residues Ala475 and Tyr489, which were important for the spike to recognize and bind to the human ACE2.

To estimate their binding affinity and the contribution for each of the critical residues, we computed the binding free energies of the SARS-CoV-2 spike to ACE2 by MM-PB/SA approach (see Materials and methods). The binding free energies for SARS-CoV spike binding to human ACE2 were also computed as a comparison. As a result, the binding free energy of SARS-CoV-2 spike to the human ACE2 was -60.64 kcal/mol on average, a little lower than that of the SARS-CoV spike binding to the human ACE2 (-58.55 kcal/mol ±8.75), indicating that a stronger affinity of SARS-CoV-2 than SARS-CoV binding to the human ACE2. This observation is in good agreement with the results from the surface plasmon resonance (SPR) study, which estimated a value of 14.7 nM for the affinity of SARS-CoV-2 spike to the J o u r n a l P r e -p r o o f human ACE2 [12] . Energy decomposition analysis indicated that the ACE2 residues Tyr41 and Tyr83 contributed most to the spike binding affinity for both SARS-CoV-2 and SARS-CoV ( Figure 2B ). Basing on the crystal structure analysis, these two tyrosine residues belong to the critical recognition groups in the binding interface of the spike and human ACE2. Besides, Met82, Gly354, and Asp355 also had substantial contributions to the binding free energy between spike protein and human ACE2.

These residues were considered to be critical and their alteration in ACE2 of other animals will greatly impact the interaction of SARS-CoV-2 spike and animals' ACE2 molecules, thus affecting the dynamic of animal infection by SARS-CoV-2.

We obtained the sequences of the ACE2 gene for twenty-five pets or wild animals from GenBank at NCBI (Additional file 4: Table S4 ), some of which were suggested as the possible intermediate host for transmission to humans. After constructing the structural models for the ACE2 of these animals, we found many critical residues were replaced that resulted in changes in the binding affinity of the spike-ACE2 complex, thus changing the dynamic of animal infection by SARS-CoV-2 (Additional file 4: Table S4 ). Interestingly, their estimated binding free energies to the SARS-CoV-2 spike protein exhibited a large variation in comparison to the human ACE2 binding to the SARS-CoV-2 spike ( Table 1 ). The ACE2 from civet, pangolin, wild mouse, Neovison vison, and snake were further characterized to evaluate these animals' infection risk.

Based on the receptor-binding complex models constructed for the SARS-CoV-2 J o u r n a l P r e -p r o o f spike binding to civet and pangolin ACE2, we calculated the corresponding binding free energies. As a result, the binding free energies of SARS-CoV-2 spike protein to civet and pangolin ACE2 receptors were -56.11±2.42 kcal/mol and -54.78±2.22 kcal/mol, showing a lower binding affinity than human ACE2 (Table 1 ). This observation indicated that the SARS-CoV-2 spike protein is weakly bound to both civet and pangolin ACE2 receptors. The weak binding of the spike protein to the civet ACE2 was also supported by the virus-infected experiments in HeLa cells [3] .

Compared with the human ACE2 molecule, the civet ACE2 had a substitution from phenylalanine to serine at position 40 in the first helix at the N-terminal lobe region. This mutation broke the pi-pi stacking interactions between the N-terminal helices, resulting in the side-chain rearrangement and loose helix-packing structures ( Figure 3A ). These structural alternations disrupted almost all the hydrogen bonds in the left recognition group, only leaving the ones formed by the civet ACE2 residue Lys352 and the viral residues Gln498 and Asn501 ( Figure 3A ). This observation was also supported by the data of residue contribution to the binding free energy. The energy contribution of the civet ACE2 residue Tyr41 to the binding free energy was increased to -2.96 kcal/mol, which was -4.51 kcal/mol in the human ACE2 receptor upon the SARS-CoV-2 spike binding (Additional file 5: Table S5 ). Additionally, the salt bridge between Lys31 and Glu35 disappeared in the civet ACE2 receptor due to the Asp-to-Glu transformation at position 30, leading to an unfavorable residue contribution of Lys31 to the binding free energy (4.32 kcal/mol for the civet ACE2 receptor and -0.98 kcal/mol in the human ACE2 receptor) (Additional file 5: Table   S5 ).

Similar to the civet ACE2 receptor, there was also an Asp-to-Glu transformation in the pangolin ACE2 (at position 38 in the left recognition region). This J o u r n a l P r e -p r o o f transformation had a significant impact on the hydrogen bond network in the left recognition region. In the pangolin ACE2, the ACE2 residue Glu38 did not form hydrogen bonds with the viral residue Tyr449, which turned to have a hydrogen-bonding interaction with the ACE2 residue Gln42 ( Figure 3B) . Besides, the viral residue Gln498, which interacted with the ACE2 residue Tyr41 in the human ACE2, also formed a hydrogen bond with the pangolin ACE2 residue Lys353. In the middle recognition region, only one hydrogen-bonding interaction was detected between the ACE2 residue Ser34 and the viral residue Gln493. Also, the His-to-Ser transformation at position 34 in the pangolin ACE2 reduced its energy contribution to the binding free energy (-2.34 kcal/mol in the pangolin ACE2 and -3.36 kcal/mol in the human ACE2) (Additional file 5: Table S5 ).

Rodents have been regarded as intermediate hosts or carriers of many human pathogens [33], and murine hepatitis virus (MHV) has been used as a model system to study coronaviruses [34] . Biological studies on alpha-and beta-coronaviruses in rodents indicated the murine origin of some beta-coronaviruses [35, 36] . We constructed the ACE2 of the wild mice from Yunnan province, i.e., mole, vole, rat, mice, and bamboo rat, and employed free energy calculations to estimate their binding affinities to the SARS-CoV-2 spike. Most mouse ACE2 receptors (including rats and mice) showed higher binding free energies than the human ACE2 (Table 1) , indicating weaker binding to the SARS-CoV-2 spike. However, the mole ACE2 had a similar binding free energy to the human ACE2 due to their high sequence identity at the binding interface to the SARS-CoV-2 spike.

Notably, the bamboo rat ACE2 showed a lower binding free energy (-65.04 J o u r n a l P r e -p r o o f kcal/mol on average) than the human ACE2 (-60.64 kcal/mol on average), indicating a stronger inclination for the SARS-CoV-2 spike binding to the bamboo rat ACE2.

Looking at their binding interface, the major differences between bamboo rat and human ACE2 molecules were caused by the His-to-Gln transformation at position 34 and the Asp-to-Glu transformation at position 38 in the N-terminal helix (Table 1) .

These changes in the bamboo rat ACE2 led to a side-chain rearrangement at the binding interface to the SARS-CoV-2 spike. As a result, the hydrogen bonding between the human ACE2 residue Asp38 and the viral residue Tyr449 disappeared in the bamboo rat ACE2 receptor, and instead, the bamboo rat ACE2 residue Gln34 had a hydrogen-bonding interaction with the viral residue Ser452 (Figure 4) . Furthermore, the bamboo rat ACE2 residue Glu35 formed two hydrogen bonds with the viral residues Gln493 and Phe488, and the aromatic sidechain of the bamboo rat ACE2 residue Tyr41 employed a pi-pi stacking interaction with the viral residue Tyr449 (Figure 4 ). These additional interactions enhanced the SARS-CoV-2 spike binding to the bamboo rat ACE2.

Due to the strong positive binding free energies (> 400 kcal/mol), the SARS-CoV-2 spike protein did not bind to Neovison vison and snake ACE2 ( Table 1 ). The unfavorable binding free energies for Neovison vison and snake ACE2 molecules were caused by their structural gaps in the N-terminal helices, which were essential for receptor recognition. According to the sequence alignment (Additional file 7: Figure S1 ), a large segment in the N-terminal region (residues 1-318 in the human ACE2 receptor) was missing in the Neovison vison ACE2. The residues 19-52 in the human ACE2 contained the important residues interacting with the receptor-binding J o u r n a l P r e -p r o o f motif of the SARS-CoV-2 spike protein, which was stabilized by the second helix in the N-terminal lobe region (residues 53-84).

The unfavorable binding free energy for the snake ACE2 was caused by the large dissimilarity to the human ACE2 in protein sequence at the N-terminal lobe region (Table 1 ). In the human ACE2, two virus-binding hot spots have been identified at the binding interface of the SARS-CoV receptor-binding domain (RBD). One centers on the ACE2 residue Lys31 (hot spot 31), forming a salt bridge with the ACE2 residue Glu35 ( Figure 2B ). Another is on the ACE2 residue Lys353 (hot spot 353), forming a salt bridge with the ACE2 residue Asp38. However, the Lys31 and Glu35 in human ACE2 for the virus-binding hot spot 31 were mutated to glutamine (Gln56) and arginine (Arg60) in the snake ACE2, which broke the significant salt bridge. Besides, the side-chain conformation in the snake ACE2 receptor had an adverse effect on the SARS-CoV-2 spike binding. For example, the snake ACE2 residues Asp64 and Glu352 showed strong van der Waals conflicts with Tyr488 and Tyr505 in the SARS-CoV-2 spike.

We modeled SARS-CoV-2 spike-ACE2 complexes and calculated binding free energies to evaluate spike protein mutations occurring during COVID-19 pandemics, and assess their binding affinities. We found 1573 missense mutations in the spike gene, of which 226 were within the RBD region that interacts with the human ACE2.

We found 74 of them in the region of interaction interface between RBD and ACE2.

Modeling the interactions between SARS-CoV-2 spike mutants and ACE2 molecules indicated that several missense mutations within the spike RBD region increased the binding affinity toward human ACE2 thus making the virus likely more contagious, 

The structural models of SARS-CoV-2 spike mutants were constructed based on the crystal structure of the wild type SARS-CoV-2 spike (PDB entry 1m0j) by MODELLER. The complexes of SARS-CoV-2 mutants and ACE2 were then modeled by fitting the SARS-CoV-2 spike mutant into the crystal structure of the SARS-CoV-2 spike-ACE2 complex (PDB entry 6m0j) [10]. These complexes were finally minimized by a short molecular dynamics simulation (about 1 nanosecond (ns)) with the protein backbone fixed, and then without any limitation for another 1 ns.

The ACE2 protein sequences for human (Homo sapiens), and animals: civet (Paguma 

The binding affinities of the SARS-CoV-2 spike to the ACE2 receptors were estimated by the MM-PB/SA method [43] . In the MM-PB/SA method, the binding free energy of the SARS-CoV-2 spike (deltaG) was calculated as the energy difference between the free energies of ACE2, SARS-CoV-2 spike, and their complex (Equation 1).

As shown in Equation (2), these free energies were estimated by summing the internal energy in the gas phase (E gas ), solvation free energy (G sol ), and a vibrational entropy term (TS conf ). The internal energy in the gas phase (E gas ) was a standard force field energy computed from Equation (3) by the strain energies from covalent bonds (E bond and E ang ), torsion angles (E tor ), van der Waals (E vdw ), and electrostatic energies (E ele ).

The solvation free energy (G sol ) was calculated by electrostatic (G ele ) and non-polar energies (G non-polar , Equation (4) . The binding interface of the SARS-CoV-2 spike protein to the bamboo rat ACE2 receptor. Top panel: sequence alignment between the interface residues in human and bamboo rat ACE2 receptors. Bottom panel: hydrogen bonding interactions in the recognition groups. The hydrogen bonds in the recognition groups are labeled as gray dash lines: right group, Gln24-Asn487, Lys31-Asn487 and Glu76-Tyr489; middle group, Gln34-Ser494, Glu35-Phe490 and Glu38-Gln493; left group, Tyr41-Gln498, and Asp324-Gly446.

A novel coronavirus outbreak of global health concern

A pneumonia outbreak associated with a new coronavirus of probable bat origin

Structure of SARS coronavirus spike SARS-CoV-2 by full-length human ACE2

Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation

SARS-CoV-2 infection, disease and transmission in domestic cats

Infection of dogs with SARS-CoV-2

From people to Panthera: Natural SARS-CoV-2 infection in tigers and lions at the Bronx Zoo

Susceptibility of ferrets, cats, dogs, and other domesticated animals to SARS-coronavirus 2

Isolation and characterization of viruses related to the SARS coronavirus from animals in southern China

Evidence for camel-to-human transmission of MERS coronavirus

Evidence for camel-to-human transmission of MERS coronavirus

Molecular Evolution of Human Coronavirus Genomes

We would like to thank Dr. Zhang Zhang for the helpful discussion concerning the

Additional file 1: Table S1 . Source of SARS-CoV-2 virus sequences.Additional file 2: Table S2 . 2782 mutations found within spike gene of SARS-CoV-2.Additional file 3: Table S3 . 1573 spike mutations found to be missense mutations.Additional file 4: Table S4 . Critical amino acid residues of ACE2 from twenty-five pets or wild animals.Additional file 5: Table S5 . ACE2 residue contributions to binding free energy (kcal/mol). Additional file 6: Table S6 . GenBank accessions for the ACE2 receptors used in the current study. Additional file 7: Figure S1 . Sequence alignment between human (Homo sapiens) andNeovison vison ACE2.