Detection of Coronaviruses in Bats in Lebanon during 2020

Bats are considered the main reservoir of coronaviruses (CoVs), and research evidence suggests the essential role of bats in the emergence of Severe Acute Respiratory Syndrome Coronaviruses (SARS-CoV) and SARS-CoV-2. SARS-CoV-like viruses have been recently detected in bats in different countries. In 2020, we conducted surveillance for CoVs among six different bat species in Lebanon. Of 622 swab specimens taken, 77 tested positive. Alpha- and Beta- CoVs were identified in samples collected from different species. Our results show that SARS-like coronaviruses circulate in bats in this region, and we provide new data on their genetic diversity. The interaction between the spike of the detected SARS-CoV-like viruses and the human angiotensin-converting enzyme 2 (hACE2) receptor could be crucial in understanding the origin of the epidemic. The 3D protein structure analysis revealed that the receptor-binding domains of the SARS-like virus identified in Lebanon bind to the hACE2 protein more efficiently than to the spike of the SARS-CoV-2 strain. The spike of the detected SARS-CoV-like viruses does not contain the recognition site of furin at the cleavage site. Thus, our study highlights the variety of bat coronaviruses in Lebanon and suggests the zoonotic potential for other SARS-CoV-like viruses.


Introduction
After rodents, bats are the second most diverse group of mammals in the world with about 1400 different species. Bats (order Chiroptera) have recently divided mainly into Yinpterochiroptera and Yangochiroptera suborders. Bats act as important reservoir hosts for a variety of neglected and non-neglected viruses and play crucial roles in the spillover of different viruses, including emerging coronaviruses (CoVs). Some of these emerging CoVs have acquired mutations to adapt to humans directly or through intermediate hosts.
Six human CoVs were characterized by the end of 2019, including CoV-OC43, CoV-229E, CoV-HKU1, CoV-NL63, Middle East Respiratory Syndrome CoV (MERS-CoV), Severe Acute Respiratory Syndrome CoV (SARS-CoV), and SARS-CoV-2. The close relationship between the six human CoVs and CoVs identified in animals supports the zoonotic origin of these viruses. The discovery of multiple SARS-CoV-related viruses in bats highlights their role in CoV transmission [1][2][3][4]. Moreover, serological evidence indicated the presence of antibodies against bat SARS-CoVs and HKU10-CoVs at a low prevalence rate among rural inhabitants in China which highlights the potential transmission of these viruses to humans [5]. Further, the COVID-19 pandemic emphasizes the significance of bat-borne viruses and the need for continued monitoring and research to understand their emergence and transmission. At early stages of the pandemic, genetic characterization of the genome of SARS-CoV-2 revealed that this newly-emerging virus was closely related to SARS-CoV and several SARS-like coronaviruses identified in bats. Later, several strains of SARS-CoV-2-related coronaviruses had been identified in pangolins and bats. Although the identification of multiple SARS-like coronaviruses strains in wildlife supports the role of wildlife in the spillover of CoVs, the evolutionary history of SARS-CoV-2 is unclear.
In Lebanon, 21 species of bats exist [6], such as Egyptian fruit bats (Rousettus aegyptiacus), horseshoe bats (Rhinolophus sp.), and Schreiber's bent-winged bat (Miniopterus schreibersii). During the last decades, bats in their natural habitats in Lebanon have faced numerous threats, including climate change, habitat loss, hunting, and disturbance of roosting sites [7]. Understanding dynamics and diversity of CoVs in bats should be considered to prevent zoonotic transmission and add to the knowledge of virological features prior to the emergence of viruses. Our active surveillance from 2013-2015 showed the presence of HKU9-related CoV in R. aegyptiacus in Lebanon [8]. In this study, we aim to characterize the prevalence of CoVs in different species of bats in Lebanon during the SARS-CoV-2 pandemic.

