Identification of Some Potential Host Cellular Receptors and Enzymes that Potentially Refine the Bovine Coronavirus (BCoV) Tissue Tropism, Replication, and Pathogenesis

1. Introduction

Virus replication occurs in many consecutive steps to generate mature viral particles in the permissive cells. The virus always hijacks the cellular machinery and directs them to synthesize their viral proteins instead of cellular one’s [1,2]. The viral entry is a crucial step in any viral infection, including BCoV [3]. The viral entry involves many host cells and viral proteins [1]. The presence of viral-specific receptors is considered an important factor in the process of viral entry into the host cells. The host cell is called permissive for certain viruses if they express the specific viral receptors and provide suitable environments for virus replication, including any auxiliary receptors, transcription, and translation factors. Coronaviruses are enveloped viruses containing positive sense RNA genomes and belong to the order Nidovirales and are classified into four genera (α, β, λ, and δ) [4]. The genus β-coronavirus includes five important human coronaviruses (the Sever Acute Respiratory Syndrome Coronavirus-1 (SARS-CoV-1), the Middel East respiratory syndrome coronavirus (MERS-CoV), the SARS-CoV-2, human coronavirus-HKU1 (HCoV-HKU1), and the human coronavirus-OC43 (HCoV-OC-43). This genus β-coronavirus includes some other important viruses affecting animals, particularly the BCoV and the equine coronavirus (ECoV) [4]. The coronavirus’s genome size ranges from 27-31 kb in length and has a unique organization. The full-length genome is flanked with two untranslated regions at the 5′ and 3′ ends. Coronaviruses are characterized by the production of a subset of sub-genomic messenger RNA (mRNA) at their 3′ end [5]. The 5′ end of the genome of most coronaviruses contains a large gene called Gene-1, which consists of two overlapping open reading frames (ORFs) with a ribosomal frameshifting between those two ORFs. However, the 3′ end of the genome is mainly occupied by the common structural proteins interspersed with some small accessory proteins. There are four major structural proteins in most coronaviruses, including the spike glycoprotein (S), the envelope (E), the membrane (M), and the nucleocapsid protein (N). Some members of the genus β-coronavirus, including BCoV and HCoV-OC43, have an additional structural protein called hemagglutinin esterase (HE); thus, their genome is a little larger in size compared to other coronaviruses (31 Kb) [6]. The S glycoprotein is a key player in all coronavirus replication. There are several proteins, including some cellular receptors, co-receptors, and cellular enzymes are involved in the BCoV/host interaction. Both the BCoV spike (BCoV/S) and (BCoV-HE) proteins play important roles during BCoV replication and pathogenesis (6). BCoV/S has potential binding to the 5-N-acetyl-9-O-acetylneuraminic acid (Neu5,9Ac2), suggesting its possible roles as BCoV receptors [7]. On the other hand, the BCoV-HE acts as a receptor-destroying enzyme during BCoV replication [7]. However, there is a lack of comprehensive understanding of the interplay of the BCoV-S/BCoV-HE with the cellular receptors during BCoV replication. The availability of specific receptors is one of the main factors that make the target cells permissive to coronavirus (CoVs) infection. Each group of CoVs recognizes certain types of receptors and may require the presence of additional auxiliary receptors to facilitate virus attachment and downstream replication. The SARS-CoV-2 uses the angiotensin-converting enzyme 2 (ACE2) as the main receptor, and the chaperone GRP78 acts as an auxiliary receptor [8,9]. It was also shown that MERS-CoV utilizes the dipeptidyl peptidase-4 (DPP4) as receptors in humans and dormitory camels [10,11]. The amino-peptidase N (APN), also called cluster of differentiation -13 (CD13), acts as a major receptor for the transmissible gastroenteritis virus that causes enteric infections in pigs [12]. The carcinoembryonic antigen cell adhesion molecule 1 (CEACAM-1) also acts as a major receptor to another coronavirus called murine hepatitis virus (MHV) [13]. Although the presence of the CoVs receptors and co/receptors is important for the success of viral replication, most coronaviruses require the presence of host cell enzymes that help in the cleavage of the CoV-S and CoV-HE proteins to initiate the process of viral infection. Usually, these host cell enzymes are highly expressed in the target tissues, which play crucial roles during the CoVs entry to the host cells, particularly the mucosal surfaces of the respiratory and enteric tracts of the affected hosts. Homology modeling is a well-established method that has been shown to produce quite accurate models for a protein sequence if an X-ray structure of a protein with a sufficient degree of sequence similarity is available [14]. The method is based on the fact that the structural conformation of a protein is more highly conserved than its amino acid sequence and that small or medium changes in sequence normally result in little variation in the 3D structures [15]. The quality of the model is directly linked to the identity between template and target sequences. As a rule, models built with over 50% sequence similarities are accurate enough for drug discovery applications [16].

During the CoVs replication, various host cell proteases, including transmembrane serine protease 2 (TMPRSS2), furin, and cathepsin-L (CTS-L), are involved in processing the S protein. However, the precise order of protease cleavage and the interactions between these host proteases are not yet fully understood [17]. It has been recently shown that the TMPRSS2 (transmembrane protease serine-2) is a serine protease that plays an important role in SARS-CoV-2 replication, particularly during viral entry to the host cells. The TMPRSS2 usually cleaves the viral spike glycoprotein, which enhances the interaction between the virus and the cellular receptors, thus facilitating its entry into the host cells [18]. The Furin is another host cell enzyme that is a subtilisin-like proprotein convertase. Furin usually cleaves the target proteins at a polybasic amino acid sequence R-X-(K/R)-R (where R is arginine, K is lysine, and X can be any amino acid). The Furin cleavage to the SARS-CoV-2-S protein enhances the pathogenicity transmissibility and increases the virus infectivity in the target host [19,20]. Although the roles of the receptors mentioned above and enzymes were intensively studied in SARS-CoV2, there is a lack of knowledge about the roles of these receptors and enzymes in BCoV infection and the replication cycle. The main goals of the current study are (1) to use the in-silico prediction and docking tools to study the roles of these proteins and enzymes in BCoV replication. (2) to confirm the molecular docking results and in silico prediction of the short-listed potential receptors and host enzymes by drafting the host gene expression profiles during the BCoV replication of some bovine host cells. The results of this study highlighted some unknown aspects of BCoV tissue tropism and pathogenesis through the identification of some potential novel receptors and enzymes involved in BCoV entry and virus replication. It will also pave the way for the development of some novel vaccines and antiviral therapies for BCoV infection in cattle.

2. Materials and Methods

2.1. BCoV Spike (BCoV/S) and Host Receptor Protein Structure Prediction

Homology modeling was used to build the predicted homology structures of proteins based on 3D confirmations of most identical template proteins. The MODELLER tool of the Biovia Discovery Studio v22.1.021297 can be used to predict the protein/protein interactions through a comparative modeling approach. It has been used to indicate the 3D structural model for query protein under study. The protein sequences of the BCoV-S, the bovine (ACE2, Furin, TMPRSS2, Cathepsin-L, NRP1, DPP4, APN, CEACAM-1, and AXL) were retrieved and downloaded in FASTA format from the NCBI (SciENcv - Home (nih.gov) and Uniprot databases (https://www.uniprot.org/).

The protein sequence of BCoV/S glycoprotein (references strain; Mebus strain (ID: P15777) was downloaded from the UniProt database. The query sequences of other bovine host proteins including the (ACE2 (ID: Q58DD0), Furin (ID: Q28193), the TMPRSS2 (ID: A2VDV7), the NRP1 (ID: E1BMX5), DPP4 (ID: P81425), the APN (ID: P79098), CEACAM-1 (ID: Q6VAN8) and AXL (ID: F1N0D3)) were retrieved from the UniProt database. The homology modeling procedure requires the alignment of the query protein sequences with template protein sequences from the BLAST method (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastp&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome). The query protein sequence aligned with 100 most identical template sequences with a maximum identity of 99% and a minimum of 30% in BLAST with BLOSSUM62 algorithm protocol. Five identical template sequences (identity: maximum 99-minimum 65%) were chosen to be aligned by load structure and alignment method. Further, the 3D-build homology model protocol was run to align the template structures and predict query sequence 3D confirmation based on template proteins.

About (1-10) different homology models per query sequence can be generated using this protocol. It also gives us the opportunity to choose the best-predicted model. In this study, we predicted the top five different homology models for each query sequence. We used the Discrete Optimized Protein Energy (DOPE) score to assess the quality of the generated five different 3D protein models. The DOPE and the PDF total energy scores were used to evaluate the accuracy of the protein models. The more negative values of the DOPE score and the PDF total energy, the more stable and accurate the predicted protein model. The best model having the lowest DOPE score was selected for further in silico computational studies [21]. The correspondence homology model is then energy minimized (add hydrogen and CharMm force field) using the energy minimization method (full minimization). All the generated 3D structures of the all-atom models were verified through the Ramachandran plot [22] and the Verify 3D model tool of the Biovia Discovery Studio v22.1.021297. The Ramachandran plot suggested the stability of structure based on amino acid residues located in the most favored, highly allowed, and allowed region or site of the plot. If the verified score from the Verify 3D method of the protein model is higher than the verified expected low Score value, then the model is of acceptable quality. The closer the verified score result is to the verified expected high Score value, the better the quality of the model.

2.2. Mapping the BCoV-Spike Protein Domains

The InterProScan is an international initiative that was conceived to streamline the efforts of the signature database providers (http://www.ebi.ac.uk/interpro/) and contains unique signature patterns that could be diagnostic for protein families, domains, or functional sites. It has important tools for the computational analysis of the functional classification of the newly determined sequences that lack biochemical characterization [23]. We used InterPro’s family, domain, and functional site definitions to annotate and classify uncharacterized sequences. The sequence-based searches are done using the InterProScan (https://www.ebi.ac.uk/interpro/about/interproscan/) and the PROSITE server (https://prosite.expasy.org/, Release 2024_03 of 29-May-2024), which combines the search methods from all mentioned databases. The Web interface allows text-based and sequence-based searches using a sequence retrieval system (SRS) [24]. The InterPro is accessible for interactive use via the EBI Web server (http://www.ebi.ac.uk/interpro), which can also be reached via each of the member databases [25].

2.3. Molecular Docking with CDOCKER

The protein structure of the BCoV-S N-terminal domain (NTD) (residues 15–298), the BCoV-HE protein structure file was downloaded in the PDB file format from the Protein Data Bank (PDB) (https://www.rcsb.org/). In the molecular docking process of BCoV protein with the ligand, the N-terminal domain (NTD) of BCoV spike protein (PDB ID:4h14) was used. To confirm the structure similarity of NTD-from BCoV/S homology modeled and BCoV/S-NTD, we analyzed/compared the amino acid sequences from the PDB database (PDB ID: 4H14). First, protein sequences were aligned, and then the homology model as template and 4H14 as target sequence’s projected 3D-structure with the alignment calculated using PyMOL molecular viewer yielded an RMSD value of 0.17. The ligand molecule, sugar molecule N-acetyl-9-O-acetylneuraminic acid (Neu5,9Ac2), the most common type of sialic acid, generally acts as the terminal sugar in cell surface glycans and polysaccharides, was taken in 3D-SDF file format from PubChem database (https://pubchem.ncbi.nlm.nih.gov/). On the other hand, heparan sulfate (HS) glycan (PubChem CID), a compound from the bovine membrane receptor heparan sulfate proteoglycan 2 (HSPG2), was taken in 3D-SDF file format from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The HS was also docked to BCoV/S with the help of the CDOCKER tool. The interaction of HS was investigated to see the interaction site in the BCoV/S protein. The proteins and ligand molecules were prepared using the built-in protein and ligand preparation methods using the Biovia Discovery studio v22.1.021297. This process typically involves several subsequent steps, including removing the water molecules, adding hydrogens, and minimizing energy to ensure accurate molecular docking. Both protein and ligand structures are manually examined and reverted to pre-reaction form. Grids with a grid space of 0.5 Å represent the receptor structure implicitly. We typically use ten docking poses in the CDOCKER tool from Biovia Discovery Studio docking protocol. Further, ligand-protein interaction binding energy (BE) calculation was done to calculate the binding energy of complexes. The ΔG determines the binding energy and stability of any given protein-ligand complex or the binding affinity of a ligand to a protein. These computational tools enable the visualization of the ligand-target interaction (molecular docking) and the identification of the compounds that bind more efficiently with the target with the highest binding energy [26]. The docking output files were generated and analyzed to identify the most energetically favorable binding poses of the Neu5,9Ac2 with the BCoV/S glycoprotein. This analysis typically involves the examination of the docking scores, the ligand-protein interactions, and the visualization of the docked complexes. To observe the exact binding affinity of the ligand to the protein, the calculation of binding energy is likely an internal step. The binding energy calculations were performed using the Biovia Discovery Studio to define the binding energy of ligands and proteins [27].

2.4. Molecular Docking to Study the Protein-Protein Interaction

The protein-protein molecular docking was performed using the ZDOCK docking tool available in BIOVIA, Discovery studio v22.1.021297. The ZDOCK method uses a fast Fourier transform to perform a 3D search of the spatial degrees of freedom between two proteins. The molecular docking approaches, which take two (or more) structures as input and predict their complex structure, are increasingly being used for this purpose [28]. In each run of the ZDOCK, we used the default angular sampling (3,600 ligand rotations) and a single 2.8 GHz 64-bit Opteron processor with 64 GB available system RAM. The protocol ran for 2000 different interaction poses of protein-protein interaction with 60 clusters. The best pose from the ZDock output files can be selected based on the binding site and largest cluster choice. The resulting poses from the ZDOCK docking method must be further refined with the Energy refinement method (E_RDock or RDock). Therefore, the ZDOCK docking result poses were further refined for the best pose with the highest docking interaction energy/binding affinity with the refinement method, RDock. This server refined protein-protein docking out of 2000 poses of the ZDOCK with multiple conformations and calculates the weighted energies with the scoring function of E_RDock [E_RDock = E_Sol + beta + E_elec2] for every binding conformation, where the E_RDock denotes the scoring function of RDOCK, E_vdw1, and E_vdw2 denotes van der Waals non-bond interaction energy of the protein complex after the first and second CHARMm minimization, E_elec1 and E_elec2 is an electrostatic energy term of the protein complex after the first and second CHARMm minimization, and E_Sol represents desolvation energy of the protein complex calculated by ACE method [28]. For protein-protein complex, the total desolvation score is simply the sum of the ACE scores of all receptor-ligand atom pairs within a distance cutoff of 6 Å. The protein-protein interactions interface was visualized by using the Biovia Discovery studio.

2.5. Bovine Cell Lines and BCoV Isolates

Bovine Pulmonary Artery Endothelial Cells (BEC) (Cat. No. CRL-1733) were obtained from the ATCC. Dr. Udeni B. R. Balasuriya, Louisiana State University, kindly provided the Madine Darby Bovine Kidney (MDBK) cells). Both cell lines were tested for Bovine Viral Diarrhea Virus (BVDV) and found to be negative. BEC cells were cultured in F12 media (ATCC; Cat. No. 30-2004), supplemented with 10% Horse Serum (HS) (Gibco; Ref. No. 26050-088) and 1% 10,000 ug/mL streptomycin and 10,000 units/mL penicillin antibiotics (Gibco; Ref. No. 15140-122). MDBK cells were cultured in Minimum Essential Medium Eagle (MEM) media (Sigma-Aldrich; Cat. No. M0200-500ML), supplemented with 10% horse serum and 1% streptomycin and penicillin antibiotics. Both cell lines were incubated at 37^oC at 5% CO₂ for subsequent culture. Bovine coronavirus (BCoV) enteric isolate ‘Mebus’ [29] was obtained from BEI resources (NIAID, NIH; Cat. No. NR-445). The BCoV respiratory isolate was kindly provided by Dr. Aspen Workman (Animal Health Genomics Research Unit, USDA, A.R.S., US Meat Animal Research) [30].

2.6. The Next Generation Sequencing (NGS)

The BEC cells were independently infected with either BCoV Enteric (BCoV/Ent) or BCoV Respiratory (BCoV/Resp) isolates at a Multiplicity of Infection (MOI) of 1. The cells were monitored for up to 4 days for the development of Cytopathic Effect (CPE). The cells were collected from sham and BCoV-infected groups, and total RNA was extracted using an RNeasy Mini Kit (QIAGEN; Lot. No. 57501521) following the manufacturer’s instructions. RNA samples were submitted to LC Sciences (LLC, 2575 West Belfort Street, Houston, TX 77054, USA) for NGS. RNA library preparation was briefly performed using Illumina’s TruSeq-small-RNA-sample preparation protocols (Illumina, San Diego, CA, USA). Quality control and quantification of the DNA libraries were conducted using an Agilent Technologies 2100 Bioanalyzer High Sensitivity DNA Chip. Single end 50bp sequencing was carried out on Illumina’s Hiseq 2500 sequencing system following the manufacturer’s instructions. Differential expression of mRNAs based on normalized deep-sequencing counts was analyzed using various statistical tests, including the Fisher exact test, Chi-squared 2X2 test, Chi-squared nXn test, the student’s t-test, or ANOVA, depending on the experimental design.

2.7. The RNAs Extraction and the Quantitative Real-Time PCR (qRT-PCR)

Following the manufacturer’s instructions, the total RNA was isolated from MDBK and BEC control and the infected cells using TRIzol LS Reagent (Invitrogen; REF. No. 10296010). RNA concentration was assessed using a Nano-Drop OneC (Thermo Fisher Scientific). cDNA was synthesized from RNAs using a high-capacity reverse transcription kit (Applied Biosystems, Lot. No. 2902953) following the manufacturer’s instructions. The qRT-PCR was performed using Power-Up SYBR Green Master Mix (Applied Biosystems, Lot. No. 2843446) on a QuantStudio3 System (Applied Biosystems). The oligos used for qRT-PCR were designed using the online software Primer3 [31]. The relative gene expression was normalized to the bovine β-actin following 2-^^Ct method [32]. The oligos used in this study are listed in Table 1.

2.8. Statistical Analysis

Results are expressed as means ± SD and were analyzed with the GraphPad Prism v9. A one-way analysis of variance (ANOVA) with Tukey’s or Dunnett’s post hoc test was used to compare multiple groups. P values < 0.05 were considered statistically significant. Statistical significance in figures is indicated as follows: (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001).

3. Results

3.1. The Homology Modelling for the BCoV-Spike and the Potential Host Cell Receptor Proteins

To perform the homology modeling, we searched for the experimentally determined template reference sequences closely related to the query sequence and shared at least 65% sequence identity. The BCoV spike protein sequence (Strain: Mebus) of 1363 amino acids (AA) in length. The sequences of the target molecules from bovine as a receptor are ACE2 (804 AA.), the host cell Furin (797 AA.), the TMPRSS2 (490 AA.), the NRP1 (924 AA), the DPP4 (765 AA.), the APN (965 AA) and the CEACAM-1 (436 AA). Generally, the homology models, with lower DOPE scores and lower PDF total energy values, are considered more structurally stable and reliable. However, if the total energy values for multiple models are very similar, the DOPE score can be used as the deciding factor based on its statistical potential. Therefore, out of five, we selected the two best models with the PDF total energy values and lowest DOPE score for each protein homology model (Table S1). Out of these five models, based on the Lowest DOPE score and PDF total energy, one of the best models that have been further considered for other parameters, such as Ramachandran plot and VERIFY 3D (Table 2), were used to assess the protein 3D model’s structural integrity. The best BCoV spike protein model has a DOPE score of −137057.54, and a PDF total energy value of 51589.75 (model 1) was chosen in the case of the BCoV. The spike protein Ramachandran plot shows that the maximum of the total amino acid residues is in the most favored regions and then is in the additional most allowed areas. A small proportion of residues are in the generously allowed regions. Based on Ramachandran plot results, all the predicted models, including BCoV/S, ACE2, NRP1, CEACAM-1, APN, DPP4, AXL, Furin, TMPRSS2, and the CTS-L, indicate a highly accurate homology model of each structure (Figure 1A–J). The modeled proteins verify that the score is closer to the expected high score, showing that the predicted homology models are of good quality (Table 2).

3.2. The Differential Genes Expression Profiles of Some Host Cell Receptors in BCoV/Ent and BCoV/Resp Infected Cells In-Vitro

The Next-generation sequencing (NGS) was performed to examine the expression profiles of several known coronavirus receptors following infection with either the BCoV/Enteric (BCoV/Ent) or the BCoV/Respiratory (BCoV/Resp) in BEC cells in-vitro. The results revealed that the (NRP1, AXL, CEACAM1, FURIN, and HSPG2) were downregulated in BCoV/Ent infected cells while they were upregulated in cells infected with the BCoV/Resp isolate compared to the sham infected group of cells (Figure 2). The expression level of ACE2 was upregulated in the BCoV/Ent infected group, while no significant change was observed in the expression of ACE2 in the case of the BCoV/Resp infected cells (Figure 2). The expression of the bovine CatapsinL was upregulated in the BCoV/Ent infected cells and downregulated in the BCoV/Resp infected cells (Figure 2). These findings suggest that BCoV isolated have unique differential gene expression profiles in the infected bovine cells, particularly the receptors under study.

3.3. Identifcation of teh BCoV Spikeglycoprotein Structural Domains

In the BCoV-spike-protein sequence, the conserved domain region was identified through the PROSITE and Interpro servers. Based on the PROSITE server and Interpro server results, the BCoV spike glycoprotein has two chains: spike protein S1 and S2. The amino acids from 629 to 1297 contain the S1/S2 cleavage region, the S2 fusion subunit of the spike (S) glycoprotein, and the HR1 and HR2 region from BCoV (Figure 3A). The protease action related to the S1/S2 cleavage region is between 758-797 amino acid residues of the BCoV spike protein structure. The chain S1 contains two domains: an N-terminal galectin-like domain (NTD) and a receptor-binding domain (S1 RBD), also referred to as the RBD site present in between the S1 C-terminal domain (CTD) (Figure 3B). The NTD of S1 is of 15-298 amino acids. However, the S1-RBD site has 343-491 amino acid residues, which are present between the CTD of the S1 protein, and it is 329-617 amino acids long [33].

3.4. Docking of the Bovine Neu5,9Ac2 with the BCoV-S Protein

It was previously reported that the Neu5,9Ac2 is a receptor for BCoV/S-NTD [34]. The BCoV/S-NTD (4h14) shows sequence and structure similarity with the BCoV-isolate BCoV/S-NTD region (Figure S1). To investigate the molecular interaction of the BCoV/S with the bovine cell surface receptor molecule, we performed molecular docking of the BCoV/S-NTD (PDB ID: 4H14) with molecule Neu5,9Ac2. The docking result gives ten poses of ligand-protein interaction complexes. Further, ligand-protein interaction binding energy (BE) calculation was done to calculate the binding energy of complexes. Binding energy (ΔG) determines the binding energy and stability of any given protein-ligand complex or the binding affinity of a ligand to a protein. The binding energy calculation result of the ten best ligand-protein complexes indicated a strong binding affinity between the Neu5,9Ac2 molecule and the BCoV/S-NTD. The interaction between Neu5,9Ac2 and BCoV/S NTD with negative binding energy (-79.55 kcal/mol) is shown as the best ligand-protein interaction pose in Figure 3. This interaction involves some key amino acid residues, including Asp-187, Gly-189, His-185, and Lys-196 (Figure 4A–C). The docking results revealed that Neu5,9Ac2 interacts with a specific pocket above the β-sandwich core in the BCoV/S-NTD.

3.5. The interaction Interface between the BCoV Spike Glycoprotein and the Bovine ACE2.

The energy minimization for the bovine ACE2 protein, the BCoV/S protein, and protein-protein dockings were performed on the ZDOCK server. The ACE2 interacted putatively with the BCoV/S-NTD of the S1 protein based on the homology-modeled BCoV/S protein. The ACE2-BCoV/S protein-protein interaction interface residues are given in Table S2 and Figure 5. The total interaction energy between contact residues of bovine ACE2 and the BCoV/S were computed and reported in units of kcal/mol. The best docking pose shows a ZDOCK interaction score of 34.52 and an E_RDock score of -6.26 kcal/mol. The molecular interaction of the ACE2 shows interaction with BCoV/S at the receptor binding domain (RBD) of the C-terminal domain (CTD) of the S1 protein chain. ACE2 interaction residues are SER19, THR20, THR21, GLU23, GLN24, THR27, GLU30, LYS81, THR 82 and other residues (Table S2). Spike S protein CTD-RBD site residues involved in interaction are TYR661, ASP662, SER663, GLY665, ASN666, SER611, THR425, ARG419, ALA527 and other residues (Table S2). The net interaction between ACE2 and the spike protein is attractive, as indicated by the negative binding energy of their interaction energy (RDOCK interaction energy) (Table 3). Specifically, the receptor-protein interaction analysis showed differences in the key residues at the interface between ACE2 and S1 protein CTD (Figure 5A–C). This might be due to the local sequence difference of BCoV/S protein among different species, suggesting that bovine ACE2 is a putative receptor for BCoV/S due to low binding energy (-6.26 kcal/mol) towards RBD or domain B of BCoV/S protein. QRT-PCR analyzed the expression profile of bovine ACE2 mRNA in the MDBK and BEC cells following the infection with either the BCoV/Ent or the BCoV/Resp isolates. In the case of the MDBK-infested cells, the ACE2 mRNA was upregulated in cells infected with both BCoV isolates. However, a more pronounced ACE2 expression was observed in the BCoV/Resp group (Figure 5D). In BEC cells, the ACE2 mRNA expression was significantly downregulated in both the BCoV/Ent and BCoV/Resp infected groups compared to the sham (Figure 5E). This indicates that ACE2 may play a distinct role in BCoV/Ent and BCoV/Resp infections, depending on the type of the infected cells. Cells.

3.6. The Proposed Model for the Interactions between the Bovine NRP1 and the BCoV/S and the Conformation of the NRP1 Expression Levels in Some Bovine Cell Lines Infected with BCoV

To test the potential roles of the bovine NRP1-as novel receptors for BCoV, the ZDock docking method of protein-protein interaction was performed to identify the interaction sites of the NRP1 with BCoV/S protein. The NRP1 potentially bids with both the NTD and RBD sites of BCoV/S. The interacting interface of amino acid residues for both proteins is given in Table S2 and Figure 5A–C. The NRP1 protein shows binding specificity with NTD and CTD site of BCoV/S with residues Leu322, Leu333, val269, Arg307, Tyr138, Gln 266, and other amino acid residues (Table S2). The residues from BCoV/S-NTD that interacted with NRP1 are Thr134, Leu194, Arg197, Phe232, Leu238, and other amino acid residues given in Table S2. The protein-protein docking between BCoV/S with the bovine NRP1 shows greater interaction, as suggested by the E_RDock score (-22.99 kcal/mol). The interaction affinity between both proteins is due to multiple non-covalent interactions (15 conventional hydrogen and 6 Pi bonds) (Figure 6A–C and Table 3). The binding of the NRP with the BCoV/S glycoprotein is highly stable, as suggested by the negative E_RDOCK binding energy values (-22.99 kcal/mol). Thus, the bovine NRP1 exhibits a higher binding affinity to the NTD and CTD or domain B of the BCoV/S. The in-vitro experiments revealed that the mRNA expression of bovine NRP1 in MDBK cells was significantly upregulated in both the BCoV/Ent and BCoV/Resp infected groups compared to the sham (Figure 6D). Similarly, NRP1 mRNA expression was significantly increased in BEC cells in both infected groups compared to the sham (Figure 6E). Notably, NRP1 expression was higher in both MDBK and BEC cell lines in the BCoV/Resp infected group.

3.7. The Proposed Model of the Interaction between the BCoV/S Proteins and the Bovine Type-I Membrane Protein Receptor Carcinoembryonic Antigen-Related Cell Adhesion Molecule 1 (CEACAM-1) and the Confirmation of the Bovine CECAM1 during BCoV Infection in Bovine Cell Lines

To test the potential roles of the bovine CEACAM-1 protein as a receptor for the BCoV, the ZDock docking method of protein-protein interaction was performed to map the interaction sites of the CEACAM-1 with BCoV/S glycoprotein. Our results show that the bovine CEACAM-1 protein binds with the BCoV/S-NTD. The interacting interface of amino acid residues for both proteins is shown in (Figure 7A–C and Table S2). The bovine CEACAM-1 protein shows high binding specificity to the BCoV/S-NTD site of Spike glycoprotein mainly through the residues (ARG143, ASN170, THR171, ASN178, LYS196, LYS196, and other amino acid residues given in Table S2). The residues of the BCoV/S-NTD that interacted with CEACAM-1 are (LEU233, PHE232, VAL229, Phe232, THR145, ASP148, and some other amino acid residues shown in Table S2). The best docking pose shows a ZDOCK interaction score of 16.8, and E_R Dock interaction energy is -6.29 kcal/mol. These parameters are mainly due to multiple non-covalent interactions, including nine conventional hydrogen bonds, 12 Pi bonds, and one salt bridge interaction between the two proteins (Figure 7A–C and Table 3). Based on the results of the E_R Dock Score of the interaction energy of the best pose, the low binding affinity of BCOV/S glycoprotein with CEACAM-1 was observed compared to the interaction of BCoV/S-NTD known receptors Neu5,9Ac2. The in-vitro experiments demonstrated that the bovine CEACAM-1 expression level in the MDBK cells was downregulated in the case of cells infected with both BCoV-infected groups. Marked prominent downregulation was observed in the BCoV/Ent infected group (Figure 7D). In the case of the BEC cells infected with BCoV, the bovine CEACAM-1 expression levels were significantly downregulated in the case of BCoV/Ent infected cells, while its expression levels were significantly upregulated in the BCoV/Resp group compared to the sham (Figure 7E).

3.8. The Proposed Model for the BCoV/S Glycoprotein Interaction with the Bovine Aminopeptidase N (APN) and the Conformation of the Expression Levels of the APN Genes during BCoV Infection in Some Bovine Cell Lines

To explore the potential roles of the bovine APN as a receptor for the BCoV, the ZDock with RDOCK docking method of the protein-protein interaction was performed to map the BCoV/S protein interaction sites with the bovine APN protein. The interaction between the two proteins suggested that the APN showed binding with the RBD of the BCoV/S glycoprotein. The interacting interface amino acid residues for both proteins are shown in Figure 8A–C. The APN protein shows a binding specificity with the BCOV/RBD, particularly the BCoV/S-CTD site, through the mapped amino acid residues (ASN397, TYR451, ARG514, LYS579, ASN583, CYS524, and with some other amino acid residues listed in Table S2). The residues of BCoV/ S1-CTD interacted with the APN are CYS826, ASP849, GLU843, SER839, GLU843, SER839, ILE512, GLY511, PRO516, ASN495 and other amino acid residues given in Table S2. The best docking pose with low-energy binding conformation shows a ZDOCK interaction score of 17.6 and an E_RDock score of -2.96 kcal/mol. The interaction between proteins is due to multiple non-covalent bonds, including 19 conventional hydrogen bonds, two Pi bonds, and one salt bridge interaction between the two proteins (Table 3). The in-vitro experiments demonstrated that the gene expression levels of the bovine APN in the MDBK cells were downregulated in the case of the BCoV/Ent-infected cells. At the same time, no notable change was observed in the BCoV/Resp infected group (Figure 8D). Conversely, the BEC cells’ bovine APN gene expression levels were significantly upregulated in both BCoV-infected groups. Notably, a more obvious upregulation was observed in the BCoV/Resp-infected cells than in the sham-infected cells (Figure 8E).

3.9. The Proposed Model for the Interaction between the BCoV/S Glycoprotein with the Bovine Dipeptidyl Peptidase 4 (DPP4) and the Confirmation of the Expression Levels of the DPP4 during BCoV Infection in Some Bovine Cell Lines

We used the ZDock docking method to map the interaction between some residues of the bovine DPP4 with BCoV/S. The bovine DPP4 protein was mapped to bind to the C-terminal domain (CTD) of the BCoV/S protein. The interacting interface of amino acid residues for both proteins are shown in (Table S2 and Figure 9A–C). The DPP4 protein showed binding specificity with the BCOV/RBD, especially in the BCoV/S-CTD site with the mapped residues (TRP184, HIS185, TRP186, LYS196, TRP401, LYS422, and some other interacted amino acid residues presented in (Table S2). The residues of the BCoV/S-CTD interacted with the bovine DPP4 protein are (ASP451, ASP500, LEU448, ASN146, THR144, GLU182, HIS185, THR188, GLY189 and other interacted amino acid residues given in the Table S2. The best docking pose shows a ZDOCK interaction score 22.44 and an E_R Dock score of -6.23 kcal/mol. These parameters might be due to multiple non-covalent interactions, including 20 conventional hydrogen bonds, nine Pi bonds, and one salt bridge interaction between the two proteins (Table 3). The gene expression levels of the bovine DPP4 in MDBK cells were significantly upregulated in the case of the BCoV/Resp infected cells, while no significant changes were observed in the BCoV/Ent infected cells compared to the sham (Figure 9D). In the case of the BEC cells, bovine DPP4 gene expression levels were significantly upregulated in the case of the BCoV/Ent infected cells, while no significant changes were observed in the BCoV/Resp infected cells compared to the sham (Figure 9E).

3.10. The Proposed Model of the Interaction between the BCoV/S Glycoprotein with the Bovine AXL Protein and the Confirmation of the Bovine AXL Protein Expression Levels in Some Bovine Cell Lines Infected with BCoV Isolates

The in silico docking procedure reveals the potential for very strong interactions between the bovine AXL protein and the BCoV/S glycoproteins. The bovine AXL protein interacted with the NTD of the BCoV/S protein. The interacting interface of the amino acid residues from both proteins is shown in (Table S2 and Figure 10A–C). The AXL protein showed strong binding affinities with the NTD in the BCoV/S with the residues (LYS524, SER583, LYS665, TYR697, ASN698, ARG703, ILE707, GLU579, CYS587, LYS694, LYS695, ARG706). The residues of the BCoV/S1-NTD interacted with the bovine AXL protein are (TYR237, VAL191, ASN198, CYS160, HIS185, GLU182, VAL181, THR188, TYR136, GLN141, GLN161, CYS193, LYS196, GLY189, GLU156 and GLY239). The best docking pose shows a ZDOCK interaction score 17.56 and an E_R Dock score of -13.68 kcal/mol. The AXL shows high binding/interaction energy with BCoV/S due to multiple non-covalent interactions, including 21 hydrogen bonds and nine Pi bonds between the two proteins (Table 3). The in-vitro experiments revealed that the gene expression of the bovine AXL in the MDBK cells was significantly upregulated in the BCoV/Resp infected cells, while there were no obvious changes observed in the AXL in the case of the BCoV/Ent infected cells compared to the sham infected cells (Figure 10D). Similarly, the bovine AXL gene expression in the BEC cells was significantly upregulated in the BCoV/Resp infected groups. At the same time, no significant changes were observed in the BCoV/Ent-infected cells compared to the sham-infected cells (Figure 10E). These findings strongly suggest that the bovine AXL is specifically upregulated during BCoV/Resp infection in both cell lines, highlighting its potential pivotal roles in the context of the BCoV/Resp isolates.

3.11. The Proposed Model of the Interaction between the Heparan Sulfate (HSPG2) and the BCoV/S Glycoprotein and the Confirmation of the HSPG2 Expression during BCoV Infection in Some Bovine Cell Lines

The docking of the heparan sulfate (HS) of the proteoglycan protein (HSPG2) interacted with BCoV/S-NTD (PDB ID: 4h14). The docking result generated ten poses of ligand-protein interaction complexes. The best pose shows the interaction as indicated by the strong binding affinity between the HS molecule and the BCoV/S, with the highest binding energy (-243.27 Kcal/mol). The HS interaction involves some key amino acid residues from BCoV/S, including (TRP184, HIS143, GLU156, ASN198, ARG197 and LYS192) (Figure 11A–C). Therefore, the docking (molecular interaction) result for the BCoV/S with the HS glycan of the HSPG2 receptor from the bovine cell surface showed high binding affinity towards the BCoV/S protein. The gene expression level of the bovine HSPG2 in the BCoV-infected MDBK cells was significantly downregulated in the case of the BCoV/Ent-infected cells, while there was no significant change in the expression level of the HSPG2 observed in the BCoV/Resp infected cells compared to the sham (Figure 11D). Similarly, in the BCoV-infected BEC cells, the bovine HSPG2 gene expression was downregulated in the case of the BCoV/Ent-infected cells, while there were no significant changes observed in the BCoV/Resp infected cells compared to the sham (Figure 11E). The findings indicate that bovine HSPG2 is specifically downregulated during BCoV/Ent infection in both cell lines, highlighting its potential role as a novel receptor for BCoV.

3.12. The Proposed Model for the Putative Interaction between The bovine Host Cell Furin Enzyme and the BCoV/S Glycoprotein and Mapping the Putative Furin Cleavage Sites in the BCoV-S Glycoproteins

The online ProP 1.0 Server software prediction result showed one pro-peptide cleavage site (Arg764, Arg765, Ser766, Arg767, Arg768, and Ala 769 (RRSRR|A)) amino acids present across the full-length genome sequence of the BCoV/S (Uniprot ID: P15777) (Figure S2). The multiple sequence alignments of the BCoV/S (Mebus) sequence with other betacoronaviruses spike protein sequences mapped the protease-specific residues at (Arg764, Arg765, Ser766, Arg767, Arg768, and Ala 769 (RRSRR/A)) (Figure 12A and Figure S3). This finding strongly suggests that the host bovine cell furin enzyme could bind to the BCoV-S glycoprotein and cleave this protein at these polybasic amino acid residues. The protein-protein interaction of the bovine host cell Furin and the BCoV/S glycoprotein suggested the potential binding of Furin at the cleavage site identified by the ZDOCK Score 16.72 and E_RDock score -13.14 kcal/mol. Among the ten refined poses from E_RDock, the best pose of the Furin binding to the S1/S2 region mapped near the protease cleaving site of polybasic amino acid residues (Arg764, Arg765, Ser766, Arg767, Arg768, and Ala 769 (RRSRR|A) at BCoV S1/S2 cleavage sites (Figure 12A). The protease Furin interacted with spike protein at the amino acid residues (Glu230, Leu227, Asp228, Asp174, Asp177, Asn192, Met189), to the spike protein S1/S2 cleavage RRSRR site amino acid residues (Val758, Asp759, Tyr760 Ser761, Thr762 Lys763, Arg764, Arg765, Ser766, and Arg767) (Figure 12A and Figure 13A–C). The protein-protein interaction interface of BCoV/S-Furin residues is shown in (Table S3). The protein-protein interface interacted residues, Lys736 from BCoV/S, form an electrostatic interaction, forming a salt bridge interaction with ASP174, ASP177, and Asp228 residues of the Furin. The Arg767 and Ser766 of the BCoV/S RRSRR cleavage site show an electrostatic interaction and hydrogen bonding to the Glu230 of the Furin (Table 3 and Table S3). The gene expression levels of the bovine Furin in the MDBK cells were upregulated in both BCoV isolates infected cells, with a significant upregulation observed in the case of BCoV/Ent infected cells (Figure 13D). Similarly, the bovine Furin gene expression in the BEC cells was significantly upregulated in both infected groups (Figure 13E). These results indicate that bovine Furin expression is upregulated during BCoV infection, regardless of the viral isolate.

3.13. The Putative Model for the TMPRSS2 Docking with the BCoV/S Glycoprotein and the Conformation of the TMPRSS2 Expression Levels during BCoV Infection in Some Bovine Cell Lines

The MSA of the BCoV/S (Mebus) strain sequences with other Betacoronaviruses spike protein sequences suggested the interactive protease-specific residues are (Arg764, Arg765, Ser766, Arg767, Arg768, and Ala 769 (RRSRR/A)) (Figure 12B). The protein-protein interaction takes place between the bovine TMPRSS2 amino acid residues (Val278, His294, Thr339, Asp415, Asn416, Trp459, Gly460, Ser461, Gly462) near the RRSRR site of BCoV/S protein residues (Asp738, Ser740, Thr741, Ser742, Ser743, Ser761, Arg764, Thr762, Lys763, Arg767, Thr771 and Ile779) (Table S3). The results of the protein-protein docking between the BCoV/S/ and TMPRSS2 using the ZDOCKER show the (E_RDock energy by -0.62 kcal/mol) due to the presence of some multiple non-covalent interactions (conventional hydrogen and Pi bond) between the two proteins (Table 3). The best pose of the TMPRSS2 binding at the active sites (substrate binding site and catalytic site) of the BCoV/S2 protein is shown in (Figure 14A–C). The expression levels of the bovine TMPRSS2 gene in the case of the BCoV-infected MDBK cells were significantly upregulated in the BCoV/Ent-infected cells, while no changes were observed in the BCoV/Resp-infected cells compared to the sham (Figure 14D). Similarly, in the case of the BCoV-infected BEC cells, the bovine TMPRSS2 gene expression levels were significantly upregulated in the case of the BCoV/Ent-infected cells, while there were no significant changes observed in the case of the BCoV/Resp infected- cells compared to the sham (Figure 14E). These findings strongly suggest that the bovine TMPRSS2 expression is upregulated, particularly during BCoV/Ent infection in the MDBK and the BEC.

3.14. A Proposed Putative Model for the Cathepsin-L (CTS-L) Docking with the BCoV/S and the Confirmation of the Cathepsin-L Expression Levels during BCoV Infection in Bovine Cells

The MSA of the BCoV/S from the Mebus strain with other Betacoronaviruses, particularly SARS-CoV-2, suggested no conserved cathepsin-L cleavage site (CS-1 and CS-2) or region present in the S chain of spike protein as SARS-CoV-2 (Figure S4A,B). The protein-protein interaction of the bovine host cell CTS-L and the BCoV/S glycoprotein suggested the potential binding site identified by the ZDOCK Score (19.84 and E_RDock score -16.45 kcal/mol) (Table 3). Among the ten refined poses from the E_RDock, the best pose is the CTS-L binding to the residues of BCoV-NTD and the downstream of the BCoV/S1/S2 region of polybasic amino acid residues (Figure 15A–C). The BCoV/S amino acid residues in this interaction through hydrogen bonding are (Thr38, Ala33, Ser35, Ile36, Gl72, Arg768, Pro723, Gln722, Phe713) (Table S3). The Pi-interaction residues from the BCoV/S side are (Ile36 and Met291). The protease CTS-L interacted with spike protein through amino acid residues (Lys267, Thr275, Asp326, Gly190, Phe186, Tyr159, Asp326, Glu142, Thr143 and Pro149) (Table S3). The Pi-interaction residues from CTS-L side are (Trp154 and Phe145. The Arg786) residues from the polybasic RRSRR site are present at the BCoV/S1/S2 junction of the BCoV/S, forming hydrogen bonding with (Lys276) residue of the bovine CTS-L. The bovine CTS-L gene expression levels in the MDBK cells were significantly downregulated in the case of the BCoV/Ent infected cells, while there were no changes observed in the BCoV/Resp infected cells compared to the sham-infected cells (Figure 15D). In the case of the BCoV-infected BEC cells, bovine CTS-L gene expression was significantly downregulated in the case of the BCoV/Resp-infected cells, while there were no changes observed in the case of the BCoV/Ent infected cells compared to the sham-infected cells (Figure 15E).

3.15. The Proposed Model of the Neu5,9Ac2 Docking with the BCoV Haemagglutinin Esterase (BCoV-HE) Protein

The X-ray crystal structure of the BCoV/HE complex with the 4, 9-O′-diacetyl sialic acid retrieved from the RCSB PDB database (PDB ID:3CL5; 1.80 Å resolution) and was chosen as the receptor for ligand Neu5,9Ac2 [35]. The docking of Neu5,9Ac2 with HE indicated a strong binding affinity between the sugar molecule and the BCoV/HE, with the highest binding energy (-122.73 Kcal/mol). The docking of Neu5,9Ac2 results suggested its interaction with a specific pocket where 4,9-O′-diacetyl sialic acid interacts with the HE. This interaction involves critical amino acid residues, including Asn-264, Glu-265, Ser-221, and Asn-236 (Figure 16A–C). Our docking (molecular interaction) result of the BCoV/S and the BCoV/HE proteins with the known BCoV receptors, Neu5,9Ac2, showed high binding affinity towards both the proteins of BCoV. However, based on binding energy, HE interaction with receptor Neu5,9Ac2 (-122.73Kcal/mol) is more stable and more significant compared to the interaction of the BCoV/S protein with the Neu5,9Ac2 (-79.55 kcal/mol). This model suggests that HE may help in the interaction of BCoV with or without spike protein to bind to its specific cellular receptor Neu5,9Ac2. The BCoV-HE could interact with its cellular receptors and may cause changes to the cell surface, which could lead to infection of the host cells.

4. Discussion

Bovine coronavirus (BCoV) is one of the most common pathogens affecting cattle of all ages. BCoV possesses dual tissue tropism in cattle, particularly in the respiratory and enteric tracts. The process of BCoV entry to the target cells is a crucial step in the replication cycle of the virus. This process requires the orchestration between several viral and host cell proteins. The coronavirus spike glycoproteins contributed substantially to other viral infections and the molecular pathogenesis of coronaviruses, including BCoV. Meanwhile, among all coronaviruses, the hemagglutinin esterase enzyme is found in only BCoV and the human coronavirus-OC43 (HCoV-OC43). The viral receptors are among the key players in the viral entry and could fine-tune the viral tissue tropism. Some host cell proteases also enhance the coronavirus replication from entry to the release of the virus from the host cells. Previous studies have shown that the 9- O-acetylated sialic acids could potentially act as receptors for the BCoV. The tropism of coronaviruses is a complicated process that requires the availability of some factors from the viral side, including some attachment proteins, particularly the spike glycoprotein. This process also involves some cellular factors, including the receptors and other transcription and translation factors. It also requires some factors from the infected host, particularly the availability of host enzymes that help activate some essential proteins [36,37,38]. BCoV possesses a dual tissue tropism in cattle. The virus mainly affects the digestive and respiratory tract of the affected animals; this pattern is called pneumoenteritis [39]. The viral tropism primarily depends on the availability of specific viral receptors, some other transcription translation factors, and some host cell enzymes [40]. Members of the family coronaviridae utilize many host cell receptors to attach to their target cells [41].

Most coronaviruses require activation through the cleavage of their spike glycoprotein by some host cell proteases to facilitate the viral infection [32,33,34,42]. Although BCoV was discovered long ago, little is still known about viral tropism, especially the roles of the host cell receptors and the enzymes in fine-tuning the viral tissue tropism, pathogenesis, and replication [35]. The current study aims to identify some novel BCoV receptors with higher affinities that bind to the S and HE proteins to initiate and facilitate viral infection and pathogenesis. The BCoV/S glycoprotein is the main viral protein involved in the process of viral attachment to the host cells [43]. The BCoV/S protein consists of two subunits (S1 and S2): the N-terminal domain (NTD) and BCoV-S1-CTD are in the S1 subunit, whereas the fusion peptide (FP) and heptad repeat (HR) domains 1 and 2 are located in the S2 subunit of the BCoV/S glycoprotein. The BCoV/S protein usually attaches to the cell membrane by interacting with viral receptors on the surface of the target cells, initiating the viral infection. Spike protein S2 mediates the fusion of the virion and cellular membranes by acting as a class I viral fusion protein. Also, it acts as a viral fusion peptide, which is unmasked following the S2 cleavage site occurring upon virus endocytosis [44,45]. The distal S1 subunit of the coronavirus spike protein is responsible for receptor binding. Either the S1-NTD or the S1-RBD at the C-terminal domain of the BCoV-S1 protein chain, or occasionally both, are involved in the binding to the host receptors [46,47]. The observed binding between Neu5,9Ac2 and the BCoV/S-NTD is particularly interesting because BCoV haemagglutinin esterase (HE) also utilizes this O-acylated sugar molecule as a substrate. While our study identified interacting residues like Asp-187, Gly-189, His-185, and Lys-196, another report [34] proposed a slightly similar set of critical residues for Neu5,9Ac2 binding (Tyr-162, Glu-182, Trp-184, and His-185).

BCoV/S and BCoV/HE proteins act synergistically and harmoniously to orchestrate the BCoV infection in the target cell [7,48]. The BCoV/S is mainly involved in the initial attachment of the virus to the host cells, while the BCoV/HE destroys the sialic acid in the cell’s surface, promoting the viral release from the host cells [7].

Most of the Betacoronaviruses, including the BCoV isolates, use the 9-O-acetyl-SA as receptors; however, during the evolution of these viruses, some isolates started to recognize other forms of sialic acid (the 4-O-acetyl-SA isoform) [49]. It was recently shown that the HE gene of the SARS-CoV-2 recognizes and binds to the 9-O-acetyl-SA in contrast to the type-II HE protein of SARS-CoV-2 uses the other isoform; the 4-O-acetyl-SA [49]. The ligand-interaction sites of the BCoV HE and the subset of coronavirus S glycoproteins evolved to recognize 9-O-acetyl-SA via hydrogen bonding [50].

In the case of most Betacoronaviruses, the RBD of the CTD or domain of the spike glycoproteins showed the highest variability within S1 subunits across various members of the coronaviruses, including Betacoronaviruses. This phenomenon allows coronaviruses to bind to various types of host cell receptors [46]. It has been proved that SARS-CoV-2 uses the ACE-2 as a valid receptor for the viral entry into the target host cells [51]. However, there are no records about the potential roles of ACE2 as receptors for BCoV. In our study, based on the interaction between bovine ACE2 and BCoV /S interaction with its CTD, the bovine cell surface ACE2 could act as a putative receptor for BCoV/S (Figure 17). This finding is based on the BCoV/S binding to the ACE2 receptor’s tendency to show lower interaction energy [52]. This claim is also supported by the nature of the interaction energy, which is consistent with the structure of the virus receptor interface [53]. Additionally, we demonstrated that the expression level of the bovine ACE2 is highly expressed in the MDBK following BCoV infection, whereas the ACE2 gene expression is downregulated in bovine endothelial cells (BEC) under the same infection conditions (Figure 5D,E).

We also investigated the potential roles of Neuropilin-1 (NRP1) as a receptor for the BCoV. Our ZDock method result suggested the high interaction affinity (-22.99 kcal/mol) of NRP1 with BCoV/S-NTD and BCoV/S-CTD region (Figure 5 and Figure 16, and Table 2). Furthermore, the gene expression studies in some bovine cell lines revealed a significant upregulation of the expression levels of the NRP1 following BCoV infection (Figure 6D,E). This raises the possibility that BCoV could potentially use the NRP receptors for viral entry into the host cell. According to some previous reports, the higher expression of NRP1 will facilitate virus-host cell interactions, especially in cells that do not express other potential BCoV receptors [54,55]. Previous research indicates that Neuropilin-1 (NRP-1) exhibits stronger binding affinity to the CTD region of S1, and this interaction stabilizes the folded conformation of the S protein [56]. Conversely, in the absence of NRP-1, the SARS-CoV-2 S protein tends to stretch and unfold.

It is well known that MHV utilizes the CEACAM1 protein for a dual function. The CEACAM-1 acts as a functional receptor for MHV in addition to enhancing the activation of the MHV-S glycoprotein by inducing some conformational changes, allowing the fusion of the virus with the target cells [57]. However, there is no data about the potential roles of CEACAM1 in the BCoV replication. Our in-silico data showing low binding affinity interaction between BCoV/S and CEACAM1 suggested that CEACAM1 interacted with low affinity to the NTD of BCoV/S protein (Figure 17). The gene expression analysis in BCoV-infected bovine cells revealed downregulation of CEACAM1 gene expression levels in the case of the BCov/Ent isolate infected cells, whereas a significant upregulation in the CEACAM1 gene expression levels was observed in the case of the BCoV/Resp infected cells (Figure 2 and Figure 7D,E). These findings indicate that CEACAM1 activation is specific to the BCoV/Resp isolate infection. The variations in the binding of coronaviruses to CEACAM1 can be directly attributed to the structural differences between their N-terminal domains, particularly those of BCoV and MHV. Unlike MHVs sugar-binding ancestors, contemporary MHVs utilize their NTDs for specific interaction with the host cell protein CEACAM1 [58].

The aminopeptidase N (APN), a widely found enzyme on cell surfaces, is involved in various cellular processes like survival, migration, blood pressure regulation, and even virus uptake [59]. Coronaviruses like the porcine delta coronavirus (PDCoV) can utilize APN as a receptor to enter cells, highlighting APN’s potential role as a viral infection gateway from different host species [60]. Aminopeptidase N from porcine functions as a receptor for the enveloped RNA virus TGEV [12]. This underscores the wide variety of membrane-bound proteins viruses exploit to infiltrate cells. The APN is a common receptor for the members of the alpha coronaviruses, particularly the HCoV-229E and the TGEV in pigs [61,62,63]. Based on our protein-protein interaction (E_RDock), the interaction energy of the best pose of APN interaction with BCoV/S shows low-affinity interaction (-2.96 kcal/mol). This result confirms that the binding interaction of BCOV/S is not strong, and it shows a putative binding with the APN. The RBD of CTD from BCoV/S is involved in interaction with Bovine APN (Figure 17). Meanwhile, the host APN gene was not expressed in MDBK cells following BCoV infection, but the limited expression of the APN in the BCoV infected BEC cells (Figure 8D,E). This pattern of APN expression suggesting a subtile expression of APN in response to BCoV infection in the BEC.

The DPP4 is a cell surface protease that exhibits exopeptidase activity and is expressed on the surface of various cell types, including those found in the human airways [64]. Researchers have compared the potential binding interactions between the receptor binding domain (RBD) of different spike variants of SARS-CoV-2 and DPP4 with the interactions observed in the experimentally determined structure of the MERS-CoV complex with DPP4 [65]. Members of the Betacoronaviruses, especially the SARS-CoV-2 and the MERS-CoV, utilize the ACE2 and DPP4 as receptors, respectively [66]. The pathogenicity of MERS-CoV is caused by the specific binding of its S1B domain or S1-CTD to the human DPP4 receptor [67]. In our study, docking of the BCoV spike with the DPP4 shows low interaction affinity of the DPP4 towards the BCoV/S (E_R Dock score -6.23 kcal/mol). However, the binding of DPP4 with BCoV/S1-CTD region or S1B domain is correlated with MERS-CoV binding of its S1B domain or S1-CTD to the human DPP4 (Figure 17). Additionally, the DPP4 gene expression level was noticed to be specific to the BCoV/Resp infected MDBK cells, while in the case of the BEC cells, it was only expressed in the BCoV/Ent infected cells (Figure 9D,E). This indicates an inconsistent pattern of bovine DPP4 expression in response to BCoV infection. This finding suggests that DPP4 is not a strong candidate receptor for the BCoV.

The AXL is a member of the TAM tyrosine kinase receptor family and appears to be the most efficacious at mediating viral entry, given their prevalent use among enveloped viruses [68]. The AXL receptor is found on the human cell membrane and contributes to the pathogenesis of SARS-CoV-2 infection [69]. The AXL has been reported to bind to the NTD of the spike protein S1 subunit and could synergistically work with the ACE2 to facilitate SARS-CoV-2 entry into the target cells [70]. In our study, the NTD of the BCoV/S is predicted to be involved in the interaction with Bovine AXL (Figure 17). The AXL shows high binding affinity with NTD of the BCoV/S due to multiple non-covalent interactions, suggesting the possibility of its role in infection with its potential synergistic effect with other receptors such as ACE2 as shown in the context of the SARS-CoV-2 infection [42]. The AXL gene expression profiles in the bovine cell lines infected with BCoV support our in-silico prediction and modeling. However, it indicates that the AXL gene expression level is upregulated specifically during the infection with the BCoV/Resp isolate (Figure 2 and Figure 10D,E). Therefore, this AXL could be a strong candidate receptor, especially in the case of BCoV/Resp infection.

The Heparan sulfate proteoglycans of the HSPG2 play important roles during the entry of many viruses to the target cells [71]. Glycan-binding may act as the initial step for cellular attachment, bringing the virus close to the epithelial cell membrane where it can interact with its protein receptor(s). The general mechanism appears conserved among coronaviruses (CoVs) such as MERS-CoV, SARS-CoV-1, and HCoV-NL63, which are related to SARS-CoV-2 [71,72]. The spike protein initially interacts with cell surface HS, resulting in an increase in the expression levels of the RBDs in the “up/open” conformation, which in turn leads to enhanced binding affinities to the bovine ACE2 receptors and subsequent stabilization of the interaction of the trimer and the virion with cells [73]. The HS captures SARS-CoV-2 virions, increasing the local viral concentration and presentation of these virions to the ACE2 receptores and to the other proteinaceous receptors as well [74]. Our in silico and molecualr docking results are showing the potential roles of the BCoV/S/NTD in the interaction with the HS of the HSPG2 (Figure 17). The HS shows high binding energy with NTD of BCoV/S due to multiple non-covalent interactions such as hydrogen and pi-bonds, suggesting the possibility of its role in infection with its synergistic effect with other receptors such as ACE2. The BCoV/S-HS interaction could interact with its other cellular receptors, such as ACE2, may cause changes to the cell surface, which could lead to infection of the host cells. However, our results showed the downregulation of the HSPG2 gene expression, especially during the BCoV/Ent infection (Figure 2 and Figure 11D,E). Therefore, HS binding inhibition could also be a novel antiviral therapeutic for inhibiting BCoV infection in cattle.

Previous studies on SARS-CoV-2, MERS-CoV, and other coronaviruses have shown that activation of the spike protein is often a complex process involving multiple cleavage events occurring at distinct sites and with the involvement of several host proteases [75,76]. For virus entry and infections, proteolytic cleavage is widely used to activate the fusion machinery of viral glycoproteins. Furin binding to RRSRR (Arginine-Arginine-Serine-Arginine-Arginine) site of BCoV spike protein, which is related to the conserved cleavage site of SARS-CoV-2 cleavable PRRAR|S residues at receptor-binding (S1) and fusion (S2) domains of the spike protein [77]. Such a motif may allow Spikes to be cut into S1 and S2 by Furin and TMPRSS2-like proteases before maturity, which provides S1 with the flexibility to change the conformation to better fit the host receptor. The arginine residues at the SARS-CoV-2 spike protein catalytic site are popped out of the closed state of the S protein to form multiple non-covalent interactions with the furin [78]. In our study, the BCoV/S interaction with the furin shows that furin binds near the spike protein S1/S2 cleavage RRSRR site specific for Furin and TMPRSS2 proteases (Figure 17). Our gene expression data strongly supports the in-silico prediction of host cell Furin targeting the BCoV-S glycoprotein; it is also shows that host Furin is highly expressed during BCoV infection, with particular upregulation of the Furin expression observed in response to BCoV/Resp isolate infection in bovine cells (Figure 2 and Figure 13D,E). Our recent study shows the host cell miRNA-16a targeting the host cell furin, limiting the BCoV replication and enhancing the host immune response [79]. The TMPRSS2 has been shown to proteolytically activate the S glycoprotein of many coronaviruses, including SARS-CoV-2, SARS-CoV-1, MERS-CoV, 229E, as well as influenza virus [80,81,82]. The TMPRSS2 triggers HKU1-mediated cell-cell fusion and viral entry and binds with high affinity to both HKU1A and HKU1B RBDs [83].

Our modeling results show that furin binds firmly with the S protein RRSRR site, which is conserved in BCoV as RRSRR at the S1/S2 site, whereas the TMPRSS2 cleaves the S protein in the lungs of SARS-CoV-2 infected person and promotes pathogenicity [84]. SARS-CoV-2-S glycoprotein (S1/S2) is cleaved by furin protease, and subsequently, the TMPRSS2 mediates the cleavage and activation of the S2 region on the S2 protein [85]. The TMPRSS2, which colocalizes with the ACE2 at the cell membrane, has been identified as the dominant proteolytic driver of the SARS-CoV-2/S protein activation that potentiates the viral infection in the aerodigestive tract [86]. The TMPRSS2 protease can be opted as a potential therapeutic target for bovine host-specific viruses. In this study, TMPRSS2 binds near the S1/S2 junction RRSRR protease cleavage-specific region for furin and TMPRSS2 proteases. The in-vitro results also demonstrated high expression levels of the TMPRSS2especaiaaly during the BCoV/Ent isolate infection in the bovine cells (Figure 14D,E). This finding strongly suggests the involvement of the TMPRSS2 in the entry of BCoV enteric isolate into the host cell.

The CTS-L is a member of the lysosomal cysteine protease family [17]. The expression level of CTS-L mRNA is higher than ACE2, FURIN, and TMPRSS2 in human lung tissues during SARS-CoV-2 replication [87]. In the context of SARS-COV-2, the S protein contains two highly conserved CTS-L cleavage sites (CS-1 and CS-2) among all known SARS-CoV-2 variants. CTS-L cleavage promoted S to adopt receptor-binding domain (RBD) “up” activated conformations, facilitating receptor-binding and membrane fusion [17]. The highly conserved CS-2 site seems to be the most essential for the life cycle of SARS-CoVs, while CS-1 is likely to play an auxiliary role in virus infection. However, with the help of multiple sequence alignment, it has been observed that the BCoV/S has no conserved region for CTS-L cleavage site as the SARS-CoV-2 sequence of the S chain. On the basis of our protein-protein interaction result, cathepsin showed binding to BCoV/S in NTD and downstream of the S1/S2 region of CTD of the S chain (Figure 17). However, our gene expression analysis of the host genes during BCoV infection in some bovine cells showed the down-regulation of the CST-L expression upon BCoV infection (Figure 15D,E). To understand the catalytic role of bovine cathepsin towards BCoV/S glycoprotein, the interaction of BCoV/S with cathepsin could be further investigated in vitro.

Therefore, there is limited knowledge about the bovine-specific virus-host interactions that determine cellular entry of SARS-CoV-2. Viruses display considerable redundancy and flexibility because they can exploit weak multivalent interactions to enhance affinity. To date, studies of BCoV entry have focused almost entirely on ACE2 and Neu5,9Ac2. In this study, we aimed to see the interaction of bovine cell surface molecules as receptors for BCoV/S, including Nue5,9Ac2, NRP1, DPP4, CEACAM-1, APN ACE2, AXL, HS of HSPG2, TMPRSS2, and Furin shown in Figure 17. The results from our study suggested that NRP1 and AXL showed greater binding affinity towards BCoV/S other than receptor Neu5,9Ac2 and HS. However, bovine HE also showed a stronger binding affinity with Neu5,9Ac2. Our in silico structural interaction data provide a blueprint for understanding the BCoV specificity for the different bovine receptors.

5. Conclusions

Our proposed in silico modeling and docking show the potential interaction between the BCoV/S protein and some known coronaviruses receptors (ACE2, DPP4, NPR1, APN, AXL, HS, and CEACAM-1) at various degrees. The NRP1-BCoV/S interaction and the AXL- BCoV/S interaction show the highest parameters of refined dock score (E_RDock score). The AXL interacts with BCoV through the NTD, while the NRP1 interacts with BCoV/S through both ends (NTD and CTD). Therefore, based on the docking result from our study, both the NRP-1, and the AXL might have some potential roles of co-receptors for BCoV entry into the host cell. Consistently, our gene expression analysis of the BCoV infection in some bovine host cells revealed that NRP-1 is significantly upregulated during BCoV/Ent and BCoV/Resp infection, whereas AXL was upregulated only during BCoV/Resp infection. However, regarding cell proteases, our data confirms that the host cellular proteases (Furin and TMPRSS2) have conserved recognition sites within the BCoV/S glycoprotein. The in-vitro experiments also showed high expression of Furin during BCoV/Ent and BCoV/Resp infection, whereas TMPRSS2 was upregulated only during BCoV/Ent infection. As SARS-CoV-2 CTS-L cleavage sites, the multiple sequence alignment of BCoV/S showed the absence of any CTS-L-specific conserved region for its catalytic activity. CTS-L showed binding affinity towards NTD of the S1 chain. The interaction of BCoV-HE shows high-affinity interaction towards the Neu5,9Ac2 receptor from bovine. The BCoV/HE protein may augment the interaction between the BCoV and Neu5,9Ac2, enhancing the virus’s attachment to the cellular receptors and the downstream BCoV replication steps. Our in-silico study predicted that the ACE2, NRP1, and AXL receptors and protease CTS-L from bovines would be potential receptors and co-receptors for bovine coronaviruses. Our results show that the use of ACE2, NRP1, AXL, and CTS-L inhibitors may provide anti-BCoV therapeutic activity. Our data strongly suggests the potential roles of the bovine ACE-2, NPR1, and AXL as potential receptors for BCoV infection. Both bovine host cell Furin and TMPRSS2 could potentially cleave the BCoV-S glycoprotein and play pivotal roles during the activation of the BCoV infection. The current study highlighted the roles of the identified potential receptors and host cell proteases in fine-tuning BCoV replication and molecular pathogenesis. However, further studies are needed to confirm the above-mentioned points.