An immunoinformatics approach for the design of a multi-epitope subunit vaccine for urogenital schistosomiasis

Retrieval of Schistosoma haematobium HSP sequence for vaccine construction

The amino acid sequences of S. haematobium HSPs were retrieved from the NCBI protein database (www.ncbi.nlm.nih.gov) and subjected to multiepitope vaccine designing. To test for the proteins ability to induce immune response within a host body, all the retrieved protein sequences were subjected to an antigenicity test using the ANIGENpro database (scratch.proteomics.ics.uci.edu). Antigenic probability of ≥0.8 was used to determine which proteins to be chosen for the next step of the multi epitope vaccine construct (Pandey, Bhatt & Prajapati, 2018).

Prediction of cytotoxic T-lymphocytes and helper T-lymphocytes epitope and immunogenicity

To get an immunogenic cytotoxic T-lymphocytes (CTL) epitope with potential of inducing cell mediated immunity and form memory cells, all the three highly antigenic HSP sequences in FASTA format were fed into the NetCTL 1.2 server (http://www.cbs.dtu.dk/services/NetCTL/) to predict the CTL epitope at threshold score of 0.75 (using the default setting). The server predicts CTL epitopes from inputted protein sequences based on the training dataset (prediction of MHC-I binding peptides, proteasomal C-terminal cleavage, and Transporter Associated with Antigen Processing (transport efficiency) and only epitopes with combined score greater than 0.75 were selected as CTL epitopes (Khatoon, Pandey & Prajapati, 2017; Pandey, Bhatt & Prajapati, 2018). The selected CTL epitopes (Scores >0.75) were further analysed in Immune Epitope Design Database (www.iedb.org) for MHC class I immunogenicity prediction. The immunogenicity score determines the probability of eliciting an immune response (higher score means greater probability of immune response and vice versa) (Pandey, Bhatt & Prajapati, 2018).

The helper T-Lymphocytes (HTL) epitopes of 15mer length for mouse allele (H2-IAb, H2-IAd and H2-IEd) were predicted for the HSPs of S. haematobium using the Immune Epitope Design Database. The output epitopes were ranked based on lower percentile rank scores that is, the lower the ranked score the higher the binding affinity for HTL receptors (Khatoon, Pandey & Prajapati, 2017). Only epitopes with higher binding affinity (≤1.5) with MHC-II were selected for the final multi-epitope vaccine construct. In order to determine the ability of the predicted epitopes in activating Th1 type immune response followed by IFN-γ production, the top 17 HTL epitopes were subjected in FASTA format to the IFN epitope server by using the Motif and SVM hybrid as approach selection and IFN-γ versus other cytokine as model of prediction.

Construction of multi-epitope subunit vaccine

To design a viable vaccine construct capable of inducing innate and adaptive immune response, the high scoring CTL and HTL epitopes from the present immunoinformatics predictions were used. These epitopes were joined together using AAY and GPGPG linkers which were added at the intra-epitope position thereby linking the CTL and HTL epitopes respectively (Khatoon, Pandey & Prajapati, 2017). Also, a TLR-4 agonist (RS-09; Sequence: APPHALS) was used adjuvant (Pandey, Bhatt & Prajapati, 2018; Shanmugam et al., 2012) and linked with the multi-epitopes (CTL and HTL) using EAAAK linker (Lee et al., 2014).

B cell epitopes prediction for S. haematobium proteins

The BCPREDS server (ailab.ist.psu.edu/bcpred) was used to predict the linear B-cell epitopes for the final vaccine construct. The 493 amino acid sequence of the final vaccine construct was inputted into the BCPREDS server in plain format followed by the selection of a fixed length epitope prediction method and length of the epitope. The default method of BCPREDS was selected as the prediction method for the amino acid epitopes (El-Manzalawy, Dobbs & Honavar, 2008; Pandey, Bhatt & Prajapati, 2018). The specificity threshold was set as default at 75% and conformational epitopes were predicted using the ElliPro server for the tertiary protein structure of the vaccine construct (Pandey, Bhatt & Prajapati, 2018).

Prediction of antigenicity, allergenicity and physiochemical parameter of vaccine construct

The potential of the vaccine construct in eliciting immunological response through binding to the B and T cell receptor was determined by first determining the antigenicity of the vaccine construct. The ANTIGENpro server (scratch.proteomics.ics.uci.edu); a sequence-based alignment free and pathogen independent predictor was used to generate the antigenicity index (Khatoon, Pandey & Prajapati, 2017).

The non-allergenic potential of the multi-epitope vaccine construct was determined using the AllerTOP v.2.0 server (www.ddg-pharmfac.net/AllerTOP). Several physiochemical parameters of the vaccine construct which include amino acid composition, the theoretical pI, instability index, in-vitro and in-vivo half-life, aliphatic index, molecular weight and grand average of hydropathicity (GRAVY) parameters were all assessed using the ProtParam server (https://web.expasy.org/protparam/) (Onile, 2014).

Prediction, refinement and validation of the tertiary structure of vaccine construct

The 3D structure of the vaccine construct was predicted using RaptorX structure prediction server (raptorx.uchicago.edu/). This server predicts secondary and tertiary structures as well as contact map, solvent accessibility, disordered regions and binding sites of protein sequences. RaptorX server also assign some confidence scores P-value, GDT (global distance test) and uGDT (un-normalized GDT), and modeling error at each residue to indicate the quality of a predicted 3D mode. In order to improve and refine the predicted 3D model of vaccine construct, the output model of RaptorX server was further modelled using GalaxyRefine server (http://galaxy.seoklab.org/) (Heo, Park & Seok, 2013). This server uses CASP10 based refinement approach to reconstruct and repack the protein side chain, followed by dynamics simulation to relax the structure. ProSA-wed (https://prosa.services.came.sbg.ac.at/prosa.php) was employed in the tertiary structure validation; also the refinement output was validated using Ramachandran plot analysis (Mordred.bioc.cam.au.ck/~rapper/rampage.php).

Disulfide engineering for vaccine stability

The Disulfide by Design v2.12 server was used for the in silico disulfide engineering of the refined model 3D structure of the vaccine construct. The refined model was first run to determine residue pairs useful for the purpose of disulphide engineering. Only residue pairs with energy (kcal/mol) less than 2.2 and Chi³ (χ³) between −87 and +97 degree were selected for mutation (Rana & Akhter, 2016). A total of 5 residue pairs were eventually selected for mutation by using the create mutate icon of the Disulfide by Design server.

Codon adaptation and in silico cloning

To determine and optimise the expression rate of the vaccine construct in a proper expression vector, the primary sequence of the vaccine construct was inputted into the Java Codon Adaptation Tool (JCAT) for codon optimization which was performed in the host Escherichia coli strain K12. Due to difference in codon usage of E. coli and native host S. haematobium from where the vaccine construct sequences originated, three different options in the JCAT server were selected which include: avoid rho-independent transcription, prokaryotes ribosome binding site and restriction enzymes cleavage sites. To determine the high-level of protein expression, the output of the JCAT which include codon adaptation Index and percentage GC content were used. E. coli pET-28(+) was used as vector to clone the adapted nucleotide sequence (provided by JCAT server) of the final vaccine construct using the restriction cloning module of SnapGene tool (Ali et al., 2017).

Molecular docking of vaccine construct with TLR4

The ClusPro 2.0 server was used for the molecular protein–protein docking to check for the binding affinity between the vaccine construct and TLR-4 receptor (Kozakov et al., 2017; Pandey, Bhatt & Prajapati, 2018). The refined sub unit multi-epitope vaccine protein and TLR-4 PDB file (4G8A) obtained from the RCSB-Protein Data Bank were used as the ligand and receptor respectively.

Retrieved S. haematobium HSP sequences for vaccine construct

To design the multi-epitope subunit vaccine for urinary schistosomiasis consisted of seven (7) HSP sequences (HSP90 (KGB37337.1), HSP11 (XP_012799478.1), HSP97 (KGB34481.1), HSP75 (KGB38609.1), HSP83 (KGB37622.1), HSP10 (KGB41273.1), HSP60 (KGB40455.1)) were retrieved from the NCBI database (https://www.ncbi.nlm.nih.gov) in FASTA format. Only three (HSP11, HSP97, HSP90) proteins had an antigenic probability of ≥0.8 after retrieved protein sequences were analyzed using ANTIGENpro for antigenicity prediction (Table 1), consewuently they were selected for the vaccine construct.

Predicted CTL, HTL epitope and immunogenicity

By using the NetCTL 1.2 server for prediction of the CTL receptor specific immunogenic epitopes, a total of 23 CTL epitopes (9mer length) were obtained for the inputted three HSP protein sequences and only 12 epitopes with high immunogenicity scores as obtained from IEDB server were chosen as the final CTL epitopes to undergo vaccine designing (Table 2).

The prediction of the HTL epitope was carried out using the IEDB MHC-II epitope prediction module, all three protein sequences were subjected to the module. The mouse alleles used for prediction were H2-1Ad, H2-1Ed and H2-1Ab. However, only 17 HTL epitopes with lowest percentile with ≤1.5 with MHC-II were selected and subjected for the final multi-epitope vaccine construct (Table 3A). Also, all the 17 HTL epitopes have the capacity of inducing IFN-γ, and this was evident from the positive score obtained from the IFN epitope server output (Table 3B).

Constructed multi-epitope subunit vaccine

The constructed multi-epitope subunit vaccine consisted of 12 CTL epitopes and 17 HTL epitopes with high binding affinity fused together with the help of AAY, GPGPG and EAAAK linkers. The adjuvant and the CTL epitopes were combined by the EAAAK linkers, the intra-CTL and intra-HTL epitopes were joined by intra-CTL AAY and intra-HTL GPGPG linker respectively to form the final vaccine constructs composed of 493 amino acid residues. For maximum ability of the candidate vaccine to elicit adequate immune response, adjuvant and CTL epitopes were joined together using an EAAAK linker (Fig. 1A).

Predicted B cell epitopes for S. haematobium proteins

A total of 14 epitopes having 0.99 and above probability score with 20mers length were selected as the linear B cell binding epitopes for the final vaccine construct (Table 4). Discontinuous epitopes of 125 amino acids were predicted from the final 3D model of vaccine construct with the score of 0.734 (Fig. 1B).

Predicted antigenicity, allergenicity and physiochemical parameter of the vaccine construct

The antigenicity test of the designed vaccine construct showed the antigenicity probability to be 0.905150, this represents the antigenic nature of the vaccine construct. An allergenicity test carried out for the vaccine construct also showed that the vaccine candidate is non-allergenic in nature.

The prediction of the physiochemical parameters of the vaccine construct was characterized by using the ProtParam server which evaluated for the following parameters; The molecular weight was 51409.68 kDa, theoretical pI 9.11, while the total number of negative and positive charge residues were 36 and 49 respectively. The estimated half-life was 4.4 h in mammalian reticulocytes in-vitro, 20 and 10 h in yeast and E. coli respectively. The extinction coefficient was found to be $36790{\mathrm{M}}^{-1}\mathrm{c}{\mathrm{m}}^{-1}$, at 280 nm measured in water assuming all pairs of cysteine residues are reduced, instability index was given to be 35.01. The value of the aliphatic index was 77.81 while the grand average of hydropathicity (GRAVY) was −0.138.

Predicted, refined and validated tertiary structure of vaccine construct

The 3D model and tertiary structure were predicted using RaptorX server (Fig. 2). The best template for the homology modeling was PDB ID: 5ja1A with 493 (100%) amino acid residues modelled as single domain with 4% disorder. The secondary structure information includes helix 31%, Beta sheet 17% and coiled structure 51%. To refine the vaccine construct, GalaxyRefine server was used. Out of all refined models, model 3 was selected as the best based on various parameters including GDT-HA (0.9239), RSMD (0.484), MolProbity (2.69), Clash Score (32.7) and Poor rotamers (1.9).

The refinement output was also validated using Ramachandran plot analysis, which revealed that 93.3% of the residue are located in most favoured regions, 5.3% in allowed regions and only 1.4% in outlier region and the quality and potential errors in the crude 3D as verified by ProsSA-web showed Z-scores of −6.01 for ProSA (Figs. 3 and 4).

Disulfide engineering for vaccine stability

Disulfide engineering was done to stabilize the final model of the vaccine construct using Disulfide by Design v2.12. The result showed a total 58 pairs of residues that are useful for the purpose of disulphide engineering (Fig. 5) but 10 mutations were created from the only five pairs of residues (GLU56-CYS249, TYR60-HIS493, VAL116-GLN242, GLN243-GLU489) that were found with energy value less than 2.2 and Chi3 value range between −87 and +97 degree.

Codon adaptation and in silico cloning

The codon representing the subunit vaccine candidate was adapted to the codons of the E. coli K 12 strain using JCAT server and showed that the optimized codon sequence has a length of 1,450 nucleotides with a codon adaptation index (CAI) of 0.95 and average GC content of 54.56%. The adapted codon sequence was later inserted into pET28(+) vector to form the restriction clone (Fig. 6). The adapted codon sequences used in insilico cloning are made up of the following 1,400 nitrogenous bases “GCTCCGCCGCACGCTCTGTCTGAAGCTGCTGCTAAACTGGACATGACCGAAGAATGGCTGTACGCTGCTTACCAGACCGAAGAAATCGACGGTACCCTGGCTGCTTACTCTATCGCTGCTGGTGAACCGACCTACGCTGCTTACTCTCAGGCTCAGCTGATCGAAGAATACGCTGCTTACAACTCTAAAAACGCTGTTGAAGAATACGCTGCTTACTTCACCGTTATCGAACAGTGCCTGTACGCTGCTTACGGTTCTGACGGTTGCTTCGTTTCTTTCGCTGCTTACTTCTCTGGTAACATCACCGGTATCTTCGCTGCTTACAAACTGGACGGTGAAGCTAACGTTTACGCTGCTTACGGTTCTAAAATCGAAAACGACCTGTACGCTGCTTACGAAGTTATCTCTCTGATCGACGAATACGCTGCTTACACCTCTTCTACCGACGGTGACCTGGTTGGTCCGGGTCCGGGTCGTAAAGTTTCTCAGATCAAAATCGTTACCTCTCGTGTTAAAGCTGGTCCGGGTCCGGGTCCGCGTAAAGTTTCTCAGATCAAAATCGTTACCTCTCGTGTTAAAGGTCCGGGTCCGGGTCAGCCGCGTAAAGTTTCTCAGATCAAAATCGTTACCTCTCGTGTTGGTCCGGGTCCGGGTCGTTCTCTGGTTCGTATCCGTAAAGTTCTGGACTCTATCGCTGCTGGTCCGGGTCCGGGTGTTTCTCAGCAGCAGGCTATGGAATCTGTTTGCCACCCGATCATCGGTCCGGGTCCGGGTGACATCGTTTCTCAGCAGCAGGCTATGGAATCTGTTTGCCACCCGGGTCCGGGTCCGGGTATCGTTTCTCAGCAGCAGGCTATGGAATCTGTTTGCCACCCGATCGGTCCGGGTCCGGGTAACCGTTCTCTGGTTCGTATCCGTAAAGTTCTGGACTCTATCGCTGGTCCGGGTCCGGGTTGCGCTTTCCAGGCTGCTATCTGCTCTCCGGCTTTCAAAGTTAAAGGTCCGGGTCCGGGTGACCAGCCGCGTAAAGTTTCTCAGATCAAAATCGTTACCTCTCGTGGTCCGGGTCCGGGTTTCAACCGTTCTCTGGTTCGTATCCGTAAAGTTCTGGACTCTATCGGTCCGGGTCCGGGTGGTTGCGCTTTCCAGGCTGCTATCTGCTCTCCGGCTTTCAAAGTTGGTCCGGGTCCGGGTAACTTCAACCGTTCTCTGGTTCGTATCCGTAAAGTTCTGGACTCTGGTCCGGGTCCGGGTCGTGGTTGCGCTTTCCAGGCTGCTATCTGCTCTCCGGCTTTCAAAGGTCCGGGTCCGGGTTACCGTGTTCTGACCACCAAAGAACTGGTTAAAGCTTTCACCCGTGGTCCGGGTCCGGGTACCCGTCTGTTCCTGGCTCAGAAAAA”.

Molecular docking of vaccine construct with TLR4

A total of 29 models were generated from the molecular protein–protein docking between the vaccine construct and the TLR4 receptor using ClusPro 2.0 server. Only one model (000.23) with the lowest energy score (−1,250.0 and highest binding affinity among other predicted docked complex fulfilled the desired criteria for best-docked complex and was eventually selected (Fig. 7).