New putative animal reservoirs of SARS-CoV-2 in Italian fauna: A bioinformatic approach.

SARS-CoV-2 is a virus belonging to the betacoronavirus family, causing fatal respiratory disease in humans, which became pandemic in 2020. Italy is one of the most affected countries by COVID-19, particularly in the northern regions. Several studies consider COVID-19 a zoonotic disease and, since Italy is the repository of a high biodiversity, SARS-CoV-2 infection in animals can be considered as a reservoir of the virus or favor the spreading between animals and humans. In this work, we analyzed the amino acid sequences of ACE2 protein of the most common domestic and wild animals present in Italy. Among the latter, we focused on ACE2 of the Chiroptera species present in Italy to identify the primary reservoir in this region. First, we reproduced in silico the Chiroptera ACE2/viral spike (S) protein interactions on the human ACE2/SARS-CoV-2 S complex model and identified the critical residues for the binding. In silico molecular docking of ACE2 belonging to Chiroptera vs SARS-CoV-2 S protein pointed to Rhinolophus ferrumequinum as a bat living in Italy, that may be a potential primary reservoir of the virus. On the other hand, a sequence similarity search on ACE2 of domestic and wild animals living in Italy pointed to domestic (horses, cats, cattle and sheep) and wild (European rabbits and grizzly bears) animal species as potential SARS-CoV-2 secondary reservoirs. Molecular docking of ACE2 belonging to these species vs S protein of Bat coronavirus (Bt-CoV/Rp3/2004) suggests that the primary reservoir Rhinolophus ferrumequinum may infect the secondary reservoirs, domestic and worldwide animals living in Italy, determining a specific risk of SARS-CoV-2 infection.

Introduction

Severe respiratory syndrome coronavirus-2 (SARS-CoV-2) is a novel, highly transmissible betacoronavirus that caused pandemic Coronavirus Disease 2019 (COVID-19) [1]. Outbreaks occurred worldwide, including Italy, where the number of infections was exceptionally high in northern regions [2]. Among betacoronaviruses, SARS-CoV-2 is the third virus causing fatal severe human respiratory disease after SARS-CoV and Middle East respiratory syndrome coronavirus (MERS-CoV) [1].

Identification of animal species that can potentially be reservoirs for the virus is strategic to controlling SARS-CoV-2 infection and preventing future outbreaks caused by a mutated form of the virus [3, 4]. From the genetic proximity of SARS-CoV-2 to RaTG13 coronavirus, a bat origin for the current COVID-19 outbreak was hypothesized. Concurrent evidence also proposed pangolins as a potential intermediate reservoir [1, 5, 6, 7, 8]. Concerning the secondary reservoir, susceptibility to the SARS-CoV-2 infection was demonstrated for several domestic animals both in silico [1, 9] and experimentally [10, 11]. SARS-CoV-2 can replicate efficiently in young cats and ferrets. Both these animals can also transmit the infection to other cats via airborne transmission in an experimental setting. Conversely, dogs have shown low susceptibility to experimental infection, while pigs, chickens, and ducks are not susceptible [6]. Besides pets, several studies have demonstrated that numerous wild animals are potential host intermediate for SARS-CoV-2 infection [12].

Mammalian biodiversity is recognized as a risk factor for zoonotic disease emergence [13, 14, 15]. COVID-19, as a zoonotic disease, represents a high treat for regions with high biodiversity; it represents a high treat for regions like Italy, which is among the European countries richest in biodiversity [16, 17] (https://www.isprambiente.gov.it/en/archive/news-and-other-events/ispra-news/year-2015/may/biodiversity-in-italy). Strategies to contain infection are being adopted following the WHO recommendation [18] (https://www.who.int/publications/i/item/clinical-management-of-covid-19). However, given the high biodiversity, great attention is on identifying mammalian species occupying the Italian territory, which may act as potential new primary or secondary reservoir of SARS-CoV-2. Lastly, ecologically speaking, the control of a new virus generation may prevent endangered species extinction and avoid a severe threat to Italy's biodiversity.

Here we present a bioinformatic investigation aimed at identifying the host range of SARS-CoV-2 most diffuse in Italy. In particular, considering that the virus originates from bats [19], we wondered i) if Chiroptera species exist in Italy as potential primary reservoir for the virus, ii) if these bats could infect animals in the wild and domestic environments, and iii) if these animals could be a secondary reservoir in Italy (Figure 1).

Figure 1

Bioinformatic workflow followed for the identification of primary and secondary reservoir in Italy.

