Designing an efficient multi-epitope vaccine displaying interactions with diverse HLA molecules for an efficient humoral and cellular immune response to prevent COVID-19 infection.

BACKGROUND: The novel SARS-CoV-2 coronavirus, the causative agent of the ongoing pandemic COVID-19 disease continues to infect people globally and has infected millions of humans worldwide. However, no effective vaccine against this virus exists.
METHOD: Using Immunoinformatics, epitopic sequences from multiple glycoproteins that play crucial role in pathogenesis were identified. Particularly, epitopes were mapped from conserved receptor-binding domain of spike protein which have been experimentally validated in SARS-CoV-1 as a promising target for vaccine development.
RESULTS: A multi-epitopic vaccine construct comprising of B-cell, CTL, HTL epitopes was developed along with fusion of adjuvant and linkers. The epitopes identified herein are reported for the first time and were predicted to be highly antigenic, stable, nonallergen, nontoxic and displayed conservation across several SARS-CoV-2 isolates from different countries. Additionally, the epitopes associated with maximum HLA alleles and population coverage analysis shows the proposed epitopes would be a relevant representative of large proportion of the world population. A reliable three-dimensional structure of the vaccine construct was developed. Consequently, docking and molecular-dynamics simulation ensured the stable interaction between vaccine and innate-immune receptor.

Introduction

The recent emergence of highly pathogenic SARS-CoV-2 coronaviruses, the causal agent of COVID-19, has created a global concern for public health [1,2]. Since its first emergence in late 2019 in China, the virus has infected more than 7 million individuals and 414,340 deaths, with exponential growth of numbers in many countries. The clinical symptoms associated with this virus infection includes fever, headache trailed by severe respiratory distress and the most severe pathology in individuals above 60 years of age and especially in those with underlying chronic illness like diabetes and hypertension [3]. Moreover, COVID-19 is highly transmissible and the reproductive number (R0) of the infection is estimated to be 2.2, which indicates that one infected individual either symptomatic or asymptomatic can transmit the virus to 2.2 other individuals in a population with a mean incubation period of 5.8 days [3]. Given the severity of the disease, which is above that of seasonal influenza or pandemic H1N1 influenza and the ability of rapid transmission of virus by asymptomatic individual triggering its ability to cause pandemic outbreak, development of a safe and effective vaccine to fight against this novel virus is urgently needed [4,5]. Despite, frequent periodic outbreaks of human coronavirus infection reported every decade in the twenty-first century like severe acute respiratory syndrome coronavirus (SARS‐CoV-1) epidemic in 2002, middle- east respiratory syndrome coronavirus (MERS‐CoV) in 2012 and now SARS‐CoV-2, till date there are no existing vaccines available in the market to tackle this severe pathogenic virus.

Comparative homology between SARS-CoV-2 and SARS-CoV-1 genomes, both belonging to genus Betacoronavirus and exhibiting genome size of approximately 30 Kb [6] showed high sequence and structure conservation, speculating that both coronaviruses shares the same receptor, ACE2 (angiotensin-converting enzyme 2) found in human lungs [7]. On this basis, the lessons learnt through the challenges faced during SARS-Co-1 vaccine trials like undesired immunopotentiation [8], waning of antibody response after short term [9] should be prioritized while designing SARS-CoV-2 vaccine strategies. Another consideration for effective SARS-CoV-2 vaccine development is the intended target population. Because, the clinical manifestation of COVID-19 disease is more severe in older individuals due to immune senescence [10], it is important to develop vaccine that protects this high-risk segment of the population with reduced immunity. Many leading multinational companies/consortium are presently working relentlessly exploring different vaccine production platforms and technologies for SARS-CoV-2 like whole virus vaccines (Johnson & Johnson, Codagenix/Serum Institute of India); subunit vaccines (Novavax, Clover biopharmaceuticals, University of Queensland/CEPI); nucleic acid vaccines (Inovio/Beijing Advaccine Biotechnology Co./CEPI, CureVac/CEPI/, Moderna/NIH/CEPI) [11]. However, current vaccine development through conventional techniques might take years and each of the strategy has its distinct advantages and disadvantages leading to serious health consequences as adverse side effects [12]. Considering the present scenario of the COVID-19 pandemic worldwide, it is therefore essential to quickly integrate different vaccine design approaches so as to minimize both time and costs in the process of target vaccine candidate identification, which can then be evaluated for efficacy and safety.

The availability of pathogen whole genomes and a pool of robust vaccinomics driven sequence analysis tools have made it possible to computationally search on a genome scale for promiscuous epitopes as potential vaccine candidates which has the ability to elicit desired immune response against the pathogenic organism [13]. Subunit vaccine formulations containing cocktail of defined T-cell and B-cell multi-epitopes which determines the antigenic part of the pathogens have the potential to stimulate both humoral and cellular immune systems within the human body. Besides nullifying the chances of pathogenicity reversal as in case of live attenuated vaccines, the major advantage of employing subunit vaccine includes its efficacy in a person with reduced immunity. Application of epitopic vaccine construct has been successfully demonstrated in several disease conditions like HIV-1 [14], cancer [15-17], and also allergic and autoimmune diseases [18].

