Mutations of the SARS-CoV-2 Spike Glycoprotein Detected in Cats and Their Effect on Its Structure and Function

The high frequency of Severe Acute Respiratory Syndrome Coronavirus-2 (SARS-CoV-2) mutations and homology of the Angiotensin-Converting Enzyme-2 (ACE2) cell receptors in various hosts help the virus transcend species barriers. In this study, we investigated the mutations of the SARS-CoV-2 spike glycoprotein detected in cats and their effect on its structure and function. Interestingly, some of these mutations are reported here in cats for the first time. Structural analysis showed seven residue substitutions in the spike glycoprotein. Four of the detected mutations are located on the spike surface, which are critical interaction points for neutralizing antibodies. Furthermore, three of the reported mutations could facilitate viral binding to the ACE2 host receptor, influence S1/S2 cleavage, destabilize the β-hairpin structure of the S2 and enhance viral infectivity. Structural modeling and phylogenic analysis of the ACE2 receptor provided an indication of the binding capacity of the virus to the specific cell receptors of different species and hosts. The presented work highlights the effects of the residue substitutions on viral evasion, infectivity and possibility of SARS-CoV-2 spillover between humans and cats. In addition, the work paves the way for in-depth molecular investigation into the relationship between SARS-CoV-2 receptor binding and host susceptibility.

