The surface glycoprotein sequence of the pneumonia virus discovered at the Wuhan seafood market (QHD43416.1 from MN908947.3 reference genome) was retrieved from NCBI (22). To scrutinize required HLA binding epitopes, a TepiTool from IEDB was used (23). A set of 12 MHC class I super-types (A^*01:01, A^*02:01,A^*03:01, A^*24:02,A^*26:01, B^*07:02, B^*08:01, B^*27:05, B^*39:01, B^*40:01, B^*58:01, and B^*15:01) were used, and the two highest-scoring epitopes (based on percentile rank and IC50 values) for each allele were selected. A percentile rank is calculated by the comparison of the peptide's predicted binding-affinity against a panel of a variety of peptides randomly selected from the Swiss-Prot. Hence, a lower percentile rank numerical value depicts better binders. Additionally, all the predicted peptides were checked for their IC50 value, and those with IC50 ≤ 500 nM were taken into account. Specific immune responses are based on CD4⁺ and CD8⁺ T cells, and protective vaccines should thus induce specific T-cell responses based on peptides represented by MHC-I and MHC-II alleles. The rationale behind prioritizing HLA binding epitopes is to ensure the specific immune response in infected macrophages.

For MHC-II-binding peptide epitopes, the seven-allele method was used. This selection is based on the median of consensus percentile ranks among the seven commonly encountered DR alleles, namely, HLA-DRB1^*03:01, HLA-DRB1^*07:01, HLA-DRB1^*15:01, HLA-DRB3^*01:01, HLA-DRB3^*02:02, HLA-DRB4^*01:01, and HLA-DRB5^*01:01. Epitopes with a median consensus percentile rank ≤ 20.0 were designated as good binders. The scrutiny of promiscuous peptide epitopes was established based on the median of the consensus percentile rank of the seven preselected alleles.

For B-cell epitope prediction, BepiPred 2.0 from Immune Epitope Database Analysis Resource (IEDB-AR) was used (24). IEDB-AR is linked to IEDB and offers computational analysis regarding both B and T cell epitope prediction and their subsequent analysis. BepiPred 2.0 works on the basis of a randomly chosen forest algorithm that has been trained on epitopes acquired from antibody–antigen models obtained from interactive protein structures (25).

Owing to the significance of spike protein, the selected epitopes were manually screened for their presence in this zone. The epitopes were further examined for antigenic potential via VaxiJen version 2.0 (26). A threshold value of 0.5 was taken into account. Non-antigenic peptides (having VaxiJen score <0.5) were discarded, while antigenic epitopes (with threshold value > 0.5) were further prioritized for their immunogenicity. The Immune epitope database (IEDB) tool for immunogenicity score calculation was used to predict immunogenicity scores for all MHC-I predicted epitopes (27). This tool is designed to predict immunogenicity of the peptide based on amino-acid position and properties. Immunogenic epitopes were then verified for their presence in IEDB database.

Shortlisted top-scoring epitopes were checked for their binding affinity with each other for determining the final sequence of the chimeric vaccine. The epitopes were analyzed using a HADDOCK web server (Guru interface) (28). Clusters representing two epitopes, which possessed the highest interaction scores, depicting their maximum interaction, were refined by removing the water molecules, which may hinder their interaction, and then having them dock to the third epitope. Likewise, evaluation of clusters with three epitopes was done. The refined and the highest-scoring cluster was docked to the fourth epitope to obtain the final sequence.

To facilitate the process of vaccine development, a flexible linker GGGGS was added between each epitope. This helps to restore protein folding by allowing interaction between different domains (16, 29). Additionally, another linker EAAAK was added at the N terminal to separate bi-functional domains. Designed vaccines were then tested with different epitopes, including Truncated Ov-ASP-1 Protein (residues 10–153) and Beta defensin (45 residues long), and constructs having higher antigenicity and that are predicted to produce high antibody titers were added with the multi epitope vaccine construct to the enhance immune response (30). Three different constructs were designed in this study, one comprising the top-scoring CD4 and CD8 epitopes lying in the S1 domain, while another is formed by taking two epitopes from the S1 domain and two from the S2 domain, representing MHC-I and MHC-II binders. Finally, the third one is formed by adding a B-cell epitope to the second one but with a different adjuvant.

