Applying high throughput and comprehensive immunoinformatics approaches to design a trivalentsubunit vaccine forinduction of immune response against human emerging coronaviruses SARS-CoV, MERS-CoV and SARS-CoV2

Coronaviruses (CoV) cause diseases such as severe acute respiratory syndrome (SARS), Middle East respiratory syndrome (MERS) and Coronavirus disease 2019 (COVID-19). Therefore, this study was conducted to contrast a trivalent subunit vaccine against SARS, MERS and COVID-19. The CTL, HTL, MHC I, and IFN-γ epitopes were predicted. Moreover, to stimulate strong helper T lymphocytes (HTLs) responses, Pan HLA DR-binding epitope (PADRE) was used. Also, for boosting immune response, β-defensin 2 was added to the construct as an adjuvant. Furthermore, TAT was applied in the vaccine to facilitate the intracellular delivery. Studies show stability fully higher expression rate. Physico-chemical factors analysis via ProtParam online server showed that an estimated half-life of the subunit vaccine was 30 hours in mammalian reticulocytes, more than 20 hours in yeast, and more than 10 hours in E. coli. Furthermore, the instability index was 39.09 that indicates a subunit vaccine is a stable protein. Also, the tendency of vaccine solubility upon overexpression in E. coli was computed 0.95% showing an admissible percentage of solubility in an overexpressed state. Outputs of the primary model with the Galaxyweb server showed that the Galaxyweb program model was applied to generate the primary 3D structure of the protein vaccine. Obtained results showed that 43.76% alpha-helix, 16.93% extended strand, and 39.31% random coil structure was in the tertiary model. Besides, the renement of modeled structure was carried out in GalaxyRene online server via molecular dynamics simulation. After the renement, applying via the GalaxyRene server resulted in modeling showed a high-quality tertiary model. The analysis of all discontinuous and continuous B-cell epitopes showed that the identied epitopes on the protein surface could interact comfortably with antibodies, and they were commonly exible. DisEMBL online server was utilized to specify unstructured regions and disorders in the subunit vaccine. Obtained results indicated that disorder regions were sited among positions 1–23, 464–478, 689–697.


Background
Coronaviruses (CoVs) are positive-strand viruses from the family Coronaviridae that are named according to their crown-like shape [1]. CoVs have caused major outbreaks in recent years [2,3]. Various human CoVs have spread globally, causing mostly mild illness, while others are highly pathogenic. The rst severe acute respiratory syndrome (SARS)-CoV emerged in china in 2001 and resulted in a global outbreak. However, this outbreak was successfully contained. In 2012, another CoV variant named Middle East respiratory syndrome (MERS)-CoV emerged in Jeddah, Saudi Arabia. However, this outbreak never perfectly diminished, and hospitals and communities outbreaks continue to occur annually [4]. Recently, a novel coronavirus, SARS-CoV-2 initiated an outbreak wich turned into a pandemic, raising global consternation. This newly surfaced virus that caused a new disease named coronavirus disease (COVID-19) has a considerably lower mortality rate. However, the rapid transmission of the virus has turned COVID-19 into a fatal pandemic, which indicates a raising need for novel treatments and preventive methods.
Efforts have been made for vaccine and drug design regarding CoVs. As a result, some vaccines are currently being tested in phase I clinical trials for SARS-CoV-2 in the USA (NCT04283461) and China (NCT04313127). Despite the efforts, outcomes of therapeutics and preventive approaches aimed at CoVs, in particular vaccines have not been adequately successful. Several aspects of CoVs may highlight the importance of preventing.Notably, various subtypes of CoVs, including some subtypes like MERS-CoVs that have caused past outbreaks, have remained relevant [5].
CoV proteins can be classi ed into two subgroups, structural and non-structural proteins (nsps). The genome of CoVs is surrounded by structural proteins comprising spike (S) protein, nucleocapsid (N), envelope (E), and membrane (M) proteins [6]. S protein exerts various functions and mediates cellular entrance of CoVs.
Following attachment of viral envelope-anchored S protein to its receptor, CoVs can enter the cytoplasm. Protease breakdown of S protein is required for virus-cell fusion and the entry of genomic RNA [6][7][8].S protein of SARS-CoV-2 is highly similar to SARS-CoV. Interestingly, polyclonal antibodies for SARS-CoV suppress spikeregulated entry of SARS-CoV-2 [9]. HCoV-NL63, a virus that only causes mild symptoms, also has S protein and recruits the angiotensin-converting enzyme 2 (ACE2) for entering the host cell, while showing very low pathogenicity. This indicates that additional elements are responsible for the pathogenicity of CoVs that have caused signi cant outbreaks. [10]. Furthermore, the N protein has been suggested as a promising candidate for cross-protective vaccine design. This is because of the high similarity of SARS-CoV, SARS-CoV-2, and MERS-CoV N proteins, while being signi cantly different from CoVs that mainly cause mild illness. [11]. The E and M proteins have been suggested as vaccine candidates in RV assessments. The E and M proteins exert crucial roles in the viral assembly of CoVs. Further, the N protein is essential for viral RNA synthesis. Deletion of E protein has inhibited the virulence of CoVs, and a number of researches have been directed at the e cacy of live attenuated vaccines that comprise recombinant SARS-CoV or MERS-CoV with a mutated E protein [12]. Also, M protein has been utilized for DNA vaccine construction against SARS-CoV [13]. Most vaccines comprise attenuated or inactivated whole virus and structural proteins. However, such vaccines could not confer complete protection and whole vaccines may cause adverse reactions. The non-structural proteins (nsp), a part of ORF1a polyprotein, may confer additional protection. Following S protein, nsp3 is the second-highest pathogenicelement of SARS-CoV-2 protein [11].
Taken together, an e cient multipotent preventive solution against various subclasses of the rapidly emerging human CoVs remains to be discovered. This is highly important because of the recent CoV outbreak, which has raised an urgent need for the development of vaccine research. For this purpose, bioinformatics analyses have been utilized to aid the vaccine design process. Hence, we utilized in silico methods to design a trivalent novel cross-protective multi-epitope vaccine against SARS-CoV, MERS-CoV, and SARS-CoV-2 via spike (S) protein, nucleocapsid (N), envelope (E), membrane (M) protein, nsp3, and nsp8 antigens.

