Decoding bovine coronavirus immune targets: an epitope informatics approach

Abstract

Bovine coronavirus (BCoV) poses a significant threat to the global cattle industry, causing both respiratory and gastrointestinal infections in cattle populations. This necessitates the development of efficacious vaccines. While several inactivated and live BCoV vaccines exist, they are predominantly limited to calves. The immunization of adult cattle is imperative for BCoV infection control, as it curtails viral transmission to calves and ameliorates the impact of enteric and respiratory ailments across all age groups within the herd. This study presents an in silico methodology for devising a multiepitope vaccine targeting BCoV. The spike glycoprotein (S) and nucleocapsid (N) proteins, which are integral elements of the BCoV structure, play pivotal roles in the viral infection cycle and immune response. We constructed a remarkably effective multiepitope vaccine candidate specifically designed to combat the BCoV population. Using immunoinformatics technology, B-cell and T-cell epitopes were predicted and linked together using linkers and adjuvants to efficiently trigger both cellular and humoral immune responses in cattle. The in silico construct was characterized, and assessment of its physicochemical properties revealed the formation of a stable vaccine construct. After 3D modeling of the vaccine construct, molecular docking revealed a stable interaction with the bovine receptor bTLR4. Moreover, the viability of the vaccine’s high expression and simple purification was demonstrated by codon optimization and in silico cloning expression into the pET28a (+) vector. By applying immunoinformatics approaches, researchers aim to better understand the immune response to bovine coronavirus, discover potential targets for intervention, and facilitate the development of diagnostic tools and vaccines to mitigate the impact of this virus on cattle health and the livestock industry. We anticipate that the design will be useful as a preventive treatment for BCoV sickness in cattle, opening the door for further laboratory studies.

Introduction

Bovine coronavirus (BCoV), a member of the genus Coronavirus and family Coronaviridae, is an important pathogen affecting cattle populations worldwide (Liu et al. 2006). The virus is associated with respiratory and gastrointestinal infections, causing significant economic losses in the livestock industry. It can infect various parts of cattle respiratory and digestive systems, including the upper lower respiratory tract, respiratory system, and intestinal tract, at various ages and has recently been recognized as a primary cause of neonatal calf diarrhea (Saif 2010; Fulton et al. 2015). Cattle can contract BCoV through two primary routes, the fecal–oral route and inhalation of aerosols (Kin et al. 2016), and commonly results in conditions such as calf diarrhea (CD), winter dysentery (WD), and bovine respiratory disease complex (BRDC) (Gong et al. 2020). BCoV can infect cattle of all ages; however, the disease severity is greater in calves. BCoV infection contributes to winter dysentery (WD) development in adult dairy cattle, causing a dramatic decrease in milk production and significant economic losses (Vlasova and Saif 2021).

Since its first discovery in the United States, BCoV has been detected on five continents, including America, Europe, Asia, Oceania, and Africa; however, the incidence rate and timing of BCoV infection vary among nations. Prior to 2000, BCoV was reported to cause intestinal symptoms in America, Europe, and Asia.

The genome of BCoV contains an enveloped positive-sense RNA with a pleomorphic structure and a size ranging from 65 to 210 nm (Clark 1993). It is characterized by a double layer of surface projections, consisting of short (hemagglutinin) and long (spike) projections. The large genome is composed of 5 important structural proteins encoded by single-stranded RNA, two of which play important roles in viral attachment and infection, namely, the spike glycoprotein (S) and nucelocapsid (N). For viral RNA detection assays, the nucleocapsid (N) is frequently the target since it is highly conserved among strains (Cho et al. 2001). Along with other CoVs, an outer-surface spike (S) glycoprotein is present on BCoV. The S protein comprises an S1 subunit, which contains the dominant neutralizing epitopes, and an S2 subunit that facilitates viral membrane fusion. Proteins are important for the viral life cycle, and they were an excellent selection for our study.

Despite recurrent outbreaks in recent years, BCoV lacks efficient mitigation and prophylactic methods. Some BCoV vaccines are available to prevent gastrointestinal disease in infants (Cho et al. 2001; Hasoksuz et al. 2002; Fulton et al. 2016). Three inactivated vaccines are available and are given to pregnant cows and calves throughout pregnancy to improve the humoral immunity of newborn calves against three different neonatal gastrointestinal diseases (Cho et al. 2001; Fulton et al. 2016). One modified live virus vaccine is provided orally to stimulate a strong immune response and prevent enteric disease in newborn calves (Cho et al. 2001; Fulton et al. 2016). Adults are known to shed the virus while being asymptomatic, which poses a challenge for the control of disease among young adults. Hence, there is a need for an effective vaccine for adults, which is currently lacking. Addressing the limitations of current BCoV vaccinations and extending their applicability to older cattle is crucial for enhancing the overall efficacy and coverage of vaccination strategies. The fact that existing vaccinations primarily target calves leaves a significant gap in addressing infections in older cattle and the broader population. There are several challenges associated with extending the application of this vaccine to older cattle. Exploring solutions to these challenges involves research into effective vaccine design, vaccine formulations, adjuvants, and delivery methods that can stimulate robust immune responses in older animals. This might demand the need for in silico design using bioinformatics tools that are based on data from numerous vaccine studies (María et al. 2017). Addressing the limitations of current BCoV vaccines and extending their applicability to older cattle requires a multifaceted approach in which the use of bioinformatics tools to design the vaccine could accelerate the process of finding an efficacious vaccine against BCoV.