The four important glycoproteins that majorly contribute to the structure of all coronaviruses are the spike protein (S), small envelope protein (E), membrane protein (M), and nucleocapsid (N) protein [9]. Particularly, the spike protein, containing putative S1 and S2 domains, mediates coronavirus entry into the host cells by binding to the ACE2 receptor which triggers a cascade of events leading to fusion of both viral and host membranes [19,20]. A recent report on the structure of SARS-CoV-2 spike showed that specifically a fragment located in the middle region of the S1 subunit defined as ‘receptor-binding domain’ (RBD; residue 317–569) of human novel coronavirus ARS-CoV-2 strain spike protein (Uniprot Id P59594) binds to the peptidase domain of ACE2 [21]. Therefore, it is proposed that the RBD region of the spike protein should be used for the development of SARS-CoV-2 vaccine, because it is not only a functional domain for receptor-binding and viral entry but also a critical neutralization determinant of coronaviruses [22-24]. Furthermore, experimental data suggests that RBD-based vaccines may induce broad protection against both human and animal SARS-CoV variants [22]. Besides spike protein, other structural glycoproteins also have diverse functional role in the viral pathogenesis, like the membrane protein mediates assembly and binding of virus particles [25], nucleocapsid protein packages the viral genome into a helical ribonucleocapsid and hence fundamentally assist in viral self-assembly and pathogenies [26,27] while the envelope protein play an important role in the assembly of the virion and in virus replication [28].

Recently, few studies have been reported on epitope identification for SARS-CoV-2 vaccine construct using immunoinformatics approach [29-32]. However, these in silico analysis are majorly on spike protein and most importantly, regardless of the importance of RBD region of spike protein in the pathogenesis of SARS-CoV-2 virus, none have focused on identification of conserved epitopes spanning the antigenic RBD region. An advantage of subunit vaccine comprising of only the RBD of the SARS-CoV spike protein is shown to minimize host immunopotentiation [8]. Consequently, this study aims to design a multi-epitope SARS-CoV-2 vaccine construct using all antigenic surface structural proteins, including the RBD-spike protein, envelope protein, membrane protein, and nucleocapsid protein. Promiscuous vaccine candidates comprising of most potential T-cell and B-cell epitopes from each of the four antigenic protein sequences have been mapped along with analysis of the epitopes conservancy, potential immunogenicity, allergenicity, toxicity, and conservancy analysis against different isolates to ensure formulation of robust widely protective SARS-CoV-2 vaccine with the ability to generate the humoral as well as cell mediated immunity. Physicochemical characterization and population coverage analysis of the final multi-epitope vaccine construct was assessed followed by homology modeling and docking of validated vaccine model with the immune receptor (TLR-3) present on lymphocytes cells. Consequently, the stability of the interaction of the vaccine with the receptor was determined by molecular dynamics simulation studies.

Methodology

Selection of SARS-CoV-2 proteins for potential epitope screening