Primary analysis of the retrieved sequences
In this study, spike (S) protein, nucleocapsid (N), envelope (E), membrane (M) protein, nsp3 and nsp8 antigens (totally 6 protein) were used for bioinformatics analysis to nd CTL, HTL, IFN-γ and MHC I binder's epitopes.
The highly conserved regions of β-defensin 2 were employed as an adjuvant because the vaccine requires a strong induction of T cell immune responses. The signal sequences were excluded from all constructed structures. The work ow of the designed work representing the overall procedures of vaccine design is shown in Fig. 1.
Cytotoxic T lymphocyte (CTL) prediction CTL epitopes were selected according to the highest speci city and lowest sensitivity for their adaptive immune receptor. Top three high score epitopes were chosen from spike S, N, E, M, nsp3 and nsp8 antigens to get a total of 18 nal CTL epitopes (Table 1). On the other hand, we utilized the Immune Epitope Database (IEDB) tool to  predict MHC-I peptide binders to human alleles HLA-B*57:01, HLA-A*26:01, HLA-A*11:01, HLA-B*35:01,  Finally, 18 epitopes were chosen from the predicted HTL epitopes according to their lowest percentile rank ( Table 2).The lowest percentile rank re ecting their good binding a nity for the MHC-II receptor.

Construction of the subunit vaccine
To design the nal vaccine construct, the top-scored epitopes of CTL and HTL were joined together using HEYGAEALERAG linker.To facilitate cell penetration of the proposed vaccine, the TAT peptide of HIV (GRKKRRQRRRPPQ) was added at the N-terminal of the construct. In addition, for increasing the immunogenicity of the vaccine, a 41 amino acid-long sequence of β-defensin 2 was added after TAT peptide as an adjuvant and then followed by the PADRE sequence (AGLFQRHGEGTKATVGEPV). The TAT peptide and βdefensin 2 were linked to the PADRE sequence throughGPGPG and EAAAK linkers, respectively. The order of the nalconstruct is summarized in Fig. 2.