A new field for creating effective multiepitope vaccines (MEVs) has recently evolved as a result of recent advancements in immunoinformatics technology and the understanding of the host immune response. This has significantly accelerated the improvement of vaccines (María et al. 2017). An effective multiepitope vaccine (MEV) should encompass promising antigenic epitopes sourced from viral proteins added to adjuvants. This combination aims to provoke an optimal protective immune response, thereby enhancing overall vaccine efficiency (Rana and Akhter 2016). Prioritized epitopes need to have human leukocyte antigen (HLA) binding patterns that work with the host’s major histocompatibility complex (MHC) molecules. Bovine leukocyte antigen (BoLA) molecules, which resemble mammalian MHC in both structure and function, are known as MHCs in cattle (Takeshima and Aida 2006). Epitopes taken from viruses mimic natural pathogenic elements, making them capable of inducing both humoral and cell-mediated immune (CMI) responses. Importantly, these epitopes carry a reduced risk of causing allergenic reactions, resembling the natural immune reaction to pathogens (Skwarczynski and Toth 2016; Tahir Ul Qamar et al. 2018). The aforementioned approach is effective enough to target different viruses, including those that affect humans and animals.

Researchers and scientists are increasingly applying this methodology to combat emerging and re-emerging infectious diseases, employing the power of immunoinformatics to predict antigenic epitopes and design effective vaccines. These methods have become significant, and vaccines developed through these methods have demonstrated effective in vivo protection and have progressed to phase I clinical trials by saving time and costs (Kar et al. 2020; Tahir ul Qamar et al. 2020). The first epitope-based vaccine was designed against N. meningitidis using bioinformatics, after which many of the vaccines were designed with promising results (Hoque et al. 2021). For example, epitope-based vaccines targeting pathogens such as the influenza virus have been developed by leveraging bioinformatics tools to identify conserved epitopes that offer broad protection against multiple strains. Similarly, vaccines against human papillomavirus (HPV) have benefited from epitope mapping using bioinformatics, leading to the design of vaccines that target specific immunogenic epitopes associated with HPV-induced cancers (Friend Tambunan and Aditya 2012). In the fight against malaria, bioinformatics has played a crucial role in predicting antigenic epitopes of Plasmodium falciparum for inclusion in vaccine candidates (Pritam et al. 2020). Moreover, bioinformatics-driven epitope prediction has advanced HIV vaccine research by identifying conserved regions capable of eliciting immune responses, paving the way for the development of vaccines that target critical epitopes (Pandey et al. 2018). These examples illustrate the transformative impact of bioinformatics in epitope-based vaccine design across a spectrum of infectious diseases.

In this research, an array of immunoinformatics tools was employed to create the initial MEV targeting BCoV, with the aim of providing effective protection against this pathogen. We successfully identified highly antigenic viral protein epitopes, such as B-cell epitopes, cytotoxic T lymphocyte (CTL) epitopes, and helper T lymphocyte (HTL) epitopes. The current vaccine incorporates all the prioritized epitopes, which are connected using appropriate linkers and adjuvants to ensure an optimal immune response. Through a computational biology approach, we conducted a comprehensive assessment of various immunological and physicochemical parameters, including stability, flexibility, and solubility. Subsequently, a 3D model of the construct was generated, refined, and subjected to thorough quality assessment and validation processes.

Furthermore, we assessed the binding affinity between the construct and bovine Toll-like receptor-4 (bTLR4) using molecular docking. The stability of this interaction and associated molecular dynamics (MD) simulation were employed to validate and confirm the binding. Finally, we proceeded to perform in silico cloning of the final construct with codon optimization in a prokaryotic expression system. This step was aimed at facilitating future large-scale production with improved translation efficiency.

Results

Protein sequence retrieval

The spike glycoprotein and nucelocapsid protein sequences of the BCoV isolate India were obtained from the NCBI database under accession numbers UZN72603.1 and UZN72609.1, respectively. The retrieved protein sequence of the nucleocapsid is 448 aa long, while the spike protein is 1,363 aa long. With MEGA 11 software, a phylogenetic study of the chosen protein was conducted using the neighboring joining-tree bootstrap technique. Any bootstrap value greater than 70–80% is regarded as a high bootstrap value, and the branch is likely to be reliably dominant. For the spike glycoprotein, analysis revealed four subgroups and one outgroup, while the nucleopcapsid protein showed two primary subgroups (Fig. 1).

Fig. 1

Phylogenetic tree analysis of A Spike glycoprotein and B nucleocapsid protein using MEGA 11 software

B-cell epitope prediction

B-cell epitopes are essential for the spread of virus resistance. With a default threshold score of 0.51, the ABCpred method for predicting linear epitopes was used to predict a total of 139 and 46 linear epitopes from spike and nucelocapsid, respectively. After confirming the immunogenicity of the epitope, 19 promising B-cell epitopes from the spike glycoprotein and nine from the nucleocapsid protein were predicted (Supplementary Table 1). The top five B-cell antigenic epitopes were selected to create a stable vaccine that was simple to obtain (Table 1).