The final sequences of the chimeric vaccine constructs were screened for its antigenic potential and solubility using ANTIGENpro and SOLpro (31). Allertop version 2.0 was used to check the probability of the construct to cause an allergic reaction (32). Sequence of the finalized vaccine candidate in FASTA format was given as an input to ExPASy server, in order to calculate various parameters like molecular weight, theoretical PI, half-life of the protein, instability index, amino acid composition, aliphatic index, and GRAVY (33).

Secondary structures of the vaccine constructs were predicted using PDBsum (34). This step was executed to better understand the structures of predicted vaccines. PDBsum is a database that is exclusively designed to show the molecules that build DNA or proteins, ligands, and metal ions along with the illustration of graphical representation of their interactions with each other. To generate 3D structures of the vaccine candidates, 3Dpro was used (31). The predicted models were then refined using Galaxy refine server (35). This server is responsible for subjecting the predicted 3D model to structural perturbations and subsequent structural relaxations. It generates five different models. All five models for each vaccine construct were screened for GDT-HA, RMSD, and poor rotamers, and the finest predicted models were taken to the next step.

The finalized models were further evaluated using ERRAT scores and Ramachandran plot analysis for verification. In order to obtain stabilized vaccine constructs, energy minimization was carried out using online YASARA server. YASARA deals with molecular-dynamics simulations of the given models in solvent, using an exclusive forcefield that has been derived from Amber, whose constraints have been improved to minimalize the impairment done to protein structure during the process of energy minimization (36).

In order to study the binding affinity of the putative vaccine candidates with immune receptors, molecular docking technique was adopted. Prioritized vaccine constructs were docked to ACE2 receptor (PDB ID: 3sci), TLR2 (PDB ID: 2Z7X), and TLR4 (PDB ID: 4G8A). Vaccine 3, having a B-cell epitope, was also checked for its interaction with a B-cell receptor (BCR) CD79 (PDB ID: 3KG5). For this protein–protein docking validation process, the HAADOCK server (Guru Interface and refinement interface) was used (28). Additionally, to obtain a graphical illustration of the interactions between vaccine and receptor, PDBsum was used (34). Moreover, in order to verify the binding affinity of our multiepitope peptide vaccines with HLA alleles, all our vaccine constructs were docked with class I and class II Superfamily alleles to reveal the interaction of epitopes with MHC alleles when combined as well. Hence, for this purpose, class I [HLA A^*02 01 (PDB ID 4U6Y), HLA B^*51 01 (PDB ID 4MJI)] and class II [HLA-DRB1^*1402 (PDB ID 6ATF)] were used; they represent broad-spectrum peptide-binding repertoires.

Population coverage of epitopes was determined using IEDB for prioritized epitopes, as it helps to determine the percentage population that can respond to the particular epitope and can elicit an immune response against it.

The finalized epitopes in Table 1 were examined for their interactive ability with one another using HADDOCK. All possible combinations of epitopes along with a flexible linker GGGGS between them were explored. For vaccine 1, the E1S1–E4S1 combination had the highest haddock refinement score.

