We proposed a novel computational framework to trace susceptible animals to SARS-CoV-2 (Figure 1). First, we aligned the amino acids sequence of the M2 region of hACE2 protein against the NCBI non-redundant protein sequences (nr) database. Based on predefined filter rules, we identified 31 species to model the structure of their ACE2 proteins and calculated the binding affinities between their ACE2 proteins and spike protein for five SARS-CoV-2 variants, respectively. It was reported that three species, ferret, cat, and tiger, were experimentally or naturally infected with SARS-CoV-2 (McAloose et al., 2020; Shi et al., 2020; Zhang et al., 2020). Thus, we used the docking scores of the three species as the threshold values to determine susceptible animals to SARS-CoV-2. In addition, we analyzed the binding sites of a few representative species to explain how spike proteins interact with their ACE2 proteins. Finally, we simulated some mutations on human ACE2 protein and constructed the mutational human ACE2 protein model. The binding affinity prediction was performed on the mutated proteins.

Yan et al. (2020) reported that SARS-CoV-2 spike protein direct bonds to an M2 region of hACE2 protein. Thus, we extracted the amino acid sequence of the M2 region of hACE2 and aligned it with the NCBI non-redundant protein sequences (nr) database. There are 407 species with alignment identity greater than or equal to 70% (Supplementary Data Sheet 1), which all belonged to the Chordata phylum and contained 326 genera.

We downloaded the ACE2 protein sequences of the 407 species from NCBI (Supplementary Data Sheet 2) and aligned the amino acids sequences of the M2 region of their ACE2 proteins. By setting tiger, cat, domestic ferret as thresholds and a few other criteria (see section “Materials and Methods”), we selected 30 species plus human (overall 31 species) for subsequent analysis (Supplementary Data Sheet 3). A phylogenetic tree was constructed based on the ACE2 proteins of these species using MEGA7 (Figure 2). The phylogenetic tree showed that these species clustered four main orders: Primates, Rodentia, Artiodactyla, and Carnivora. The amino acid sequence of Primate ACE2 protein is the closest to that of humans. The suspected host bat (Rhinolophus sinicus) and pangolin (Manis javanica) were far away from the human, indicating that there might be some intermediate host if they are a genuine original host.

We selected the amino acid sequences of the M2 region within ACE2 protein across 30 species (Supplementary Data Sheet 4) to predict their protein structures using the I-TASSER pipeline. We calculated the sequence and structural similarity of human ACE2 protein with that of the 30 species (Table 1). We performed a Spearman correlation test between the two vectors to test the consistency between sequence and structural similarity. The sequence and structure similarity metrics (including TM score, GDT-TS score, and RMSD) have a correlation of 0.379, 0.406, and −0.407 with a p-value of 0.039, 0.026, and 0.026, respectively, indicating that the two vectors are correlated thus statistically consistent. The top three species with similar structures to human ACE2 were Vulpes vulpes (TM score: 0.9831), Pan paniscus (TM score: 0.983), and Ursus arctos horribilis (TM score: 0.9828). However, the top three species ranked by sequence similarity are Gorilla gorilla gorilla (99.48%), Pan paniscus (99.48%), and Pan troglodytes (99.48%). Interestingly, Rhinolophus sinicus only has a sequence identity of 80.28% with the human in the M2 region of ACE2 protein, which ranks the last in the 30 species. However, its TM score is 0.9529, ranking 15th. These results indicated that variations in more amino acids in the protein sequence do not necessarily lead to more remarkable changes in structure.

The HADDOCK score estimates the stability of the complex of ACE2 protein binding with viral spike protein. A complex with HADDOCK scores less than 0 is considered stable. We showed the HADDOCK scores of the ACE2 protein from 31 studied species docking with spike protein in SARS-CoV-2 (Wuhan) and Omicron in Table 2. Equus caballus ranked first in the list, followed by Oryctolagus cuniculus for SARS-CoV-2 (Wuhan). In contrast, Panthera tigris altaica ranked first, followed by Mus musculus for Omicron. Previous studies have suggested that SARS-CoV-2 could infect Homo sapiens, Panthera tigris altaica, Mustela putorius furo, Felis catus and Rhinolophus sinicus rank 20, 10, 15, 5, and 6, respectively in the list of Wuhan. Thus, we suspect that the top 20 species are potentially susceptible animals for SARS-CoV-2 Wuhan. The 20 species include wild animals, husbandry animals, and pets, indicating that the infection of SARS-CoV-2 is unbiased across these three categories of animals.