Analogously to many viruses [20], infection of SARS-CoV-2 occurs through the binding of viral surface glycoprotein, spike (S) protein, to the angiotensin-converting enzyme 2 (ACE2) [1, 21]. ACE2 is expressed in several cell types, including type II pneumocytes, myocardial cells, cholangiocytes, enterocytes, and oral mucosal epithelium [22]. According to a structural biology view, infection of SARS-CoV-2 depends on an optimal interaction of S protein with the ACE2 receptor. Therefore, species where ACE2 can efficaciously interact with SARS-CoV-2 S protein can be considered valuable reservoirs of the virus.

In this context, starting from the study of the recently solved human ACE2 (hACE2)/S protein complex [9] and considering the sequence similarity of ACE2 in Chiroptera and domestic and wild animal species living in Italy, we followed two bioinformatics approaches to identify new potential primary and secondary reservoirs of SARS-CoV-2 occupying Italian territory (Figure 1).

In silico molecular docking of ACE2 belonging to Chiroptera vs SARS-CoV-2 S protein pointed to Rhinolophus ferrumequinum as a bat living in Italy, that may be a potential primary reservoir of the virus. On the other hand, a sequence similarity search on ACE2 of domestic and wild animals living in Italy indicated 10 species as potential SARS-CoV-2 secondary reservoirs. Molecular docking of ACE2 belonging to these species vs S protein of Bat coronavirus (Bt-CoV/Rp3/2004) suggested that the primary reservoir Rhinolophus ferrumequinum may infect the secondary reservoirs, domestic and wild animals living in Italy, determining a specific risk of SARS-CoV-2 infection.

Methods

Construction of homology models

Homology models of residues 1–358 of Chiroptera ACE2 proteins and bat coronavirus spike protein were built using ExPASy SWISS-MODEL [23, 24]. The templates were chosen according to i) sequence similarity; ii) global model quality estimate (GMQE) value [24]. GMQE is a quality estimation method combining properties from the target–template alignment and the template search method. GMQE value is expressed as a number between 0 and 1: the higher the value, the higher the reliability; iii) qualitative model energy analysis (QMEAN) value [25, 26]. QMEAN is a composite estimator based on different geometrical properties that provide global and local absolute quality estimates based on one single model. The PDB models generated were used for the molecular docking calculation.

Molecular docking