Table 1 List of the top five ABCpred linear B-cell epitopes of the BCoV spike glycoprotein and nucleocapsid

Cytotoxic T-lymphocyte (CTL) epitope prediction

Considering the strong tendency to connect with particular BoLA alleles, the CTL epitopes were predicted. The FASTA sequences of the spike and nucelopasid proteins were uploaded to NetMHCpan 4.1 to select the dominant BoLA allele. The epitopes were selected based on their high prediction score and low percentile rank against the BoLA alleles. To develop a successful vaccine, the epitopes containing the greatest binding scores for each targeted allele were used. Interestingly, the spike glycoprotein interacts more with the selected BoLA allele than with the nucelocapid allele, which is associated with only five BoLA alleles (Table 2).

Table 2 List of highly conserved, antigenic, nonallergenic and nontoxic MHC-I epitopes of the BCoV spike and nucleocapsid proteins

Helper T-lymphocyte (HTL) epitope prediction

The MHC class II epitope was identified by utilizing NetMHCIIpan 2.1 with selective BoLA-DRB3 alleles. By considering an IC50 value of 50 nM, the lowest percentile rank score and greatest prediction score of the 15-mer MHCII epitopes were selected. The selected epitopes demonstrated strong binding affinities with the various subtypes of the BoLA DRB3 allele. Using the five most potent BoLA-DRB3 alleles (BoLA-DRB3_1501, BoLA-DRB3_0101, BoLA-DRB3_1101, BoLA-DRB3_14011, and BoLA-DRB3_1201), 11 strongly bound (threshold < 50.00) epitopes were obtained, of which nine distinct epitopes were identified from the spike protein and two from the nucelocapsid protein (Table 3).

Table 3 List of antigenic, nonallergenic and nontoxic MHC class II epitopes of the BCoV spike glycoprotein and nucleocapsid protein obtained from the NetMHCIIpan 2.1 tool

Multiepitope vaccine design

To create a vaccine, the identified epitopes were used. Linkers were used to connect the 10 predicted linear B-cell epitopes using KK, 12 CTL epitopes using AAY, and 11 HTL epitopes using GPGPG. At the N-terminus, the adjuvant β-defensin 2 (UniProt ID: P85150), which is connected through the EAAAK linker, was employed. The final vaccine, which included linkers and adjuvants, measured 615 peptides in length, as shown in Fig. 2.

Fig. 2

Multiple epitope vaccine constructs against BCoV are depicted graphically. A 615 amino acid long vaccine construct consisting of an adjuvant at the N-terminus is linked with the multiepitope sequence through the EAAAK linker (green). The BCE, CTL and HTL epitopes are fused with the support of KK (blue), AAY (dark pink) and GPGPG (dark green) linkers, respectively

Physiochemical properties of the constructed vaccine

The vaccine was found to be nonallergenic, nontoxic, and highly antigenic, with an antigenic score of 0.52 at a 0.4% threshold according to the VaxiJen server. The physicochemical properties of the constructed vaccine were predicted using the ExPASy ProtParam server. The molecular weight of the construct, which was found to be 63.78 kDa, reflects its good antigenicity and ease of purification. The basic nature of the peptide is indicated by its pI value of 9.90. At 0.1% absorption, the extinction coefficient was calculated to be 53.220, considering that all cysteine residues were reduced. The half-life of the protein was determined to be 100 h in human reticulocytes, > 20 h in yeast, and > 10 h in Escherichia coli when measured in vivo, suggesting its capacity for long-lasting exposure and immune system stimulation of the host. Furthermore, the construct’s stability was confirmed with an instability index of 11.99. The strong thermostability and hydrophilicity, as indicated by the GRAVY (grand average of hydropathic) index of -0.043 and the aliphatic index of 84.49, respectively, led to enhanced interactions in the polar environment of the body. Taken together, these findings revealed that this structure is a strong candidate for vaccination. The protein was soluble upon overexpression, according to the SOLpro website, with a probability of 0.80. The outcomes of every anticipated physiochemical property are shown in Table 4. Overall, the results suggested that this construct could be a potential vaccine candidate.

Table 4 Physiochemical characteristics of the vaccine

Secondary modeling, refinement and validation of the vaccine construct

Using the transform-restrained Rosetta tool, the resulting 3D structure of the vaccine antigenic peptide was designed. The top five scoring models were downloaded and refined via the GalaxyRefiner online tool. The GalaxyRefine server developed five models based on the square root of deviation (RMSD) and MolProbity technique. Model 1 was selected due to its greater docking performance relative to the Ramachandran model and because it had the highest score (95.8%). GalaxyRefine estimated a clash score of 16.1. MolProbity (2.015), GDT-HA (0.9721), and RMSD (0.356) were among the additional parameters that were calculated. A model with a lower Z score was considered to be of higher quality because it reflects the model’s overall quality. The first model’s Z score is -6.79, while the refined model’s Z score is -6.99, indicating that the improved model is not much different from the first model (Fig. 3).

Fig. 3

Demonstration of the 3D model’s structural refinement, quality evaluation, and validation of the vaccine construct; A The tertiary structure of the improved construct shows a helix, strand, and random coil; B The ProsA Z score (-6.99); C A Ramachandran plot of the improved model (95.8%)