There is no significant correlation between ACE2 sequence similarity and ACE2-RBD structural similarity. The similarity of M2(ACE2)-RBD(SARS-CoV-2) structure between humans and other species was evaluated using RMSD and TM score, and the result is shown in Table 3. Surprisedly, compared with M2(ACE2)-RBD(SARS-CoV-2 Wuhan) structural similarity between human beings with other 31 species, M2(ACE2)-RBD(Omicron) decreased for most of the 31species.

HADDOCK score of ACE2 protein docking with spike protein of Omicron decreased sharply for most 31 species studies. The top-eight species were Manis javanica, Mus musculus, Panthera tigris altaica, Ursus arctos horribilis, Mustela putorius furo, Mesocricetus auratus, and Homo sapiens (Table 2). In addition, a phylogenetic tree for 7 SARS-CoV-2 variants was constructed with spike protein sequences (Figure 3), showing an intriguing evolutionary relationship between Omicron with other variants that evolved in human patients. We speculated that the progenitor of Omicron rapidly evolved by accumulation of mutations conductive to infecting that host and then jumped into human beings.

Zhou et al. (2020) reported a cryo-electron microscopy structure of the complex, in which human ACE2 protein binds with the receptor-binding domain (RBD) of the viral spike protein. The structure indicated that human ACE2 protein interfaced with the RBD of viral spike protein by seven residues: Q24, D30, H34, Y41, Q42, K353, and R357 of ACE2 interact with Q474, K417, Y453, N501, Q498, N501, and T500 of RBD.

At the N terminus of cat ACE2 protein (Figure 4A), we found that the ACE2 protein binds with RBD by three H-bonds. The THR739, TYR780, and ASN770 of cat ACE2 interacted with ASP428, ASP389, and LEU518 of RBD (Figure 4B). These H-bonds were less than 3 Å in the distance and provided a strong binding force.

The HADDOCK score of Mustela putorius furo was larger than that of Felis catus, and its ACE2 protein had six H-bonds with spike protein (Figure 4C). SER459 of RBD bond with GLU799 (2.8 Å) at the N terminus of the ACE2 of Mustela putorius furo. At the ACE2 loop, THR740 (2.334 Å), and ILE741 (2.532 Å) of ACE2 interfaced with LYS417 of the spike protein. LYS417, GLN474, and TYR489 of RBD loop bond with ILE741 (2.609 Å), ILE793 (2.626 Å), and PRO135 (2.638 Å) of ACE2.

Molecular docking between 31 ACE2 proteins and RBD of the Spike protein of SARS-CoV-2 were conducted for five strains (including Wuhan, Gamma, Delta, and Omicron), respectively. Of 31 species, 20 species exhibited the lowest value of HADDOCK score of ACE2 docking with Spike protein of Omicron and nine species for Delta variants (Figure 5), indicating that binding affinity of spike protein with ACE2 of host elevated with the accumulation of mutations in SARS-CoV-2 variants. For the Delta strain, the species with the strongest binding affinity to the Spike protein was Macaca mulatta. For the Omicron strain, the strongest binding affinity is Panthera tigris altaica.

The structural similarity of M2(ACE2)-RBD(SARS-CoV-2) between humans and other species were compared by the TMalign using metrics TM score and RMSD, and the results were shown in Supplementary Figures 1, 2, respectively. Surprisingly, compared with M2(ACE2)-RBD(SARS-CoV-2) structural similarity in the wild-type Wuhan strain, the similarity decreased in variants for most species of 31 animals (Supplementary Figures 1, 2). However, the most striking finding to emerge from our study was that structural similarity of M2(ACE2)-RBD(ARS-CoV-2) between Homo sapients and Octodon degus and was elevated for Delta virus compared with that of Wuhan strain. Furthermore, the structural similarity of M2(ACE2)-RBD(ARS-CoV-2) between Homo sapients and Bos taurus was elevated for the Omicron virus compared with that of the Wuhan strain.

To test the effect of mutations of ACE2 on the binding affinity between hACE2 and SARS-CoV-2, we chose 11 residues to mutate and modeled the protein structure of the mutated ACE2 proteins. These residues were selected based on the frequency of mutations in other species (compared to humans). Table 4 shows the result of a few mutational ACE2 proteins docking with spike protein. The HADDOCK score of the original human ACE2 protein was −100.6 ± 17.5, which is less than the HADDOCK score of these mutated ACE2 proteins. This result indicates that the mutant ACE2 protein has a stronger binding capacity to the spike protein. Especially, the mutation Q340R has the lowest HADDOCK score of -139.2, which is the most stable in binding. It is known that residue 340 is not a binding site, suggesting that mutations in the non-binding sites can also affect the binding of ACE2 protein and spike protein. In the virus infection process from wild hosts to humans, SARS-CoV-2 may gradually adapt to the human ACE2 protein.