Population coverage
The predictionof MHC-I and MHC-IIbased coverage of the applied epitopes in the constructed vaccine accomplished by IEDB analysis resources for the world population.As SARS-CoV-2 is global pandemic, hence, worldwide analysis has been accomplished. The world population coverage of MHC-I and MHC-II found to be 88.77%, and 97.91%, respectively, that showed in Fig. 3.
Optimized codon and mRNA structure of the chimeric vaccine We use the JCAT online server to obtain the codon-optimized DNA sequence of the proposed vaccine to verify for cloning and expressing in E. coli strain K12 as a host. The optimized codon sequences had 2091 nucleotides long, GC-content was reached to 52.9 % and its codon adaptation index was raised to 0.98. The graphical illustration of the optimized gene was shown in Fig. 4A. In addition, the mfold server was used to get free energy related to the whole mRNA structure and its 5′ end. Its details are shown in Table 3. As shown in Fig.  3B, C, and D, the minimum free energy of the predicted secondary RNA structure was ΔG = -595.30 kcal/mol at 37 °C without hairpin nor pseudoknot at the 5′ side.
Evaluation of allergenicity, antigenicity, physiochemical parameters and solubility of the vaccine construct Unveiling allergenicity of the proposed vaccine was a crucial consideration to avoid allergic reactions. The resulting output of the Algpred online server determined the non-allergenic behavior of the vaccine constructs with a score of -0.59 (Threshold=-0.4).Likewise, the antigenicity of the construct was predicted by the VaxiJen v2.0 online server, and it con rmed that the protein was a probable antigen with an overall prediction score of antigen 0.4996.Furthermore, the physicochemical parameters were evaluated, and the molecular weight of vaccine protein was found to be 74.8 kDa which is good to support the antigenic nature of the vaccine construct.The predicted instability index of the proposed vaccine was 39.09, which is below 40 hence, classifying the protein as stable.The isoelectric point (pI) of the construct was 9.30 and it shows the basic nature of the vaccine.
The estimated half-life (in-vitro) in mammalian reticulocytes was found to be 30 h, while invivo studies in Yeast and E. coli more than 20 and 10 h, respectively. In addition, the aliphatic index was 80.26, which determines the thermostable nature of the proposed vaccine, and the grand average of hydropathicity obtained to be 0.063, which proves the hydrophilic nature of the constructed vaccine.The value for the solubility tendency of the proposed vaccine upon overexpression in E. coli host was 0.55.
Secondary and tertiary structure of subunit vaccine PSIPRED server was employed to predict the secondary structure of the proposed vaccine (Fig.5). This server showed that the vaccine is made of 43.76% alpha-helix, 16.93% extended strand, and 39.31% random coil.The tertiary structure was predicted as 3 domain structures using the galaxyWEB web server which is a templatebased protein structure modeling web server that resulted in a modeled protein structure (Fig. 6). The 1FD3_C, 5MKE_A, and 5MKF_A were used as the best-suited template to model the tertiary structure. All 697 amino acid residues were modeled, while the 70-92, 681-688, and 690-697 amino acids of the residues were unreliable local regions.
Re nement and validation of the 3D structure The generated 3D structure was submitted for the re nement process by the GalaxyRe ne online server to enhance the quality of predicted 3D modeled structure beyond the accuracy.Furtheremore, the re ned model from GalaxyRe ne was analyzed by the ramachandran plot.Validation of 3D structure by ramachandran plot before re nement demonstrated 540(92%), 37 (6.3%), and 9(1.5%) of residues were placed in favored, additionally allowed and disallowed regions, respectively (Fig. 7A). In the re ned model, 550 (93.7%), 26 (4.4%), and 7 (1.2%) of residues were placed in favored, additionally allowed and disallowed regions, respectively (Fig.  7B).

Prediction of intrinsic protein disorder
DisEMBL online server identi ed and annotated disordered regions. Amino acids in the input sequence are considered disordered by remark-465 de nition in the following regions: 1-23, 464-478, 689-697.Based on IUPred results, amino acids were considered disordered when the con dence score is higher than 0.5. Disordered regions were also presented in Fig. 8.

Predicted B-cell epitopes
Ellipro software and the Discotope online server were utilized to predict the linear B-cell and discontinuous B cell epitopes, respectively. The linear B-cell and discontinuous B cell epitopes on the 3D structure are shown in  Table 4.
Vaccine protein disul de engineering Disul de engineering leads to an increase in vaccine stability with mutating the residues present in high mobility areas with cysteine. 76 residue pairs were found to have the ability of disul de bond formation. But, when all these residue pairs were more analyzed on the other parameters like energy on chi3 value and B-factor, only 1 pairs of residues provided the disul de bond formation. Those residue pairs were GLY360-SER368 and ASP580-GLY601 (Fig. 10). These two residues were mutated to the cysteine. The energy value analyzed for the residue screening was less than 2.2 and the value of chi3 was between -87 to +97.