Protein–protein docking analyses of a multiepitope-based vaccine against TLR4

Molecular docking was used to determine the ability of the improved vaccine construct to bind the bovine TLR4 (bTLR4) immune receptor and to determine whether this immune activation was effective. The ClusPro server produced 26 standard outputs among these, and the model with the lowest binding energy score was chosen because it represents good binding affinity. Model 7 was chosen as the best-docked complex because it has the lowest energy criterion, which is -1592.2. According to these findings, this MEV is a great possible vaccination candidate (Fig. 4).

Fig. 4

Docking structure of MEV constructed (red) with the bovine receptor TLR4 (blue) visualized through the Discovery Studio tool

Molecular dynamics simulation of the VRR complex

For the analysis of protein‒receptor (vaccine-TLR4) docking, the server iMODS was used. The highest-ranked model chosen for assessing functional mobility showed that TLR4 and vaccine constructs moved closer to each other, indicating robust and stable binding, as shown in Fig. 5A. According to the deformability investigation, the distortions of the docked complex were greatly reduced compared to those of the monomeric TLR4 protein (Fig. 5B). Furthermore, B-factor analysis revealed slight atomic aberrations in the docked complex, as shown in Fig. 5C. The combined and particular variances of the complex are represented as green and blue bars, respectively, in the variance analysis bar plot, which is inversely proportional to the eigenvalue (Fig. 5D). In addition, 1.341878e-07 was identified as the eigenvalue of the docked complex. Here, the docked complex had significantly greater eigenvalues, which indicates that the complex is more stable (Fig. 5E). The interaction between residue pairs of proteins in a complex is depicted by covariance matrix analysis, where uncorrelated and anti-correlated motions are represented by white, red, and blue colors, respectively (Fig. 5F). Elastic network analysis was used to evaluate the stiffness of the protein complexes. Greater protein stiffness is indicated by the darker gray spots at certain portions shown in Fig. 5G. According to an elastic network model, the atoms of docked protein molecules are connected by “springs” with varying strengths (the stiffer springs are represented by darker grays indicating stiffer regions). TLR4 complexes and vaccine constructs appear to be stable based on the results of the iMODS simulation. The molecular dynamics simulation results suggest that our vaccine model is stable.

Fig. 5

The molecular dynamics simulations of the Vaccine and TLR4 docked complex (A); Main chain deformity simulation (B); B-factor values calculated by normal mode analysis, quantifying the uncertainty of each atom (C); Variance (D); The Eigenvalue of the docked complex (E); Covariance map, red: correlated, white: uncorrelated and blue: anti-correlated (F); The elastic network of the model, describes relation between the atoms and springs (darker gray) (G)

Codon adaptation and in silico cloning

The E. coli expression system is essential for effective vaccine expression in the in silico cloning process. To maximize protein expression, the Java codon adaptation tool (JCat) was used for codon optimization. The optimized codon sequence has a length of 1,845 nucleotides. It has a codon adaptation index (CAI) of 0.93 (0.8–1.0) and an average GC content of 53.92% (30–70%), which points to a high likelihood of the final vaccine being effectively expressed in the E. coli host. Using the EMBOSS backtranseq server, the sequence of amino acids was subsequently determined and translated into the nucleotide sequence. To ensure complementation in the direction of vector translation, this optimized sequence was reversed, and restriction sites were inserted into the 5’ and 3’ ends. The restricted sequence was successfully ligated into the pET28a (+) vector, resulting in a 5,842 bp clone. SnapGene software was used to construct and visualize the cloned map (Fig. 6).

Fig. 6

The MEV-BCoV construct was in silico cloned and inserted into the pET28a (+) expression vector. At the Nhel and BstEII restriction sites, the vaccine sequence was codon-optimized and inserted into the vector. The vaccine construct is represented by the red area, and the vector backbone is represented by the remaining area

Discussion

Bovine coronavirus (BCoV) is the primary cause of respiratory and gastrointestinal sickness in cattle, which results in large economic losses in the global beef and dairy cow industries. Cattle can develop BCoV, which is known to cause calf diarrhea (CD), winter dysentery (WD), and respiratory diseases in cattle (BRDC), mostly through the fecal-oral route and aerosol inhalation (He et al. 2016; Kin et al. 2016; Geng et al. 2023). The available BoCV immunizations are permitted to prevent intestinal disease only in newborns (Cho et al. 2001). One live virus immunization containing BCoV and three Vaccines that have been rendered inactive can be obtained for oral administration to only newborn calves (Awadelkareem and Hamdoun 2022). Since BCoV is believed to cause respiratory illness in bovine populations of various age groups (calves and adults) and Winter dysentery affects mature cattle and is closely related to human coronaviruses and the number of animal ailments, researchers and cattle farmers have paid special attention to it (Saif 2010). Therefore, as a top priority, to curtail BCoV outbreaks, there is a need for rapid and efficient vaccination programs for its prevention and control in calves and adult cattle.