The binding affinity of the E1S1–E4 S1 combination with the other two epitopes was determined to find the best combination of three epitopes. E1 S1 –E4 S1 –E3 S1 was thus formed. Finally, the vaccine construct obtained after combination analysis was E1 S1–E4 S1–E3 S1–E2 S1 (Table 3). Similarly, for vaccine 2, E1- S1 –E4 S1 was the first combination, and it was followed by E1S1 –E4 S1–E2 S2 and E1 S1–E4 S1–E2 S2 –E1 S2. Each probable combination lined up for putative vaccine design along with their corresponding HADDOCK scores is present in Table 3. Moreover, truncated Ov-Asp1 (IVVAVTGYNCPGGKLTALERKKIVGQNNKYRSDLINGKLKNRNGTYMPRGKNMLELTWDCKLESSAQRWANQCIFGHSPRQQREGVGENVYAYWSSVSVEGLKKTAGTDAGKSWWSKLPKLYENNPSNNMTWKVAGQGVLHFTQ) was attached to the N terminal of both the putative vaccines using another linker, EAAAK. The finalized vaccines together with the linkers and adjuvant were 212 amino acids long. Ov-ASP-1 reportedly has ability to activate antigen-processing cells (APCs) which define its good adjuvanticity for a number of vaccines and antigens (30). They are thus added in vaccine constructs to improve the efficacy of these new generation subunit vaccines.

In order to ensure both cell and humoral mediated responses, a potent B-cell epitope was added to vaccine 2 based on the best docking scores predicting the combination pattern of epitopes. Vaccine 3 was created with an order; E4 S1–E5 S1–E2 S2–E1 S2–E1 S1 and the corresponding docking score are enlisted in Table 3. For comparison purposes, another adjuvant beta-defensin (GIINTLQKYYCRVRGGRCAVLSCLPKEEQIGKCSTRGRKCCRRKK) was added to this combination. Beta defensin has previously been reported as a potent adjuvant when conjugated with MERS-CoV antigens (38). Vaccines containing defensins as adjuvants have been shown, both in vivo and in vitro, to activate the primary innate antiviral immune response and mediate other immunomodulatory activities against a number of viruses, including coronaviruses (38, 39). Vaccine 3, after addition of this adjuvant at the N terminal along with EAAAAK and GGGGS linkers, consisted of 143 residues. The final combination of epitopes of all three vaccine constructs have been shown in Figure 1.

Various physiochemical properties were examined for both the constructs. The molecular weight of vaccine 1 is 23235.26 g/mol while the theoretical pI is 9.50, depicting the basic nature of the peptide construct. The instability index II showed that the construct is stable with a score of 24.79. The GRAVY (GRand AVerage of hydropathY) index was calculated to be −0.479, validating the hydrophilic nature of the construct that can form interactions with surrounding water molecules. The aliphatic index 67.12 illustrated that the construct is thermostable in nature.

Vaccine 2 has a molecular weight of 23013.07 g/mol, and its theoretical pI is 9.33. Hence, this construct was also found to be basic in nature. Likewise, instability analysis showed that the protein is stable with a score of 24.50. The GRAVY index testified the hydrophilic nature of this construct as well (−0.492). The thermostable nature of the construct was established by the value of aliphatic index, 62.97. The predicted values of antigenicity for both the vaccines were found to be 0.883591 and 0.946425, respectively. This ensured highly antigenic nature of the constructs. Similarly, the solubility upon overexpression was predicted to be 0.864955 and 0.951926. Furthermore, both vaccine constructs designed in this study were designated as non-allergenic by AllergenPro.

Vaccine 3 has a molecular weight of 15084.28 g/mol and its theoretical pI is 9.25. Therefore, this vaccine construct was also found to be basic in nature. The GRAVY index testified the hydrophilic nature of this construct as well (−0.253). The thermostable nature of the construct was established by the value of aliphatic index, 55.87. The predicted values of antigenicity for this particular the vaccine was 0.883570. This ensured highly antigenic nature of the construct. Similarly, the solubility upon overexpression was predicted to be 0.806206. Furthermore, like both the previous vaccine constructs designed in this study, this vaccine was also found to be non-allergenic by AllergenPro.