Molecular docking of vaccine construct with the TLR-3
Protein-protein docking was done with the help of the Hdock online server between the nal vaccine construct as a ligand and the TLR-3 as a receptor. We have procured 20 models and the model with the lowest binding energy was selected. The Docking score for the best acceptable model was -326.86, whereas the ligand RMSD (Å) was 205.96.PatchDock server was used to rea rm the docking between the nal vaccine construct and TLR-3. Top 20 models were provided by the server and the model with the top-ranked and lowest binding energy was chosen.The obtained output encompasses the score (17080), Area (2493.80), ACE (408.17), and transformation (-2.83 -0.57 -1.83 -49.29 3.75 160.90). The docked complex achieved from the server has been shown in Fig. 11A and B.

Molecular dynamics simulation results
MD simulation was used to analyze the physical movements of the atoms and molecules of the proposed vaccine. At the initial step of MD simulation, energy minimization was done and several parameters including pressure and temperature were measured.Regular mode investigation allowed the demonstration of the stability of proteins and the large scale mobility. The iMODS server performed such an investigation based on internal coordinates of the receptor-ligand complex (Fig. 12A). In 3D model, the direction of each residue was given through arrows and the length of the line represented the extent of mobility. The proposed vaccine protein and TLR3 receptor were oriented towards each other. The eigenvalue found for the complex was 4.953620e-08 (Fig. 12B). The B-factor values deduced from normal mode investigation was analogous to root mean square (Fig. 12C). Hinges in the chain indicated the feasible deformability of the complex calculated through the complication of each residue (Fig. 12D). The variance related to each normal mode was inversely related to the eigenvalue. Covariance matrix described the occlusion between pairs of residues was correlated, uncorrelated or anti-correlated motions were represented with red, white, and blue colors, respectively (Fig. 12  E).The result also provided an elastic network model (Fig. 12F) that identi ed the pairs of atom connected by springs. Each dot in the diagram was indicated one spring between the corresponding pair of atoms and colored based on the grade of rigidity.

In silico cloning
To examine the expression of the recommended vaccine into a vector, cloning was done by utilizing the SnapGene server.As the expression system of the human and E. coli K12 differ from each other and the codon requires the adaptation as per the host expression system, the codon sequence that was optimized by JCat was reverse transcribed and modi ed as per the E. coli K12. The restriction enzym sitesNot I and Xho I were introduced at the N and C-terminal, respectively. The adapted codon was entered into the pET21b (+) vector and a recombinant plasmid with the length of 7514 base pair was obtained as displayed in Fig. 13.
Immune simulation C-ImmSim server results revealed a noticeable increase in the generation of immune responses. The humoral immune response was identi ed by high levels of IgM along with IgG + IgM (Fig. 14A). Estimation of the B cell population uncovered an increase in B memory cells and B-cell isotype IgM with a corresponding reduction in antigen concentration (Fig. 14B). In addition, the high response was observed in the TH1 and TC cell populations with corresponding memory development (Fig. 10C, D). This pro le identi ed the development of immune memory and consequently enhanced clearance of the antigen upon subsequent exposures. In contrast, regulatory T cells were estimated in low level (Fig. 14E). Moreover, these predictions were consistent with the induced level of IFN-γ produced after immunization with the proposed vaccine (Fig. 14F).

