key: cord-0880064-g2bt8mke
authors: Khan, Salman; Shaker, Bilal; Ahmad, Sajjad; Abbasi, Sumra Wajid; Arshad, Muhammad; Haleem, Abdul; Ismail, Saba; Zaib, Anita; Sajjad, Wasim
title: Towards a novel peptide vaccine for Middle East respiratory syndrome coronavirus and its possible use against pandemic COVID-19
date: 2020-11-06
journal: J Mol Liq
DOI: 10.1016/j.molliq.2020.114706
sha: fe1d341f65e9cc146037636febea9c7915bf1a83
doc_id: 880064
cord_uid: g2bt8mke

Middle East respiratory syndrome coronavirus (MERS-CoV) is an emerging health concern due to its high mortality rate of 35%. At present, no vaccine is available to protect against MERS-CoV infections. Therefore, an in silico search for potential antigenic epitopes in the non-redundant proteome of MERS-CoV was performed herein. First, a subtractive proteome-based approach was employed to look for the surface exposed and host non-homologous proteins. Following, immunoinformatics analysis was performed to predict antigenic B and T cell epitopes that were used in the design of a multi-epitopes peptide. Molecular docking study was carried out to predict vaccine construct affinity of binding to Toll-like receptor 3 (TLR3) and understand its binding conformation to extract ideas about its processing by the host immune system. We identified membrane protein, envelope small membrane protein, non-structural protein ORF3, non-structural protein ORF5, and spike glycoprotein as potential candidates for subunit vaccine designing. The designed multi-epitope peptide then linked to β-defensin adjuvant is showing high antigenicity. Further, the sequence of the designed vaccine construct is optimized for maximum expression in the Escherichia coli expression system. A rich pattern of hydrogen and hydrophobic interactions of the construct was observed with the TLR3 allowing stable binding of the construct at the docked site as predicted by the molecular dynamics simulation and MM-PBSA binding energies. We expect that the panel of subunit vaccine candidates and the designed vaccine construct could be highly effective in immunizing populations from infections caused by MERS-CoV and could possible applied on the current pandemic COVID-19.

subunit vaccine in vaccine formulation is also not an appropriate choice to have as they harbor several different antigenic determinants, all of which are not desired and some may have detrimental effects on protective immunity induction [19] . This rationale has led to create an interest in designing a peptide vaccine for MERS-CoV that comprise of epitopes required for desired cellular and humoral immune response [20] . Peptide vaccines are considered sufficient enough to induce appropriate T and B-cell immunity, while at the same time excluding the risk of reactogenic and allergenic responses [21] . Moreover, peptide vaccine can be implicated for induction of broad-spectrum immunity against multiple variant of the pathogen [22] . On the contrary, because of the small size peptide vaccine is weakly immunogenic and require a carrier molecule for adjuvanting and chemical stability [19] . Several peptide vaccines are under process of development for different viral diseases, such as hepatitis C virus (HCV) [23] , human papilloma virus (HPV) [24] , human immunodeficiency virus (HIV) [25] , and influenza [26] etc. A multi-epitope vaccine contains a series of overlapping epitopes that can elicit effective cytotoxic T cells, T helper cells and B cells responses against viral pathogen [27] . Unlike, subunit or single epitope-based vaccines, multi-epitope vaccines have a distinctive concept design with additional properties. These properties include: (i) the chances of adverse effects or pathological responses are less because of reduced unwanted components, the introduction of adjuvant capacity enhances chances of long-lasting and immunogenicity responses, (iii) as the final formation has multiple epitope the spectra of targeted viruses is expanded, (iv) they are capable of producing strong cellular and humoral immunity simultaneously, (v) the epitopes in multi-epitope peptide can be recognized by T cell receptors (TCRs) of various T-cell subsets [19] .

Keeping in view the broad applicability of multi-epitope vaccine, in this work we explored the complete proteome of MERS-CoV exoproteome and secretome for cytotoxic T lymphocyte (CTL) and helper T lymphocyte (HTL) epitopes containing B-cell epitopes. Additionally, molecular docking, molecular dynamics simulation and MM-PBSA binding free energies were performed to underpin multi-epitope peptide binding pose and interactions with TLR3 receptor proteins.