Although the use of vaccines to minimize the effects of calf respiratory illnesses in dairy and beef cattle is widespread, strong scientific support is lacking. It has been difficult to create effective immunization methods (Murray et al. 2016). Additionally, traditional vaccine development processes are costly, time-consuming, and labor-intensive. The production of effective in silico vaccines with less in vitro testing is made possible by the computationally assisted approach of next-generation vaccinology (Pyasi et al. 2021). The most immuno-dominant epitopes produced by infectious agents that mimic natural infections can now be predicted. Furthermore, the addition of adjuvants and linkers helps to improve immunogenicity. As a result, it offers the opportunity to engineer and modify epitopes to improve the stability and effectiveness of the designed multiepitope vaccine construct, which has been demonstrated to be beneficial (Chauhan et al. 2019). Interestingly, these vaccines outperform monovalent vaccines because they have the natural ability to rapidly and strongly elicit innate, humoral, and cellular responses (Amanna and Slifka 2011). Several studies on multiepitope vaccines that used immunoinformatics methods to fight infection have shown encouraging results (Almofti et al. 2021).

For efficient MEV-BCoV design, we specifically selected BCoV structural proteins (spike and nucleocapsid) owing to their regulatory roles in virus infectivity and pathogenicity (Fulton et al. 2015). The structural protein sequences of BCoV were obtained from NCBI to predict potential B- and T-cell epitopes. Phylogenetic tree analysis was also performed. Compared to the nucleocapsid protein, which has been shown to have a poor bootstrap value, the spike glycoprotein was found to have high bootstrap values and excellent branch reliability, suggesting that the predicted epitopes can be considered strong vaccine candidates in these regions. B cellcells initiate a humoral immune response that destroys viruses and creates memory to protect against subsequent exposure, although this response often occurs only partially and worsens with time (Bacchetta et al. 2005). Instead, the CMI response elicited by CTL and HTL precisely controls the spread of pathogens by either destroying infected cells or secreting antiviral cytokines that promote lifelong immunity (Arpin et al. 1995). Epitopes that can induce both B-cell and T-cell immunity are known to be good candidates for vaccines (Kumar Verma et al. 2015). One study, however, suggested designing MEV-BCoV using two distinct structural proteins, the spike protein and hemagglutinin esterase, which play important roles in immunological defense (Awadelkareem and Hamdoun 2022). Our studies, in contrast, focused on making such predictions using the spike protein and nucelocapsid of structural proteins. Consequently, each epitope type was used in the construction of the vaccine. Our findings predicted the highest-scoring linear B-cell epitopes for each examined protein. robust affinity epitopes for binding with experimentally confirmed alleles are typically an excellent choice for use in the construction of MEV constructs (Sohail et al. 2019). We selected BoLA alleles due to their prevalence among bovine species. To identify the CTL and HTL epitopes, a thorough analysis of all susceptible dominant BoLA alleles of class I/II molecules was carried out. The highest-ranked epitopes (CTL and HTL), determined by the established threshold as confirmed by BoLA class-I/II, were selected for construction. Furthermore, a number of spike glycoproteins and nucleocapsid epitopes interact with various BoLA alleles, indicating the possibility of broader immune responses to diverse BCoV strains.

The successful completion of the requisite parameters was followed by adjuvant and linker attachment of the promising epitopes. Adjuvant-β defensin-2 is connected to the N-terminus of the vaccination design because it improves effectiveness, stability, and long-term survival. It functions as the initial line of defense. in the fight against the numerous infections that affect dairy cattle (Gurao et al. 2017). At the infection site, it binds to its corresponding receptors, activating both developing dendritic cells and naive T cells (Mackenzie-Dyck et al. 2014; Gurao et al. 2017). As linkers improve the expression, folding, stability folding, and stability of separate domains, linkers were added as essential elements in the vaccine construct (Pyasi et al. 2021). The main purpose of vaccines is to effectively elicit an immune response while having little to no effect on the host. Later, the developed vaccine showed strong antigenicity, nonallergenicity, and nontoxicity, along with strong solubility and other physiochemical features. A thorough structural assessment employing 3D refinement and analysis through the Ramachandran plot demonstrated a stable and high-quality model. Subsequently, to predict the stability and binding affinity of the interacting ligand and receptor complex (vaccine construct and bovine-TLR4 receptor), molecular docking and molecular dynamics (MD) simulations were performed. experiments were performed. Since the actual receptor for BCoV is still unknown, common virus infection TLR receptors, including TLR7 and TLR8, can be considered. However, a key receptor involved in recognizing the viral glycoprotein (spike glycoprotein) and initiating the immune response, bovine TLR4 (bTLR4), is selected as a receptor against the MEV-BCoV vaccine construct. TLR4 is involved in viral infections, becomes activated in the innate response to nonbacterial microbial pathogens in vivo and is considered to be involved in the sensing of SARS-CoV-2 infection (Jung and Lee 2021). TLR4 plays a vital role as a receptor in triggering a proinflammatory response, responding to both viral and noninfectious stimuli (Molteni et al. 2016; Ghosh et al. 2021a). The docking score indicated a significantly strong binding affinity, showing stable interactions within the ligand‒receptor complex. This discovery was confirmed and validated through molecular dynamics simulations. The MD simulation results of the normal mode analysis utilizing iMODS showed that upon molecular binding, both TLR4 and the vaccine candidates stably interacted with each other (Hayward et al. 1997). This finding suggested that the developed vaccine activated TLR4, potentially leading to enhanced immune responses in the host.

