In-silico design of envelope based multi-epitope vaccine candidate against Kyasanur forest disease virus

Kyasanur forest disease virus (KFDV) causing tick-borne hemorrhagic fever which was earlier endemic to western Ghats, southern India, it is now encroaching into new geographic regions, but there is no approved medicine or effective vaccine against this deadly disease. In this study, we did in-silico design of multi-epitope subunit vaccine for KFDV. B-cell and T-cell epitopes were predicted from conserved regions of KFDV envelope protein and two vaccine candidates (VC1 and VC2) were constructed, those were found to be non-allergic and possess good antigenic properties, also gives cross-protection against Alkhurma hemorrhagic fever virus. The 3D structures of vaccine candidates were built and validated. Docking analysis of vaccine candidates with toll-like receptor-2 (TLR-2) by Cluspro and PatchDock revealed strong affinity between VC1 and TLR2. Ligplot tool was identified the intermolecular hydrogen bonds between vaccine candidates and TLR-2, iMOD server confirmed the stability of the docking complexes. JCAT sever ensured cloning efficiency of both vaccine constructs and in-silico cloning into pET30a (+) vector by SnapGene showed successful translation of epitope region. IMMSIM server was identified increased immunological responses. Finally, multi-epitope vaccine candidates were designed and validated their efficiency, it may pave the way for up-coming vaccine and diagnostic kit development.


Results
Conserved region of KFDV envelope protein. The multiple sequence alignment analysis (Supplementary File 1) reveals that except few amino acids, the entire KFDV envelope protein sequences remain conserve. Phylogenetic tree was constructed from the alignment file, shows separate cluster of 2012 KFDV envelope protein, rest of the sequence were almost the same (Fig. 1A). Based on the conserve score (Supplementary File 2) conservation graph has been drawn, value 1.0 on the Y-axis represents more conserved (Fig. 1B). Amino acids, which got > 0.90 conserve score was considered as a conserved region and taken into further study, conserved amino acid sequence and their position were given in Table 1.

Conserved regions of KFDV envelope protein containing T-cell and B-cell epitopes.
Here we found that the conserved region contains epitopes for T-cell and B-cell targets. Totally 35 epitopes were obtained  . IEDB-MHC-I top 10 percentile epitopes were also predicted in MAPPP tool, they are listed in Table 2. Combined predictor of proteasomal cleavage, TAP transport, and MHC process has produced 34 epitopes (Supplementary File 5), and top 10 ranked epitopes were listed in Table 3. There are 20 MHC class II epitopes (Supplementary File 6) were identified by IEDB-MHC-II tools, whereas MHC2Pred tool recognized 504 epitopes (Supplementary File 7) on the conserved region of the envelope protein. Based on IEDB = MHC-II percentile > 0.10 (FILT-1) and MHC2Pred epitopes which scored > 0.60 as well having overlap sequence with FILT-1, totally19 epitopes were selected as MHC II targets (Table 4). B-cell liner epitopes prediction tools ABCpred, BCpred and SVMTriP were identified 15, 13 and 4 epitopes respectively (Table 5).