Molecular docking was performed using HDOCK server (http://hdock.phys.hust.edu.cn/), which calculates protein-protein interaction through a hybrid algorithm of template-based and template-free docking. The crystal structure of S protein receptor binding domain (RBD) (PDB ID: 6M0J), corresponding to the residues 437–508 of S protein, was docked against the homology models of bat ACE2 [27]. The calculation was carried out imposing as constraints 8 Å distance between S protein A475 and ACE2 E23 and 5 Å distance between S protein N501 and ACE2 D355. The results were analyzed according to i) docking score value (kcal/mol), generated by HDOCK scoring function [28]; ii) ligand root-mean-square deviation (RMSD) of atomic positions value (Å). RMSD measures the distance between the docked pose and a model of the same ligand predicted by template-based homology modelling generated by HDOCK (the lower the value, the more similar the poses).

Molecular dynamics

The structures resulting from molecular docking simulation were subjected to 1 ns molecular dynamic (MD) simulations using GROMACS [29, 30], (Gromos96 53a6 force field) [31]. The structure was immersed in explicit water using the SPC model [32]. The protein was solvated; the system was neutralized by adding Na+ ions, energy minimized, and equilibrated using NVT and NPT runs. The temperature and pressure were kept constant at 300 K and 1.01325 bar using the Berendsen weak coupling-method [33]. The results were used for an MD simulation using Particle Mesh Ewald for long-range electrostatics under NPT conditions. Coordinates were saved every 100 ps. Trajectory files containing the structural coordinates of receptor-ligand complex sampled every 100 ps were fitted in the box and converted in PDB coordinates using the trjconv tool of GROMACS Package. The structures were visualized with Maestro by Schrödinger [34] (see Supplementary Figures S1 and S7).

Analysis of ACE2 sequences of human and Chiroptera species present in Italy

Amino acid sequences of Chiroptera ACE2 living in Italy were retrieved from UniProt database. The sequences of Rhinolophus ferrumequinum (UniProt ID: E2DHI2 and B6ZGN7) and Myotis daubentonii (UniProt ID: E2DHI8), were aligned to hACE2 sequence (UniProt ID: Q9BYF1) using Clustal Omega [35, 36]. Residues 1–358, including all the amino acids interacting with S protein RBD, were considered for alignment.

Analysis of ACE2 sequences of wild and domestic animal species present in Italy

Sequence similarity search over domestic and wild animal ACE2 sequences was carried out using as template 1–358 ACE2 residues of Rhinolophus ferrumequinum and Myotis daubentonii. ACE2 sequences of species not living in Italian territory were excluded. ACE2 sequences of domestic and wild animals with distribution areas on the Italian territory were selected if presenting sequence homology not less than 70%.

Analysis of human SARS-CoV-2 S and spike protein of Bat coronavirus Rp3

Sequence alignment was carried out using the RBD (residues 339–490) of bat coronavirus Rp3/2004 S protein (UniProt ID: Q3I5J5), and residues 336–518 of SARS-CoV-2 S protein (UniProt ID: P0DTC2). The alignment was carried out using Clustal Omega [35, 36] (see Table 1).

Table 1

ACE2 sequences of wild and domestic species present in Italy with their UniProt codes and the percentage of identity with human and ACE2 of Chiroptera species present in Italy.

Wild and domestic species present in Italy	UniProt ID	% Identity with hACE2	% Identity with M. daubentonii ACE2	% Identity with R. ferrumequinum (E2DHI2) ACE2	% Identity with R. ferrumequinum (B6ZGN7) ACE2
Equus caballus (Horse)	F6V9L3	84.72	82.40	84.36	84.92
Felis catus (Cat)	Q56H28	83.06	80.44	82.12	82.68
Oryctolagus cuniculus (Rabbit)	G1TEF4	83.90	81.84	79.33	79.89
Ursus arctos horribilis (Grizzly bear)	A0A3Q7TE16	82.22	78.49	80.16	80.73
Canis lupus familiaris (Dog)	A0A5F4BS93	80.56	77.65	78.49	79.05
Vulpes vulpes (Red fox)	A0A3Q7RAT9	80.27	77.93	79.05	79.61
Mus musculus (Mouse)	Q8R0I0	80.83	78.61	78.71	79.33
Sus scrofa (Pig)	K7GLM4	80.55	77.34	77.93	78.49
Bostaurus (Cattle)	Q58DD0	78.67	76.19	78.49	79.05
Ovis aries (Sheep)	W5PSB6	79.10	75.98	78.27	78.83

Results

Based on the evidence that SARS-CoV-2 S derives from the bat [19], the first step of our study was aimed to identify Chiroptera living in Italy that may be new reservoirs of SARS-CoV-2 and, as such, responsible for new outbreaks in Italy.

hACE2 receptor is essential for infection of SARS-CoV-2 in humans. The crystal structure of SARS-CoV-2 S receptor bound to ACE2 has been recently solved [9]. In a preliminary step, we analyzed hACE2/SARS-CoV-2 S complex and identified in hACE2 the residues K31, H34, E37, D38, Q42, K353, and N330 as those critical for the interaction with S protein; among them, K31 and K353, revealed to be critical for S recognition [1].

The study of the binding mode of SARS-CoV-2 S protein in complex with ACE2 of Chiroptera may offer important information on the possibility that these species could be a reservoir for the virus. Using ExPASy SWISS-MODEL and having as a template the crystallographic coordinates of hACE2 (PDB ID: 6M0J) [9] we calculated the homology models of the ACE2 proteins presenting not less than 70% sequence homology with hACE2 (Figure 2B). These were subjected to molecular docking against SARS-CoV-2 S protein, imposing as restraints the distances with K31 and K353, that are the most critical residues for S recognition. The data resulting from molecular docking indicates similar binding stability for all the ACE2/SARS-CoV-2 S protein complexes (see Supplementary Table S2). Analysis of the structural models evidences the conservation of the residues within the binding site, highlighting a critical similarity of hACE2 with the ACE2 of the Chiroptera species. In particular, we focused on ACE2 of Myotis daubentonii and Rhinolophus ferrumequinum. These are the solely ACE2 sequences retrieved from the UniProt database that belong to Chiroptera species living in Italy. The homology models of ACE2 belonging to Chiroptera species in complex with S protein were subjected to 1 ns molecular dynamics simulations in water using protein S/hACE2 complex (PDB ID: 6M0J) as control.

Figure 2

(A) Rhinolophus ferrumequinum and Myotis daubentonii ACE2 in complex with SARS-CoV-2 S RBD. The structures are superimposed to hACE2 structure in complex with SARS-CoV-2 S RBD. The data, extracted from the last step of molecular dynamics, are relative to 1–358 ACE2 residues. (B) Multiple sequence alignment including hACE2 and ACE2 of Rhinolophus ferrumequinum and Myotis daubentonii. Residues involved in the interaction with S protein are shown: the residue numbers in red represent the residues essential for hACE2/SARS-CoV-2 S binding. Conserved residues are highlighted in blue, while similar residues are highlighted in green. Bold black residue numbers indicate amino acids essential for the stability of hACE2/SARS-CoV-2 S interaction.

Analysis of molecular dynamics considering the variation of RMSD during the time indicated each complex reaching a conformational steady-state, as observed from the presence of RMSD plateau (see Supplementary Figure S1) [35]. In Figure 2A the structural coordinates of Myotis daubentonii and Rhinolophus ferrumequinum ACE2 derived from the last step of molecular dynamics are superimposed to the crystal structure of hACE2. The poses are very similar (RMSD on backbone heavy atoms: 5.83 Å) and the residues defining the protein-protein interaction framework are almost conserved (see Supplementary Tables S3–S7, Figures S2–S6).

Figure 2B shows the sequence alignment of 1–358 residues of hACE2 and ACE2 of Myotis daubentonii and Rhinolophus ferrumequinum (Omega Clustal). 78.77% sequence identity is observable for Myotis daubentonii, 77.93% for Rhinolophus ferrumequinum (UniProt ID: E2DHI2), and 81.24% for Rhinolophus ferrumequinum (UniProt ID: B6ZGN7). The identity refers to all the residues ranging from 1 to 358. Most of the residues essential for the binding with S protein are conserved. In particular, Rhinolophus ferrumequinum preserves the residues E37, Q42, N330, and K353, with D38 (Homo sapiens) replaced by N38 (Rhinolophus ferrumequinum), which conserves the physical-chemical properties of the amino acid. Myotis daubentonii preserved H34, E37, D38, and K353, with K31 and Q42 (Homo sapiens) replaced by N31 and E42 (Myotis daubentonii), respectively. Considering these results together with the data deriving from in silico calculations, Rhinolophus ferrumequinum is to be considered a potential primary animal reservoir within the Italian distribution area [37].

Once we identified Rhinolophus ferrumequinum as a potential primary reservoir of SARS-CoV-2 in Italy [19], the objective of our study was to identify a hypothetical secondary animal reservoir in the Italian territory. If confirmed by experimental data, this information is critical for predicting a possible contagion line in the wild and domestic fauna present in Italy.

Residues 1–358 of ACE2 sequences of Rhinolophus ferrumequinum (UniProt ID: E2DHI2 and B6ZGN7), Myotis daubentonii (UniProt ID: E2DHI8) and Homo sapiens (UniProt ID: Q9BYF1) were used as a template in a similarity search (UniProt database) to identify ACE2 sequences of wild and domestic animals present in the Italian region. As a result, ACE2 sequences of 10 wild and domestic animal species present in Italy were selected, showing no less than 70% identity on 1–358 residues (Figure 3).

Figure 3

Multiple sequence alignments including ACE2 amino acid sequences of Rhinolophus ferrumequinum (UniProt ID: E2DHI2 and B6ZGN7) (A–B), Myotis daubentonii (UniProt ID: E2DHI8) (C), with ACE2 amino acid sequences of wild and domestic animal species present in Italy. The residues in blue are preserved, the residues in green have similar chemical characteristics.

Figure 3B indicates that the binding sites of ACE2 sequences in domestic and wild animal species show the highest number of conserved residues if considered in comparison to ACE2 of Rhinolophus ferrumequinum: the binding site of Rhinolophus ferrumequinum ACE2/SARS-CoV-2 S complex has >80% sequence identity with the binding site of domestic and wild animal species present in Italy.

To evaluate the possibility that the selected wild and domestic animals may be the second reservoir of a mutated form of coronavirus infection, we simulated the interaction between ACE2 of domestic and wild animal species present in Italy (Figure 3B) and Bat-CoV spike protein.

Preliminarily, homology models of the SARS-like bat coronavirus (Bat-CoV/Rp3/2004) S protein and ACE2 proteins of domestic and wild animals were built (see Supplementary Table S8). After, molecular docking calculations were carried out imposing intermolecular distance restraints to preserve the binding poses: 5 Å distance was imposed between ACE2 E37 and Bat-CoV S protein Y477. 1 ns molecular dynamics performed on all the structural models revealed conformational stability consistent with RMSD plateau in RMSD vs time plot. (see Supplementary Figure S7, Table S9).

In Figure 4A, the structure of Rhinolophus ferrumequinum ACE2 in complex with Bt-CoV/Rp3/2004 S RBD is superimposed to ACE2 structures (residues 1–358) of domestic and wild animal species as derived from the last step of 1 ns molecular dynamics. The ACE2 structures are very similar (RMSD on backbone heavy atoms: 4.76 Å), suggesting that the BtCoV/Rp3/2004 S protein is potentially adapt to undertake efficacious interaction with ACE2 of domestic and wild animals. In Figure 4B multiple sequence alignments indicate that ACE2 residues involved in the interaction with BtCoV/Rp3/2004 S protein are conserved in all domestic and wild species except for Canis lupus familiaris and Vulpes vulpes.

Figure 4

(A) The structure of Rhinolophus ferrumequinum ACE2 in complex with Bt-CoV/Rp3/2004 S RBD is superimposed to ACE2 (residues 1–358) of domestic and wild animal species as derived from the last step of 1 ns molecular dynamics ribbon representation. (B) Residues of ACE2 (belonging to Italian domestic and wild animals) interacting with Bt-CoV/Rp3/2004 S RBD, compared to residues of hACE2 interacting with S protein RBD. Residues highlighted in blue are preserved; residues in green have similar chemical characteristics.

Discussion

Italy is among the European countries the wealthiest for biodiversity [17]. (https://www.isprambiente.gov.it/en/archive/news-and-other-events/ispra-news/year-2015/may/biodiversity-in-italy). As mammal biodiversity has been recognized as risk factor for zoonotic disease emergence, the Italian region is exposed to a high risk of SARS-CoV-2 pandemic outbreak and the relative risk of endangered animal extinction.

As a great deal of evidence indicates that SARS-CoV-2 originated from the bat, the present work aimed to identify the Chiroptera species living in Italy that could be the primary reservoir in this region.

Accordingly, performing in silico calculation based on multiple sequence alignment and homology modeling, molecular docking, and molecular dynamics, we identified Rhinolophus ferrumequinum as the best Chiroptera candidate to be the potential primary reservoir for the SARS-CoV-2 virus in Italy. The sequence alignment of Rhinolophus ferrumequinum ACE2 with hACE2 sequence indicated that most of the residues involved in SARS-CoV-2 S binding are conserved. Moreover, Rhinolophus ferrumequinum ACE2 shows the highest identity with the ACE2 proteins of wild and domestic species present in Italy (Figure 3B), indicating a potential role as a human-animals bridge infection in wild and domestic environments.

A low number of domestic animals have been found infected by SARS-CoV-2, despite the high viral circulation during the current pandemic. This may depend on a still inadequate screening of domestic animals or a low human-to-animals transmission rate in natural conditions. The low transmission rate might depend on an innate or acquired resistance of domestic animals to the virus or in the differences in the structure of ACE2 proteins.

Our data show that a conserved pattern of residues is present in ACE2 proteins of many animal species, in particular horses (Equus caballus), cats (Felis catus), cattle (Bos taurus), and sheep (Ovis aries), among the domestic animals, European rabbits (Oryctolagus cuniculus) and grizzly bear (Ursus arctos horribilis) among the wild animals. Considering these data, animal species present in Italy could be considered SARS-CoV-2 secondary reservoirs, or else could behave as host intermediate if infected by bat coronaviruses. On the contrary, ACE2 proteins of animals like dogs (Canis lupus familiaris) and red fox (Vulpes vulpes) that do not exhibit this residue pattern, are not eligible as SARS-CoV-2 host intermediates. These results agree with recently published data that exclude pigs as infection transmitters and recognize dogs as low susceptible to the experimental infection [6]. On the other hand, cats (and other animals) may be a silent intermediate host of SARS-CoV-2, because infected cats have never shown any visible symptoms that might be recognized by their owners [10]. The high conservation rate of ACE2 residues may suggest that also farm animals such as cattle and sheep may be endangered by a viral transmission, and finally, the horse can be a potential intermediate host transmitting SARS-CoV-2 to humans.

Taken together, our data point to Rhinolophus ferrumequinum as a potential primary reservoir for SARS-CoV-2 in Italy, and domestic and wild animal species can act as the second reservoir. However, the similarity of ACE2/spike protein complexes among these species does not explain the low human-to-animals transmission rate of SARS-CoV-2 infection. An innate or acquired resistance of domestic animals to the virus potentially accounts for the rare infection of domestic or wild animals occurring during the pandemic.

Declarations

Author contribution statement

Michela Buonocore, Carmen Marino: Performed the experiments; Wrote the paper.

Manuela Grimaldi, Angelo Santoro, Mohammad Firoznezhad: Analyzed and interpreted the data.

Orlando Paciello: Conceived and designed the experiments.

Francesco Prisco, Anna Maria D'Ursi: Conceived and designed the experiments; Wrote the paper.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.