Discussion
Design novel vaccines for preventing the ever-rising global burden of disease is an important issue that cannot be ignored. Coronaviruses (CoV) are a wide family of zoonotic viruses that reason illness ranging from the common cold to more severe diseases such as Severe Acute Respiratory Syndrome (SARS-CoV) Middle East Respiratory Syndrome (MERS-CoV) and SARS-CoV-2 [14,15]. While many coronavirus vaccine studies targeting various structural proteins were conducted, most of these endeavors nally stopped soon after the outbreak of SARS and MERS [11].
In the present study, we purposed to design an epitope-based peptide vaccine with structural and nonstructural proteins against the more severe form of human coronavirus SARS-CoV-2, SARS-CoV, and MERS-CoV via immunoinformatics strategies utilizing various in silico tools. Advancement in genomics has alternated the procedures and concepts to vaccine candidate design and chosen [16]. Reverse vaccinology, a novel method to merge immunogenetics and immunogenomics via immunoinformatic has been used tremendously to introduce novel vaccine [17,18]. Thus, in this study, we developed a novel trivalent subunit vaccinevaccine from immune-dominant peptides of spike (S) protein, nucleocapsid (N), envelope (E), membrane (M) protein, nsp3, and nsp8 antigens. In our suggested construct, the CTL, MHC-I, MHC-II, epitopes along with IFN-γ containing epitopes were predicted and used in our designed vaccine. Cytotoxic CD8 T lymphocytes (CTL) limit the spread of infectious pathogens via recognizing and killing infected cells or secreting particular antiviral cytokines [19,20]. Thus, T-cell epitope-based vaccination is a unique role of eliciting a potent immune response against infectious pathogens such as viruses [21]. Helper T lymphocyte (HTL) responses act an important role in the induction of both humoral and cellular immune responses. Therefore, HTL epitopes are believably to be a key part of prophylactic and therapeutic vaccines [22]. Several pieces of research introduced IFN-γ as the effective cytokine of both the adaptive and natural immune system eliciting antiviral, immune regulatory as well as antitumor activities. The secretion of IFN-gamma is the key part of the Th1 response critical for the rein of pathogens such as viruses [23]. Furthermore, once the development of the novel multi-epitope vaccine begins, there will be great concern about its immunogenicity. Thus, adjuvants were established and attached to vaccine construction to enhance its e cacy. There are some complications with the most routine aluminum adjuvants such as anthrax, macrophagic myofasciitis, and post-vaccination phenomena [24]. Recently, the role of TLRs in the relation of innate and adaptive immunity for vaccine design has been highlighted. TLRs are a subgroup of pattern recognition receptors on antigen-presenting cells (APCs), which activating via an appropriate TLR agonist will elicit cellular and humoral reactions as well as innate responses [25,26]. There are multiple ways to induce TLRs. β-defensin 2 which induce signal transduction via TLR-4 regulates adaptive immunity during infection [27]. Self-destructive signaling to delete activated APCs and thus modifying the adaptive immunity is also possible through utilizing β-defensin 2 [28].
In our vaccine construction human β-defensin 2 serves as an adjuvant that promotes the expression of primary immune-inducing and antiviral compounds level [29]. To increase the quality of vaccine delivery, a cellpenetrating peptide, namely TAT, from HIV was engaged to involve T cells via APCs [30]. PADRE T helper epitope can bind to a prominent number of HLA-DR types with an a nity spanning from intermediate to high [31]. Also, the YQVNNLEEI epitope located in a conserved region of EBOV NP added to the C terminal of the vaccine construct. Some studies indicated that this epitope has a protective CTL epitope against COVID-19 [32].
Different parts of the vaccine structure were connected through AAY, GPGPG, HEYGAEALERAG, EAAAK, and GGGS linkers. These sequences exert a critical role in the 3D structure and subsequent exposure of candidates [33]. The AAY linker was used to join HTL epitopes. AAY enhances two subsequent epitopes cleavage [34]. EAAAK linker was used to link PADRE and adjuvant to decrease the interaction with other vaccine domains and cause more effective separation. We used the GGGS linker to connect CTL epitopes. GGGS linker containing sequences are powerful stimulants of the immune response [35,36]. We utilized the GPGPG linker to connect TAT and adjuvant. Through polypeptide and DNA vaccination, GPGPG can act as a prominent HTL inducer.
This spacer, however, as an assistant for epitope presentation by HLA-II, inhibits generation of junction [37,38]. The HEYGAEALERAG linker provides the speci c cleavage target for both the lysosomal and proteasomal degradation system. The twelve-residue peptide consists of all ve appropriate cleavage sites, A7-L8, A5-E6, Y3-G4, R10-A11, and L8-E9, that are determined for eukaryotic proteasome complexes [39]. For joining the CTL epitopes to HTL epitopes and HTL epitopes to the PADRE and YQVNNLEEI epitope to the histidine chain, the HEYGAEALERAG linker was used.
Allergenicity is one of the main problems in vaccine development. Recently, most vaccines stimulate the immune system into an allergic reaction [40]. based on the WHO/FAO, if a sequence has a recognize of at least six contiguous amino acids over a window of 80 amino acids (0.35% sequence identity) to identi ed allergens, it is considered to be potentially allergenic. In this study, AlgPred online server was applied for the allergenicity assessment of the designed vaccine. The result of the allergenicity assessment showed that the vaccine isn't allergen. Thus, the subunit chimeric vaccine was constructed with 697 residues in length. Restriction sites, Instability elements, and all the cis-acting sites were omitted from the construct which signi cantly interfering with cloning. Moreover, codon optimization was performed to facilitate the high-level expression of the multiepitope vaccine in E. coli and improve translational and transcriptional output. mRNA secondary structure is a serious parameter in the expression of subunit vaccines. The outputs from mRNA prediction by Mfold online server showed that the mRNA had enough stability for impressive translation in the host. Studies show that higher stability result fully leads to higher expression rate. Physico-chemical factors analysis via ProtParam online server showed that an estimated half-life of the subunit vaccine was 30 hours in mammalian reticulocytes, more than 20 hours in yeast, and more than 10 hours in E. coli. Furthermore, the instability index was 39.09 that indicates a subunit vaccine is a stable protein. Also, the tendency of vaccine solubility upon overexpression in E. coli was computed 0.95% showing an admissible percentage of solubility in an overexpressed state. Outputs of the primary model with the Galaxyweb server showed that the Galaxyweb program model was applied to generate the primary 3D structure of the protein vaccine. Obtained results showed that 43.76% alpha-helix, 16.93% extended strand, and 39.31% random coil structure was in the tertiary model. Besides, the re nement of modeled structure was carried out in GalaxyRe ne online server via molecular dynamics simulation. After the re nement, applying via the GalaxyRe ne server resulted in modeling showed a high-quality tertiary model. The analysis of all discontinuous and continuous B-cell epitopes showed that the identi ed epitopes on the protein surface could interact comfortably with antibodies, and they were commonly exible. DisEMBL online server was utilized to specify unstructured regions and disorders in the subunit vaccine. Obtained results indicated that disorder regions were sited among positions 1-23, 464-478, 689-697.
An intrinsically disordered protein (IDP) is classically described as a protein with absenting a xed or wellstructured 3D fold. IDPs include a spectrum of states from wholly unstructured to partially structured and contain (pre-) molten globules, random coils, and large multi-domain proteins connected via exible linkers.
Disordered areas are mostly functional. Many disordered segments fold on binding to their biological targets and some others organize exible linkers causing an assembly of the macromolecular array [41]. Then, the docking process between TLR3 and the designed vaccine was accomplished and the foremost docked model was selected and investigated via MD simulation to compute the stability of the vaccine-protein complex.
Outputs of MD simulation support that the subunit vaccine can interact with TLR3 protein be ttingly. The interactions between the two molecules were modi ed during the MD simulation.
The infectious form of the attenuated vaccine may reverse or inactivation process may fail to repress the virulence property, In the case of conventional vaccines [42,43]. Furthermore, some recent clinical trial studies enhanced safety concerns while corroborating on conventional method for vaccine development [44]. Some when the cellular lines employed for the generation of the vaccine may be contaminated via microorganisms that do not produce any cytopathic effect while undergoing replication along with the vaccine virus [45,46].
These contaminations perhaps proved dangerous for the host, when the quality control process failures too. However, the recombinant multi-epitope vaccine does not carry any risk of the emergence of virulence property as it contains the antigenic sections only (not the entire virus itself). Moreover, several instruments have been employed in the current study to accredit the non-allergic and nontoxic behavior of the designed vaccine construct. However, the predicted bioinformatics results were based on the various investigation of sequences and different immune servers. We offer future wet lab-based analyses employed in vivo models for the experimental validation of our subunit vaccine.