Selected epitopes display cross-protection against Alkhumra hemorrhagic fever virus.
Uniportkb-human BLAST analysis showed similarity with four human proteins such as shish-7, Isoform 2 of Ribonuclease T2, Ribonuclease T2 and Extra-cellular ribonuclease with E-value of 2e−1, 7.7e−1, 2.1e0 and 2.1e0 respectively (Supplementary File 8). These similarity hits can be ignored, as E-value is ≥ 0.1, because as rule of thumb an E-value should be < 10-4 to assure the homology. Whereas UniProtKB BLAST revealed that selected epitopes were sharing the highest similarity with Alkhumra Hemorrhagic Fever Virus (AHFV) 20 polyprotein, E-value was 1.7e−43 (Supplementary File 9). So, the proposed vaccine candidates also confer immunity against AHFV.
Multi-epitope vaccine candidate sequences their allergenicity and antigenicity. Totally, two multiple-epitope vaccine candidates were designed from selected high ranked, multi-immunogenic and over lapping epitope sequences. Vaccine construct 1 (VC1) composed with an adjuvant protein β defensin, 6 MHC-I epitope and 11 B-cell epitopes. Vaccine constructs 2 (VC2) composed of VC1 and 4 MHC-II epitopes. Complete information of vaccine construct has given in Table 6. Vaccine candidates VC1 and VC2 were found to be non-allergic in behavior, VC1 got a higher antigenic score (0.6667) and identified as a better antigen, while VC2 antigenic score was 0.5835. 3D-strucure of designed vaccine candidates. Phyre2 server has built three-dimensional structure of multi-epitope vaccine candidates. The cryo-em structure of tick-borne encephalitis virus complex was used as template for VC1 and VC2 tertiary structure, and both the model got 100% score in terms of model confidence.  www.nature.com/scientificreports/ Initial refinement of vaccine candidates VC1 and VC2 on ModRefiner have produced refine structure with TMscore of 0.9507 and 0.9622 respectively. Further refinement of VC1 and VC2 were done with GalaxyRefine server, based on quality score the best refined structure were chosen and named as VC1R and VC2R respectively. Those were used as ligand molecules in docking study ( Fig. 2A, B), quality scores of ligands were given in Table 7.
The Ramachandran plot analysis of modeled vaccine candidates revealed that 97.4% of VC1R protein residues were in favored region and 2.6% of residues in allowed region (Fig. 2C). Similarly, 95.9% of VC2R protein residues were in favored region and 4.1% of the residues in allowed region (Fig. 2D). None of the amino acid residues were fallen in outlier region, which indicates the good quality of protein structure. ProSA-web analysis of yielded Z-score of − 2.74 and − 2.66 for vaccine candidates of VC1R and VC2R respectively. (Fig. 2E, F).