The secondary structure of vaccine 1 includes six helices, 35 beta turns, seven gamma turns, and nine helix-helix interactions. The secondary structure of vaccine 2 has eight helices, 22 beta turns, 12 gamma turns, and nine helix–helix interactions. For vaccine 3, secondary structure consisted of two beta strand, one hairpin, one sheet, four helices, 23 beta turns, 23 gamma turns, and one helix–helix interaction. Helix–helix interaction presents facts about different pairs of helices, interacting with each other with the vicinity of the protein structure, whereas beta turns depict four consecutive residues. These four residues are represented by i, i + 1, i + 2, and i + 3. This is possible when the measured distance between the alpha Carbon atom of the first residue (i) and alpha carbon atom of the fourth residue (i + 3) is <7 Å plus the two residues between them are not helical. A gamma turn comprises of three residues i, i + 1, and i + 2. This is possible when a hydrogen bond is present between the two residues (i.e., i and i + 2). Moreover, the phi angle and the psi angle of the second residue i.e., i + 1 lies within a range of 40 degrees in one of the next two cases: (1) classic [phi i + 1(75), psi i + 1(−64)] or (2) inverse [phi i + 1(−79), psi i + 1(−69)].

The 3Dpro tool, which works on the basis of ab initio method for predicting tertiary structure, was used to predict three dimensional structures of proposed vaccine constructs. This strategy was adopted due to the lack of fine homolog proteins that could be exploited for homology modeling. The obtained models were then refined via several structure perturbations and subsequent structure relaxations using GlaxyRefine server. The obtained best models are shown in Figure 2. The ERRAT score for 3D models of three vaccines were calculated as 74.1379, 67.5676, and 74.2574, respectively. While Ramachandran plot analysis showed 97.1% residues in favored region for vaccine 1, 98.1% residues in the favored region for vaccine 2 and 86.5% for vaccine 3 (Figure 2). These analyses authenticated the reliability and stability of the predicted structures.

Energy minimization by a YASARA server was performed. For vaccine 1, the YASARA force field was applied to 2,032 atoms. A total of 5,282 water molecules were found. The initial energy was −68794.3kJ/mol (Z score −1.90), which was minimized to −97974.1 kJ/mol (−1.93). For vaccine 2, the YASARA force field was applied to 2,006 atoms while the water molecules were 5,208. Initial energy was −66687.0 kJ/mol (Z score −2.08); however, the final energy was 101214.7 kJ/mol (Z score −1.47). For vaccine 3, the YASARA force field was applied to 2,093 atoms while the water molecules were 4,185. Initial energy was −53609.9 kJ/mol (Z score −3.35); however, the final energy was −60374.6 kJ/mol (Z score −3.16).

SARS-CoV spike protein has been studied previously for its exceptional binding affinity with human ACE-2. It should be noted that, structurally, SARS-CoV-2 and SARS-CoV spike proteins are highly homologous in nature, sharing 76.5% identical amino acids. Atomic level studies between SARS-CoV and ACE-2 show promising interactions between the two, and therefore, owing to the structural and sequence similarity, it is anticipated that an ACE-2 blocker might be handy in curbing SAR-CoV-2 (40). For vaccine 1 and the ACE-2 receptor, therefore, docking was carried out. A HADDOCK server clustered 36 probable structures into seven different clusters, which represented a total of 18.0 % of the water-refined models. The top-scoring cluster had a score of 39.8 +/– 29.1 and a Z score of −1.6. Similarly, for vaccine 2 and the ACE-2 receptor, HADDOCK clustered 18 structures in three clusters, which represented 9.0% of the water-refined models. The top-scoring cluster had a value of 0.3 +/– 9.8 and a Z score of −1.3. Likewise, for vaccine 3, HADDOCK clustered 22 structures into five clusters, which depicted 11% of the water refined models generated by HADDOCK. Here, the best cluster had a score of 147.5 +/– 15.0 and a Z score of −1.2 (Figure 3).

