Anti-COVID-19 multi-epitope vaccine designs employing global viral genome sequences

Retrieval of sequences and multiple alignment

GISAID (available at https://www.epicov.org/epi3/cfrontend) and NGDC databases were mined to retrieve updated complete genomes of SARS-CoV-2. The retrieved genomes were then aligned using the Multiple Alignment using Fast Fourier Transform tool (MAFFT) (available at: https://www.ebi.ac.uk/Tools/msa/mafft/) (Katoh et al., 2002). Default parameters were used encompassing 200 PAM/κ = 2 scoring matrix and a gap penalty of 1.53. These aligned sequences were visualized and analyzed in Unipro UGENE v.34 software (Okonechnikov et al., 2012) to generate a consensus sequence of all genomes and later, nBLAST was performed to analyze its diversity with Reference sequence NC_045512.2. The ORFs from this sequence were retrieved and subjected to further analysis of predicted proteins.

CTL epitope mapping

CTLpred (available at: http://crdd.osdd.net/raghava/ctlpred/) was used for prediction of CD8 T cell epitopes (Bhasin & Raghava, 2004). Consensus approach was used that predicted CTL epitopes based on both Support Vector Machine (SVM) and Artificial Neural Network (ANN) methods using default thresholds (0.36 for SVM method and 0.51 for ANN method). For further analyses top 100 peptides were selected. These predicted epitopes were also verified by HLApred (available at: http://crdd.osdd.net/raghava/hlapred/) and ProPred1 (available at: http://crdd.osdd.net/raghava/propred1/) (Singh & Raghava, 2003).

Antigenicity

All CTL epitopes were screened for antigenicity using VaxiJen (available at: http://www.ddg-pharmfac.net/vaxijen/) (Doytchinova & Flower, 2007). VaxiJen is an online server that predicts the antigenicity of peptides using the physicochemical properties of the proteins irrespective of their lengths and the need of alignments. We used VaxiJen server for our CTL epitopes using the default threshold 0.5 for viruses.

Immunogenicity

After screening for the antigenic potential of CTL epitopes, epitopes having immunogenic potential were screened using the IEDB analysis resource for MHC class I Immunogenicity (available at: http://tools.iedb.org/immunogenicity/) (Calis et al., 2013). This tool works best for 9-mer epitopes but larger lengths can be employed too. Only those epitopes that showed a strong binding affinity with HLA super family alleles were shortlisted. Default threshold was used for prediction of results. Epitopes having percentage rank less than 0.5 were qualified as strong binders.

Strong binding affinity with HLA alleles and population coverage analysis

CTL Epitopes that were both antigenic and immunogenic were screened for their binding affinity with HLA superfamily alleles using Net MHC 4.0 (available at: http://www.cbs.dtu.dk/services/NetMHC/) (Nielsen et al., 2003; Andreatta & Nielsen, 2016). The tool predicts results on the basis of an ANN that has been trained using 81 different HLA alleles. MHC binding affinity of shortlisted epitopes was checked against HLA superfamily alleles (i.e., HLA-A0101, HLA-A0201, HLA-A0301, HLA-A2402, HLA-A2601, HLA-B0702, HLA-B0801, HLA-B2705, HLA-B3901, HLA-B4001, HLA-B5801, HLA-B1501). HLA supertype representative are clusters of HLA alleles based on their overlapping peptide binding specificity, hence can cover 80% of population in a less computationally exhaustive environment. The concept has been effectively utilized in prediction of promiscuously recognized epitopes against SARS, HCV, vaccinia virus and HBV (Frahm et al., 2007; Sidney et al., 2008). Furthermore, population coverage analysis was performed using IEDB resource server. The tool predicts population coverage of epitopes based upon the restriction database (Bui et al., 2006).

HTL epitope mapping and population coverage analysis

The ORFs that didn’t yielded prioritized CTL epitopes were subjected to Helper T Lymphocytes (HTL) epitope (CD4 T cell epitopes) prediction by the HLApred server. This was done to cover the whole proteome of the virus. All alleles of MHC-II were selected, and experimental as well as predicted binders were included in the analysis. A default threshold was used, and all predicted epitopes that shared no identity with humans were selected. The predicted epitopes were screened for their antigenic potential using VaxiJen as discussed earlier. The antigenic epitopes were further screened for their IFN-gamma induction potential using IFNepitope using hybrid approach that involves both SVM and Motif based approach. The prediction model chosen for analysis was IFN-gamma vs non IFN-gamma (available at: http://crdd.osdd.net/raghava/ifnepitope/scan.php) (Dhanda, Vir & Raghava, 2013). The epitopes screened positive for IFN-gamma were checked for overlapped B cell epitopes predicted by ABCpred (available at: http://crdd.osdd.net/raghava/abcpred/) using a default threshold (Saha & Raghava, 2006b). Lastly, the shortlisted epitopes were checked for their strong binding affinity with all the HLA DR alleles present in Net MHCII server (available at: http://www.cbs.dtu.dk/services/NetMHCII/). These alleles were DRB1_0101, DRB1_0103, DRB1_0301, DRB1_0401, DRB1_0402, DRB1_0403, DRB1_0404, DRB1_0405, DRB1_0701, DRB1_0801, DRB1_0802, DRB1_0901, DRB1_1001, DRB1_1101, DRB1_1201, DRB1_1301, DRB1_1302, DRB1_1501, DRB1_1602, DRB3_0101, DRB3_0202, DRB3_0301, DRB4_0101, DRB4_0103, DRB5_0101. HLA-DR alleles are present on the surface of antigen-presenting cells and present the epitopes to helper T-cells initiating an immune response (Smith-Garvin, Koretzky & Jordan, 2009). HLA DRB alleles were reported to be involved in viral clearance, hence all the available HLA DRB alleles in server were selected for analysis (Barrett, Ryan & Crowe, 1999; Keicho et al., 2009). Lastly, Population coverage analysis was performed using IEDB population coverage analysis tool.

Multi-epitope vaccine design and construction

Epitopes were arranged as per their arrangement in the reference sequence, the CTL epitopes were linked together using a flexible linker GGGGS as recommended by Chen, Zaro & Shen (2013), and HTL epitopes with GPGPG linker to construct a multi-epitope (Saadi, Karkhah & Nouri, 2017). Moreover, a second multi-epitope construct was also designed where β-defensin, an adjuvant, was added at the N terminus of multi-epitope with EAAAK linker (Shey et al., 2019).

Sequence-based physio-chemical properties of the multi-epitope design

Different physicochemical properties of multi-epitopes were predicted using ProtParam tool including prediction of in vivo half-life, theoretical isoelectric point (pI), molecular weight, instability index, aliphatic index and Grand Average of Hydropathy (GRAVY) index of peptides (Wilkins et al., 1999; Gasteiger et al., 2005). The degradation of peptides was based on amino acid at N-terminus hence N-end rule in Protparam was employed to predict in vivo half-life of both peptides. Instability index with a threshold of 40 was used in order to determine the stability of the protein in a cellular environment. The aliphatic and GRAVY index of both peptides was calculated based on their amino acid profiles.

Allergenicity of multi-epitope

Both multi-epitope constructs (with and without adjuvant) were screened for their allergenicity using two online servers; AlgPred (available at: http://crdd.osdd.net/raghava/algpred/) (Saha & Raghava, 2006a) and Allergen FP (available at: http://ddg-pharmfac.net/AllergenFP/) (Dimitrov et al., 2014). In AlgPred, all four approaches were used including mapping of IgE epitopes and PID approach (allergenicity prediction is based on similarity with experimental IgE epitopes), BLAST search on allergen representative peptides (makes prediction by performing BLAST search against 2890 allergen representative peptides), MEME/MAST motif approach (predicts allergenic motifs using MAST) and hybrid approach. Allergen FP predicts the allergenicity of peptide using an alignment-free method by implementing a four-step algorithm and compared by Tanimoto coefficient.

Antigenicity of multi-epitope design and modeling

The antigenicity profiles of both multi-epitope constructs were predicted again to confirm if the epitopes were antigenic or not. We used two servers this time; VaxiJen (available at: http://www.ddg-pharmfac.net/vaxijen/) (Doytchinova & Flower, 2007) and AntigenPro (Magnan et al., 2010). Scratch Protein Predictor server available at http://scratch.proteomics.ics.uci.edu/ hosted the tool Antigen Pro for antigenicity estimation. AntigenPro, like VaxiJen, predicts antigenicity using physicochemical properties, irrespective of the length of peptides as well as its alignment.

Secondary structure prediction

For evaluation of further structural characteristics of designed multi-epitopes, PDBsum (available at http://www.ebi.ac.uk/thornton-srv/databases/pdbsum/) (Laskowski et al., 2018) was used. PDBsum provides insight into the unique structural aspects of proteins, peptides and their ligands.

Tertiary structure prediction and refinement

Tertiary structures of both vaccines were predicted by Scratch Protein Predictor-hosted 3Dpro tool (available at http://scratch.proteomics.ics.uci.edu/) (Cheng et al., 2005) and refined by Galaxy refine (available at http://galaxy.seoklab.org/cgi-bin/submit.cgi?type=REFINE) (Ko et al., 2012). The model showing highest Ramachandran favored residues and minimum poor rotamers were selected. The refined structures were checked for their stability and flexibility using molecular dynamics (MD) simulation studies. Structural flexibility of a protein/peptide is important for its molecular recognition and its function, so a coarse-grained protein model implemented in a webserver CABS-Flex 2.0 (available at: http://biocomp.chem.uw.edu.pl/CABSflex/) (Kuriata et al., 2018) was used for near-native dynamics of both vaccines. Default distance restrained parameters were used, number of cycles and cycles between trajectory frames were raised to 100. Temperature of simulation was also kept default.

Molecular docking with Toll-Like Receptor 8

The models were then checked for their interaction with Toll-Like Receptor 8 (TLR8). Active and passive residues for these analyses were predicted using CPORT (De Vries & Bonvin, 2011). Molecular docking was performed by using guru level interface of HADDOCK 2.2 (available at: http://haddock.science.uu.nl/services/HADDOCK/haddockserver-guru.html) (Van Zundert et al., 2016) with default parameters and a representative structure from the top-ranked docked cluster with minimum HADDOCK score was refined using Refinement Interface (available at: http://haddock.science.uu.nl/services/HADDOCK/haddockserver-refinement.html) (Van Zundert et al., 2016). Moreover, NMA analysis of refined complexes was also performed using iMOD in order to determine their deformation potential (Kremer, Mastronarde & McIntosh, 1996) (available at: https://bio3d.colorado.edu/imod/paper/). The interacting residues between each complex were predicted using PDBsum analysis (Laskowski et al., 2018).

Codon optimization and in silico cloning in Escherichia coli

Codon optimization of both multi-epitopes was performed using JCAT (available at: http://www.jcat.de/) (Grote et al., 2005) and codons were adapted as per E. coli K12 strain. Default parameters were used and plasmid UC19 was used for in silico cloning of both multi-epitopes. Moreover, a His6 tag was added at the both ends of sequences for the purification of multi-epitope proteins.

Retrieval of sequences and multiple sequence alignments

The recruited 475 complete SARS-CoV-2 genomes showed a high level of conservancy upon multiple sequence alignments and in the phylogenetic tree (Fig. 1). The consensus sequence was found to be 99% identical, with 10 gaps and no mismatches with Ref seq NC_045512.2 from Wuhan. There were 11 ORFs available in Refseq NC_045512.2 for protein coding sequence (ORF1ab (YP_009724389.1), Spike Glycoprotein (YP_009724390.1), ORF3a (YP_009724391.1), Envelope Protein (YP_009724392.1), Membrane Glycoprotein (YP_009724393.1), ORF6 (YP_009724393.1), ORF7a (YP_009724393.1), ORF7b (YP_009725318.1), ORF8 (YP_009724396.1), ORF9 (YP_009724397.2) and ORF10 (YP_009725255.1)). These ORFs were subjected to epitope mapping.

CTL epitope mapping

A total of 13 CTL epitopes were shortlisted for designing multi-epitopes. These shortlisted epitopes were predicted to be antigenic (as determined by VaxiJen tool), immunogenic (as determined by MHC-I IEDB Immunogenicity tool), and had potential to bind with HLA superfamily alleles (as predicted by NetMHC) (Table 1). All the screened epitopes had 88.42% population coverage. The details of epitopes predicted from each ORFs filtered at each stage along with their properties have been provided in Table S1.

HTL epitope mapping

A total of 10 HTL epitopes were scrutinized for multi-epitopes. These shortlisted epitopes were predicted to lack homology with Human proteome (as determined by HLApred), had overlapping B cell epitopes (as determined by ABCpred), showed potential to induce IFN-gamma response (as determined by IFNepitope server) and demonstrated strong binding affinity with HLA DR alleles (as predicted by NetMHCII) (Table 2). All the epitopes had 89.70% population coverage (as predicted by IEDB population coverage analysis tool). The detail of step wise scrutiny and all the initially predicted epitopes is provided in Table S2.

Multi-epitope structure design and modeling

Two multi-epitope constructs were designed with and without β-defensin adjuvant to analyze and compare whether the designed multi-epitope is efficacious only in presence of adjuvant or has the ability to elicit immune response solely. Cov-I-Vac was the multi-epitope construct without adjuvant and Cov-II-Vac was with adjuvant (Fig. 2).

Sequence-based physio-chemical properties of multi-epitope

Physio-chemical analysis of Cov-I-Vac and Cov-II-Vac revealed that both constructs are reasonably stable (38.48 and 39.91 respectively), have molecular weight <50 kDa, lower GRAVY indexes (i.e., 0.394 and 0.249, respectively), with a pI of 9.97 and 10.03, respectively. The N terminus of Cov-I-Vac has Threonine so its half-life in mammalian reticulocytes was predicted to be 7.5 h (if tested in vitro) using Protparam. Meanwhile, Cov-II-Vac, because of Glycine at N terminus, has a half-life of >30 h in mammalian reticulocytes, as predicted using N-end rule. However, no difference in in vivo half-life was predicted, and was found to be >20 h in yeast and >10 h in E. coli.

Antigenicity, allergenicity and immunogenicity of Cov-I-Vac and Cov-II-Vac

The antigenicity potential was found to be retained in both multi-epitopes and was estimated as 0.6153 and 0.5947 for Cov-I-Vac and Cov-II-Vac, respectively. Both constructs were reported as antigenic by AntigenPro as well with an antigenicity score of 0.516406 and 0.386605.

Moreover, both vaccines were soluble upon overexpression in E. coli. The safety profile also showed that constructs do not exhibit any experimentally reported allergens and thus projected to be non-allergens.

Secondary structure of Cov-I-Vac and Cov-II-Vac

The secondary structure prediction of Cov-I-Vac shows that it has a sheet and beta hairpin between two strands (Ala160-Il161, Gly164-Gly165), six helices, 42 beta turns and two gamma turns (Fig. 3A). The secondary structure of Cov-II-Vac shows that it has a sheet and two beta hairpins among three strands (Val193-Thr198, Gly201-Val208, Gly216-Gly217), 10 helices, five helix-helix interactions, and 37 beta turns (Fig. 3B).

Tertiary structure prediction and refinement

The tertiary structure of Cov-I-Vac had 10.2% poor rotamers with 90.5% Rama favored residues; however, the refined structure selected for further analysis had an RMSD of 0.541 with 0% poor rotamers and 93% Rama favored residues. Similarly, the initial structure of Cov-II-Vac had 5.9% poor rotamers with 89% Rama favored residues, but the refined shortlisted structure with RMSD 0.548 had 92.9% Rama favored residues and 0% poor rotamers.

CABS-flex analysis

The trajectory of 10 models was generated after simulation (Fig. 4). The Root Mean Square Fluctuation (RMSF) of this simulation was found to be between 3.9 and 5.8 Å.

The trajectory of Cov-II-Vac after fast simulation

Trajectory analysis of fast simulation of Cov-II-Vac was analyzed and the RMSF of the structure was found to be 0.75–7.25 Å. The maximum fluctuation was observed at residue 227 while minimum fluctuation was residue 39–59 (Fig. 5).

The active and passive residues of Cov-I-Vac, Cov-II-Vac and TLR8 are provided in Table S3. The molecular docking of Cov-I-Vac with TLR8 resulted in the top-ranked cluster with a HADDOCK score of −40.8 +/− 5.5. Low values of HADDOCK score, especially those that are negative, usually imply significantly high interaction between proteins. HADDOCK refined the representative structure of this complex and clustered the resulting 20 structures into one complex that represented 100% of the refined water models generated. After undergoing molecular refinements, the HADDOCK score of the Cov-I-Vac-TLR8 complex improved significantly to reach −251.3 +/− 2.3. The details of the refined interaction along with statistical parameters are provided in Table 3. Cov-II-Vac also showed significant interaction with TLR8 that was −83.5 +/− 6.3 which was a higher interaction score than Cov-I-Vac. However, after refinements, the interaction of Cov-II-Vac with TLR8 was lower compared to Cov-I-Vac-TLR8 interaction, that is, −217.4 +/− 3.5. The details of HADDOCK score with all its parameters are mentioned in Table 3.

For an in depth analysis of protein–protein interactions, PDBSum analysis revealed that a total of 36 residues of Cov-I-Vac interact with 43 residues of TLR8, with one salt bridge, 15 hydrogen bonds and 210 non bonded interactions. The detail of the interaction is illustrated in Fig. 6 and atomic level interaction is provided in Table S4. Similarly, 36 residues of Cov-II-Vac interacted with 41 residues of TLR8, prominent among them were residues of multi-epitopes of 12 hydrogen bonds and 168 non bonded contacts. The interacting residues of the complex are illustrated in Fig. 7 while distances of detailed atomic interaction are provided in Table S5.

NMA analysis

NMA analysis of both multi-epitopes showed that Cov-I-Vac has higher stability than Cov-II-Vac, as the Eigenvalue, the energy required to deform the structure, was higher for Cov-I-Vas as shown in Fig. 8.

Codon optimization and E. coli expression

Codon optimization results of both vaccines Cov-I-Vac and Cov-II-Vac (with GC content 58.99% and 58.49%, respectively) were within the optimum limits (Pandey, Bhatt & Prajapati, 2018). The CAI value was predicted to be 1.0, which projects high level expression of our designed vaccine constructs in E. coli K12 strain. The total length of the Cov-I-Vac clone was 3.6 kbp while Cov- II-Vac clone was 3.8 kbp. Both sequences were designed to be added to the plasmid between restriction sites AfiIII-pciI and BspQI-sapI (Fig. 9).