The major structural glycoproteins that majorly contribute to the structure of SARS-CoV2 and play vital roles in viral pathogenesis were selected for the purpose of vaccine development. The sequences of spike protein (UniProt Id: P59594), envelope protein (UniProt Id: P59637), membrane protein (UniProt Id: P59596), and nucleocapsid protein (UniProt Id: P59595) were retrieved from UniProt database (http://www.uniprot.org/) in FASTA format. Further, in order to avoid cross- reactivity and generation of autoimmune disorder, sequence similarity search of each of the above selected proteins were done to find out the homology with the human proteome by using NCBI-Protein-Protein Blast (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastp&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome).

Mapping of B-cell epitopes

The humoral immune response specific B-cell epitopes were predicted for the aforementioned SARS-CoV-2 proteins by using the ABCpred server (http://crdd.osdd.net/raghava/abcpred/) at the threshold of 0.51 and window length of 10. The server with 65.93% accuracy, utilizes feed forward and recurrent neural network algorithm based on training data set of 700 B-cell epitopes and equal number of non-B-cell epitopes from SwissProt database [33].

Mapping of T-cell epitopes

Cytotoxic T-Lymphocyte (CTL)Epitopes Mapping

The NetMHCpan 4.1 Server (http://www.cbs.dtu.dk/services/NetMHCpan-4.1/) was used to screen the antigenic CTL epitopes for the selected SARS-CoV-2 coronavirus proteins. NetMHCpan server predicts binding of peptides to any MHC molecule of known sequence using artificial neural networks. The binding affinity data covers major twelve supertypes of human, including five major supertypes of HLA-A locus, A1, A2, A3, A24, and A26, and the HLA-B locus of seven major supertypes, B7, B8, B27, B39, B44, B58, and B62 [34] (Supplementary Table S1). For the present study representative HLA alleles belonging to the major twelve HLA supertypes were selected. In addition to NetMHCpan 4.0 Server, NetCTL 2.0 (http://www.cbs.dtu.dk/services/NetCTL/) and IEDB (Immune Epitope Database And Analysis Resource) tool ‘Proteasomal cleavage/TAP transport/MHC Class I combined predictor’ (http://tools.immuneepitope.org/mhci/) servers. The total score generated by NetCTL 2.0 for each CTL epitope gives a combined scoring which includes the proteasomal cleavage of the C-terminal end of the peptide, TAP (N-terminal) transport efficiency, and MHC Class I affinity [35]. IEDB predicts the antigenic CTL epitopes and their binding to MHC Class I alleles based on the IC (50) (nM) score [36]. Epitopes with IC50 values are predicted as high (<50 nM), intermediate (<500 nM) and of least affinity (<5000 nM) for different epitopes. The top three high scoring CTL epitopes that were consensually predicted by all the tools, NetMHCpan 4.1, NetCTL, and IEDB tools for each protein were selected for further analysis.

Helper T-Lymphocyte (HTL) epitopes mapping

The NetMHCIIpan-4.0 server (http://www.cbs.dtu.dk/services/NetMHCIIpan-4.0/) was used to predict the immunogenic HTL epitopes [34]. The server predicts peptide binding to any MHC II molecule of known sequence using Artificial Neural Network. It is trained on an extensive dataset of HLA alleles, covering the three human MHC class II isotypes HLA-DR, HLA-DQ, HLA-DP, as well as the mouse molecules (H-2). All HLA alleles belonging to DRB1, DRB3, DRB4, DRB5, DP, and DQ were selected for the present epitope prediction (Supplementary Table S1). The high scoring epitopes predicted by NetMHCIIpan-4.0 server were also validated using the IEDB tool ‘MHC II Binding predictions’ (http://tools.iedb.org/mhcii/) by selecting IEDB recommended alleles, consisting of a reference panel of 27 HLA alleles (Supplementary Table S1). The epitopes that showed strong binding to maximum number of HLA alleles by both the tools were finally selected for further analysis.

Characterization of screened epitopes

Epitope conservation analysis

Conservancy analysis of all the shortlisted CTL, HTL and B-cell epitopes was performed by ‘Epitope Conservancy Analysis tool of IEDB (http://tools.iedb.org/conservancy/).The tool computes the percentage of epitope sequence conservancy amongst a given set of protein sequences at a given identity level. The higher the percentage of different SARS-CoV-2 strains carrying the screened epitope in their protein sequences, the greater would be its effectiveness as a widely protective vaccine candidate.

Toxicity prediction analysis

The ToxinPred tool (https://webs.iiitd.edu.in/raghava/toxinpred/index.html) that employs support vector machine method [37] was used for the prediction of possible toxicity of all the screened CTL, HTL and B-cell epitopes. The tool allows identifying highly toxic or nontoxic short peptide sequences along with analysis of hydrophobicity, hydropathicity, hydrophilicity, charge, and molecular weight of the predicted peptides.

Antigenicity analysis

The antigenicity of the screened epitopes was calculated using VaxiJen 2.0 [http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html]. VaxiJen uses an alignment-free approach to identify antigens directly [38]. Instead of concentrating on predicting antigenic regions, the method uses statistical models based on discriminating positive and negative sets of bacterial, viral, and tumor antigen datasets to predict whole-protein antigenicity. VaxiJen has shown high prediction accuracy of up to 89%.

Allergenicity analysis

Allergens are proteins or glycoproteins recognized by immunoglobulin E (IgE), which is produced by the immune system in allergic individuals. It is said that a protein is considered an allergen if it has a region or peptides identical to a human IgE epitope [39]. Allergenicity of the epitopes was checked by the AllerTOPV 2.0 tool [https://www.ddg-pharmfac.net/AllerTOP/index.html]. The tool has been developed based on the robust k-nearest neighbors (kNN) method and predicts allergenicity/nonallergenicity with 86.7% sensitivity, 90.7% specificity, and 88.7% accuracy [40].

Criteria for selection of final multi-epitopes for vaccine construct

For each of the four proteins viz., spike, membrane, nucleocapsid, and envelope, one HTL epitope, one B-cell epitopes and 3 CTL epitopes (that were consensually predicted as high scoring epitopes binding to maximum number of HLA alleles by the NetMHCpan 4.1, NetCTL as well as IEDB tool) were selected for final vaccine development based on the following criteria (1) highest antigenicity score (2) high score by the predicted algorithm, (3) promiscuous peptides, which means epitopes binding to numerous histocompatibility alleles of the world population, (4) nonallergenicity and nontoxic (5) conservancy above 80% across multiple isolates and strains.

Vaccine development

For devising the multi-epitope subunit vaccine, the final selected B-cell, CTL, and HTL epitopes of corresponding spike, envelope, membrane, and nucleocapsid proteins were combined together in the same sequential manner to obtain the final vaccine construct. The N-terminal of the designed vaccine sequence was adjuvanted with β-defensin, which is known to exhibit antimicrobial and antiviral activity and the ability to enhance protective immune response [41]. Additional sequences that play a vital role in the chimeric vaccine construct are linkers, acting as spacers that facilitate protein folding and stability to achieve a proper epitope display for binding to MHC Class I and MHC Class II molecules [42]. In this study, the helix forming linker ‘EAAAK’ was added to join the adjuvant and the first B-cell epitopes, while ‘KK,’ ‘AAY,’ and ‘GPGPG’ linkers were added to join together all B-cell epitopes, CTL epitopes and HTL epitopes, respectively. These linkers were appended to attain the efficient separation of individual epitopes in vivo [43].

Population coverage analysis of vaccine construct

The predicted epitopes constituting the vaccine construct with the corresponding HLA alleles (Class I and Class II) were submitted to the population coverage analysis tool of the IEDB database (http://tools.iedb.org/population/) by keeping the default parameters on (109 countries covering 16 different geographical regions) (Supplementary Table S2). Population coverage analysis tool calculates the fraction of individuals predicted to respond to a given set of epitopes with known MHC restrictions. For individual epitope, the tool computes the following: (i) projected population coverage, (ii) average number of epitope hits/HLA combinations recognized by the populations, and (iii) minimum number of epitope hits/HLA combinations recognized by 90% of the population (PC90). These calculations are made on the basis of HLA genotypic frequencies assuming nonlinkage disequilibrium between HLA loci [44].

Physicochemical features and secondary structure prediction of the vaccine construct

Physicochemical properties of the vaccine construct comprising of amino acid composition, molecular weight, number of positive/negatively charged amino acids, instability index, aliphatic index, GRAVY, hydropathy plot, and solubility was analyzed using the ProtParam tool (https://web.expasy.org/protparam/) [45]. Instability index was calculated to estimate the stability of proteins in a test tube. The aliphatic index shows the relative volume occupied by aliphatic side chains and GRAVY is the average of hydropathy values of amino acids in a protein [45].The secondary structure of the predicted vaccine was carried out by online tool PSIPRED (http://bioinf.cs.ucl.ac.uk/) and the percentage of helix, strand, and coils were determined.

Tertiary structure prediction and validation of the vaccine construct

The tertiary structure of vaccine construct was generated using Robetta server [46]. The best 3D model developed was analyzed further using various structure validation tools. PROCHECK checks the geometrical and stereochemical constraints of the protein structure by inspecting the accuracy of the dihedral angles (i.e. u and ψ) in the Ramachandran plot [47]. VERIFY3D analyzses the compatibility of the 3D protein model with its amino acid sequence [48] while ERRAT examines the nonbonded atomic interactions [49]. A higher score from VERIFY3D and ERRAT indicates a more reliable structure for further analysis.

Analysis of molecular interaction between vaccine candidate and immune cell receptor

Molecular docking between the receptor Toll-like receptor 3 (TLR-3) and vaccine was performed with the help of the ClusPro server (https://cluspro.org), which is a widely used tool for protein–protein docking [50]. Subsequently, molecular dynamics (MD) simulations was performed to investigate the dynamical motion and stability of the docked vaccine-receptor complex with the help of GROMACS version2018.4 molecular dynamics package [51,52] using Amber ff99SB-ILDN force field in TIP3P water model. The system net charge neutralization was done by adding required Na+ and Cl- ions. After electro-neutralization, the solvated model was subjected to energy minimization using the steepest descent algorithm until the maximum force of 1000 kJ mol-1 nm-1 had been reached. Following energy minimization, using an NVT ensemble at 271 K temperature equilibration was carried out via the V-rescale method and the pressure was equilibrated using an NPT ensemble at 1 bar with the Parinello–Rahman algorithm for 2 ns. After equilibration, final production MD simulation was carried out for 50 ns. The particle mesh Ewald method was used to treat long-range electrostatic interactions. The short-range electrostatic and van der Waals interactions were calculated by specifying a cutoff 12-Å distance. The bonds involving hydrogen atoms were constrained using the linear-constraint-solving algorithm. The periodic boundary conditions were applied during production MD simulation. Trajectory data harvested from production MD was used to compute the dynamics stability of the system. The root mean square deviations (RMSD) for protein backbone and root mean square fluctuation (RMSF) for Cα-atoms were generated using gmx rms and gmx rmsf utility tools. The plots were rendered using grace program (http://plasma-gate.weizmann.ac.il/Grace/).

Results

Selection of B-cell epitopes

Recent findings have demonstrated that neutralizing antibodies nAbs targeting RBD regions of spike protein is potential for prophylaxis and treatment of human coronavirus infection including SARS, MERS and also COVID-19 other by inhibiting receptor binding and viral entry into the host cell [53-56]. Therefore, in this study, B-cell epitopes from the specific RBD region of spike protein as well as from other antigenic proteins, viz., nucleocapsid, envelope, and membrane proteins were predicted using ABCPred server. The high scoring peptides were then checked for toxicity and allergenicity using ToxinPred and AllerTOP V 2.0 tool, respectively. Out of many putative epitopes, only few peptide sequences for each protein were found to be nonallergen and nontoxic to human, and were considered for further analysis (Supplementary Table S3). Finally, the best B-cell epitope from each of the antigenic protein, with highest antigenicity score as predicted by VaxiJen server was used for vaccine construction (Table 1). Only one epitope DISNVPFSPD was predicted to fall within the RBD region (residue 317–569) of the spike protein and the same was also found to possess highest antigenicity, thereby indicating as a potential vaccine candidate. Likewise, epitopes VFLLVTLAIL, IGFLFLAWIM, and GDGKMKELSP were selected from Small Envelope Protein, membrane protein and nucleocapsid protein, respectively. The molecular weight, charge, hydrophobicity, hydropathicity, and hydrophilicity values of the selected epitopes are given in Table 1.

Table 1.

Predicted B-cell epitopes from SARS-CoV-2 proteins and their corresponding immunogenic properties to design multi-epitope vaccine construct.

Uniprot_ID	B-Cell Epitope	Position	Score	Conservancy	Toxicity	Hydrophobicity	Hydropathicity	Hydrophilicity	Charge	Mol wt.
Spike Protein (P59594)	DISNVPFSPD	454	0.74	80.00%	Nontoxin	−0.09	−0.38	0.1	−2	1090.28
Small Envelope Protein (P59637)	VFLLVTLAIL	25	0.83	100.00%	Nontoxin	0.46	3.2	−1.54	0	1101.6
Membrane Protein (P59596)	IGFLFLAWIM	23	0.77	80.00%	Nontoxin	0.48	2.46	−1.74	0	1210.7
Nucleocapsid Protein(P59595)	GDGKMKELSP	98	0.72	90.00%	Nontoxin	−0.28	−1.23	0.92	0	1061.36

Selection of CTL epitopes

CTL epitopes from the four major SARS-CoV-2 antigenic proteins were identified using the NetMHCpan 4.1 server, NetCTL server and the IEDB tool. Nontoxin, nonallergen, ten overlapping high scoring epitopes predicted by both the algorithms were selected for each protein for antigenicity analysis (Supplementary Table S4). Finally based on highest antigenicity score, three epitopes were selected for each protein. Notably CTL epitopes predicted for spike protein were VLSFELLNA, TSTGNYNYK, and VVVLSFELL; envelope protein epitopes were RLCAYCCNI, VVFLLVTLA, LAILTALRL; membrane protein epitopes were NLVIGFLFL, GTITVEELK, RINWVTGGI; and for nucleocapsid protein CTL epitopes were KLDDKDPQF, ELSPRWYFY, and DPQFKDNVI (Table 2). The affinity between identified CTL epitopes and their corresponding major histocompatibility complex (MHC) class I alleles alongwith their immunogenic properties is shown in Table 2.

Table 2.

Predicted CTL epitopes from SARS-CoV-2 proteins to design multi-epitope vaccine construct with their corresponding MHC Class I alleles and their immunogenic properties.

Proteins	Epitope	High affinity HLA Alleles binding to predicted epitopes as analyzed using the NetMHCpan 4.1 server	NetCTL Score	NetMHCpan 4.1 Score	Position	Conservancy	Toxicity	Hydrophobicity	Hydropathicity	Hydrophilicity	Charge	Mol wt.
Spike Protein (P59594)	VLSFELLNA	HLA-A*02:01, HLA-A*01:01, HLA-A*03:01, HLA-A*24:02, HLA-A*26:01, HLA-B*07:02, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	0.8431	0.7462	498–506	88.89%	Nontoxin	0.14	1.08	−0.64	−1	1102.44
	TSTGNYNYK	HLA-A*03:01, HLA-A*02:01, HLA-A*01:01, HLA-A*24:02, HLA-A*26:01, HLA-B*07:02, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	0.8456	1.0140	431–439	80.00%	Nontoxin	−0.28	−1.74	−0.4	1	1210.4
	VVVLSFELL	HLA-A*24:02, HLA-A*02:01, HLA-A*01:01, HLA-A*03:01, HLA-A*26:01, HLA-B*07:02, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	0.5299	3.1896	496–504	100.00%	Nontoxin	0.23	1.9	−0.89	−1	1132.52
Small Envelope Protein (P59637)	RLCAYCCNI	HLA-A*02:01, HLA-A*01:01, HLA-A*03:01, HLA-A*24:02, HLA-A*26:01, HLA-B*07:02, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	0.8192	6.1136	38–46	100.00%	Nontoxin	−0.02	1.25	−0.77	1	1157.56
	VVFLLVTLA	HLA-A*02:01, HLA-A*01:01, HLA-A*03:01, HLA-A*24:02, HLA-A*26:01, HLA-B*07:02, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	0.3465	3.9037	24–32	100.00%	Nontoxin	0.46	3.24	−1.51	0	1087.57
	LAILTALRL	HLA-B*58:01, HLA-A*02:01, HLA-A*01:01, HLA-A*03:01, HLA-A*24:02, HLA-A*26:01, HLA-B*07:02, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*15:01	0.5001	1.6850	31–39	100.00%	Nontoxin	0.15	2.06	−0.84	1	1086.55
Membrane Protein (P59596)	NLVIGFLFL	HLA-A*02:01, HLA-A*01:01, HLA-A*03:01, HLA-A*24:02, HLA-A*26:01, HLA-B*07:02, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	0.7903	1.0027	20–28	100.00%	Nontoxin	0.38	2.36	−1.4	0	1106.53
	GTITVEELK	HLA-A*03:01, HLA-A*02:01, HLA-A*01:01, HLA-A*24:02, HLA-A*26:01, HLA-B*07:02, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	0.6722	0.7941	5–13	100.00%	Nontoxin	−0.14	−0.37	0.33	−1	1117.42
	RINWVTGGI	HLA-B*07:02, HLA-A*02:01, HLA-A*01:01, HLA-A*03:01, HLA-A*24:02, HLA-A*26:01, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	0.5665	4.0169	71–79	88.89%	Nontoxin	0.04	0.46	−0.62	1	1086.41
Nucleocapsid Protein(P59595)	KLDDKDPQF	HLA-A*02:01, HLA-A*01:01, HLA-A*03:01, HLA-A*24:02, HLA-A*26:01, HLA-B*07:02, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	0.5527	0.4409	339–347	88.89%	Non-Toxin	−0.51	−2.07	1.39	0	1233.52
	ELSPRWYFY	HLA-A*26:01, HLA-A*02:01, HLA-A*01:01, HLA-A*03:01, HLA-A*24:02, HLA-B*07:02, HLA-B*08:01, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	1.0943	0.0760	104–112	88.89%	Nontoxin	−0.11	−0.86	−0.83	0	1423.72
	DPQFKDNVI	HLA-B*08:01, HLA-A*02:01, HLA-A*01:01, HLA-A*03:01, HLA-A*24:02, HLA-A*26:01, HLA-B*07:02, HLA-B*27:05, HLA-B*39:01, HLA-B*40:01, HLA-B*58:01, HLA-B*15:01	0.4606	0.6519	344–352	80.00%	Nontoxin	−0.15	−0.42	0.18	−1	1188.49

Selection of HTL epitopes

The HTL epitopes were predicted using the NetMHCIIpan-4.0 server and the IEDB server as described in the methodology section. Out of many epitopes identified (Supplementary Table S5), based on antigenicity score and binding affinity to maximum number of HLA alleles as predicted by both the tools, the final four HTL epitopes selected each for spike, envelope, membrane, and nucleocapsid proteins were VVLSFELLNAPATVC, LLFLAFVVFLLVTLA, IGFLFLAWIMLLQFA, and PRWYFYYLGTGPEAS, respectively (Table 3).

Table 3.

Predicted HTL epitopes from SARS-CoV-2 proteins to design multi-epitope vaccine construct with their corresponding MHC Class II alleles and their immunogenic properties.

Proteins	MHC II Epitope	High affinity HLA Alleles binding to predicted epitopes as analyzed using the NetMHCIIpan-4.0 server	Pos	IC50 value	NetMHCIIpan-4.0 Score	Conservancy	Toxicity	Hydrophobicity	Hydropathicity	Hydrophilicity	Charge	Mol wt.
Spike Protein (P59594)	VVLSFELLNAPATVC	DRB1*01:12, DRB1*01:15, DRB1*01:17, DRB1*01:24, DRB1*01:25, DRB1*01:27, DRB1*03:31, DRB1*03:46, DRB1*03:48, DRB1*04:27, DRB1*04:41, DRB1*07:01, DRB1*09:02, DRB1*11:35, DRB1*11:65, DRB1*11:78, DRB1*12:03, DRB1*12:18, DRB1*13:10, DRB1*13:12, DRB1*14:52, DRB1*14:53, DRB1*15:26, DRB1*15:44, DRB1*16:14, DRB3*01:08, DRB3*01:14, DRB3*02:04, DRB3*02:23, DRB3*02:25, DRB4*01:04, DRB4*01:07, DRB4*01:08, DRB5*01:02, DRB5*01:03, DRB5*02:05, DPA1*02:04 DPB1*48:01, DPA1*02:02 DPB1*118:01, DPA1*04:01DPB1*69:01, DQA1*01:04DQB1*03:07, DQA1*03:03 DQB1*06:35, DQA1*04:04DQB1*05:01	497–511	4	0.50	93.33%	Nontoxin	0.17	1.52	−0.75	−1	1576.1
Small Envelope Protein (P59637)	LLFLAFVVFLLVTLA	DRB1*01:03, DRB1*01:15, DRB1*01:19, DRB1*01:24, DRB1*01:27, DRB1*01:35, DRB1*03:38, DRB1*03:46, DRB1*03:48, DRB1*04:29, DRB1*04:41, DRB1*07:05, DRB1*09:02, DRB1*11:35, DRB1*11:65, DRB1*11:78, DRB1*12:07, DRB1*12:12, DRB1*13:10, DRB1*13:12, DRB1*14:52, DRB1*14:53, DRB1*15:21, DRB1*15:44, DRB1*16:16, DRB3*01:05, DRB3*01:12, DRB3*02:04, DRB3*02:12, DRB3*02:23, DRB3*02:25, DRB3*03:03, DRB4*01:01, DRB4*01:03, DRB4*01:06, DRB4*01:07, DRB4*01:08, DRB5*01:02, DRB5*01:08 N, DRB5*02:02, DRB5*02:03, DPA1*02:04 DPB1*48:01, DPA1*02:02 DPB1*118:01, DPA1*04:01DPB1*69:01, DQA1*01:04DQB1*03:07, DQA1*03:03 DQB1*06:35, DQA1*04:04DQB1*05:01	18–32	7	0.40	100.00%	Nontoxin	0.46	3.11	−1.61	0	1679.4
Membrane Protein (P59596)	IGFLFLAWIMLLQFA	DRB1*01:07, DRB1*01:13, DRB1*01:19, DRB1*01:26, DRB1*01:32, DRB1*01:35, DRB1*03:43, DRB1*03:46, DRB1*03:54, DRB1*04:03, DRB1*04:04, DRB1*04:16, DRB1*04:39, DRB1*04:47, DRB1*07:03, DRB1*08:06, DRB1*08:07, DRB1*08:34, DRB1*09:07, DRB1*10:02, DRB1*11:45, DRB1*11:59, DRB1*11:74, DRB1*12:07, DRB1*12:16, DRB1*13:13, DRB1*13:33, DRB1*14:42, DRB1*14:53, DRB1*15:21, DRB1*15:44, DRB1*16:16, DRB3*01:05, DRB3*01:12, DRB3*02:04, DRB3*02:12, DRB3*02:23, DRB3*02:24, DRB3*03:01, DRB3*03:03, DRB4*01:01, DRB4*01:03, DRB4*01:06, DRB4*01:07, DRB4*01:08, DRB5*01:02, DRB5*01:08 N, DRB5*01:14, DRB5*02:03, DRB5*02:04, DPA1*02:04 DPB1*48:01, DPA1*02:02 DPB1*118:01, DPA1*04:01DPB1*69:01, DQA1*01:04DQB1*03:07, DQA1*03:03 DQB1*06:35, DQA1*04:04DQB1*05:01	23–37	17	0.50	86.67%	Nontoxin	0.4	2.22	−1.59	0	1783.49
Nucleocapsid Protein(P59595)	PRWYFYYLGTGPEAS	DRB1*01:06, DRB1*01:11, DRB1*01:24, DRB1*01:29, DRB1*01:30, DRB1*01:36, DRB1*03:38, DRB1*03:40, DRB1*03:49, DRB1*040, DRB1*04:16, DRB1*04:34, DRB1*07:03, DRB1*07:11, DRB1*08:05, DRB1*08:21, DRB1*09:03, DRB1*09:09, DRB1*10:01, DRB1*10:02, DRB1*11:70, DRB1*11:83, DRB1*12:03, DRB1*12:18, DRB1*13:10, DRB1*13:12, DRB1*14:52, DRB1*14:53, DRB1*15:26, DRB1*15:44, DRB1*16:14, DRB3*01:13, DRB3*01:14, DRB3*02:07, DRB3*02:16, DRB3*02:23, DRB3*03:03, DRB4*01:04, DRB4*01:06, DRB4*01:07, DRB5*01:02, DRB5*01:01, DRB5*01:05, DRB5*01:11, DRB5*02:04, DPA1*03:02 DPB1*46:01, DPA1*02:02 DPB1*118:01, DPA1*04:01DPB1*69:01, DQA1*01:04DQB1*03:07, DQA1*03:03 DQB1*06:35, DQA1*04:04DQB1*05:01	107–121	44	0.50	93.33%	Nontoxin	−0.05	−0.66	−0.61	0	1807.2

Epitope conservation analysis

Genomes of 96 recently sequenced SARS-CoV-2 isolates from different origin were received from NCBI and their corresponding spike, membrane, envelope and nucleocapsid protein sequences (Supplementary Table S6). Consequently, each of the shortlisted CTL, HTL, and B-cell epitopes were analyzed against the corresponding homologous protein sequences using the Conservancy Analysis tool of IEDB. It was observed that all selected epitopes exhibited conservancy of more than 80% in all isolates from different counties (Tables 1–3) while eight antigenic epitopes out of total twenty showed complete 100% conservation.

Designing a multi-epitope chimeric vaccine construct

A major challenge for current vaccine development is the fact that many new subunit vaccines are poorly immunogenic and elicit insufficient immune responses for protective immunity. Adjuvants are therefore needed in vaccine formulations to enhance the immune response and for proper delivery of vaccine antigens to antigen-presenting cells [42]. Hence, β-defensin (TLR-3 agonist) was added at the N-terminal region of the vaccine sequence to induce humoral and cellular immunity [57]. Additional sequences that play an important role in designing of a vaccine construct are linkers which acts as spacers for facilitating folding in order to achieve a proper epitope display [58]. Finally, the designed chimeric vaccine construct was composed of adjuvant, 4 B-cell epitopes, 12 CTL epitopes, 4 HTL epitopes linked together in the same order utilizing EAAAK, KK, AAY, and GPGPG linkers, respectively. The schematic illustration of the multi-epitope vaccine construct is shown in Figure 1.

Figure 1.

Structural arrangement of the final vaccine construct.

MHC molecules are extremely polymorphic in their peptide binding regions, and as a result, they exhibit widely varying binding specificity. Therefore, population coverage analysis is crucial for design of applicable epitope-based subunit vaccines that bind to several alleles of HLA supertypes [59]. In order to calculate the population coverage of the predicted putative epitopes, the predicted putative epitopic core sequences along with their binding MHC Class I and MHC Class II alleles were submitted to the population coverage analysis tool housed at the Immuno Epitope Database.

The IEDB results corroborated that with the administration of designed vaccine construct, 98.81% of the world population would be covered, which encompasses 109 countries covering 16 different geographical regions (Figure 2; Supplementary Table S2). It was noted that the population coverage rate exceeded 90% in almost all countries across worldwide (Figure 2; Supplementary Table S2).

Figure 2.

Population coverage analysis of the final multi-epitope vaccine construct across world as predicted by the population coverage analysis tool of the IEDB database (http://tools.iedb.org/population/).

The physicochemical characteristics, such as the molecular weight, stability, hydrophilicity, charge of the vaccine construct is crucial for their uptake, processing and antigen presentation, and the concomitant induction of immunity [42]. The molecular weight of the final vaccine construct was calculated to be 37.57 kDa that ensures it has good antigenicity. The vaccine candidate was found to be basic in nature due to its 9.26 theoretical pI value. It is composed of 19 negatively charged amino acid residues and 33 positively charged amino acid residues. It is observed that charge affects the physical stability of the formulation and the degree of antigen association. Specifically, cationic particles generally interact more efficiently with antigen presenting cells than neutral and anionic particles due to electrostatic binding to cell-surface heparane sulfate proteoglycans, thereby enhancing antigen acquisition [42]. The Aliphatic index value was found to be 109.33 which defined the high thermal stability of the vaccine in experimental set up. The estimated half-life of designed vaccine construct was predicted to be 30 hours in mammalian reticulocytes, in vitro and greater than 20 hours and 10 hours in yeast and Escherichia coli, in vivo, respectively. The vaccine construct was predicted to be highly stable with Instability index value less than 40 (31.12 Instability Index) (Supplementary Table S7).

From the PSIPRED program output, it was found that secondary structure of the multi-epitope vaccine constitutes 44.44% alpha helix, 13.45% extended beta sheet, and 42.1% coils as shown in Supplementary Figure S1.

Three-dimensional structure of the multi-epitope final vaccine was modeled by using Robetta server. The final structure developed has been shown in Figure 3(a). The quality of the predicted 3D model was assessed using various structure validation tools. The calculated Ramachandran plot using PROCHECK revealed that all residues fall in allowed regions, with 92.8%, 7.2% in most favored region and additional allowed region, respectively and no single amino acid present in disallowed region. Further high scoring values, 82.7% in VERIFY3D and 89.12% in ERRAT conclusively demonstrated the accuracy and reliability of the three-dimensional structure of the multi epitopic vaccine construct for further studies (Supplementary Figure S2).

Figure 3.

Molecular three-dimensional structure of (a) modeled multi-epitope vaccine construct (b) docked vaccine-immune receptor (TLR-3) complex. (c) Surface View of Docked Vaccine-Immune receptor. (d) Interacting Residues of Docked Vaccine-Immune receptor.

It has been demonstrated that TLR-3 contributes to a protective innate immune response to severe Acute Respiratory Syndrome Coronavirus infection [60]. Therefore, the final vaccine construct was docked with TLR-3 (PDB Id: 2A0Z) using Cluspro server. The ribbon representation and surface view of the docked complex with lowest energy score of −919.8 as viewed using Pymol is shown in Figure 3(b,c). High binding affinity was observed between the vaccine construct and the immune receptor as displayed by multiple strong hydrogen bond (H-bond) interactions with less distance between the interacting residues. The residues of docked vaccine-TLR3 complex showing H-bond interactions were MET1-HIS60, SER117-PHE644, ARG329-HIS682, TRP330-LEU690, PHE12-TYR283, SER117-ASN678, THR223-PRO681, ALA302-THR26, ARG329-TYR683, THR300-THR26, and TYR162-LYS41 with distance of 2.39 Å, 1.69 Å, 2.54 Å, 2.08 Å, 1.79 Å, 2.62 Å, 2.42 Å, 2.49 Å, 2.48 Å, 2.64 Å, and 2.26 Å, respectively (Figure 3(d)). The stability of the docked complex was determined by performing molecular dynamics simulation of 50 ns using GROMACS. The stable RMSD (Figure 4) after 40 ns affirmed the protein-receptor complex attained equilibrium and tend to display a stable trend till 50 ns. While the RMSF plot of the vaccine showed more fluctuations than receptor (TLR3) indicating that the structure of the vaccine made more movements to refine the interaction with receptor in order to elicit higher immune response.

Figure 4.

Root mean square deviation (RMSD) and root mean square fluctuation (RMSF) analysis of protein backbone and side chain residues of MD simulated vaccine construct.

Discussion

SARS-CoV-2 continues to infect people globally, even under the current public health emergency with the concomitant urgency to identify appropriate therapeutic targets to prevent and treat COVID-19 infection and control its spread [61,62]. Although SARS-CoV-1 and MERS-CoV coronaviruses outbreak originated almost a decade before, only a small number of SARS-CoV-1 and MERS-CoV reached to phase I clinical trials before funding dried up and eradication of the virus was possible through nonpharmaceutical interventions and preventive measures [63]. Nevertheless, few of the advances made in developing vaccines and therapeutics for SARS-CoV and MERS-CoV could be exploited for countering current SARS-CoV-2 [11,63]. Furthermore, besides identification of vaccine candidates, it is equally important to understand the key steps needed before bringing a SARS-CoV-2 vaccine into clinical trials. Firstly, the vaccine construct is tested in appropriate animal models to confirm it is protective. In case of SARS-CoV-2, it is seen that the virus does not grow in wild-type mice and only confers mild disease in transgenic animals. One alternative approach to validate the vaccine is by in vitro neutralization assays using serum from vaccinated animal. Secondly, toxicology analysis of vaccine in animal models, e.g. in rabbits has to be performed before commencing the clinical development trials [11,63].

Earlier experiments have stated that developing effective neutralizing antibodies generated by B-cell epitopes have greater potency to neutralize coronaviruses infection [53]. Besides the humoral response, the role of cell mediated immune response mediated by T- cells in combating viral infections has been known to be equally important for protecting against re-infection. Whilst neutralizing antibodies can prevent viral entry, the body also requires SARS-CoV specific CD4 + T helper (HTL) cells for the development of these specific antibodies. Similarly, CD8+ cytotoxic T cells (CTL) are essential for the recognition and killing of infected cells [64,65]. Experiments on SARS-CoV virus-infected mice models have suggested that inefficient immune activation and a poor virus-specific T cell response underlay severe infection [66,67].

The four major structural glycoproteins namely spike protein, envelope protein, membrane protein and nucleocapsid protein is known to play vital roles in coronavirus pathogenies and greater potency to fight against coronaviruses infection. Amongst the SARS-CoV structural proteins, the spike protein has been found to be primarily responsible for eliciting neutralizing antibodies with its major immunodominant epitopes found specifically in the receptor binding region of the protein which is responsible for interaction with host receptor, ACE2 [68-70]. Besides the spike glycoprotein, persistently high levels of neutralizing antibodies for nucleocapsid protein and T cell responses were also seen in the SARS-recovered individuals after two years of post-infection [71,72]. Similarly, Tsao et al. showed several epitopes in the nucleocapsid protein they could induce specific CTL responses in transgenic mice models. In another study conducted by Li et al. showed that 11% of their SARS subjects gave positive T-cell responses against epitopes present in the nucleocapsid protein [73]. Additionally, animal studies using DNA vaccines suggest that the membrane protein of coronaviruses may also induce T-cell response [74,75]. These comprehensive studies demonstrate that cocktail of immunogenic epitopes from spike, membrane, envelope, and nucleocapsid proteins could be ideal for vaccine development against the novel SARS-CoV-2 coronavirus strain.

The search for potential immunogenic epitopes to be used as vaccine candidates by conventional techniques is time consuming and too expensive. Now a days, immunologists are frequently using in silico approaches to identify the most promiscuous antigenic epitopes which has high affinity to maximum HLA alleles before being subjected to wet-lab confirmatory analysis and thus accelerating the experimental planning in development of epitope-based vaccines. In the present study, using immunoinformatic approaches, a vaccine construct comprising of a cocktail of potential T-cell and B-cell epitopes from multiple antigenic proteins from SARS-CoV-2 have been predicted which would elicit both cellular and humoral immune responses to combat COVID-19 disease. Schematic representation of the workflow for the development of multi-epitope vaccine is shown in Figure 5. The selected epitopes were found to be highly antigenic, nonallergenic and nontoxic for human application. Further the promiscuous peptide sequences showed high conservation against different isolates to ensure formulation of robust widely protective SARS-CoV-2 vaccine. BLAST screening revealed that all the epitopes were nonself, since they lack any significant similarity with human proteome. The physicochemical characteristics, such as the molecular weight, stability, hydrophilicity, of the vaccine construct was optimum for their uptake and evoking immunity. Additionally, the selected epitopes were observed to promiscuously bind with high affinity to maximum HLA alleles, and indicating that the vaccine construct to exhibit approximately 97.97% of the world population coverage. The epitopic sequences were linked together with adjuvant and linkers to enhance immunogenicity and further the three-dimensional structure of the vaccine construct was developed using Robetta server. Various structure validation tools like Procheck, ERRAT, and VERIFY3D confirmed the reliability and accuracy of the developed structure. Efficacy of this vaccine was explored by molecular docking and dynamic simulation with TLR-3 immune receptor. Molecular docking and simulation analysis confirmed the strong and stable binding between vaccine and receptors with multiple hydrogen bond interactions.

Figure 5.

Schematic representation of the workflow for the development of multi-epitope vaccine against SARS-CoV-2 infection.

Conclusion

The present work successfully identified novel peptides of virulent, essential and antigenic proteins, which can act as promising targets for vaccine development against SARS-CoV-2. The methodology employed utilizing immunoinformatic approach can be applied as framework for future characterization of the B-cell and T-cell epitopes of emerging novel coronavirus strains as well as any other pathogenic microorganisms, and which would be of great help for advancing the peptide-based vaccine approach. However, the immune stimulation potential of the predicted multi-epitopic chimeric vaccine construct is needed to be validated in animal models for their efficacy and safety. In the present scenario of not having any direct antiviral agent and vaccines against SARS-CoV-2, strict implementation of social distancing measures and hygiene safety precautionary measures is of utmost importance to check the further spread and control of this virus.

Expert Opinion

The promiscuous epitopes capable of inducing humoral and cellular immune response can be considered as potential vaccine candidates against SARS-CoV-2.

Click here for additional data file.