Due to inconsistencies in mRNA codons, gene expression differs across various hosts, highlighting the essential role of codon optimization in achieving elevated modes of communication (Ali et al. 2017). The codon-optimized vaccine design must have a suitable GC content and CAI value in the E. coli expression vector for optimized protein expression levels (Chen 2012). Based on earlier research, E. coli is the system that is most recommended for large-scale generation of recombinant proteins (Pei et al. 2005). The targeted vaccine candidate was successfully cloned in silico into the pET28a (+) cloning vector of the E. coli K12 expression host after the codon was optimized.

The experimental validation of the predicted multiepitope subunit vaccine is crucial for confirming its effectiveness and safety. Experimental validation is necessary to ensure that the predicted epitopes can induce the desired immune response and to assess the immunogenicity, safety, and efficacy of the vaccine. This validation often involves techniques such as peptide binding assays, T-cell proliferation assays, and animal studies to assess the ability of the predicted epitopes to elicit a specific immune response, identify conserved epitopes, and assess potential cross-reactivity and immunodominance. In the pursuit of identifying a promising BCoV subunit vaccine candidate, our study employed an immunoinformatic approach to predict potential epitopes. These vaccines offer advantages over monovalent vaccines due to their inherent potency in eliciting innate, humoral, and cellular immune responses (Amanna and Slifka 2011; Tahir ul Qamar et al. 2020). According to the findings of our study, the multiepitope vaccine under design may be subjected to in vitro and in vivo experimental assessments in the hopes of creating a vaccine against BCoV.

Although the current study resulted in the development of an MEV-BCoV vaccine candidate through in silico methods and provided a starting step for a strong foundation for experimental studies, it is essential to recognize that BCoV, an RNA virus, undergoes frequent mutational changes. These adaptations allow it to account for antigenic variation under field conditions. Previous research has highlighted the advantages of similar in silico approaches in creating effective vaccines for various infectious diseases, including SARS-CoV-2 (Naz et al. 2020), Chandipura virus (Pavitrakar et al. 2022), Nipah virus (Majee et al. 2021), Ebola virus (Ullah et al. 2020), and Zika virus (Kumar Pandey et al. 2018). However, some researchers have performed experimental validation to confirm their in silico findings and have claimed that this approach could facilitate the development of an effective and time- and cost-efficient vaccine. Mustafa (2013) validated the in silico method of Mycobacterium tuberculosis (TB) by measuring antigen-specific cellular and humoral immune responses in vitro using peripheral blood mononuclear cells and sera from TB patients and BCG-vaccinated healthy subjects. He found that the bioinformatics approach facilitated the identification of novel candidates for TB diagnosis and vaccination (Mustafa 2013). In another study, Khalili et al. 2017 confirmed the previously constructed epitopes of HBV using computational methods. They expressed the protein in E. coli and tested the designed antigen using a chemiluminescent immunoassay. They stated that using bioinformatics tools enabled the rational design of multiepitope antigens in a more economical, intelligent and knowledge-based way. They also suggested that their results could serve as preliminary evidence that computational predictions could be applied as initial steps of biological studies and their subsequent experimental conditions (Khalili et al. 2017). Similarly, in 2019, Shruthi et al. performed in silico identification and wet laboratory validation of cryptic B-cell epitopes in ZnT8 by measuring ZnT8-specific isotypes (IgG, IgM and IgA) in the sera of normal glucose-tolerant (NGT), type 1 diabetic (T1DM) and type 2 diabetic (T2DM) patients by indirect ELISA. Their results revealed significantly decreased levels of IgG and IgA isotypes in T1DM patients without complications. They identified novel cryptic B-cell epitopes in the ZnT8 autoantigen against which naturally occurring autoantibody levels were found to be reduced in diabetes (Shruthi et al. 2019). The evidence collectively indicates that the integration of immunoinformatics streamlines the optimization and validation of experimental procedures, thereby improving the efficiency and accuracy of subsequent research efforts.

Leveraging immunoinformatics technology provides a deeper understanding of host immune responses and streamlines both time and cost. However, it is crucial to acknowledge the inherent limitations, such as algorithm variability and length constraints (Gazi et al. 2016). For instance, the ABCpred tool predicts linear B-cell epitopes utilizing a recurrent neural network (RNN), but it has several limitations, such as fixed-length patterns and sequence accuracy (Saha and Raghava 2006). Similarly, the NetMHCpan 4.1 and metMHCIIpan 2.1 tool predictions, trained on binding affinity (BA) data, may have limitations in terms of prediction performance since they only model a single event of peptide–MHC binding, neglecting other biological features involved in the process (Reynisson et al. 2021). This limitation may hinder the prediction of vaccine candidates and thus open the door for future experimental/laboratory studies.

Conclusion