Identification of binding energy between multi-epitope vaccine candidates and TLR2. Dock-
ing analysis by Cluspro and PatchDock has revealed that both vaccine candidates (VC1 and VC2) have strong binding capacity with TLR2. Cluspro protein-protein docking result page displayed 10 best-docked confirmations, the top ranked model from balanced coefficient module was downloaded in PDB format and used for LigPlot and iMODS analysis. TLR2 with VC1 exhibited lowest energy value of − 9473.6 kcal/mol, whereas VC2 showed − 876.4 kcal/mol (Fig. 3).
PatchDock server produced transformation file (posture of 3D-rotational angle) then it was rescored and global energy was obtained by FireDock. The best docking posture was ranked based on the global energy, top ranked docking score, energy contributed by hydrogen bond and Van der waals forces were listed Table 8. The highest global energy of TLR2 with VC1 and VC2 was found to be − 55.71 and − 49.10 respectively.
Interaction and stability of TLR2 and Vaccine candidates. Interaction of TLR2 and Vaccine candidates were visualized by LigPlot tool. Hydrogen bonds were highlighted in dash line (Fig. 4A), totally seven hydrogen bonds were found between TLR2 with VC1. The receptor's residue Asn294 showed highest binding affinity with 2.49 Å distance and least binding energy was found towards Phe325 with 3.18 Å distance. VC2 construct also showed identical interaction with TLR2 residues (Fig. 4B).
Molecular dynamics simulation study was done by iMODS server to check the stability and physical movement of atoms in docking complex. Simulation were performed in normal mode analysis and obtained results of VC1-TLR2 and VC2-TLR2 docking complexes were depicted in Figs. 5 and 6 respectively. The deformability graph of both complexes was illustrated in Figs. 5B and 6B, hinges are highlighting the deformability regions in the complex. The B-factor (Figs. 5C, 6C) calculates root mean square value and represents the uncertainty of each atom in the docking complex. Eigenvalues of VC1-TLR2 and VC2-TLR2 docking complexes were found to be 1.135892 × 10 −5 and 6.807877 × 10 −6 respectively (Figs. 5D, 6D). The variance matrix graph of residues displayed in Figs. 5E and 6E. Covariance matrix indicates (Figs. 5F, 6F) coupling between pairs of residue experience correlated (red), uncorrelated (white) and anti-correlated (blue) motions. The elastic of docking complexes was shown in Fig. 5G, it explains the relation between the atoms (darker gray). iMODS simulation results are suggest that vaccine constructs and TLR2 complexes are stable.
Codon optimization and in silico cloning of KFD-VC1 and KFD-VC2. Java Codon Adaptation Tool was used to check codon optimism of vaccine candidates in Escherichia coli (strain k12) expression system. It revealed that VC1 and VC2 multi-epitope vaccine construct composed of 927 and 1050 nucleotides respectively. The Codon Adaptation Index (CAI) was observed, 1.0 for VC1 and VC2 was 0.98. The GC content of VC1 found to be 51.13% and VC2 was 52.19%. The obtained values indicate that both vaccine candidates were having clon-  www.nature.com/scientificreports/ ing efficiency. SnapGene tool was used for in-silico insertion of adapted codon sequence of vaccine construct into pET30a (+) expression vector between XhoI and NedI restriction sites and clones were obtained successfully (Fig. 7).
Immune response simulation. The immunogenic profile of constructed vaccine candidates was obtained from C-IMMSIM server. It found that our vaccine candidates could able to elicit both humoral and cellular mediated immune responses (Figs. 8,9). Simulation results showed that secondary and tertiary immune responses is higher than primary response. Antigenic molecules were found to be cleared off after three doses of vaccination, on the other hand B and T memory cell population got increased at the maximum of 500 cells/ mm 3 and 2000 cells/mm 3 respectively, this feature makes our construct VC1 and VC2 as suitable KFDV vaccine candidates. This study also found that elevated level of cytokines such as IFN-γ and IL-2, which important for inhibition of viral replication and T-cell mediated immunity. Increased trend of IgM and IgG antibody titer was observed after third injection, at the same time antigen level were decreased. Above observed immune elicit properties ensured that vaccine constructs will be effective in human subjects.