Sample Collection and Processing
Between 27 June and 4 August 2020, we captured 311 bats from 6 different species using mist nets in mountain caves located in Akkar (n = 200), Zgharta (n = 46), Zahle (n = 62) and Saida (n = 3) (Figure 1). These included greater horseshoe bats (Rhinolophus ferrumequinum, n = 109), lesser horseshoe bats (Rhinolophus hipposideros, n = 49), Schreiber's bent-winged bats (Miniopterus schreibersii, n = 51), mouse-eared bats (Myotis sp., n = 76), long-fingered bats (Myotis capaccinii, n = 6), and Geoffroy's bats (Myotis emarginatus, n = 20). The bats were identified using external morphology and confirmed through DNA barcoding. Oral and rectal swabs were collected from each bat. Immediately after, the bats were released. The swabs were placed in a viral transport medium and transported to the laboratory on ice. The samples were then stored at −80 • C for further processing. The protocols for bat capture and sampling were approved by the Institutional Animal Care and Use Committee at St. Jude Children's Research Hospital, USA (No. 3130).
Viral RNA was extracted from the pooled oral and rectal swabs using the QIAamp viral RNA extraction kit (Qiagen, Hilden, Germany). The extracted RNA was subjected to RT-PCR using the One-Step RT-PCR kit (Qiagen, Hilden, Germany). The PCR cycler conditions for the amplification were 50 • C for 30 min followed by 95 • C for 15 min, then 50 cycles of 94 • C for 30 s (denaturation), 50 • C for 30 s (annealing), 72 • C for 60 s (extension), then 72 • C for 10 min (final extension). PCR products of 441bp were visualized on a 1% agarose gel electrophoresis, and the targeted amplicons were purified using the QIAquick gel extraction kit (QIAgen) and sequenced at the Macrogen DNA sequencing service. The sequences were assembled using SeqMan DNA Lasergene 15 software (DNASTAR, Madison, WI, USA) and submitted to GenBank (MW880969 to MW881013). Complete spike sequences of 3 representative samples of β-CoVs were amplified using primers designed based on BM48-31/BGR/2008 (NC 014470.1) (Table S2) and using the One-Step RT-PCR kit (Qiagen, Hilden, Germany). The PCR cycler conditions for the amplification were 50 • C for 60 min followed by 95 • C for 15 min, then 50 cycles of 94 • C for 30 s (denaturation), 50 • C for 30 s (annealing), 72 • C for 4 min (extension), then 72 • C for 10 min (final extension). Amplicons were then purified using the QIAquick gel extraction kit and sequenced at The Hartwell Center for Biotechnology, St. Jude Children's Research Hospital. Sequences were submitted to GenBank.

Phylogenetic Analysis
A BLASTN analysis was conducted, and closely related sequences were downloaded from GenBank. Previously identified CoVs in Lebanon and full genome sequences of identified CoVs in bats were included. The final dataset was aligned using MAFFT v7.305b [10], manually trimmed, and phylogenetic trees were constructed using the maximum likelihood method with Kimura's 2-parameter distance model and 1000 bootstrap using MEGA X [11].