The complete reference proteome of MERS-CoV (Proteome ID: UP000139997) was retrieved from Swiss-Prot reviewed universal protein knowledgebase database (Uniprot) [28] and subjected to subtractive proteomics pipeline for identification of potential vaccine candidates [29] . First in the process, proteins subcellular localization was predicted using an online tools: CELLO2GO [30] tool with E-value set at 0.001 and Virus-PLoc server (http://www.csbio.sjtu.edu.cn/bioinf/virus/). Proteins having localization in extracellular matrix and membrane were considered and analyzed further along the framework. All the screened proteins were evaluated for redundancy in Cluster Database at High Identity with Tolerance J o u r n a l P r e -p r o o f Journal Pre-proof (CD-Hit) with sequence identity threshold set at 50%. [31] Next, comparative homology analysis of the filtered proteins with the human proteome (H. sapiens, taxid: 9606) was performed with parameters of E-value: 0.005, bit score: ≥ 100 and sequence identity: ≤ 30% in online BlASTp tool of the National Center for Biotechnological Information (NCBI) [32, 33] . To avoid any homology of the protein sequence with the beneficial bacteria of human gut, a subsequent BlASTp was run against Lactobacillus rhamnosus (taxid: 47715) taking the identity sequence cut-off ≤ 30%, bit score of 100 and E-value of 0.005 [34, 35] .

The CTL epitopes for the proteins were predicted using NetCTL 1.2 server (http://www.cbs.dtu.dk/services/NetCTL/). NetCTL 1.2 server combines prediction of Cterminal cleavage, binding affinity for reference set of MHC class I which is restricted to 12 MHC class I super type, and TAP transport efficiency. The prediction of MHC peptide binding, proteasome cleavage event, TAP transport efficiency is done using neural networks, NetChop neural networks, and weight matrix based methods, respectively. Predicted CTL epitopes with overall combine score of 1 were selected. CTL epitopes from the pathogen proteins will elicit cell mediated immunity for inhibiting the growth and development of the pathogen and produce memory cells for remembering the pathogen in future encounters [36, 37] . The predicted CTL epitopes were then subjected to VaxiJen to evaluate antigenicity of the epitopes [38] . Those with values greater > 0.4 were considered antigenic and considered only. The response from Helper T-lymphocyte (HTL) epitopes constitute the major part of cellular immunity and aid in producing various cytokines and immune cell to clear the pathogen [39] . The HTL epitopes induce both CTL and HTL immune responses by secreting different lymphokines and could be crucial part of prophylactic and immunotherapeutic vaccine [40] . HTL epitopes interacting with a reference set of HTL alleles were predicted using IEDB MHC-II epitope prediction module keeping the available parameters as default [41] . All predicted epitopes further subjected to allergenicity, toxicity and IFN-gamma inducing peptides prediction by using online AllergenFP v.1.0 [42] , ToxinPred [43] , and IFNepitope server [44] , respectively. The predicted epitopes were ranked based percentile rank score: lower score describe higher affinity binding of the epitope for HTL binding. The predicted HTL epitopes for its affinity to activate T1 type immune response were cross-validated using IFN epitope server [44] .

The shortlisted conserved epitopes were then linked to each other through a flexible linker AAY. β-defensin was used as an adjuvant to recruit naive T-cells and immature dendritic cells to the infection site [45] . The CTL-epitopes peptide is linked to β-defensin at the N-terminal using EAAAK linker while the CTL and HTL epitopes were linked to each other AAY and GPGPG linker, respectively. The structure of the construct was predicted using Raptor X [46] and ab ignition I-tasser [47] . The physicochemical parameters of the vaccine construct were evaluated using ProtParam server [48] . The generated structures of the multi-epitope peptide were then utilized in PDBsum to create Ramachandran plot for the structures [49] .The antigenicity and allergenicity were predicted by employing VaxiJen [38] and AlgPred servers [50] , respectively.

Identification of B-cell epitopes is vital in vaccine designing considering the fact that it leads to production of B lymphocytes that differentiate into antibody producing plasma cells and memory cells [51] . B-cell epitopes of 20 amino acids long in the input protein sequences were identified using BCPREDS [52] . The epitope specificity threshold was allowed as default by 75%.

Reverse translational analysis was carried out using Codon Usage Wrangler Tool (https://www.mrc-lmb.cam.ac.uk/ms/methods/codon.html) to generate optimized cDNA. The cDNA sequence was then used in GeneScript to compute GC content and codon adaptation index (CAI) score. Ideally, the GC content of a sequence should be in between 30% to 70%, while the CAI score of 1 is considered good.

Multi-epitope peptide vaccine sequence was then used in molecular docking studies in PatchDock [53] . TLR3 (PDB ID, 1ZIW) was used as receptor proteins in the process. The best predate complexes were subsequently subjected to FireDock for refinement [54] . The affinity of the vaccine construct for TLR3 was also validated by ClusPro 2.0 docking software [55] .

The vaccine construct immunogenic potential was elucidated using C-immune server [56] . The simulation parameters were treated as: 10 number of antigen dose, random seed (12345), volume of simulation 100, number of simulation steps,1000, HLA alleles used (MHC class I A0101 allele, B MHC class I B0702, DR MHC class II DRB1_0101 allele), and time step of injection employed is 1.

MD simulation was run for 120 nanoseconds on the most suitable docked solution of TLR3-MEPVC. AMBER16 was used to execute this MD simulation [57] . The ff14SB force field was used to parameterize the complex [58] . The complex was solvated into a TIP3P water box. Na+ ions were added to the system to counterbalance its charge. Hydrogen atoms of the system were minimized for 500 steps, solvation box and carbon alpha atoms for 1000 steps, and all non-heavy atoms for 300 steps. Later on, in an attempt to manage the system temperature, it was heated to 300 K (NVT) for 20 picoseconds through langevin dynamics [59] where a 2 femtosecond time step long restraint of 5 kcal/mol-A 2 on carbon alpha atoms was permitted. System was then unrestrained for 100 picoseconds while being equilibrated. NPT ensemble was put in practice for 50 picoseconds to manage the pressure on the system. Finally, simulation run of 120 nanoseconds was achieved at a time step of 2 femtoseconds. AMBER CPPTRAJ [60] was then used to examine the trajectories for to evaluate complex dynamics and stability.

MMPBSA [55] .py package was used to predict binding free energies of the system [61] . A total of 100 frames were derived from simulation trajectories to be employed in MMGBSA (Molecular Mechanics Generalized Born Surface Area) and MMPBSA (Molecular Mechanics Poisson-Boltzmann Surface Area) [62] . The examination of the free energy difference between solvated phase and gas phase of the system is the ultimate goal of this assay. The determination of net binding free energy was carried out as,

The binding free energy calculation was done on 100 frames extracted from complete length of MD simulation as well as from last 20-ns to thoroughly validate intermolecular affinity and stability. Entropy of the system was estimated using a bash script following the protocol provided by Duan et al., 2016 [63] .

The complete MERS-CoV proteome contains 13 proteins in total, namely: Replicase polyprotein 1ab, Replicas polyprotein 1a, Membrane protein, Envelope small membrane protein, Nonstructural protein ORF4b, Non-structural protein ORF3, Non-structural protein ORF5, Nonstructural protein ORF4a, Nucleoprotein, and Spike glycoprotein. These proteins were first investigated in subcellular localization of subtractive proteomics to screen only those candidates that are exposed for direct interaction with the host cells. Such proteins identification represent an important goal as they are involved in virulence and pathogenesis and their vital importance from stimulating immune responses [64] . Only proteins that are localized to either extracellular (exoproteome) or outer membrane (secretome) were opted in subcellular localization filter. Five proteins: Membrane protein, Envelope small membrane protein, Non-structural protein ORF3, Non-structural protein ORF5, Spike glycoprotein were shortlisted in this analysis as highlighted in bold in Table 1 . Next in the framework, these proteins were subjected non-redundancy check.

In this, non-redundant dataset was only considered. All the five proteins mentioned above were non-redundant which means that each protein represented in the dataset is only singly repeated [35] . The comparative homology check with the humans revealed no hit, indicating highly nonsimilarity between these virus proteins and the human proteome. This is quite vital in avoiding the auto-immune responses and designing a safe vaccine [65, 66] . Lastly, a comparative homology analysis between the MERS-CoV shortlisted proteins and the beneficial gut bacterium, Lactobacillus rhamnosus, was performed to eliminate the risk of accidental inhibition of good bacteria in the host. Bacteria of such type play profound role in improving functions of the host immune system, aid in digestion of food for extracting vital nutirents, and prevent the growth of hostile infection causing bacteria etc [34] . All the five proteins were noticed to show no homology to the said representative gut bacterium proteome allowing the safe use of the filtered proteins epitopes. 

Cytotoxic T cell lymphocytes (also known as Killer T cells) are important component of immune system for tumor surveillance and against intracellular pathogens including bacteria, protozoan and viruses [36] . When these cells recognized its antigen, it become activated and follow three step mechanisms to kill the infected or malignant cells [67] . First, it secretes cytokines primarily IFN-γ and TNF-α which have antiviral and antitumor effects [68] . Secondly, it produces and release cytotoxic granules. These granules have perforin that make the target cell membrane porous and granzymes which being serine proteases hydrolyze intracellular proteins, block viral proteins synthesis and eventually leads to apoptosis [67] . All the five proteins shortlisted in the previous step was analyzed for 9mer cytotoxic T-cell receptor epitopes. Briefly, 5 CTL epitopes were obtained for membrane protein, 3 for envelop small membrane protein, 3 for non-structural protein ORF3, 4 for non-structural protein ORF7 and 38 for the spike protein. Among these predicted epitopes, we found only 4 antigenic epitopes for membrane protein, 0 for envelop small membrane protein, 1 for non-structural protein ORF3 and ORF5, and 17 for spike protein.

Only antigenic epitopes were considered while non-antigenic were discarded. The selected epitopes are also non-allergic, non-toxic and are IFN-gamma inducing peptide. Predicted CTL epitopes for the screened 5 subunit vaccine proteins are shown in Table 2 . 

T helper lymphocytes or CD4+ cells are important player of the immune system by instigating and shaping adaptive immune system. T helper lymphocytes recognize peptides displayed on the MHC class II molecules. They release cytokines that aid in activating several other immune cells. They are essential for B cell antibody class switching, enhancing bactericidal activity of macrophages, and growth and activation of cytotoxic T cells [51, 69] . The five proteins were subjected to MHC-II epitope prediction of 15mer length as shown in Table 3 . Only epitopes of lowest percentile score were selected and as such epitopes with percentile score < 3 were opted only. Using this criterion, we predicted 3 HTL epitopes for membrane protein, 9 for envelope small membrane protein and non-structural protein ORF3, 25 for non-structural protein ORF5 and 31 for the spike protein. Antigenicity evaluation revealed 0 epitopes for membrane protein, envelope small membrane protein and non-structural protein ORF5 4 and 3 respectively, 6 for non-structural protein ORF3 and 17 for spike proteins. The selected epitopes besides being antigenic are also non-allergic, non-toxic and are IFN-gamma inducing peptide. 

A multi-epitope peptide of 874 residues containing sequence conserved epitopes of both T cell immunity (CTL and HTL) and B cell immunity was designed. This design was based on the fact that individual peptide epitope is weakly immunogenic, therefore, fusing the screened epitopes to make a multi-epitope peptide could potentially enhance their immunogenicity capacity. Further, linking this peptide to an adjuvant can further increase the chance of good immune presentation [19] . In total, 53 epitopes were predicted for cytotoxic T cell immunity among which 19 epitopes were found antigenic and only considered for multi-epitope vaccine construct designing. Similarly, 77 HTL epitopes were predicted among which 30 were antigenic and considered. The CTL epitopes were joined together using AAY linkers. These linkers possess the cleavage sites J o u r n a l P r e -p r o o f for proteasome in the mammalian cells thus aiding in the production of natural epitopes and stopping the generation of junctional epitopes and by this way enhancing the presentation of the epitopes presentation to the host immune system [70] . The HTL epitopes were interconnected to each other through GPGPG linker for adequate presentation and stimulation of T cell responses. The EAAAK linker was employed to link the peptide of CTL epitopes to β-defensin adjuvant at the N-terminal and HTL epitopes at the C-terminal. β-defensin is a potent immune adjuvant and functions by stimulating the production of lymphokines which subsequently cellular immunity and production of antigen specific immunoglobulin [45] . The multi-epitope vaccine construct is schematically presented in Fig.1 . 

To ensure the non-allergic responses of the vaccine construct upon administration, the allergenicity of the construct was evaluated [71] . The prediction was based on SVM that securitized the amino acids composition of the input sequence for allergic sequences. The J o u r n a l P r e -p r o o f construct sequence was declared as non-allergen with score of -0.80559987 (cut-off = -0.4). The percent of positive predictive value deduced was 0 % in contrast to negative predictive value that was equal to 0 %. Similarly, the antigenicity of the construct was re-evaluated to certain its immune system evoking potential. The overall antigen prediction score for the vaccine construct is 0.4744 that indicate the probable antigenic behavior of the designed chimeric peptide.

The physiochemical properties of the designed construct were predicted to guide the experimentalists in synthesis and understanding the formulation of the vaccine. The molecular weight estimated for the construct is 140 kDa that explain two things: (i) first the medium size of the construct will allow its convenient purification and secondly, small and medium size chimeric proteins harbor appropriate number of epitopes for specific and targeted immune responses against the pathogen thus limiting the chances of allergic and adverse reactions [72] . The vaccines composed of whole organism or a single large protein have large number of antigenic determinants all of which are probably not required as they could led to non-specific immunity [19] . The theoretical pI computed for the vaccine construct is 8.05 indicating its basic nature. The instability index of the protein is estimated as 27.07 and classify it as stable. This instability index tells an estimate of protein stability in a test tube and is based on statistical analysis of 32 stable and 12 unstable proteins [73] . It has been determined that stable nature of proteins is due to the occurrence of certain dipeptides that differ significantly from non-stable. The protein instability index calculated by Protparam assigned a weightage of instability to different 400 dipeptides. Proteins with instability index of higher than 40 are categorized as unstable. The aliphatic index represents relative volume occupied by a protein aliphatic side chains involving amino acid residues: leucine, isoleucine, alanine and valine [74] . It is regarded as a key factor conferring thermostability to globular proteins. The aliphatic index for the query vaccine peptide is 106.10 presenting the high thermostability of the construct. The estimated half-life of the construct in mammalian reticulocytes (in-vitro), yeast (in vivo) and Escherichia coli (in vivo) is 5.5 hours, 3 and 2 minutes, respectively. This half-life tells the predicted time for the protein to be disappeared after its synthesis. The grand average of hydropathicity (GRAVY) is 0.664. Lower the value implies greater hydrophilicity of the protein and has more probability of interacting with water residues. GRAVY is the cumulative hydrophilicity score of all the amino acids divided by total number of protein sequence.

modelled. The best structure was modelled using the A chain of the 4KNH template by Raptor-X with p-value of 3.18e-04. The 3D model of the predicted structure is shown in Fig.2.   Fig.2 . Predicted 3D structure of in silico designed multi-epitope peptide vaccine construct.

Ramachandran diagram or Ramachandran plot is a method of visualizing Phi or φ (the bonds between nitrogen and carbon alpha) and Psi or ψ (the bond between carbon and carbon alpha) of a given polypeptide [75] . The Phi and Psi are also known as Ramachandran or dihederal angles and hold control local structural importance vital for protein folding. Ramachandran plot is a vital assessment factor for evaluating three dimensional structure of a protein. This plot can be divided into four Quadrants (from I to IV) that give important information about secondary structure of a protein. Quadrant I represent conformations that are left handed alpha helices. Quadrant II mostly accommodates beta strands and is the region for most favorable conformations and biggest among all of the regions. Quadrant III is the following biggest region after II and is the place for right handed alpha helices. Lastly, Quadrant IV is the disfavored region and is the area for steric clash conformations. The Ramachandran plot statistics are given in Table 4 . It can be observed that majority of the protein residues i.e. 527 that make 80.5% of the total protein residues are in the most favored regions. This means that most of the residues are not in correct and usual geometry. In count to these regions, additionally allowed and generously allowed regions covers 122 residues (18.6%) and 4 (0.6%), respectively. Only 2 residues (0.3%) were found in the disallowed regions having conformation φ around 0-180 degrees and ψ around -180 to 0 degrees. The number of glycine and proline residues are 123 and 170, respectively. Proline is frequently found to be present in the outlier region of the Ramachandran plot and there is a reason to explain this. The cyclic side chain of proline face J o u r n a l P r e -p r o o f constrains upon rotation by its inclusion in pyrrolidine ring thus Phi values are confined around -60. Glycine, on the contrast, lacks beta carbon atom and therefore least steric compared to other amino acids. This support the large area distribution of glycine on the Ramachandran plot and allows conformation forbidden to other amino acids: allowing their occurrence in protein turn regions where others would be sterically hindered. The different parameters of G-factor of the vaccine construct are also provided in the Table 4 . G-factor provides information about geometry of proteins that how unusual it is. Values below -0.5 are considered unusual, while those having values less -1.0 are highly unusual. The structure was immediately submitted to GalaxyRefine server to rebuilds the side chains and run side chains repacking followed by structure relaxation [76] . The Ramachandran for the refined vaccine construct is provided in S- Fig.2 . It is quite clear that the structure quality is improved by increasing the percentage of residues in the most favored region. Also, the G-factor value of the structure is improved to 0.00. Table 4 . Structure assessment of multi-epitope peptide vaccine construct. Overall average -0.31/0.00

Disulphide engineering is a key approach in biotechnology to increase thermal stability of protein under study and assist in understanding its dynamics [77] . The introduction of disulphide bonds in the protein is believed to lower the conformational entropy and at the same time increases the free energy of the denatured state [78] . Thus both these actions lead to overall stability of protein fold conformation. Total of 70 pair of residues were found that could be manipulated for the purpose of increasing protein stability. However, only pair of residues with combined energy of > 2 kcal/mol were opted. The complete set of residues with potential to be mutated is listed in Table 5 . The wild and mutated multi-epitope vaccine construct are provided in Fig.3 . 

The codon adaptation tool is a simple way to adjust the codon usage of the input sequence to the selected prokaryotic and eukaryotic system [79, 80] . The difference in codon usage of the host from that organism from which the desired sequence stems from plays a major role in minor expression of that sequence. Codon adaptation of the vaccine construct primary sequence to E. coli 

Molecular docking of the designed vaccine construct with TLR3 as a receptor was accomplished using PatchDock, which generated more than 20 complexes of TLR3 and vaccine construct conformations. The top 10 best predicted solutions were subjected to the fast interaction refinement in molecular docking (firedock) to refine the predicted solutions based the global energy. Lower the global energy of a solution implies the stable formation of complex and vice versa. As such solution 3 with global energy of -24.02 kcal/mol, attractive van der Waals energy of -38.57 kcal/mol, repulsive van der Waals energy of 28.88 kcal/mol, atomic contact energy of 11.49 kcal/mol and hydrogen bonding energy of -4.00 kcal/mol was considered the best predicted complex ( Table 6 ). The docked conformation of the TLR3 and construct with respect to each is shown in Fig.5 . The top complex conformation of the TLR3-MEVC was also understand via ClusPro that complement good binding of the molecules with lowest energy of -1268.5 kcal/mol. J o u r n a l P r e -p r o o f 

The vaccine construct appears to produce robust IgM antibody response even in the absence of recurrent antigen exposure. The IgM antibody titer seen reach close to 700000 scales (Fig.6 ). Among the cytokines and interleukins, high level interferon gamma was reported to the antigen compared to others (Fig.6) . Additionally, appropriate and decent cellular activity with their respective memory from the host immune system were revealed along with proliferating dendritic cells and enhanced activity from macrophages (S- Fig.4 ). Different humoral immunity isotypes can be easily noticed, an indicator of isotype switching and memory cell formation (S- Fig.5 ). This immune profile leads to a conclusion of a good agreement on natural development of immunity against the virus and protection against the virus. 

The multi-epitope peptide vaccine construct complex was seen stably docked with the TLR3 receptor and the epitopes are exposed throughout the simulation time for efficient and easy processing by the host immune system. The dynamics of the complex was evaluated through two structural parameters. First, root mean square deviation (rmsd) was investigated that measures the distance between the carbon alpha atoms of superimposed proteins [81] . The average rmsd estimated for the system is 5.2 Å. The system is reported with some major rmsd peaks that after visual inspection found because of the moving loops of the system. This may be an approach towards complex stability as both the interacting molecules are stably docked with each other's. Secondly, we evaluted root mean square fluctuations (rmsf) for the system. RMSF sheds light on the fluctuating residues from there mean position [82] . The average rmsf for the system is 4.6 Å. As mentioned earlier, these variation are because of fluctuating loops of the system which is not altering intermolecular conformations and does not affect overall stability. The rmsd and rmsf plots of the system can be seen in Fig.7 . Additionally, we performed hydrogen bond analysis for J o u r n a l P r e -p r o o f the complex to decipher the strength of intermolecular interactions during simulation time. As can be seen in Fig.8 , the molecules are involved in good number of hydrogen bonds throughout the simulation time, enhancing stability of the complex. 

The docking prediction for the vaccine construct conformation with respect to the TLR3 receptor was validated by performing computationally cost effective and sufficiently reliable MMMGBSA and MMPBSA approaches. Both these approaches were run over trajectories of the MD simulations. It is evident from the Table 7 that robust energies are noticed for the system and the values are quite low indicating a good agreement on the docking findings. The intermolecular ineractions are dominated by both van der Waals and electrostatic energies with energy value of -188.44 kcal/mol and -473.80 kcal/mol, respectively. Summing, both these mentioned values are contributing factor to the significantly low gas phase energy of -626.25 kcal/mol. Opposed to this, non-favorable energy of solvation was noticed for the system (564.3822 kcal/mol in MM-GBSA and 527.0717 kcal/mol in MM-PBSA). The net binding energy of the system in MM-GBSA is -61.8682 kcal/mol whereas in MM-PBSA this energy is -J o u r n a l P r e -p r o o f 99.1787 kcal/mol. Both these net values support highly stable system and excellent intermolecular binding affinity between the vaccine construct and TLR3 receptor. The net entropy concluded for the system is 57.76 kcal/mol which again suggesting stable nature of the vaccine at the docked site. Additionally, to evaluate intermolecular affinity towards end of the simulation time, MMPB/GBSA of last 20-ns was performed (S- Table 1 ). The net MMGBSA and MMPBSA in this time period are -77.4 kcal/mol and -107.4 kcal/mol, respectively. 

This study shows that membrane protein, envelope small membrane protein, non-structural protein ORF3, non-structural protein ORF5 and spike glycoprotein could act as potential subunit vaccine candidates for designing a subunit-based vaccine against MERS-CoV. Similarly, these proteins were mapped for novel set of antigenic B and T cell epitopes that could be used in the design of a multi-epitope peptide vaccine construct to circumvent the limitations of subunit vaccines. Further, we have shown that the designed multi-epitope vaccine construct docking affinity for TLR3 immune receptor, thus maximizing its presentation to the human immune system and ensuring the evoking potential of both humoral and cellular immune responses. The prediction of docking was validated by molecular dynamics simulation and binding free energies that validated complex high stability and intermolecular affinity. Lastly, the use of the subunit vaccine candidates and multi-epitope vaccine construct in experimental in vivo model is highly 

Global status of Middle East respiratory syndrome coronavirus in dromedary camels: a systematic review

Middle East respiratory syndrome: emergence of a pathogenic human coronavirus

Diab, others, Middle East respiratory syndrome coronavirus in dromedary camels: an outbreak investigation

Boudjelal, others, High prevalence of MERS-CoV infection in camel workers in Saudi Arabia

Middle East respiratory syndrome coronavirus: another zoonotic betacoronavirus causing SARSlike disease

MERS-CoV spillover at the camelhuman interface

MERS epidemiological investigation to detect potential mode of transmission in the 178th MERS

Human--dromedary camel interactions and the risk of acquiring zoonotic Middle East respiratory syndrome coronavirus infection

Temporal dynamics of Middle East respiratory syndrome coronavirus in the Arabian Peninsula

Treatment strategies for Middle East respiratory syndrome coronavirus

Middle East Respiratory Syndrome Vaccine Candidates: Cautious Optimism

An updated roadmap for MERS-CoV research and product development: focus on diagnostics

Searching for animal models and potential target species for emerging pathogens: Experience gained from Middle East respiratory syndrome (MERS) coronavirus, One Heal

Prospects for a MERS-CoV spike vaccine

others, Characterization of novel monoclonal antibodies against MERS-coronavirus spike protein

MERS coronavirus envelope protein has a single transmembrane domain that forms pentameric ion channels

Middle East respiratory syndrome coronavirus M protein suppresses type I interferon expression through the inhibition of TBK1-dependent phosphorylation of IRF3

Peptide vaccine: progress and challenges

Peptide-based synthetic vaccines

Design and development of synthetic peptide vaccines: past, present and future

Computer aided epitope against Lassa virus

Development of the schedule for multiple parallel "difficult" Peptide synthesis on pins

Granadillo, others, Safety and immunogenicity of a human papillomavirus peptide vaccine (CIGB-228) in women with high-grade cervical intraepithelial neoplasia: first-in-human, proof-of-concept trial

McElrath, others, Evolution of human immunodeficiency virus type 1 cytotoxic T-lymphocyte epitopes: fitness-balanced escape

Heterosubtypic protective immunity against influenza A virus induced by fusion peptide of the hemagglutinin in comparison to ectodomain of M2 protein

Multi-epitope vaccines: a promising strategy against tumors and viral infections

UniProt: a hub for protein information

TiD: Standalone software for mining putative drug targets from bacterial proteome

CELLO2GO: a web server for protein subCELlular LOcalization prediction with functional gene ontology annotation

CD-HIT: accelerated for clustering the nextgeneration sequencing data

Basic local alignment search tool

Database resources of the national center for biotechnology information

From probiotics to psychobiotics: Live beneficial bacteria which act on the brain-gut axis

Subtractive genome analysis for in silico identification and characterization of novel drug targets in Streptococcus pneumonia strain JJA

CD8+ Cytotoxic-T-Lymphocyte Breadth Could Facilitate Early Immune Detection of Immunodeficiency Virus-Derived Epitopes with Limited Expression Levels

Large-scale validation of methods for cytotoxic T-lymphocyte epitope prediction

VaxiJen: a server for prediction of protective antigens, tumour antigens and subunit vaccines

CD4+ T cells: differentiation and functions

Development of a multi-epitope peptide vaccine inducing robust T cell responses against brucellosis using immunoinformatics based approaches

The immune epitope database (IEDB): 2018 update

AllergenFP: allergenicity prediction by descriptor fingerprints

Consortium, others, In silico approach for predicting toxicity of peptides and proteins

Designing of interferon-gamma inducing MHC class-II binders

Avian antimicrobial peptides: the defense role of $β$-defensins

Template-based protein structure modeling using the RaptorX web server

I-TASSER server for protein 3D structure prediction

PDBsum: summaries and analyses of PDB structures

AlgPred: prediction of allergenic proteins and mapping of IgE epitopes

Fundamentals and methods for T-and B-cell epitope prediction

BepiPred-2.0: improving sequencebased B-cell epitope prediction using conformational epitopes

PatchDock and SymmDock: servers for rigid and symmetric docking

FireDock: fast interaction refinement in molecular docking

The ClusPro web server for protein--protein docking

Computational immunology meets bioinformatics: the use of prediction tools for molecular binding in the simulation of the immune system

The FF14SB force field

Langevin stabilization of molecular dynamics

CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data

MMPBSA.py: An efficient program for end-state free energy calculations

The MM/PBSA and MM/GBSA methods to estimate ligandbinding affinities

Interaction entropy: a new paradigm for highly efficient and reliable computation of protein--ligand binding free energy

Exoproteome and secretome derived broad spectrum novel drug and vaccine candidates in Vibrio cholerae targeted by Piper betel derived compounds

Barh, others, Pan-genome analysis of human gastric pathogen H. pylori: comparative genomics and pathogenomics approaches to identify regions associated with prediction of potential core therapeutic targets

Identification of putative vaccine candidates against Helicobacter pylori exploiting exoproteome and secretome: a reverse vaccinology based approach

The CD8 T cell response to respiratory virus infections

A review on T cell epitopes identified using prediction and cell-mediated immune models for mycobacterium tuberculosis and Bordetella pertussis

Helper T cells and lymphocyte activation

Design, structure prediction and molecular dynamics simulation of a fusion construct containing malaria pre-erythrocytic vaccine candidate, PfCelTOS, and human interleukin 2 as adjuvant

Vaccine allergies

Pangenome and immuno-proteomics analysis of Acinetobacter baumannii strains revealed the core peptide vaccine targets

Correlation between stability of a protein and its dipeptide composition: a novel approach for predicting in vivo stability of a protein from its primary sequence

Thermostability and aliphatic index of globular proteins

Objectively judging the quality of a protein structure from a Ramachandran plot

GalaxyRefine: protein structure refinement driven by side-chain repacking

Protein disulfide engineering

Disulfide by Design 2.0: a web-based tool for disulfide engineering in proteins

Codon usage: nature's roadmap to expression and folding of proteins

JCat: a novel tool to adapt codon usage of a target gene to its potential expression host

Significance of root-mean-square deviation in comparing three-dimensional structures of globular proteins

Determination of ensemble-average pairwise root mean-square deviation from experimental B-factors

Original Draft Preparation, Muhammad Arshad, Abdul Haleem, Saba Ismail, Anita Zaib: Data Curation

The authors declare no conflict of interest. Fig.1 . Predicted B cell epitope region in the vaccine construct. Fig.2 . Ramachandran plot for the 3D structure of multi-epitope peptide vaccine construct. Fig.3 . Codon optimization score of the vaccine construct. Fig.4 . T cell immunity profile against the vaccine construct antigen. Fig.5 . B cell immunity profile against the vaccine construct antigen. Table 1 . Binding free energies for the TLR3-MEPVC complex.

The authors declare that they have no conflict of interest.J o u r n a l P r e -p r o o f