The SARS CoV-2 spike glycoprotein initiates a viral infection through its attachment to angiotensin-converting enzyme 2 (ACE2) cell receptors, which are distributed widely in different organs, including oro/nasopharyngeal mucosae, the lungs, gastrointestinal tract (GIT) mucosae, the brain, kidneys and the liver (Hamming et al., 2004). The spike glycoprotein is a trimeric fusion molecule (as illustrated in Figure 1), and its corresponding gene consists of 3,822 nucleotides (21,384) that encode 1,273 amino acids. The spike glycoprotein is composed of two main subunits, an amino (N)-terminal S1 subunit and a carboxyl (C)-terminal S2 subunit, that are cleaved at the furin cleavage site (S1/S2 cleavage region) . The detailed molecular structure of the spike gene with detailed regions and locations is illustrated in Figure 1 and tabulated in Supplementary Table S1, according to the Wuhan strain (Reference Sequence NC_045512.2, GeneID: 43740568) (Guruprasad, 2021;NCBI, 2021). The S1 subunit interacts with the ACE2 receptor subsequently promoting the viral infection (Hoffmann et al., 2020), and this interaction has an impact on both cross-species and human-to-human transmissions (Lu et al., 2015;Wan et al., 2020). The S1 subunit switches between a standing up receptor-binding domain (RBD) (e.g., one or two RBD-up conformers), amenable to ACE2 binding, and closed-down positions, undergoing spontaneous conformational transitions between ensembles of closed and open receptor accessible forms Wang et al., 2020). In contrast with the buried S2 subunit, the S1 subunit domains are located on the surface of the spike glycoprotein, effectively protecting the fusion apparatus . Upon proteolytic activation at the S1/S2 region, the S1 subunit is dissociated from the S2 . Following the latter dissociation event, the S2 subunit A B FIGURE 1 | The structure of the spike glycoprotein and its coding gene. (A) shows a schematic surface view of the spike ectodomain and its constituting domains and regions. (B) illustrates the structure of the native spike glycoprotein in its pre-fusion and post-fusion states as determined by Cryo-EM (PDB ID 6xr8 and 6xra, respectively) (Cai et al., 2020) in which the spike trimer is composed of three interacting chains or protomers depicted here in distinct colors. (C) demonstrates the detailed structure of the spike gene with the locations of the different regions shown according to the Wuhan strain (Reference Sequence NC_045512.2, Gene ID: 43740568) (Guruprasad, 2021;NCBI, 2021). The spike glycoprotein is composed of 2 subunits (S1 and S2). The S1 subunit includes the N-terminal domain (NTD), the receptor binding domain (RBD) which carries the receptor-binding motif (RBM), and two structurally conserved subdomains (the C-terminal domains 1 (CTD1) and the C-terminal domains 2 (CTD2). The S2 subunit contains the N-terminal hydrophobic fusion peptide (FP), the fusion peptide proximal region (FPPR), the heptad repeat 1 (HR1) motif, the central helix region (CH), the b-hairpin region, the connector domain (CD)and the heptad repeat 2 (HR2) motif. These are followed by the transmembrane region (TM) and the intracellular regions (IC) (Guruprasad, 2021;NCBI, 2021). The figure was assembled with Biorender website (Biorender, 2021).
undergoes a series of enormous structural rearrangements in order to orchestrate the fusion of cellular and viral membranes (Duan et al., 2020;Shang et al., 2020;Wan et al., 2020;Wang et al., 2020). SARS-CoV-2 has evolved at a rate of approximately 1.5-3.3 × 10 −3 /per site/year (Chaw et al., 2020). Consequently, residue substitutions may affect viral infectivity, pathogenesis, transmission and viral virulence (Shafique et al., 2020). The high mutation rate and recombination frequency is a significant dilemma that may have a role in cross-species transmission (Elaswad et al., 2020). SARS-CoV-2, which has a bat precursor, evokes many questions regarding its natural origin, its transmissibility to different animal hosts and the possible incidence of animal viral reservoirs (Elaswad et al., 2020). Many experiments have investigated the susceptibility of different species to SARS-CoV-2 infection, including nonhuman primates, ferrets, cats, dogs, minks, tigers, mice and pigs (Hobbs and Reid, 2021). Furthermore, a case of zoonotic transmission from minks to farm workmen in Denmark was investigated (Prince et al., 2021).
The cross-species transmission of coronaviruses and their ability to jump species barriers have been reported. Moreover, many research papers have investigated the factors that may strengthen the capability of coronaviruses to overcome species barriers and increase their virulence (Elaswad et al., 2020). These features are illustrated not only by the high mutation rates of coronaviruses (Su et al., 2016), but also by the large RNA genome of these viruses that allows the occurrence of recombination and mutation, resulting in the emergence of new coronaviruses . Despite the importance of these features, the interaction ability of coronaviruses with different ACE2 receptors has become the dominant and most important virulence factor of these viruses (Bolles et al., 2011). According to the fifth report of the Office International des Épizooties (OIE) regarding the SARS-CoV-2 situation up until the end of September 2021, 584 outbreaks in animals (mostly in minks, cats and dogs) were reported globally, affecting 12 species in 30 countries (as demonstrated in Figure 2) (OIE, 2021). Therefore, the OIE strongly encourages the reporting of SARS-CoV-2 occurrence in animals (OIE, 2021).
Furthermore, in our recent study, we reported SARS-CoV-2 spillover between humans and pet animals in Egypt in 30.3% of cats (10/33) and 24.2% of dogs (8/33), as evidenced by real-time reverse transcriptase polymerase chain reaction (rRT-PCR) (Mervat et al., 2021). Moreover, thrombocytopenia, lymphocytopenia and high levels of ferritin, C-reactive protein and D-dimers were reported in hematology and clinical studies (Mervat et al., 2021).
In order to further investigate such interspecies transmission, in this study, we conducted a full spike gene sequencing of SARS-CoV-2 detected in four cats, and their treating veterinarian. Samples selection was based on the severity of clinical signs. Sequencing was combined with computational analysis, including structural visualization and homology modeling of the glycoprotein to investigate the mutations and their effect on its structure and function. Structural modelling of the ACE2 receptor in different species was conducted to study the host susceptibility to SARS-CoV-2 infection.

Ethical Approval
The study was conducted according to ethical guidelines and approved by the Ethics Committee of Animal Health Research Institute (protocol code 11429 and date of approval 12/2020).

Sample Collection and Preparation
Swabs from oro/nasopharyngeal mucosa were collected from four cats, (three of these cats had severe respiratory signs and one of them were a symptomatic), that were in contact with COVID-19 positive owners. Additionally, a sample was taken from a veterinarian who had been in contact with these cats in his clinic and he was asymptomatic. These samples were collected in December 2020, March 2021 and July 2021. The swabs were deposited into sterile tubes containing 2 ml of autoclaved phosphate buffered saline (PBS). The latter samples were collected independently from different houses and represent independent outbreaks. The possibility of SARS-CoV-2 spillover was based on infection's chronology, symptoms onset and case history, as owners were infected with SARS-CoV-2 before symptoms appear on their animals, and hence they likely have transmitted the virus to their companion animals.

Nucleic Acid Extraction and Molecular Diagnosis
According to the manufacturer's instructions, viral ribonucleic acid (RNA) was extracted automatically using a MagPurix extractor (Labgene scientific, Switzerland). This was followed by a SARS-CoV-2 screening of the extracted RNA with rRT-PCR using a transgenbiotech multiplex SARS-CoV-2 diagnostic kit (Cat no. DV101) that detects two genes (nucleocapsid and ORF1ab genes). The cycling protocol of the rRT-PCR was composed of 10 minutes at 50°C and 30 seconds at 95°C, followed by 40 cycles of 5 seconds at 95°C then 30 seconds at 60°C.

Full-Length Spike Gene Amplification
The full-length spike gene of the spike glycoprotein was amplified for the rRT-PCR positive samples using conventional RT-PCR utilizing spike gene-specific primers, as reported previously (Ren et al., 2020) and tabulated in (Supplementary  Table S2). Briefly, the RT-PCR reaction was set up using an RT-PCR kit (PrimeScript, Takara, California, USA). The cycling conditions were as follows: a cycle of 30 minutes at 50°C for the activation of reverse transcriptase (RT); then 2 minutes at 95°C for the deactivation of the RT; followed by 40 cycles of 95°C for 30 seconds, 50°C for 30 seconds and 72°C for 1.5 minutes; with a final extension cycle of 72°C for 10 minutes.

Sequencing and Phylogenetics
According to the manufacturer's recommendations, positive DNA bands with the expected sizes were cut and purified using a gel extraction kit (QIAquick, Qiagen). The BigDye Terminator 3.1 Cycle Sequencing Kit was used to sequence the purified products (Applied Biosystems, USA). Next, the Centri-Sep Purification Kit (Applied Biosystem) was used to purify the sequencing reactions. The purified sequence products were then decoded using the 3500 Genetic Analyzer (Applied Biosystems). The BioEdit Sequence Alignment Editor 7.2.5 (Hall, 1999) was used to prepare protein multiple sequence alignments. MEGA7 software was used to create a midpoint rooted maximum likelihood phylogenetic tree with 1000 bootstrap replicates to verify the tree (Kumar et al., 2016). This tree depends on the Jones-Taylor-Thornton (JTT) matrix-based model (Jones et al., 1992) and the tree is illustrated with the highest likelihood log value (-11894.02) in which related taxa and closely related branches are clustered together. The initial tree for the heuristic search was automatically generated with application of the neighbor join and BioNJ algorithms to a matrix of pairwise distances predicted using the JTT model, then the topology with the highest log likelihood value was chosen. The branch lengths were calculated by the number of residue substitutions per site, and the tree is depicted to scale. MEGA7 software was used to display the evolutionary analysis (Kumar et al., 2016).

Structural Visualization and Obtaining a Homology Model of the Spike Mutant
Cartoon illustrations were generated using UCSF Chimera molecular visualization software (Pettersen et al.,2004). Atomic coordinates from PDB (ID 6xr8 and 6xra (Cai et al., 2020) were used to represent the native structures for the pre-fusion and post-fusion (S2) spike glycoprotein ectodomain trimer, respectively; whereas atomic coordinates from PDB ID 7krq and 7krr  were used to represent the (D614G) closed and one RBD-up open conformers, respectively. The homology model of the spike glycoprotein chains (i.e., trimer) for the sequence presented in this study (hCoV-19/Cat-Kitten/ Egypt/AHRI/July/2021) was generated using the MODELLER software (Eswar et al., 2003), with the (D614G) model (PDB ID 7krr, which is in the open state with one RBD-up) as a template. The sequences were aligned using the MUSCLE program (Edgar, 2004), and the columns corresponding to missing residues in the template model were removed before submission to the MODELLER program.

Homology Modeling of the ACE2 Receptors in Various Species and Their Phylogenetic Analysis
To be able to compare the receptors from various animal species, a BLASTP (Altschul et al., 1990) search was performed using the reference sequence for human ACE2 receptor against the NCBI Protein Reference Sequence database (Refseq_protein). Over 1200 significant hits were found. Sequences of the best 249 hits (i.e., those with the highest sequence identity) were retrieved and aligned using the COBALT multiple sequence alignment tool (Papadopoulos and Agarwala, 2007) through the NCBI server. Throughout this paper, all residue numbers for the receptor refer to the human ACE2 isoform 1 protein reference sequence (accession: NP_001358344.1). The alignment blocks of interest were extracted using the extractalign tool in the EMBOSS package (Rice et al., 2000). The alignment blocks for the main residues (or regions) that were in contact with the viral RBD in the structure of the RBD-ACE2 complex were used to generate sequence logos using the WebLogo Server (Schneider and Stephens, 1990;Crooks et al., 2004), namely alignment blocks incorporating residues 30-41, residues 82-84 and residues 353-357 were used. The residues coding for the signal peptide (residues 1-17) were removed from the alignment. The alignment blocks of residues 19-90 and of 324-393 were joint, and were then used for the construction of the ACE2 phylogenetic tree.
To obtain further structural insights into the differential binding capacity of the various host species receptors with the unmutated viral RBD, we built homology models of the receptor-RBD complexes of five distinct species. Specifically, models from cat, dog, ermine, ferret and chicken ACE2 receptor sequences were built. Homology models were built (also using MODELLER software) for the peptidase domains (PD domain, residues 18-614) of the ACE2 sequences of each species with a fixed COVID-19 RBD domain. For simplicity, glycosylation sites, which are far from the receptorbinding interface, as well as bound water molecules in the template structures were ignored. However, the Zinc atoms were included in the built homology models as they may be essential for stabilizing the structures. The receptor-RBD complexes (PDB IDs 6LZG and 6M17) were used as templates. A model for the human ACE2 isoform 1 protein reference sequence was also built and used as a control. For each species-specific receptor-RBD complex, the model with the best score using MODELLER was used in the subsequent steps. Hydrogen atoms were added to each model assuming a neutral pH = 7.0, and the structure was then minimized using all atom Amber14 force fields ("amber14-all") (Maier et al., 2015) in vacuum to remove any clashes. The energy minimization was performed using OpenMM Python API (Eastman et al). Isosurface representations and electrostatic (Coulombic) potential coloring of the isosurfaces were performed using Chimera (Pettersen et al., 2004).
Next, a phylogenetic tree was constructed for all residues at the binding interface between the ACE2 and the viral RBD (19-90 and 324-393) in different species. MEGA7 software was used to create a midpoint rooted maximum likelihood tree with 1000 bootstrap replicates to verify this tree (Kumar et al., 2016). A maximum likelihood technique depending on the JTT matrixbased model was used (Jones et al., 1992) with the highest likelihood log value (-1606.65). The related taxa and closed branches were clustered together. The initial tree for the heuristic search was automatically generated with the application of the neighbor join and BioNJ algorithms to a matrix of pairwise distances predicted using the JTT model, then the topology with the highest log likelihood value was chosen. The branch lengths were calculated by the number of residue substitutions per site, and the tree was depicted to scale. MEGA7 software was used to display the evolutionary analysis (Kumar et al., 2016).

Real-Time RT-PCR (rRT-PCR)
Based on rRT-PCR molecular diagnosis, the four cat samples and the veterinarian sample were SARS-CoV-2 positive (as illustrated in Supplementary Table S3) with cycle threshold (CT) values of nucleocapsid (N) and ORF1ab genes.

The Overall Distribution of Mutations and Impact on the Structure and Stability of the Spike Protein
The overall distribution of the mutations mapped onto the predicted structure (homology model) of the studied sequence ( Figure 4) showed that four residue substitutions were located on the viral surface ( Figure 4A), which were (W152R), (L452R), (A570V) and (A899S). Furthermore, the (D614G) mutation, which is found downstream of the receptor-binding domain (RBD), leads to a loop misplacement that results in structural rearrangements .
More specifically, in the closed trimer state, the D614 structure was found to have disordered loops between residues 620 and 640 (designated as 630 loops) within its monomers. Upon mutation, the G614 showed ordered loops instead, as a consequence of the elimination of a salt bridge between D614 of CTD2 and K854 in the FPPR. It was concluded that such a mutation leads to the stabilization of the closed trimer state. On the contrary, the one RBD-up state requires at least one of the loops to be disordered and the second one to be partially disordered. The latter state precedes the formation of two and three RNA-up states, which are the precursors of the S2 postcleavage state. The (D614G) mutation, therefore, slows down the formation of one RBD-up state, as well as the shedding of the S1 before the post-cleavage state formation. The accumulation of FIGURE 3 | Phylogentic analysis of the SARS-CoV-2 sequences. The rooted maximum likelihood tree shows the relationship between sequenced samples and current circulating variants in humans and animals in which the related taxa and closed branches are clustered together illustrating SARS-CoV-2 evolution. The blue circles illustrate the sequenced SARS-CoV-2 that were detected in cats in the present study while the red circles display the detected SARS-CoV-2 in cats worldwide from 2020 till January 2022. The green triangles depict circulating SARS-CoV-2 strains in Egypt within the same time frame. The red square refers to SARS-CoV-2 Wuhan-Hu-1 human reference strain.
the one RBD-up state, but not the other RBD-up states, is thought to explain the increased infectivity, but at the same time, the reduced overall binding affinity with the receptor. As a consequence, the formation propensity of a post-cleavage state is decreased in viral strains with a (D614G) mutation .
The (A570V) and the (S1051Y) mutations are of particular interest as they are located at the interface of different chains and domains. Supplementary Figure S1 illustrates that the residue 570 is surrounded by several other hydrophobic residues in both native and predicted structures. The increase in hydrophobicity of V570 compared with A570 suggests that this mutation is likely to contribute to the stabilization of this residue in place.
The impact of (S1051Y), which is located in the b-hairpin region of subunit 2 (S2), is illustrated in Figures 4C-E. The pre-cleavage structures of the native sequence, in both the closed and open conformations, show that the residue S1051 is in contact with two histidine residues, namely H1048 and H1064, as illustrated in Figures 4C-E. In the closed pre-cleavage structure, the distances between the hydroxyl oxygen of the serine S1051 to the nearest nitrogen atoms on the rings of H1048 and H1064 are 2.71 and 2.77 Å, respectively. One of the two histidines, H1048, is on the same side of the b-sheet strand as S1051, whereas the other (H1064) is on the other strand of the b-sheet hairpin. The serine residue is therefore bridging the interaction between the two histidines, effectively stabilizing the hairpin and the overall structure of the b-sheet. In contrast, the post-cleavage structure shows that one of the two histidines (H1064) moves away from the serine residue, whereas 1048 remains in close proximity. The distances between the serine oxygen and the closest ring nitrogen atoms of the histidines are 2.82 and 4.4 Å for H1048 and H1064, respectively. This is due to a shift that is seen of the relative orientation of the two strands of the bhairpin upon the unfolding of the helices in the post-cleavage conformation. Such large conformational changes increase the mechanical stress on the b-hairpin structure as well as on the nearby disulphide bridge. The predicted model of the (S1051Y) mutation shows that the histidine interactions are diminished in the pre-cleaved conformer, which may contribute to the destabilization of the hairpin and the overall stability of the b-sheet. On the other hand, these hydrogen bonds break may contribute to lowering the energetic barrier to the large conformational changes required for the post-cleavage structure to be formed.

Structural Modeling of the ACE2 Receptor in Various Species and Their Phylogenetic Analysis Deciphers Host Susceptibility
As the mutations that were found in the cat samples have been previously reported in humans, there must be other factors that determine the host range of the virus, including viral entry into the host cell. To further understand the factors that control host range and host susceptibility, we compared the sequences of the host receptor, ACE2 (angiotensin-converting enzyme 2), in different hosts. The ACE2 receptor is a highly conserved protein that performs a specific enzymatic function in mammals. We performed an initial alignment of the ACE2 protein isoform 1 from humans with that of susceptible companion animal hosts, (namely, cats and dogs) and with a species that have not been reported to be infected with SARS-CoV-2 (i.e., an unsusceptible control species). For the latter purpose, we used the ACE2 protein from chickens.
In this initial comparison, we focused on the residues that have been shown to have direct contact with the virus RBD domain (e.g., from PDB IDs 6LZG and 6M17). Table 2 shows a comparison between the key residues that form such polar and hydrophobic contacts. One can ascertain that the number of residues substituted in such a way that alters the charge distribution of the binding interface in dogs (3 out of the 10 residues are substituted) is more than that in cats (2 out of the 10 residues are substituted), and that such substitutions are the most in chickens (7 out of the 10 residues are substituted). We also noticed that there seems to be no clear pattern or preference in these substitutions, but rather an overall change in the charge distribution was observed. More specifically, a polar to hydrophobic residue substitution Q24L is seen in cats and dogs. Whereas in chickens, a polar to a charged residue substitution Q24E is seen. Residues D31 in humans or E31 in cats and dogs, are substituted with A31 in chickens. The polar and positively charged residues Q32 and K31, respectively, which are conserved in humans, cats and dogs are substituted with negatively charged residues (K31E and Q32E) in chickens. The hydrophobic methionine residue M82, is substituted with a polar residue in cats and dogs, and even with a charged residue in chicken. Finally, the tyrosine residue Y83, which is conserved in human cats and dogs, is substituted with a more hydrophobic residue (phenylalanine) in chickens. Additionally, the alignment showed several residues on the surface of the receptor in close proximity to those that form direct contacts with the RBD, which undergo substitutions that also alter the binding interface charge distribution in unsusceptible species compared with susceptible ones. It is possible that these neighboring residues may affect the overall binding affinity of the viral RBD with the host receptor.
For instance, some substitutions may introduce an overall A B D E C FIGURE 4 | The overall distribution of the spike protein mutations observed in the studied sequences and the impact of (S1051Y) mutation on its structure and function. (A) shows the overall distribution of the mutations, found in the sequenced samples in this study, mapped onto the predicted structure by homology modelling. The mutated residues in one monomer are depicted in red and the monomer with RBD-up is depicted in dark blue. The isosurface of the trimer is shown on the left and a rotated view (by 180•) of the latter is shown in the center, while the atomic details are shown on the right. (B) illustrates the predicted structure of the sequence presented in this study, depicted in light blue, superposed onto the template used to generate it (depicted in gold, derived from PDB ID 7krr). (C-E) illustrate molecular details of the impact of (S1051Y) on the structural stability and function of spike protein.
In the presence of the native serine 1051 (PDB ID 7krr), the residue interacts with two histidine residues (His 1048 and His 1064), as illustrated in (C) (top view, showing hydrogen bonds depicted by the cyan lines) . Such interactions are likely to stabilize the beta-sheet within the native conformer. These interactions are likely to be abolished by the mutation (S1051Y), as shown in the predicted model by homology modelling (E). (D) shows the interaction of the serine and the two histidine residues in the unmutated postcleavage structure [PDB ID 6xra (Cai et al., 2020)]. The figure illustrates that post-cleavage structural changes in the protein lead to bond break between the serine residue and His 1064.
negative charge (e.g., the substitution of E35 in humans, cats and dogs, to R35 in chickens), increase hydrophobicity or even introduce stearic barriers. These observations suggest that the overall residue composition and associated charge distribution on the binding interface may be taken as an indicator (or at least a contributing factor) of the susceptibility of the host. The latter also suggests that the homology of the 10 sites on the ACE2 that make direct contact with the RBD, between humans and other hosts, may not alone be an efficient indicator of the host susceptibility, and a better indicator would be to compare all residues that are at or near the binding interface.
To study this on a wider range of animal species, we performed a BLASTP search against the protein reference sequence database and retrieved sequences from different animal species with the highest sequence identity compared with the human ACE2 peptidase domain sequence (see methods for details). The multiple sequence alignment of the latter containing 249 sequences was used for further analysis. Indeed, sequence logo representations ( Figure 5) of the alignment show that certain residues of those that form direct contact with the viral RBD, along with neighboring residues, are conserved among many animal species. More specifically, residues F32, N33, A36, E37, D355, F356, K353 and R357 (marked with black and blue arrows in Figure 5) are highly conserved in the analysed sequences of the different animal species and are present at the receptor-RBD binding interface. Residues that are in direct contact with the viral RBD (mostly through their side chains) are marked with blue or red arrows. K353 and R357 are highly conserved and are in direct contact FIGURE 5 | Conservation of residues involved in direct contact between the receptor and the viral RBD. The figure shows sequence logo representation of alignments of residues 30-41, residues 82-84 and residues 353-357 of the host receptor (ACE2) from different species. The alignments involved 249 sequences that had the highest identity with the human ACE2 isoform 1 sequence. Highly conserved residues in this region are marked with black arrows. Conserved and nonconserved residues that are involved in direct contact with the viral RBD are marked with blue and red arrows, respectively.

Human
Cat Dog Chicken The table shows key residues involved in the binding of the human ACE2 with the viral RBD. The structure showed key residues that form a set of polar contacts (as well as one hydrophobic contact) with the viral RDB, hence stabilizing the complex. Substitutions that increase the hydrophobicity and polarity of the binding surface with respect to the human receptor are depicted in gray and orange fonts, respectively. Whereas substitutions that increase negative and positive charge are depicted in red and blue fonts, respectively.
with the viral RBD. Whereas, residues D30, K31, H34, Y41, M82 and Y83 (marked with red arrows), along with two other residues (Q24 and Q42), are also directly involved in the binding, but show moderate to high variations among the different animal species. Hence the extent of substitutions in the latter (variable) residues can provide a preliminary indication of how similar is the binding capacity of the different animal receptors compared to human. Neighbouring highly variable residues at the binding interface may, however, such as E35, D38 and G354 although are not in direct contact with the viral RBD in the human receptor, may have substitutions that would enhance or interfere with binding in other species. These findings confirm the previous notion that considered that all residues at or near the binding interface would form a better indicator for the likelihood of viral RBD binding with a specific host. Based on the latter arguments, and for the purpose of finding a simple but quantitative indicator of virus entry and susceptibility based on comparative sequence analysis with the human sequence, we constructed a phylogenetic tree (Figure 6) using the alignment blocks of all the residues at the binding interface between the ACE2 and the viral RBD. Residues 19-90 and 324-393 (shown in orange color in Figure 6A) were selected for this purpose. The phylogenetic analysis and identity matrix of all the aligned residues at the binding interface between the ACE2 and the viral RBD revealed that this region in humans, gorillas and chimpanzees are closely related to each other and are identical (100%). The latter region identity between humans and cats (83.8%) is higher than that of humans and dogs (79.5%). While the identity is 78.1% in ermines (minks) and 77.4% in ferrets, it is quite low in chickens (64.7%), as presented in Figure 6B and Supplementary Table S5.
To visualize and investigate the impact of mutations at the binding interface, homology models were built for the receptor-RBD complex, by fixing the sequence of the viral RBD and varying the receptor sequence according to the species (Figure 7). Models were built using the PDB IDs 6LZG  and 6M17 (Yan et al., 2020) as templates. Figure 7A shows the isosurface representations of the obtained models (Supplementary Figure S2) for the receptor PD of the selected species, namely human, cat, dog, ferret, ermine and chicken. The binding interface is colored in different colors, whereas the remaining part of the receptor domain is colored in dim gray. An analysis of these models revealed that, on the one hand, the overall shape and structure of the binding interface region of the cat and the dog are more similar to that of the human when compared with their counter regions in the ferret and the A B FIGURE 7 | Isosurface representations and electrostatic (Coulombic) maps of the predicted ACE2 receptor structures for various animal hosts. The figure shows the isosurface representations of the predicted structures of the ACE2 (A) and the corresponding Coulombic coloring of the surface representations (B). The region that is likely to affect the binding (residues 19-90 and 324-393) is depicted in distinct colors whereas the rest of the PD domain is depicted in dim grey. The viral receptor is hidden to show the overall shape and surface charge of the binding surface on the ACE2 receptor. For clarity, the human receptor-RBD complex is shown on the left of both atomic and its isosurface representations (panels A and B, respectively), in the same orientation as of the other models. In the latter, the RBD is in blue, the region having residues 19-90 and 324-393 is depicted in orange and the remainder of the PD domain is in dim grey. In panel B, the positively charged, negatively charged, and hydrophobic (neutral) regions are depicted in blue, red and white, respectively. chicken. These larger variations in the latter (less susceptible) species may introduce stearic effects that may reduce or hinder the binding. On the other hand, the electrostatic (Coulombic) colored isosurfaces ( Figure 7B) show the appearance of several hydrophobic regions within the binding region (shown in white), as well as an increase in the positive charges (shown in blue), in the cat and the dog compared with that of the human receptor. The latter positively charged regions appear to be more in the dog than in the cat. In contrast, the surface of the ferret and the chicken receptors appear to be more negatively charged. Similar conclusions can be drawn on the analysis of polar contacts and potential hydrogen bonding in the same models, i.e., between the RBD and receptors from the different species (Supplementary Figure S2). It is clear that the number of polar contacts in cats and dogs are less than in humans. Out of the studied species, the chicken has the least contacts. The ermine and the ferret had a different distribution of predicted hydrogen bonds when compared with humans.

DISCUSSION
Coronaviruses infect a wide variety of hosts (Shafique et al., 2020). In the present study, we aimed to investigate the mutations of the SARS-CoV-2 spike glycoprotein detected in cats and their effect on its structure and function by sequencing the full spike gene, followed by bioinformatics analysis. Furthermore, we highlighted the effect of residues substitution on viral infectivity, evolution and evasion from immune response. We found that cats with a history of direct contact with COVID-19-positive owners were also infected with COVID-19, revealing the possibility of SARS-CoV-2 transmission between humans and cats.
Generally, the high substitutions and mutations in coronaviruses revealed that the virus is adapting to new hosts (Duffy, 2018). The significance of amino acid substitution depends on the switching of properties and classes; for instance, a spike glycoprotein switching from leucine to methionine in L452M may have no significance as both are hydrophobic and neutral. Similarly, switching from tyrosine to phenylalanine in the Y453F mutation may have no significance as both are aromatic and hydrophobic. In contrast, substitution from glutamine to histidine in the Q498H mutation could significantly change the protein and side chain structure because they have different features, as histidine has a positive charge while glutamine has a neutral charge (Elaswad et al., 2020). Therefore, significant changes in a spike glycoprotein could change its viral infectivity, pathogenicity and virulence, and could increase cross-species transmission .
In the present study, SARS-CoV-2 incidence in cats was detected from four clinical samples, and was consistent with molecular results reported in previous studies (Newman et al., 2020). Furthermore, previous experimental studies have shown the clinical susceptibility of cats to SARS-CoV-2 (Hamer et al., 2021), which may be due to the similarity in ACE2 receptors between feline species and humans (Sailleau et al., 2020;Patterson et al., 2020;Hamer et al., 2021;Ruiz-Arrondo et al., 2021). Moreover, cats shed SARS-CoV-2 horizontally after infection for up to 6 days (Hamer et al., 2021), and the virus can persist in cats post-infection for up to 32 days (Bosco-Lauth et al., 2020;Shi et al., 2020). The high positive incidence in our study can be attributable to the fact that all tested cat samples came from homes with COVID-19-positive human cases.
SARS-CoV-2 transmission between humans and cats was reported in this study as the results of the phylogenetic analysis revealed a similarity between SARS-CoV-2 detected in cats (OK144251, OK144252 and OK144253 to OK144255) and in human cases reported in the same months (EPI_ISL1315064, EPI_ISL1936365 and EPI_ISL3274157), as presented in Figure 3 and Supplementary Table S4. Additionally, the veterinarian who was in contact with the infected cats was asymptomatic; however, he tested positive for COVID-19 and his sequence result was identical to that of the young cat and her mother. Moreover, sequence similarity between the SARS-CoV-2 spike gene detected in the young cat and her mother confirms that the virus is the same. Thus, there is a potential support of SARS-CoV-2 transmission between humans and cats, but further studies are needed to investigate the direction of the viral transmission as the study's design has not determined the infection direction whether animal-to-human or human-toanimal or animal-to-animal viral infection occurred.
The present study revealed changes in the spike glycoproteins and the unique amino acid substitutions of the cat samples in correspondence with the reference Wuhan isolate, as illustrated in Table 1, Figure 4 and Supplementary Figure S1. Seven critical amino acid substitutions were observed. The diverse and frequent single residue mutations that emerged independently and repeatedly suggests their possible adaptive role (Kubik et al., 2021). Different spike regions are subjected to distinct evolutionary changes as mutations at receptor-binding domains (RBDs) result in both changes in the binding affinity to ACE2 receptors as well as in viral evasion from neutralizing antibodies, whereas mutations in N-terminal domains (NTDs) are mostly related to the potential escape from neutralizing antibodies (Kubik et al., 2021). Interestingly, some of these mutations are unique to the Wuhan-Hu-1 reference isolate and are reported in cat samples for the first time in this study, as they were not previously detected in cat samples in Hong Kong, France, Italy or the USA. The latter mutations are (W152R), (A570V), (Q677H), (A899S) and (S1051Y). In this study, we focused on the spike gene as it the most variable part of the viral genome and it plays crucial role in the receptor binding. It is also the protein that carries the receptor binding domain (RBD) which is the main point of entry to the host (Fung and Liu, 2019). The spike is also the main protein that exhibit mutations due to immune escape and generally have a major immunological role (Harvey et al., 2021).
Structural visualization and homology models revealed that four residue substitutions were located on the viral surface, which may have appeared as viral antigenic evolution to evade neutralizing antibodies. These mutations, (W152R), (L452R), (A570V) and (A899S), have weakened critical interaction points with the neutralizing antibodies and are associated with low sensitivity to mAbs and vaccine sera (Kubik et al., 2021). The observed (W152R) mutation is critical and is frequently observed across geographical and phylogenetic contexts; thus, it is a hotspot for mutations (Kubik et al., 2021). The substitution of tryptophan to leucine (W152L) or arginine (W152R) is the most common mutation at this location based on the SOPHiA DDM database. In addition, the substitution of tryptophan to cysteine (W152C) or glycine (W152G) is rare (Sophiagenetics, 2021). During the second wave of the SARS-CoV-2 pandemic, between December 2020 and April 2021, the (W152R) mutation increased steadily and rapidly, with one of the most notable adaptive mutations (N501Y) in the receptor-binding domain (RBD) coexisting (Kubik et al., 2021). Similarly, (W152R) substitution affects viral evasion from neutralizing antibodies and escape from immunity, as demonstrated in structural studies in which tryptophan (W) removal has weakened critical interaction points with antibodies (Kubik et al., 2021). In particular, a residue in the antibody chain interacts with W152 in a hydrophobic pi-stacking interaction, wrapping W152 inside one of the complementarity-determining region (CDR) loops in the variable chains of antibodies (Kubik et al., 2021). Indeed, W152 is a key position that contributes to the latter interaction. As a result, changes in W152 impair not only the pi-stacking interaction, but also the pocket created by the antibody's CDR loop, resulting in a significant reduction in binding affinity (Kubik et al., 2021). Furthermore, (W152R) substitution suggests viral evolution to maintain effective binding to the ACE2 receptor (Kubik et al., 2021). Therefore, monitoring the individual mutations that occur outside the receptor-binding domain (RBD) and inside the N-terminal domain (NTD) is critical.
The second mutation observed in this study was (L452R) in the receptor-binding motif (RBM), as illustrated in Supplementary Figures S2B, C. This mutation was first reported in January 2021 in California, USA in the 20C/S:452R variant, and its frequency has increased over time (Zhang W et al., 2021). Since the first observation, novel variants of SARS-CoV-2, including this amino acid substitution, have been seen in the Delta/21A/B.1.617.2, Kappa/21B/B.1.617.1 and Epsilon/21C/ B.1.427.9 viruses with increased virulence and infectivity (Lauring and Hodcroft, 2021). Additionally, in this mutation, a factor of a 3-to 6-fold reduction in antibody neutralization titers in vaccinated individuals was observed (McCallum et al., 2021). Furthermore, this substitution reduced the neutralizing activity of monoclonal antibodies against RBD by 14-to 35-fold (McCallum et al., 2021). Moreover, a 2-fold increase in viral shedding was reported in the involved variants (Deng et al., 2021). It was also found that (L452R) could strengthen ACE2 receptor-binding affinity  and increase viral infectivity . A mutation in the C-terminal domain 1 (CTD1) at position 570 was also observed. In addition, (A570V) was detected in North America and Asia (Guruprasad, 2021). Recently, biophysical and Cryo-EM studies demonstrated that some residue substitutions in the spike have affected the conformational dynamic equilibrium between the open RBD-up and closed RBD-down forms of the spike (Henderson et al., 2020). In particular, the authors found that certain residue substitutions in the regions of contact between domains could affect the equilibrium between the open and closed states. For example, the combined mutations in the A570, T572, F855 and N856 positions of SD1 and S2 inter-protomer regions have led to the dominance of the open RBD-up receptoraccessible form of the spike. Therefore, these structural rearrangements have resulted in the enhancement of the binding affinity of the spike to the ACE2 host receptor (Henderson et al., 2020). More specifically, the substitution of A570 with a more hydrophobic one, such as L570, has aided in such an equilibrium shift. Hence, the hydrophobic (A570V) mutation (Supplementary Figure S1) observed in our samples may contribute to a similar equilibrium shift. The fourth mutation on the surface was (A899S) in the fusion-peptide proximal region (FPPR). This mutation was recently observed in Jordan (Alanagreh et al., 2021), Saudi Arabia (Alghoribi et al., 2021) and Italy (Castelli et al., 2021). This substitution is recorded in lineage C.36 and is associated with high transmissibility and low sensitivity to mAbs and vaccine sera (Castelli et al., 2021).
The fifth mutation in our study (D614G) was found downstream receptor-binding domain in the C-terminal domains 2 (CTD2). Such a mutation became dominant during the emergence of the COVID-19 pandemic (Korber et al., 2020a;Yurkovetskiy et al., 2020). It is not only able to enhance the viral/ host interaction on the cellular level (Ortega et al., 2020), but it is also involved in the upward increase in SARS-CoV-2 pathogenicity, increasing viral transmission and evasion from immune response and vaccine (Korber et al., 2020b). This mutation appears to enhance the spread of SARS-CoV-2 and increases the viral infectivity that resulted from increasing the stability of the spike trimer. Certainly, in such a mutation, there is loop misplacement, structural rearrangements and interaction between downstream receptor-binding domains, which act as wedges within the spike trimer that strengthen the stability and prevent pre-maturation of the spike, and improve the strong assembly of the virus, resulting in the enhancement of the viral binding with the ACE2 receptor, as presented in Figure 4. This could facilitate viral entry, infectivity and virulence. It is unsurprising that all newly emerged variants contain this mutation and that they are more infectious than the Wuhan reference virus . The (D614G) substitution was reported to improve viral replication in human lung epithelial cells and primary human airway tissues by improving viral infectivity and stability (Plante et al., 2021). Furthermore, this mutation lowered the effectiveness of COVID-19 vaccinations in clinical studies (Plante et al., 2021).
The (Q677H) mutation in the S1/S2 cleavage site was detected in both the two cat samples and the veterinarian sample in this study. The (Q677H) mutation was firstly identified in clade 20C in the USA (Pater et al., 2021). It was also detected in the (Eta/ 21D/B.1.525) variant . As the position of such a mutation is found in a critical polybasic region of the furin-binding site and it could influence S1/S2 cleavage, SARS-CoV-2 infectivity and pathogenicity could be subsequently affected . The seventh mutation observed in this study is (S1051Y) in the b-hairpin region of subunit 2 (S2). The three-dimensional organization of HR1 (residues 910 −985), CH regions (residues 986−1035), b-hairpin (1035-1068) and UH (residues 736−781) is a significant feature of the spike glycoprotein structure (Walls et al., 2016). Structural modeling and analysis suggest that the (S1051Y) mutation ( Figures 4C-E) structure. Indeed, this mutation could play a role in viral entry during the replication cycle. The latter hypothesis requires further experimental investigation.
Despite the multiple mutations found in the sequences presented in this study, all these mutations had been previously reported in humans. Hence, although some of the latter mutations may be important for the virus strain survival within cats, these mutations are not unique to this particular host, and it is likely that a cross-infection between humans and cats has occurred. These findings raise several questions, including: What factors control the SARS-CoV-2 host ranges? Do variations in the host receptor play a role in the susceptibility of the host to SARS-CoV-2 infection? Can we obtain insights into the susceptibility of a species by comparing the species receptor with the human receptor? To answer these questions, we analyzed the ACE2 receptors from various animal species using sequence and structural bioinformatics tools.
Obviously, many host factors, may play a crucial role in susceptibility to SARS-CoV-2 infection. One of the possible factors likely to influence the virus host range is the variation in the viral receptor itself from one species to another. In this work, we showed that simple sequence comparisons ( Table 2 and Figure 5), phylogenetic analysis ( Figure 6) and the structural modeling (Figure 7) of residues at the binding interface between ACE2 receptors from different hosts and the viral RBD may be used as an initial effective indicator of the susceptibility of the host to SARS-CoV-2 cross-infection. More specifically, the amino acid composition of this region (residues 19-90 and 324-393) may provide an indication of the binding capacity of the virus to the specific species host cell, and hence a determinant for the susceptibility of the host. This study, therefore, reveals that a simple phylogenetic analysis and identity matrix of the aligned residues could provide a quantitative measure of the likelihood of virus entry into the cells of various hosts.
The phylogenetic analysis and identity matrix of the aligned residues at the binding interface between the ACE2 and the viral RBD ( Figure 6B and Supplementary Table S5) revealed that SARS-CoV-2 has the ability to interact with a broad range of ACE2 cell receptors in different species, which could increase its virulence and could enhance its inter-species transmission (Su et al., 2016) (Bolles et al., 2011). The sequence identity of the latter region between humans and cats (83.8%) is higher than that between humans and dogs (79.5%), which suggests that the potential recognition of the RBD-spike glycoprotein in cat ACE receptors is stronger than that in dogs, and these results are in accordance with a previous study in China (Luan et al., 2020). This high sequence identity in cats compared to human could explain the relatively high SARS-CoV-2 infection rate in the feline species (Hamer et al., 2021). Hence, regular monitoring of companion animals for SARS-CoV-2 infection is highly recommended. The similarity percentage in ermines (minks) (78.1%) is lower than that in cats and dogs, however the reported infection rate of SARS-CoV-2 in minks is high, according to the fifth report of the OIE (2021). This may be an indicator that there are other effective factors in minks which make this species more susceptible to SARS-CoV-2 than others, and these factors need further investigation. The sequence identity was 74.6% in greater horseshoe bats and 78.8% in Chinese pangolins. Previous studies have suggested that SARS-CoV-2 from bat CoV RaTG13 and pangolin isolates, have sequence identities of 96.3% and 89%, respectively. These studies concluded that, bats could be the natural origin of SARS-CoV-2 (Bosco-Lauth et al., 2020;Shi et al., 2020), and although pangolins cannot transmit SARS-CoV-2 directly due to the percentage of divergence, they may play an intermediate role (Calvet et al., 2021). Indeed, there is a need for further investigation to discover other factors that may point to bats being the potential natural host for SARS-CoV-2, rather than the receptor homology and binding affinity. Overall, the presented work paves the way for binding affinity predictions and free energy calculations through molecular dynamics studies, and opens the gates for further investigation into species susceptibility to SARS-CoV-2 infection.

CONCLUSIONS
This study displays critical mutations in the spike glycoprotein of SARS-CoV-2 detected in cats and in their treating veterinarian, through full spike gene sequencing, structural visualization and homology modeling. The analysis of the detected mutations revealed potential evidence for spillover of SARS-CoV-2 between humans and cats. Furthermore, the reported mutations could affect the viral evasion from immune response, the viral infectivity and pathogenicity. Hence, regular complete genome sequencing of SARS-CoV-2 in humans and nearby companion animals, as well as studying the recombination events in such viruses, are recommended in order to track newly emerging mutations and investigate their effects on viral adaptation, transmission, and their impact on therapeutics and vaccine design.

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.

ETHICS STATEMENT
The study was conducted according to ethical guidelines and approved by the Ethics Committee of the Animal Health Research Institute (protocol code 11429 and date of approval 12/2020). The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the owners for the participation of their animals in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.