The amino acids (AA) sequence of the M2 region of hACE2 (1R42) was downloaded from the PDB database (PDBid: 1R42). The amino acid sequence was aligned with the NCBI nr database by blastp (Shiryev et al., 2007; Camacho et al., 2009). BLAST parameters were chosen as follows: e-value was set to 10; the cost of gap-open and gap-extend was set to 1; the maximum number of alignment results was set to 1E8. BLAST results with an e-value more than 0.01 and a similarity of less than 70% were filtered out.

Thirty-one species, including humans and 30 other species, were selected for further analysis according to four rules. Firstly, the similarity between the M2 region of human ACE2 protein and the given species was calculated; a species more similar than Rhinolophus sinicus was selected. We selected Rhinolophus sinicus as a borderline species since its infected coronavirus has 96% similhasy to SARS-CoV-2 at the whole-genome level, and it is not quite similar to humans. Secondly, six animals with phylogenetically identical to the human ACE2 gene were selected, including Chlorocebus sabaeus, Gorilla gorilla gorilla, Macaca mulatta, Pan paniscus, Pan troglodytes, and Pongo abelii. Thirdly, four animals, including Felis catus, Mustela putorius furo, Panthera tigris altaica, and Manis javanica, were selected since they are reported to be infected with SARS-CoV-2 or coronavirus with high similarity to SARS-CoV-2 (Shi et al., 2020). The selected 31 animal species cover wild animals, pet, and husbandry animals.

The ACE2 protein sequences of the 31 selected species were aligned by ClusterX (Thompson et al., 2002). Phylogenetic analysis was conducted using the maximum likelihood method based on the JTT matrix-based model (Jones et al., 1992). One thousand bootstrap replicates were used to calculate node support by MEGA7 (Kumar et al., 2016). The phylogenetic tree was displayed and annotated by the online tool iTOL (Letunic and Bork, 2019).

We constructed the phylogenetic tree for S proteins of 7 SARS-CoV-2 variants, including Wuhan, alpha, beta, gamma, kappa, delta, and Omicron, using the same approach as the ACE2 phylogenetic tree. Except for Omicron, the Spike protein sequences were retrieved from protein data bank (PDB) for six SARS-CoV-2 variants. The genome sequence for Omicron was downloaded from the GISAID database (Shu and McCauley, 2017), and the S protein sequence was annotated by local BLAST².

The ACE2 protein structures of different species were modeled by the I-TASSR analysis, which consists of three steps, including (1) threading, (2) structural assembly, and (3) model selection and refinement (Roy et al., 2010; Yang et al., 2015). In the first step, LOMETS threaded the query protein sequence using a library of non-redundant structures to identify structural templates. In the second step, the topology of a full-length model was constructed by reassembling the continuously aligned fragments excised from templates. The structure of unaligned regions was built from scratch by ab initio folding. In the third step, the simulation of fragment assembly was utilized again as the starting selected cluster centroids. The final structural model was generated by building all-atom models from Cα traces using the optimized hydrogen-bonding networks. The similarity of human ACE2 protein structure and 30 species’ was calculated by TMscore program implemented in I-TASSER software using metrics TM score, GDT-TS score, and RMSD (Yang et al., 2015).

The docking score of the HADDOCK server (van Zundert et al., 2016) was used to predict the binding affinity between SARS-CoV-2 spike protein and the ACE2 protein of a particular species. A few hACE2 residues, including Q24, D30, H34, Y41, Q42, K353, and R357, were set as the activated sites of the ACE2 receptor in HADDOCK since these residues were reported to effectively bind to the viral spike protein (Yan et al., 2020). In addition, the activated sites of spike protein included K417, Y453, Q474, Q498, T500, and N501, respectively.

Molecular docking between 31 ACE2 proteins and spike protein of five SARS-CoV-2 were conducted, respectively. The structural similarity of M2(ACE2)-RBD(SARS-CoV-2) between humans and other species were compared by the TMalign program implemented in I-TASSER software using metrics TM score and RMSD (Yang et al., 2015).

To predict the potential effect of mutants in hACE2, we mutate a few residues in hACE2, which have single nucleotide polymorphism (SNP) in 50% of species other than humans. Then, the Python library MODELER was used to model these mutated hACE2 proteins. Again, HADDOCK (van Zundert et al., 2016) was used to predict the binding affinity between the mutated/un-mutated hACE2 proteins.