Homology Modeling, Protein-Protein Docking, and Molecular Dynamic Simulation
The 3D structure of the spike protein of the Lebanon SARS-CoV-2-like virus was predicted through homology modeling, with the 3D structure of SARS-CoV-2 spike glycoprotein obtained from PDB (PDB ID: 5WRG, 4.30 Å [12] using the MOE software. The validity of the model was confirmed using the Maestro software, which resulted in an RMSD of 2.380 Å compared to the 5WRG structure. A Ramachandran plot of the novel structure was generated and is shown in Figure S2. Protein-protein docking was performed using HADDOCK [13,14]. The ACE2 receptor in complex with the spike protein (PDB ID: 6M0J) [15] was prepared by removing the spike chain and heteroatoms using the MOE software. The structures of ACE2 along with the wide and novel spike proteins were uploaded to the HADDOCK server, and the docking poses were analyzed. Only the poses that docked around the receptor binding domain (RBD) of the spike protein were considered for further analysis using PyMOL (PyMOL Molecular Graphics System, Version 1.2r3pre, Schrödinger, LLC) and UCSF Chimerax 1.5 [16]. Finally, the docked poses were subjected to a 500 ns molecular dynamic simulation to evaluate the stability and interaction between the ACE2 receptor and spike. The details of MD simulation are provided in the Supplementary Materials.

Statistical Analysis
Statistical analyses were performed using GraphPad Prism 9 (GraphPad, San Diego, CA, USA). The association between the prevalence of CoVs among different variables was analyzed using the Chi-square test. Statistical significance was considered at p-values less than 0.05.

Diversity of CoVs in Lebanon Based on Partial Sequence of RdRp
Whole genome sequencing was attempted but was not successful. Blast analysis of 42 partial sequences retrieved using Sanger sequencing showed 28 sequences were from Rhinolophus ferrumequinum Betacoronaviruses, and the remaining 14 were related to Alphacoronaviruses (Table S3). All Betacoronaviruses detected in this study formed a single phylogenetic cluster and were SARS-like. They were most closely related to bat coronavirus BM48-31/BGR/2008 and bat coronavirus Khosta-1, which were identified in Rhinolophus blasii and Rhinolophus ferrumequinum from Bulgaria and Russia in 2008 and 2020, respectively (Figure 2a). However, Alphacoronaviruses in Lebanon formed three divergent lineages (Figure 2a). The majority (n = 9) were most closely related to BR89-55/GBR/2008 collected from Miniopterus schreibersii in Germany (bootstrap, 99%), and one sequence (Bat 51A) was most closely related to the Miniopterus bat CoV HKU 8 (Bat-CoV HKU8) with a bootstrap of 91%, whereas the other lineage comprising four sequences was related to the Porcine epidemic diarrhea virus (PEDV).

Genetic and Structural Characterization of the Spike Glycoprotein of Lebanon SARS-Like CoVs and the Potential Significance of Binding to Human ACE2
To better understand the risk that Lebanon SARS-like CoVs pose, we sequenced the complete spike protein genes of three representative samples. They were highly similar throughout the gene (Figure 3). These spike proteins were most closely related to the bat CoV Khosta-1 strain collected during 2020 in Russia and had a 94% similarity. A comparative analysis of the spike proteins of Lebanon Betacoroviruses and SARS-CoV, which showed five amino acid deletions at the S1/S2 junction (RQQ) in Lebanon Betacoronaviruses indicated the absence of a furin cleavage site (Figure 3). Among the five amino acids known to be critical for binding the S glycoprotein to the ACE2 receptor, (501N, 505Y, 455L, 486F, and 493Q), the 505Y and 455L amino acids were found among the Lebanon bat CoV.
Homology modeling of the chain A sequence of spike protein (Figures 4a,b and S3  Docking results from the HADDOCK were derived. Since minimum HADDOCK scores depicted more binding affinity of the two interacting proteins, best poses of the docking were obtained from the given results. The novel complex demonstrated significantly greater binding affinity and stability than the control complex, which had a lower HAD-DOCK. Z-score indicated how many standard deviations from the average the resultant complex was located in terms of score (the more negative, the better). Docking scores of both complexes along with other parameters are demonstrated in Table 2. Next, the interactions between the spike protein (wild and novel) and the ACE2 were analyzed. The H-bond interactions, along with their residues and distances, are reported in Table 3. In contrast to the wild-type spike, which formed 14 H-bonds with the ACE2 receptor (Table 3), 8 H-bonds were formed between ACE2 and chain C, while 6 were formed with chain A. On the other hand, the novel protein spike was able to form 22 H-bonds with the ACE2 receptor: 13 H-bonds were formed via interaction with chain C, while the other nine H-bonds were formed with chain B. The interactions of both wild and novel proteins are presented in Figures S3-S8. To validate the results obtained from docking, dynamic molecular simulation was implemented to study the stability of the wild-ACE2 and novel spike-ACE2 complexes. Our focus was only on the interactions between the ACE2 and receptor binding domain. Hence, only the RBD was considered for the MD calculations, and both complexes were subject to a 500 ns simulation. The root mean square deviation of the Cα of the proteins was plotted as a function of simulation time to measure the impact of the simulation on the stability of the SARS-CoV-2 spike-ACE2 complexes. As seen in Figure 5, the wild SARS-CoV-2 spike-ACE2 complex fluctuated at 9.50 Å with respect to its original Cα atoms at the start of the simulation time. The complex further oscillated at around 3.00 Å until around 300 ns of the simulation time, and, at around 350 ns, the complex reached a plateau and fluctuated at around 1.00 Å until the end of the simulation. On the other hand, the novel SARS-CoV-2 spike-ACE2 complex showed a similar pattern but reached stability earlier at around 180 ns of the simulation time. The complex fluctuated at around 9.00 Å at the beginning of the simulation and reached stability at 180 ns, and the complex fluctuated at around 1.00 Å toward the end of the simulation. The H-bond interactions were monitored during the simulation, and snapshots of the spike-ACE2 were taken at each 50 ns ( Table 4). As seen in Table 4, the complex holds around 15-22 H-bonds during the simulation. The residues involved in the formation of these interactions are reported in detail in Table S4.

Discussion
Horseshoe bats (Rhinolophidae) are considered potential reservoir hosts for several viruses that can cause disease in humans and animals. They are widely distributed in Asia, Europe, and North Africa. In East and Southeast Asia, SARS-CoV-like viruses have been found in multiple rhinolophids [1,3,17]. Recently, SARS-like CoVs have been identified in Rhinolophus sp. in Russia during 2020 [18]. However, information about the prevalence and genetic characteristics of CoVs in Middle Eastern bats is limited. To address this, we conducted a cross-sectional surveillance for CoVs among six bat species in Lebanon during 2020. Not surprisingly, we identified several Alpha-and Betacoronaviruses in samples collected from different species of Microchiroptera in Lebanon. However, the zoonotic origin of the newly-emerging coronaviruses (SARS-CoV, SARS-CoV-2, and MERS-CoV) remains unclear. Moreover, there are genetic differences among detected viruses in humans and animals, thus supporting the presence of intermediate hosts or missing genetic information. However, the limitations of this study include the lack of virus isolation in culture and whole genome sequences, as well as in vitro data on binding and serology.
The closely genetically-related β-CoVs detected in bats with SARS-CoV-2 (NC_045512) were RaTG13 (in Rhinolophus affinis), and RpYN06 (in Rhinolophus pusillus) viruses in China with sequence identities 96.10 and 94.48%, respectively, across full genomes [2]. However, the spike gene has a much lower sequence identity, suggestive of a genomic recombination event, thus making RpYN06 the second closest relative of SARS-CoV-2 identified to date after the bat CoV RaTG13, hence highlighting the potentially complex evolutionary history of SARS-CoV-2 [2,4]. Other SARS-CoV-2-related CoVs such as RmYN02, STT182 and STT200 were detected in R. malayanus and R. shameli bats in Asian countries [1,17].
However, for emerging zoonotic spillovers, humans exposed to the viruses should be infected through direct contact with the infected reservoir or intermediate host of coronaviruses. Other factors such as environmental factors play a role in emerging zoonotic viruses. In Lebanon, there is almost no direct contact between humans and bats, and the intensity of the reservoir bat-human interface is low, which decreases the potential for spillover risk. Globally, extensive studies of the genetic diversity of coronavirus in bats could bridge our knowledge gaps associated with the emergence of SARS-CoV-2.
In line with previous studies [2], analysis of the 3D protein structure revealed that the receptor-binding domains of the SARS-like virus identified in Lebanon bind more efficiently to the hACE2 protein than to the spike of the SARS-CoV-2 strain. Risk assessment of viruses that are potentially infectious to humans should be prioritized by national, regional, and international agencies. In summary, our results suggest a diverse range of CoVs in bats, with Betacoronaviruses and Alphacoronaviruses co-existing in the same species. The high prevalence of CoVs in Rhinolophus ferrumequinum, in particular, underscores the importance of continued surveillance and investigation of this species in understanding CoV emergence. Furthermore, our results underline the importance of continued monitoring of CoV in bats, especially in SARS-like Betacoronaviruses that have repeatedly emerged in humans.
Surveillance for coronaviruses in bats is critical as an alarm system for emerging and remerging viruses through early detection, monitoring transmission dynamics, identifying high-risk species or regions, tracking genetic diversity, and informing public health policies.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/pathogens12070876/s1, Figure S1. Graphic illustration of the selected conserved regions that was created via a web-based WebLogo application (http://weblogo. threeplusone.com/create.cgi, accessed on 20 May 2023); Figure S2. The Ramachandran plot of the novel SAR-CoV-2 spike; Figure S3. Validation of the homology modelling results of the novel SAR-CoV-2 spike; Figure S4. Sequence alignment of the wild and the novel spike proteins; Figure S5. Wild-ACE (green) interaction, Figure S6. Novel-ACE2 (green) interaction; Figure S7. Wild-ACE2 interactions; Figure S8. Novel-ACE2 interactions; Table S1. Representative α-, β-, γ-, and δ-CoV sequences retrieved from National Center for Biotechnology Information (NCBI) that was used for degenerated primers design; Table S2. Primers used for sequencing of the spike; Table S3. Nucleotide sequence identities between partial RdRp sequences detected in bats in Lebanon and nearest published sequences in GenBank based on BLASTN analysis; Table S4. Interactions between the novel spike and ACE2 receptor at different simulation times.