Discussion
This genome era enabled with bulk inflow of genomics and proteomics information of almost all clinically important organisms. It facilitates us to choosing the targets for drug design, identifying new strains or drug resistance, diagnosis kit development and personalized medicine etc. Similarity, computational approaches simplified vaccine development process, vaccine informatics is a fast-growing field where in-silico vaccine design ventures are possible. Previous report on KFDV immune response has suggested to target both T cell and B cells mediated     www.nature.com/scientificreports/ immune response for successful KFDV vaccine development 21 . We have applied few immune informatic tools to identify the epitope which elicit both cellular and humoral immune responses. In our study, KFDV envelope protein sequence was chosen as target for epitope prediction. The envelope protein of flaviviruses is major protective antigen and consists of three domains EI, EII and EIII. Especially EIII domain contains specific and sub-complex specific neutralizing epitopes, moreover it is easy to express by recombinant techniques 22 . Another challenge is genetic diversity of virus that resulted in immune escape, alwaysrecombinant vaccines design prefers to target conserved antigenic regions of virus, it might accelerates better immune response 23 . The KFDV vaccinated individuals were reported for disease occurrence, it may because of variations in circulating virus 24 . Therefore, we picked the conserved region of KFDV envelope protein sequence, to overcome the genetic diversity of viral strains.  Here we have identified both of them by using immuno-informatics tools. Epitopes or viral antigens arise the specific immune responses in the body, which induce adaptive immune responses of T-cells mediated cellular immunity and B-cells/antibody mediated humoral immunity 25,26 . All three MHC-I epitope prediction tools identified KTILTLGDY and KTAE-HLPKAW. Also, MHC II epitope KAWQVHRDWFEDLSL overlap with MHC class I epitope KAWQVHRDW, it may serve as target for both CD8 + and CD4 + cells. B-cell epitopes KRPPTDSGHD, TVVLELDKTAEH, DDI-HQENPAKTR, AANESHSNRKTASF, TGSKPCRIPVRAVA were identified by any of two-prediction tools out of three.
The high ranked, multi-immunogenic and overlapped epitopes were linked together and multi-epitope construct was obtained. Typical purified proteins are not inherently immunogenic, it requires adjuvant to accelerate the innate immune system and enhance vaccine potency 27 . We have added β-defensin at N terminal end of the multi-epitope vaccine construct. β-defensin is known to induce lymphokines production which promotes T-cell mediated cellular immunity and antigen-specific Ig production 28 .
Totally, we made two multi-epitopes vaccines construct against KFDV, 3D structure, and refinement of the vaccine candidates were achieved. The Ramachandran plots of built structures showed that none of the residues were in outlier region, > 95% of the amino acids fallen in favored region. It clearly indicates that modeled structures of vaccine candidates are good in quality.
TLRs plays vital role in innate immunity, especially they detect the virus and activates innate immunity followed by adaptive immunity. Several reports confirmed that TLR2 act as a host sensor to identifying viral envelop proteins and subsequently activates the innate immune system [29][30][31] . Our docking studies suggest that designed vaccine construct have the binding capacity with TLR2, both Cluspro and PatchDock have identified VC1 as potential binder with TLR2. Hydrogen bonds plays viral role in protein-receptor interaction, and it determines the stability of complex 32,33 . Intramolecular interactions between vaccine construct and TLR2 was visualized by LigPlot + . Immune simulation study explains the efficiency of vaccine construct to elicit the immune response. Induction of memory B-cells and T-cells is one of the criteria to be a successful vaccine candidate [34][35][36] . Each injection of VC1 and VC2 increased the level of memory B-cells and T-cells and this population level is sustaining after the third injection (Figs. 7, 8). Also, the elevated level of IFN-γ and IL2 is exhibiting the capacity of vaccine constructs to establish the antiviral state, in general IFN-γ involves in antiviral replication and act as main effector molecule in cell mediated immunity. Escherichia coli is most preferable choice of the host to obtain more quantity of recombinant vaccine protein 37 . Java Codon Adaptation tool displayed codon adaptability index

Conclusion
The present study has made an attempt to design the multi-epitope vaccine against KFDV by using immuneinformatics tools. To our knowledge this is first report for identifying T-cell and B-cell targeting epitopes of KFDV. The designed chimeric vaccine peptide could elicit immune response, but still it need to be tested on in-vitro and in-vivo models. Interestingly, our constructed vaccine has cross-protection effect against AHFV. This study has given foresight for development of new prophylaxis for KFDV control in India, and gives the directions in selecting epitopes for KFDV vaccine development.