Bovine coronavirus is associated with high morbidity and mortality. This study used bioinformatics to design a multiepitope subunit vaccine for BCoV. Certainly, a key attribute of an effective vaccine lies in its capacity to confer lifelong immunity, guarding against recurrent episodes of an infection. The purpose of this investigation was to create a subunit vaccine for bovine coronavirus (BCoV) by employing an integrated immunoinformatics approach aimed at targeting multiple epitopes. The resulting MEV-BCoV model demonstrated several features that have the potential to trigger both cellular and humoral immune responses, suggesting a promising contribution to the development of a BCoV vaccine. The designed construct meets the necessary criteria for antigenicity, allergenicity, toxicity, and various physicochemical parameters and aligns optimally with the required standards. Molecular docking and dynamic simulation revealed a notably robust binding affinity for TLR4, suggesting that TLR4 is stable within the physiological pH range. Furthermore, meticulous codon optimization and in silico restriction cloning were carried out to guarantee efficient expression in the widely used E. coli K12 strain. It is essential to highlight that our findings are solely based on computer-aided technology, setting the stage for subsequent in vivo and in vitro assessments aimed at validating the reliability, effectiveness, and safety of vaccine constructs.

Methods

Protein sequence retrieval

The complete nasal protein sequences of the BCoV spike glycoprotein (S) and nucleocapsid (N) proteins were retrieved from the National Centre for Biotechnology Information (NCBI) website (https://www.ncbi.nlm.nih.gov/) and saved in FASTA file format. MEGA 11 software was used to carry out phylogenetic analysis for both of the chosen proteins. A sequence with a similarity of more than 98% was found using the NCBI BLAST program.

B-cell epitope prediction

B-cell epitopes play a significant role in the initiation of humoral or antibody-mediated immunity. These epitopes can be found on the outside of viral antigens (Safavi et al. 2020). To identify potential B-cell epitopes in prioritized proteins, the complete FASTA sequence of the obtained protein was uploaded to the ABCPred service (http://crdd.osdd.net/raghava/abcpred/) with the default threshold value of 0.51 (Saha and Raghava 2006). The same software has been utilized by additional researchers to forecast B-cell epitopes (Ahmad et al. 2022; Aziz et al. 2022; Li et al. 2023). The ABCpred server forecasts linear B-cell epitope regions in an antigen sequence using an artificial neural network. This server will help in identifying epitope areas that will aid in the selection of candidates for synthetic vaccines (Saha and Raghava 2006).

Cytotoxic T-lymphocyte (CTL) epitope prediction

The first stage in eliciting an immune response to viral infections involves the presentation of viral antigens by MHC-I to CTLs (Pyasi et al. 2021). By uploading the FASTA sequence to the NetMHCpan 4.1 server (http://www.cbs.dtu.dk/services/NetMHCpan/), cytotoxic (CD8+) T-cell epitopes for every target protein were predicted. Using artificial neural networks (ANNs), the NetMHCpan-4.1 server forecasts peptide binding to any MHC protein with a specified sequence (Fisch et al. 2021). For other bovine diseases, this server is frequently used for the prediction of CTL epitopes (Connelley et al. 2022; Pathak et al. 2022). The dominant BoLA alleles (BoLA-HD6, BoLA-JSP.1, BoLA-T2c, BoLA-T2b, BoLA-T2a, BoLA-D18.4, BoLA-AW10, and BoLA-T5) were selected for our research because they represented Bos taurus, Bos taurus indicus, and hybrid bovine species. Additionally, nucleotide-level polymorphisms are more pronounced in cattle and buffalo (Santos Junior et al. 2020; Ghosh et al. 2021a, b; Yılmaz Çolak 2021; Patra et al. 2023). The FASTA sequences of the spike and nucleocapsid proteins were submitted to the server with binding thresholds of 0.5% and 2%, respectively, as the standard criteria for categorizing strong and weak binders. For the prediction, a score was established along with a binding score (Ysrafil et al. 2022).

Helper T-lymphocyte (HTL) epitope prediction

HTL plays a substantial role in activating both humoral and cellular immune responses, as MHC-II peptides originating from foreign proteins derived from the extracellular milieu have been identified. HTL epitopes, therefore, play an essential role in developing immunotherapeutic vaccines. The outcomes were predicted against the highly dominant BoLA-DRB3 allele (BoLA-DRB3*1501, BoLA-DRB3*0101, BoLA-DRB3*1101, BoLA-DRB3*14011, BoLA-DRB3*1201) (Fisch et al. 2021) using the NetMHCIIpan 2.1 server (https://services.healthtech.dtu.dk/services/NetMHCIIpan-2.1/), which has a high affinity for TH-cell activation (Singh et al. 2018). Using artificial neural networks (ANNs), the NetMHCIIpan server forecasts the peptide’s attachment to more than 500 HLA-DR alleles. A collection of 2,00,000 randomly selected natural peptides provided the predicted values, which are reported as nM IC50 values and as%-Rank (Nielsen et al. 2008; Pyasi et al. 2021; Pathak et al. 2022).

Prediction of antigenicity, allergenicity, toxicity and stability

VaxiJen v2.0 (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html) was utilized to test the antigenicity of the chosen epitopes. This tool is based on autocrosscovariance (ACC) translation of protein sequences into uniform vectors of primary amino acid characteristics (Doytchinova and Flower 2007). Using the ToxinPred (http://crdd.osdd.net/raghava/toxinpred/) and ExpasyProtparm tools (https://web.expasy.org/protparam/), the toxicity and immunogenicity of particular epitopes were assessed. The allergenicity was then predicted using AllerTop v.2.0 (https://www.ddg-pharmfac.net/AllerTOP/) (Dimitrov et al. 2014).

Chimeric subunit vaccine construction

A single peptide chain was constructed by combining the detected epitopes with the aid of particular peptide linkers. Peptide linkers are crucial for the folding of proteins, versatility, and division of functional domains, all of which contribute to a more stable protein structure. The screened B cell, MHC-I, and MHC-II epitopes from the target protein were joined together through linkers to generate a vaccine against several epitope sequences. B cells were linked together using the KK linker. CTL epitopes were linked together using the AAY linker, and the HTL epitopes were connected to the GPGPG linker. β-Defensin 2 was chosen as a supplementary agent. and linked through a linker at the EAAAK at the N-terminus, a vaccination sequence (Ali et al. 2017). Adjuvants in epitope-based vaccines provide several advantages, including significantly extending the vaccine’s long-term memory and helping elderly individuals with a blocked immune system regenerate (Aasim et al. 2022).

Physiochemical analysis of the constructed vaccine

To prevent unwanted immunological reactions, a vaccination needs to be allergy-free. Using the VaxiJen 2.0 online server, the antigenicity of the vaccination protein with the adjuvant was calculated. Using the AllerTop v.2.0 tool, the allergy-inducing potential of the final vaccine and each of its components was evaluated. Epitope toxicity was evaluated utilizing the ToxinPred service. The constructed solubility of the vaccine was predicted utilizing the SOLpro tool. (https://scratch.proteomics.ics.uci.edu/). SOLpro is an SVM-based approach for predicting the solubility of a protein sequence, with an estimated overall accuracy of over 74% based on tenfold cross-validation (Magnan et al. 2009). The physicochemical properties, including molecular weight, number of polypeptides, theoretical isoelectric point (PI), half-life, instability index, aliphatic index, and hydropathicity grand average of hydropathy (GRAVY), were calculated using the ExpasyProtParam online tool (Wilkins et al. n.d.; Gao et al. 2021).

Secondary modeling, refinement and validation of the vaccine construct

Using the trRosetta online tool (http://yanglab.nankai.edu.cn/trRosetta/), a secondary component of our final vaccine was predicted. The trRosetta algorithm predicts protein structures quickly and precisely. The structural organization of proteins was constructed using constrained Rosetta and direct energy minimizations (Du et al. 2021). Using the GalaxyRefine website (https://galaxy.seoklab.org/cgi-bin/submit.cgi?type=REFINE), further structure refinement was carried out. GalaxyRefiner refines the predicted structure by relaxing and repacking the side chain of the structure (Heo et al. 2013). The model’s refinement was evaluated using the GDT-HA score, RMSD score, MolProbity score, clash score, and Ramachandran’s plot score. Furthermore, the ProSA webserver validated the resulting structure by calculating the total model quality score (Rawal et al. 2021).

Molecular docking

To effectively trigger an immunological response from the host, a vaccine must interact with immune receptors in an effective manner. Therefore, the computational molecular docking approach is used to predict microscopic interactions among complex interacting macromolecules. As a key receptor for recognizing viral peptide structures that initiate the immune response, bovine TLR4 (bTLR4) was selected as a receptor for the MEV-BCoV constructed vaccine (Vaure and Liu 2014). However, since the Protein Data Bank could not provide the crystal structure of bTLR4, 3D modeling of the structure was performed utilizing the trRosetta web server using the sequence obtained from the UniProt database. Using the online protein‒protein docking service Cluspro2.0 (https://cluspro.bu.edu/home.php), molecular docking was performed to calculate the interaction between the constructed MEV and bTLR4 (Kozakov et al. 2017). Multiple docked complex models are produced in the output, each with a different estimated electrostatic interaction value and the lowest Gibbs free energy rating. Later, the PyMol tool was utilized for visualization of the docked complexes (Yuan et al. 2017).

Molecular dynamic simulation

Molecular dynamics simulations were performed with the iMODS server (https://imods.iqfr.csic.es/). iMODS performs a critical study of the structure by modifying the complex’s force field in relation to various time intervals (López-Blanco et al. 2014). The iMODS service calculates the protein’s internal coordinates using normal mode analysis (NMA) to assess the protein’s stability. By uploading the docked molecules, the server results in various graphs, which include the main-chain deformability plot, B-factor values, eigenvalue, covariance matrix, residues, and atom index model that serve as representations of the protein’s stability (Kalita et al. 2020).

Codon adaptation and in-silico cloning

For a vector to express the desired gene, the codons used by the host vector must be aligned with those of the desired gene so that high levels of activity are achieved during purification. The higher expression rates may result from the codon adaptation technique adapted to E. coli K12; for the host organism E. coli K12, the approach was utilized to boost the expression of the main order of the subunit vaccination protein, which was then sent to the JAVA Codon Adaptation Tool. Therefore, the Java Adaptation Tool (JCat) (http://www.jcat.de/) was used based on codon adaptation index (CAI) values. Rho-independent transcription termination and prokaryotic ribosome binding were avoided. To validate high expression, a protein must exhibit a desirable GC content of 30–70%, and the CAI should be between 0.8 and 1.0. Using JCAT, the protein sequence of the vaccine construct was reverse-transcribed. The SnapGene tool (https://www.snapgene.com/snapgene-viewer/) was further used to clone the cDNA sequence into the pET28a (+) vector (Pyasi et al. 2021).