TLR2 and TLR4 are well-studied Toll-Like Receptors that identify both structural and non-structural proteins of the virus and subsequent cytokine production and inflammation. They are present on the surface of cells and are triggered by viral glycoproteins. TLR agonists have the potential to initiate an immune response and actively participate in viral clearance (41). The prioritized vaccine constructs were therefore also explored for their interaction with Toll-like receptors TLR2 and TLR4. Vaccine 1 and TLR2 interaction revealed 40 structures in a total of six clusters that represented 20.0% of the water-refined models. The model with the highest score, −4.2 +/– 20.8 had a Z value of −1.2. Likewise, for vaccine 2 and TLR2, HADDOCK clustered 80 structures in 12 clusters, which represented 40.0% of the water-refined models. Here, the highest-scoring model had a score of −23.7 +/– 12.1 with a Z-value of −1.3. For vaccine 3 and TLR2, HADDOCK clustered 136 structures in 10 clusters, which represented 68.0% of the water-refined models. The highest scoring model had a score of −16.7 +/– 14.0 with a Z-value of −1.8.

Moreover, HADDOCK clustered 157 structures in 13 clusters to determine vaccine 1 and TLR4 interaction, which represented 78.5% of the water-refined models. The top-scoring model had a score of 37.9 +/– 7.8 and a Z-value of −2.2, whereas the interaction of vaccine 2 and TLR4 is determined by 47 structures in nine cluster(s), which represents 23.5% of the water-refined models. The top-scoring model had a score of −16.8 +/– 23.4 (Z-value −1.6). Similarly, HADDOCK clustered 93 structures in eight clusters to determine vaccine 3 and TLR4 interaction, which represented 46.5 % of the water-refined models. The top-scoring model had a score of 23.3 +/– 5.7 and a Z-value of −1.3.

Models from top clusters were refined using HADDOCK refinement interface. This server was used to cluster 20 structures, obtained via HADDOCK, into one cluster. This final cluster symbolized 100% of water-refined models that were generated by HADDOCK. The statistics observed in interactions of vaccine 1, vaccine 2, and vaccine 3 from their refined clusters can be seen in Table 4, and complexes are shown in Figure 4.

The PDBsum analysis of vaccine 1 with ACE2 showed 18 hydrogen bonds and one salt bridge. Additionally, 42 interface residues of vaccine 1, representing an interface area of 2,502 (A²), were found while the corresponding ACE2 had 45 interface residues covering an area of 2,319 (A²). For vaccine 2 and ACE2, there were two salt bridges and seven hydrogen bonds predicted by PDBsum. Additionally, 28 and 38 residues from vaccine 2 and ACE2 interacted with each other covering an area of 1,997 and 1,832, respectively. Likewise, for vaccine 3 there was one salt bridge and 13 hydrogen bonds predicted by PDBsum. Additionally, 27 and 22 residues from vaccine 3 and ACE2 interacted with each other, covering an area of 1,228 and 1,271, respectively.

An interaction analysis of vaccine 1 with the TLR2 interacting complex via PDBsum exhibited 19 hydrogen bonds and one salt bridge. Furthermore, 35 interface residues of vaccine 1, representing an interface area of 1,795 (A²), were found while a corresponding TLR2 had 36 interface residues encompassing an area of 1,880 (A²). For vaccine 2 and TLR2, there were two salt bridges and 14 hydrogen bonds predicted by PDBsum. Additionally, 23 and 25 residues from vaccine 2 and TLR2 interacted with each other, covering an area of 1,362 and 1,443, respectively. Lastly, for vaccine 3 and TLR2 PDBsum, 17 hydrogen bonds and five salt bridges were found. Furthermore, 25 interface residues of vaccine 3, representing an interface area of 1,104 (A²), were found while corresponding a TLR2 had 21 interface residues, encompassing an area of 1,194 (A²).

Similarly, the interaction of vaccine 1 with TLR4 exhibited eight hydrogen bonds and 18 interface residues of vaccine 1, representing an interface area of 1,171 (A²) while a corresponding TLR4 had 19 interface residues, encompassing an area of 1,146 (A²). For vaccine 2 and TLR4, there were three salt bridges and 20 hydrogen bonds predicted by PDBsum. Additionally, 33 and 34 residues from vaccine 2 and TLR4 interacted with each other, covering an area of 1,763 and 1,745, respectively. In case of vaccine 3 and TLR4, 15 hydrogen bonds and three salt bridges were found, while 34 interface residues of vaccine 3 represented an interface area of 1,397 (A²) and a corresponding TLR4 had 27 interface residues, encompassing an area of 1,396 (A²).