Material and methods
Retrieval of amino acid sequence and conservation analysis. Totally  Multiple sequence alignment (MSA) of retrieved sequences were performed in CLCworkbench (https:// www. qiage nbioi nform atics. com/ produ cts/ clc-main-workb ench/), alignment file (.aln) was generated from all retrieved envelop protein sequences, neighbor-joining method was used to the phylogenetic tree and conservation score of each amino acids were obtained.

Prediction of T-cell and B-cell epitopes.
Eight different immunoinformatic tools were used to predict T-cell (MHC class I and MHC class II) and B-cell epitope regions on KFDV envelope protein. Three different tools (IEDB-MHC-I, IEDB-combined and MAPPP) were used to obtain MHC class I epitope, proteasomal cleavage score, Transporter of Antigenic Peptide (TAP) and MHC scores. MHC class II epitope was prediction by IEDB-MHC-II and MHC2Pred tools. B-cell epitopes were identified by ABCpred, BCpred, and SVMTriP tools, the detailed information on tools algorithm, URL site and threshold values were in Table 9. www.nature.com/scientificreports/ Self-antigen and cross-protection analysis. Uniport-BLAST tool (https:// www. unipr ot. org/ blast/) was used to checked all top ranked T-cell and B-cell epitope sequences for self-antigen and cross-protection analysis. Self-antigen or similarity of human proteins with predicted epitopes were identified by Uniport-BLAST search against uniportkb_human databases. And UniProtKB database BLAST search was used to identify the cross protection of selected epitopes with other pathogenic organisms.
Construction of multi-epitope vaccine candidate. Based on high score, multi-immunogenic and overlapping sequence of T-cell and B-cell epitopes were conjugated to construct the multi-epitope vaccine candidates. Construct of vaccine has started with adjuvant β defensin (ACK99045.1) peptide sequence, in order to augment the immunogenicity of the vaccine candidate. Adjuvant was connected to MHC-I epitopes with EAAAK linker, every individual epitopes of MHC-I were linked by GGGS and all MHC-II epitopes assembled with GPGPG linker. Whereas B-cell epitopes were combined by KK linker. Usage of the linker between two epitopes will increase the immunogenicity and function of the vaccine constructs 54 .
Allergenicity and antigenicity of the vaccine candidate. Algpred (https:// webs. iiitd. edu. in/ ragha va/ algpr ed/ submi ssion. html) used to predict the allergenicity of vaccine candidate based on the similarity of known epitope. The vaccine candidate sequence was uploaded into the server, the IgE epitope and PID-BLAST search on allergen representative peptide algorithms were chosen for allergen prediction 55   www.nature.com/scientificreports/ ClusPro 2.0 server (https:// clusp ro. org) 66 and PatchDock (https:// bioin fo3d. cs. tau. ac. il/ Patch Dock/ php. php) 67 were used for this purpose. ClusPro 2.0 server was used to predict the binding energy between TLR-2 and vaccine construct, experiment was performed by uploading PDB files of receptors and ligands into server and submitted with default parameters. PatchDock was also performed at default settings, top 1000 transformed docking position generated from PatchDock server was re-scored by FireDock server (http:// bioin fo3d. cs. tau. ac. il/ FireD ock/ php. php) 68,69 and global energy of each vaccine candidate with TLR2 was obtained.
Ligplot analysis and molecular dynamic simulation of Docking complex. Interaction of TLR2 with vaccine candidates were assessed by Ligplot, docking complexes resulted from cluspro analysis were saved as PDB format and uploaded into Ligplot tool 70 and intermolecular interactions such as hydrogen bonds, hydrophobic contact was obtained in the form of 2D representation.
The molecular dynamics simulation was performed to check the stability of docking complexes. The iMODS (http:// imods. chaco nlab. org/) 71-73 web-server was used to calculate B-factor (disorder of atoms), structural deformability and eigenvalue.
In silico cloning optimization of KFD multi-epitope vaccine candidates. The Java Codon Adaptation Tool (http:// www. prodo ric. de/ JCat) was used for reverse translation and codon optimization of designed vaccine candidates 74 . Codon optimization was executed in order to express the KFD multi epitope vaccine construct in E. coli (strain K12). The output files were checked for Codon Adaptation Index (CAI) (> 0.8-1.0), and GC content (30-70%), to ensure transcription and translational efficiency of designed vaccine construct. Resulted optimized codon sequence of vaccine candidates were introduced with XhoI and NdeI restriction sties at N-terminal and C-terminal ends respectively. SnapGene (www. snapg ene. com) was utilized to insert the vaccine sequence in expression vector pET-30a(+), between XhoI and NdeI cloning site, final clones of vaccine candidates VC1 and VC2 were obtained.
Immune response simulation. Immune response profile of the constructed vaccine candidates was evaluated by C-IMMSIM server 75 . Experiment was performed by uploading PDB files of vaccine constructs, random speed, simulation volume and simulation steps were set as 12,345, 50 and 1100 respectively. Three insilco injections were given at the time steps of 1, 84 and 168 respectively (1-time step is equal to 8 h in real life), with no LPS and maintained minimum 30 days of time interval between two injections.