Conclusion
In the current study, we work to develop an e cient trivalent subunit vaccinevaccine against CoVs via different in silico servers. Immunoinformatics analysis indicated that this vaccine was able to induce both humoral and cellular immune responses effectively due to the presence of several epitopes from a combination of antigens, PADRE, and adjuvant. Thus, this subunit vaccine can conceivably be utilized for prophylactic usages against emerging human pathogenic coronaviruses SARS-CoV, MERS-CoV, and SARS-CoV2.

Cytotoxic T lymphocytes (CTL) and IFN-γ inducing epitopes prediction
Choosing CTLs epitopes is one of the most important elements of vaccine development, because CTLs kill virus-infected or damaged cells by recognizing the epitope presented by MHC-I molecule. After all, the initial step to stimulate the immune system is presenting an antigen on MHC-I. CTLPred (http://crdd.osdd.net/raghava/ctlpred/) is a freely accessible web server for the prediction of CTL epitopes. This server is based on the arti cial neural network (ANN), support vector machine (SVM), the consensus and, the combined prediction methods. In addition, IFN-γ plays an essential role in the innate and adaptive immune response by stimulating natural killer cells and macrophages and provides an elevated response to MHC antigens. IFN-γ epitope server (http://crdd.osdd.net/raghava/ifnepitope/scan.php) was also used to predict interferon-positive sequences in the selected antigens.
Helper T-lymphocyte (HTL) epitope prediction IEDB prediction tool (http://tools.immuneepitope.org/mhcii/) was employed to predict helper T cell (HTL) epitopes among the selectedantigens and each peptide was 15mer in length. The epitopes were predicted for human HLA-II allelesas target species. We used the IEDB recommended method to achieve the best possible result.
Trivalent subunit multi-epitope vaccine construct and gene optimization The N-terminal sequence of the subunit vaccine starts with TAT peptide of HIV as cell-penetrating peptides (CPPs) and followed by the adjuvant (β-defensin 2) as adjuvant. In the following, PADRE was utilized as helper T lymphocyte epitopes. Then, HTL and CTL selected epitopes were added sequentially for each SARS-CoV-2, SARS-CoV, and MERS-CoV protein, respectively. The AAY, GPGPG, EAAAK, KK, and GGGS amino acids sequences were used as suitable linkers to attach selected epitopes to achieve the multi-epitope vaccine consists of an antigenic part of a pathogenic organism to elicit an immunological response. Finally, the DNA sequence of the subunit vaccine and its adaptation to E. coli codon usage was obtained via the JCAT server at http://www.jcat.de.
Population coverage IEDB population coverage tool (http://tools.iedb.org/population/) was used to gure out the scale of the population coverage of the nal construct within worldwide population. Epitope/MHC restriction data were loaded on IEDB and, calculations were set to be done for MHC-I and MHC-II separate epitopes.

Physico-chemical parameters analysis
The codon usage was analyzed by Java Codon Adaptation Tool (JCat) carry out the codon optimization.The Mfold server (http://unafold.rna.albany.edu/)was used to predict the second structure of RNA. Mfold online server predicts structures and the true positive base pairs with minimum ΔG with its thermodynamic methods.
In addition, the physical and chemical parameters that are associated with the vaccine were discovred. Some physicochemical characteristics of the designed structure, including amino acid composition, instability index, molecular weight, PI value, estimated half-life in vitro and in vivo, aliphatic index, extinction coe cient, and grand average of hydropathicity (GRAVY) were evaluated using ProtParam online server at https://web.expasy.org/protparam/.

Prediction of solubility, antigenicity and allergenicity
SOLpro online server (http://scratch.proteomics.ics.uci.edu) was utilized to predict the propensity of a protein to be soluble upon overexpression in Escherichia Coli. This server uses a two-stage SVM architecture method to evaluate the protein solubility based on multiple representations of the primary sequence. In addition, prediction of allergenicity with high precision was discovered by AlgPred at http://crdd.osdd.net/raghava/algpred/. Furthermore, antigenicity of the proposed vaccine was assessed using another online server VaxiJen v2.0 (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html). The accuracy of this server varies from 70% to 89% depends on the organism.

Prediction of secondary and tertiary structure
The secondary structure of the predicted vaccine was carried out by using GOR IV and the PSIPRED online servers at http://expasy.org/tool/gor4.html and http://bioinf.cs. ucl.ac.uk/psipred/,respectively.The percentage of the helix (H), strand (S), and coils (C) were determined.Homology modeling of the nal protein construct was done by GalaxyTBM at http://galaxy.seoklab.org/index.html. GalaxyTBM is a template-based tertiary structure prediction online server and it predicts protein structure from the sequence when experimental structures of homologous proteins are available as templates and re nes loop or terminus regions by ab initio modeling.
Re nement and validation of the 3D modeled structure The tertiary structure of the ultimate vaccine construct received from the RaptorX server depends on the degree of similarity between available template structure options and the target protein.Therefore, the GalaxyRe ne server (http://galaxy.seoklab.org/cgi-bin/submit.cgi?type=REFINE) was used to improve the modeled structure beyond the accuracy.The CASP10 assessedthe GalaxyRe ne server and reported that the GalaxyRe ne is one of the best-performing methods to improve the local structure quality. In addition, SAVES v5.0 server at http://servicesn.mbi.ucla.edu/PROCHECK/ was employed to validate the re ned model and the Ramachandran plot of this structure was designed.
Prediction of intrinsic protein disorder DisEMBL 1.5 (http://dis.embl.de/) and IUPred (http://iupred.enzim.hu/pred.php) were utilized to predict intrinsic protein disorder and unstructured regions within a protein. DisEMBL is a public online server for predicting disorder in proteins. Also, DISOPRED3 was employed to predict intrinsic protein disorder and unstructured regions in our proposed vaccine.

Prediction of conformational and linear B-cell epitopes
We utilize Discotope 2.0 server at http://www.cbs.dtu.dk/services/DiscoTope/ to predict conformational epitopes for the modeled vaccine. It's using the calculation of surface accessibility estimated in terms of contact numbers and a new epitope propensity amino acid score to predict discontinuous B cell epitopes from the tertiary structure. The nal scores are determined by combining the propensity scores of residues in spatial proximity and contact numbers.In addition, ElliPro online server at http://tools.iedb.org/ellipro/ was employed to predict linear and discontinuous B-cell epitopes. ElliPro classi es the conformational epitopes according to their protrusion index values and clusters them based on distance R (Å), while the Discotope method utilizes a calculation of surface accessibility and epitope propensity amino acid score. Finally, conformational epitopes with high-score were selected.

Disul de bond engineering
Improving the stability of the modeled vaccine was an important step that was completed by utilizing a logical approach to disul de bond formation. Disul de binds provide notable stability to protein structure and wellde ned geometric conformation. Disul de engineering was done by Disul de by design v2.0 (DbD2) to design novel disul de bonds in the vaccine protein construct [47].

Molecular docking between vaccine construct and TLR-3
Molecular Docking was done to perceive the interaction between the receptor (TLR-3) and the ligand ( nal vaccine construct) at their stable form and to check the binding a nity between them.In humans, the TLR-3 receptor stimulates the immune system by activating the NF-κB signaling pathway followed by the production of different cytokines [48]. The protein-protein docking was performed by the Hdock online server at http://hdock.phys.hust.edu.cn/.HDOCK is a docking server that e ciently integrates various components: sequence search, template selection, and model building. The predictive power of the HDOCK server can be also improved by ranking rst of the template-based model.Moreover, the center of extremely populated clusters and the lowest binding energy of the receptor-ligand complex were predicted.After several rotations of ligand, the docked complex was selected based on the lowest binding energy.Moreover, re-a rmation was performed by the PatchDock online server at https://bioinfo3d.cs.tau.ac.il/PatchDock/ to evaluate the interaction between the designed vaccine construct and the receptor.This server gives the output based on three algorithms: Filtering and scoring, molecular shape representation, and surface patch matching.The output composed of rank, global energy, atomic contact energy, transformations, and PDB le complex.

Molecular dynamics simulation study
The molecular dynamic simulation was employed to check the physical basis of the structure and function of the biological macromolecules. The iMODS server was used to explain the collective motion of proteins by analysis of normal modes (NMA) in internal coordinates. Essential dynamics extracts the correlated motions of proteins to understand the motions that are most essential to the activity of the protein and it is a strong tool and alternative to the costly atomistic simulation that can be contrasted to the normal modes of proteins to specify their stability. This server predicted the direction and magnitude of the immanent motions of the complex in terms of eigenvalues, deformability, covariance, and B-factors [49]. The structural dynamics of the ligand-receptor complex was analyzed.
In silico cloning The optimized sequence of proposed vaccine was reversly trascripted and to ensure the restriction cloning, XhoI and Not I restriction site were added at the N-and C-terminal of sequense, respectively. Then, we used the Snap-Gene restriction cloning module to insert the optimized sequence into pET21b (+) vector between de ned restriction sites.
Immune simulation The C-ImmSim online server (http://150.146.2.1/C-IMMSIM/index.php) was used to determine the immune response pro le of the constructed vaccine [50]. C-ImmSim simulates three sections that represent three different anatomical regions in mammals including the bone marrow, tertiary lymphatic organ, and the thymus. This is an important step for the simulation of the immune response because it determines immunogenicity.
The binding of the epitope, which is the immunogenic part of an invading pathogen, together with activation and cooperation from T helper cells, is required to trigger an immune response in the affected host.

Declarations
Ethics approval This study was approved by the Ethical Committee of Babol University of Medical Sciences

Consent for publication
All authors agree to publish this study Availability of data and material The datasets applied and analyzed during the current study are available from the corresponding author on reasonable request.