For the interaction analysis of vaccine 3 and BCR (CD79), the HADDOCK server clustered 140 probable structures into 13 different clusters, which represented a total of 70% of the water-refined models. The top scoring cluster had a score of −43.6 +/– 16.0 and a Z score of −1.7. Models from top clusters were refined using HADDOCK refinement interface. This server was used to cluster 20 structures, obtained via HADDOCK, into one cluster. This final cluster symbolized 100% of water-refined models that were generated by HADDOCK. The statistics observed in interactions of vaccine 3 and BCR from its particular refined clusters can be seen in Table 4. Pdbsum analysis showed that 26 and 18 residues from vaccine 3 and BCR interacted with each other covering an interface area (A2) of 1,181 and 1,205, respectively. They formed one salt bridge and 10 hydrogen bonds.

For interaction analysis of vaccine 1 and HLA A allele, the HADDOCK server clustered 118 probable structures into 17 different clusters, which represented a total of 59.0 % of the water-refined models. The top scoring cluster had a score of −26.5 +/– 2.7 and a Z score of −2.5. Similarly, for vaccine 2 and the HLA A allele, HADDOCK clustered 97 structures in 17 clusters, which represented 48.5 % of the water-refined models. The top scoring cluster had a value −57.5 +/– 12.8 and a Z score of −2.3. Likewise, for vaccine 3 HADDOCK clustered 187 structures into three clusters, which depicted 93.5% of the water refined models generated by HADDOCK. Here the best cluster had a score of −34.7 +/– 1.9 and a Z score of −1.1. For vaccine 1 and HLA B allele, 115 probable structures were clustered by HADDOCK into 15 different clusters, which represented a total of 57.5 % of the water-refined models. The top scoring cluster had a score of −57.5 +/– 12.8 and a Z score of −2.3. Similarly, for vaccine 2 and the HLA B allele, HADDOCK clustered 84 structures into nine clusters, which represented 42% of the water-refined models. The top scoring cluster had score of −18.7 +/– 8.7 and Z score of −1.6. Likewise, for vaccine 3 HADDOCK clustered 168 structures into 10 clusters, which depicted 84% of the water refined models were generated. Here, the best cluster had a score of −41.2 +/– 18.7 and a Z score of −2.1.

Furthermore, for vaccine 1 and the HLA DRB1 allele docking, the HADDOCK server clustered 67 probable structures into 10 different clusters, which represented a total of 33.5 % of the water-refined models. The top-scoring cluster had a score of −27.8 +/– 6.0 and a Z score of −2.3. Similarly, for vaccine 2 and HLA DRB1 allele, HADDOCK clustered 64 structures in 11 clusters, which represented 32% of the water-refined models. The top-scoring cluster had score of −24.8 +/– 25.6 and Z score of −1.7. Likewise, for vaccine 3, HADDOCK clustered 93 structures into 13 clusters, which depicted 46.5% of the water refined models generated by HADDOCK. Here the best cluster had a score of −37.1 +/– 11.8 and a Z score of −1.5. Models from top clusters were refined using HADDOCK refinement interface. This server was used to cluster 20 structures, obtained via HADDOCK, into one cluster. This final cluster symbolized 100% of water-refined models that were generated by HADDOCK. The statistics observed in interactions of vaccine 1, vaccine 2, and vaccine 3 from their particular refined clusters can be seen in Table S4.

Epitope population coverage was checked by IEDB population coverage tool. Resultantly, all epitopes had a combined Class I and Class two average coverage score of 94%. This step was performed by using the entire world population datasets and the MHC restricted alleles used in this case were (A^*01:01, A^*02:01, A^*03:01, A^*24:02, A^*26:01, B^*07:02, B^*08:01, B^*27:05, B^*39:01, B^*40:01, B^*58:01, B^*15:01, HLA-DRB1^*03:01, HLA-DRB1^*07:01, and HLA-DRB1^*15:01).