Bioinformatics analysis of epitope-based vaccine design against the novel SARS-CoV-2

Background An outbreak of infection caused by SARS-CoV-2 recently has brought a great challenge to public health. Rapid identification of immune epitopes would be an efficient way to screen the candidates for vaccine development at the time of pandemic. This study aimed to predict the protective epitopes with bioinformatics methods and resources for vaccine development. Methods The genome sequence and protein sequences of SARS-CoV-2 were retrieved from the National Center for Biotechnology Information (NCBI) database. ABCpred and BepiPred servers were utilized for sequential B-cell epitope analysis. Discontinuous B-cell epitopes were predicted via DiscoTope 2.0 program. IEDB server was utilized for HLA-1 and HLA-2 binding peptides computation. Surface accessibility, antigenicity, and other important features of forecasted epitopes were characterized for immunogen potential evaluation. Results A total of 63 sequential B-cell epitopes on spike protein were predicted and 4 peptides (Spike315–324, Spike333–338, Spike648–663, Spike1064–1079) exhibited high antigenicity score and good surface accessibility. Ten residues within spike protein (Gly496, Glu498, Pro499, Thr500, Leu1141, Gln1142, Pro1143, Glu1144, Leu1145, Asp1146) are forecasted as components of discontinuous B-cell epitopes. The bioinformatics analysis of HLA binding peptides within nucleocapsid protein produced 81 and 64 peptides being able to bind MHC class I and MHC class II molecules respectively. The peptides (Nucleocapsid66–75, Nucleocapsid104–112) were predicted to bind a wide spectrum of both HLA-1 and HLA-2 molecules. Conclusions B-cell epitopes on spike protein and T-cell epitopes within nucleocapsid protein were identified and recommended for developing a protective vaccine against SARS-CoV-2.


Background
An outbreak of infection caused by a novel coronavirus has spread worldwide rapidly [1,2]. The World Health Organization (WHO) named the disease COVID-19, short for "coronavirus disease 2019". As of 27 May 2020, the WHO reported a total of 5 488 825 COVID-19 cases and 349 095 deaths globally [3], which brought a great challenge to public health worldwide. Therefore, it is imminent to prevent and control this infectious disease.
The pathogen causing the new type of pneumonia was named as severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2, also called 2019-nCoV) by the Coronaviridae Study Group (CSG) of the International Committee on Taxonomy of Viruses [4]. Its genome sequence has been released and reported by Chinese scientists and submitted to the GenBank database on 12 January 2020 [5]. Like severe acute respiratory syndrome associated coronavirus (SARS-CoV) and Middle East respiratory syndrome coronavirus (MERS-CoV), the other two viruses that caused severe epidemic problems in recent years, SARS-CoV-2 also belongs to β-coronaviruses family. Bats are proved to be their natural host [6]. At present, insufficient knowledge of the latency and contagiosity of SARS-CoV-2 increased the uncertainty of virus persistency. Specific therapeutic agents targeting the virus are currently not available. Vaccination is still the most economic and effective approach to prevent virus infection. The selection and design of protective immunogens against pathogens is a major challenge in vaccine development, especially for the newly emerging pathogens [7,8]. Traditional methods based on lab experiments could not meet the needs of the pressing situation in the event of an outbreak [9]. Bioinformatics is an interdisciplinary field specialized in organizing, storing, and processing large amounts of data generated from biological experiments. Accumulation of largescale immunological data gave rise to the field known as immunoinformatics, which provides insights into the mechanisms of immune function. As the genome and protein sequence information of SARS-CoV-2 is available, characteristics of the virus, as well as the epitopes presented in the pathogen, could be predicted by in silico analysis, which will greatly speed up the vaccine development [10][11][12].
The purpose of this study is to predict B-cell epitopes on spike protein and T-cell epitopes within nucleocapsid protein of SARS-CoV-2 by applying the bioinformatics methods and immunoinformatic tools. The step-by-step procedure of in silico analysis is depicted in Fig. 1. Epitopes information presented by this work may aid in developing a promising vaccine against SARS-CoV-2.

Data retrieval and sequence alignment
Protein sequences of SARS-CoV-2 were retrieved from the NCBI database (YP_009724390, YP_009724397). Clustal Omega is the current standard version of the Clustal family which was widely used for biological sequences alignment. The program uses seeded guide trees and Hidden Markov model engine to generate Fig. 1 Study workflow. Suitable proteins of SARS-CoV-2 were selected at the first step for epitope prediction. The second step comprised of Band T-cell epitope analysis with bioinformatics approaches. Epitope evaluation was followed and appropriate ones were chosen for vaccine design alignments [13]. In this study, alignment of protein sequences was performed on the EMBL-EBI server with Clustal Omega program. Conserved domains within predicted polypeptides of SARS-CoV-2 were analyzed with CD-search on the NCBI website.

Linear B-cell epitope prediction
ABCpred [14] and BepiPred [15] servers were employed for B-cell epitope forecast. We used a threshold value of 0.85 to achieve a sensitivity between 95.5% and 99.5% for epitope prediction on ABCpred server [14]. The length of linear B-cell epitopes normally varies from 5 to 30 residues. In this study, we used the default window length of 16 to obtain the maximum accuracy of prediction [14]. Predicted epitopes were highlighted as sphere in SARS-CoV-2 spike protein structure viewed by the pymol molecular graphics system [16]. Surface accessibility of predicted peptides was evaluated with the recently resolved protein structure [17]. We utilized vaxijen2.0 server to analyze the antigenicity of chosen epitopes [18].

Discontinuous B-cell epitope prediction
Prediction of discontinuous epitopes on spike protein (PDB ID: 6VSB chain B) was conducted via DiscoTope 2.0 server [19]. The parameter was set at − 1.0 which indicates 85% specificity and 30% sensitivity. This method is based on surface accessibility, residue statistics, and spatial information in a compiled data set of discontinuous epitopes discovered by X-ray crystallography of antigen/antibody complex structure. The contact number, propensity score, and disctope score for each amino acid are provided for conformation-based epitope prediction. Pymol was employed to illustrate the position of predicted epitopes on the 3D structure of SARS-CoV-2 spike protein recently resolved [17].

T-cell epitope prediction
We used the free online service provided by IEDB to forecast T-cell epitopes within nucleocapsid protein binding to HLA-1 [20] or HLA-2 [21] molecule. A relatively small pool of HLA alleles covering the majority of the population, over 97 and 99% for class I and class II respectively, were chosen in the analysis [22,23]. The sequences were given in plain format and the top 50% scoring peptides were retained for further analysis.

Protein coding features of SARS-CoV-2 genome
A map of the predicted open reading frames (ORFs) is depicted in Supplementary Figure 1 based on the genome sequence of the virus Wuhan-Hu-1 isolate (NCBI reference sequence number: NC_045512.2). The genomic structure of SARS-CoV-2 shares characteristics that are also found in other coronaviruses including SARS-CoV, MERS-CoV, and HCoV-NL63. All these coronaviruses contain recognizable ORFs including the replicase (ORF1ab polyprotein), surface glycoprotein (spike protein), envelope protein, membrane glycoprotein, nucleocapsid protein, and several non-structure proteins (NSP). The conserved domains of proteins encoded by the SARS-CoV-2 genome are summarized in Supplementary Table 1. Spike protein mediates the specific binding of the virion to the receptor on the host cell membrane. The overall structure of spike protein is outside the virus particle [17]. Thus, it is an ideal target for B-cell epitope screening. Compared to spike protein, nucleocapsid protein is more conserved in selected coronaviruses (Fig. 2). Though unable to induce humoral immunity, nucleocapsid protein in SARS-CoV and MERS-CoV has been experimentally tested as a robust immunogen to induce cytotoxic T-lymphocyte (CTL)mediated response [26,27], which suggests nucleocapsid protein in SARS-CoV-2 could be a good candidate for T-cell epitope prediction.

Sequence analysis of spike protein and nucleocapsid protein in selected coronaviruses
To better understand the characteristic of SARS-CoV-2, we compared its protein sequences with other selected coronaviruses. All protein sequences were downloaded from NCBI database with accession IDs shown in Fig. 2a. The total sequence identity and phylogenetic tree results were presented in Fig. 2b, c. Consistent with a recently published study [28], we found that both spike and nucleocapside proteins in SARS-CoV-2 are more closely related to that of SARS-CoV. The protein domains of spike and nucleocapsid proteins (Fig. 2d) were depicted based on previous studies on SARS-CoV [29,30] and the protein alignment result in the current study (supplementary files S1 and S2). The amino acid sequence identity result confirmed a high similarity between SARS-CoV-2 and SARS-CoV. As anticipated, the nucleocapsid protein is more conserved among selected coronaviruses compared to spike protein.

B-cell epitopes recognition
The full-length sequence of spike protein was scanned for putative sequential B-cell epitopes by two types of bioinformatics programs. A total of 28 non-overlapping peptides were identified by ABCpred server with the threshold set at 0.85 (Supplementary Table 2). For sequential B-cell epitopes prediction on BepiPred-2.0 server, a threshold value of 0.5 was applied and 35 peptides were predicted (Supplementary Table 3). Antigenicity was calculated by Vaxijen 2.0 server and peptides with the highest antigenicity scores were selected (Tables 1 and 2). The structure of SARS-CoV-2 spike protein was resolved recently with Cryo-electron microscopy (cryo-EM) [17], which could greatly facilitate the process of vaccine development. Predicted epitopes in Tables 1 and 2 were highlighted as sphere in monomer structure of spike protein viewed with pymol (Supplementary Figure 2). While most epitopes predicted were exposed on the surface of spike monomer, only epitopes Spike 315-324 (TSNFRVQP TE), Spike 333-338 (TNLCPF), Spike 648-663 (GCLIGAEHVN NSYECD), Spike 1064-1079 (HVTYVPAQEKNFTTAP) displayed good surface accessibility in spike trimer ( Fig. 3 and Supplementary File: B-cell-epitope-animation.ppt), the pattern more likely exists in nature. Conformationbased B-cell epitopes were computed on DiscoTope 2.0 server [19]. A threshold value of − 1.0 was chosen for the computation, which corresponds to a specificity of 85% and a sensitivity of 30%. The contact number, propensity score, and disctope score for each amino acid that passed the threshold were presented in Table 3. The position of these residues was viewed with pymol and highlighted as sphere (Fig. 4). Processing with a combination of B-cell epitope scanning and peptide analysis forecasted 4 potent linear epitopes and 10 residues involved in discontinuous epitopes formation.

T-cell epitopes recognition
In our study, the IEDB server was utilized following prediction methods recommended (a combination of ANN, SMM, CombLib, and NetMHCpan EL methods for HLA-1 binding prediction, and a combination of NNalign, SMM-align, CombLib, Sturniolo, and NetMHCIIpan methods for HLA-2 binding prediction). For HLA-1 binding peptide prediction, the top 50% scoring peptides were retained for further analysis. A total of 81 nonrepetitive peptides with ANN_IC50 value not higher than 500, indicative of stronger than medium binding affinity, were identified (Supplementary Table 4). Six peptides with the highest antigenicity scores by vaxijen 2.0 were chosen for next step processing. In this step, we screened all HLA-1 molecules being able to bind these peptides (Table 4). A similar strategy was applied for HLA-2 binding peptides prediction on the IEDB server and 64 peptides were identified as HLA-2 binding sequences (Supplementary Table 5). Six peptides with the highest antigenicity scores were selected for HLA-2 molecule screening and the result was presented in Table 5. In the selected peptides pool for HLA binding, Nucleocapsid 104-112 (LSPRWYFYY) was predicted as both HLA class-I and class II binding peptides.  Tables 4 and 5. A partially overlapping region was found in the CTL epitope Nucleocapsid 66-74 (FPRG QGVPI) and the helper T-lymphocyte (Th) epitope Nucleocapsid 67-75 (PRGQGVPIN), which suggests the sequence containing these two epitopes may initiate both CD4 + and CD8 + dependent immune response.

Selected T-cell epitopes feature profiling and evaluation
Peptide stability, mutation analysis, toxicity, allergenicity, hydro and physiochemical features were calculated and the results were presented in Supplementary Table 6. While no peptide listed is toxic, a majority of them are potentially allergenic. To forecast the probability of an immune response induced by the predicted HLA-1 binding peptides, the class I immunogenicity test was performed and the scores were presented in Table 6. A higher score indicates a higher potential of immune response induction.

Peptides selected and multi-epitope based vaccine design
To induce humoral and cellular immune response simultaneously, five peptides that contain four linear B-cell epitopes and three T-cell epitopes (Fig. 5a) were selected for vaccine development. To facilitate the processing of the T-cell epitopes, selected peptides from nucleocapsid protein were extended five amino acid residues at both N-and C-terminus as compared to the predicted epitopes. These peptides were scanned on the IEDB database. We found that five sequences presented in the vaccine construct were identical to the experimentally verified epitopes on SARS-CoV. These determined peptides displayed a strong or medium binding affinity to a series of MHC molecules (Supplementary Table 7). These epitopes likely possess cross-protective effects against SARS-CoV-2 as well. As shown in Fig. 5b, peptides selected in this study were joined together by using GPGPG and (GGGGS) 2 linkers. The Pan DR epitope (PADRE), a universal Th epitope that activates CD4 + cells [31], was introduced at the N terminus of the vaccine construct to enhance helper T cell activity. The vaccine construct can be generated as previously reported [32,33].

Discussion
Different from the less harmful human coronaviruses continuously circulating in the human population,  Compared to traditional vaccine development, potent epitopes can be predicted via bioinformatics analysis, which makes the vaccine design straightforward and fast [34]. As the majority of spike protein is exposed outside the virion, it could be an ideal target to search for B-cell epitopes. Spike proteins in MERS-CoV and SARS-CoV have been shown to induce robust immune response [35,36]. Baruah and Bose's study also indicated that epitopes on spike protein could be promising candidates for SARS-CoV-2 vaccine development [37]. In another study, the experimentally determined epitopes derived from SARS-CoV were screened and those harboring the same sequence in SARS-CoV-2 were identified [38]. Similar strategy was followed in Grifoni's study [39]. Though the SARS-CoV-2 spike protein displayed 77.38% sequence identity toward that of SARS-CoV, most of the antibodies against SARS-CoV spike protein showed poor cross-reactivity with that of SARS-CoV-2 [40], which indicates the spike proteins have significantly variable structure. Thus, we applied distinct strategy and evaluated the surface accessibility of predicted B-cell epitopes to exclude those buried inside the protein. Though we utilized similar bioinformatics tools and resources, the peptides we selected for vaccine development are not alike. In our analysis of SARS-CoV-2 spike protein, four peptides identified from multi-step screening displayed excellent surface accessibility. The peptide Spike 333-338 produces a random coil structure and has a high antigenicity score. Noticeably, this peptide sits in the receptorbinding domain (RBD) of spike protein of SARS-CoV-2, which has been proved to mediate the binding to ACE2 on the epithelial membrane of human lungs [41]. Likely, antibodies recognizing this epitope could also neutralize the virus and prevent infection. The epitope Spike 648-663 predicted in this study is 12 amino acids upstream of the furin cleavage site which is critical for SARS-CoV-2 biogenesis [42]. We expect that an antibody recognizing the epitope Spike 648-663 could have expanded therapeutic function.
Besides B-cell epitope prediction on spike protein, we selected nucleocapsid protein as a target protein for Tcell epitope computation for the following reasons. First,   nucleocapsid protein in SARS-CoV and MERS-CoV has been experimentally tested as a robust immunogen to induce cytotoxic T lymphocyte response [26,27]; Second, nucleocapsid protein is the predominant protein expressed in the virion during the early stage of infection [43]; Third, nucleocapsid protein was detected in the majority of SARS-CoV infected patients as early as day 1 of infection [44,45], which suggests nucleocapsid protein-based vaccine may evoke T-cell dependent immune response timely. Eighty-one and sixty-four peptides within nucleocapsid protein were predicted to bind HLA-1 and HLA-2 molecules respectively. The peptide Nucleocapsid 104-112 showed an immunogenicity score of 0.36, suggesting a capability to initiate a strong immune response. Moreover, the peptides Nucleocapsid 104-112 and Nucleocapsid 66-75 harbor experimentally determined epitopes that bind to a number of HLA-1 molecules. In this study, Nucleocapsid 104-112 and Nucleocapsid 66-75 are forecasted as multi-epitopes that bind to both HLA-1 and HLA-2 molecules, suggesting a potency of both CTL and Th-mediated immune response initiation. Though these peptides were predicted to be non-toxic, we noticed that a large number of these peptides could be allergens. Thus, special attention should be paid to potential allergic reactions during the pre-clinical and clinical trials. As nucleocapsid protein is highly conserved between SARS-CoV-2 and SARS-CoV, available information on nucleocapsid protein-based vaccine against SARS-CoV could be helpful. The final vaccine construct comprising of CTL, Th, and B-cell epitopes is predicted to initiate protective humoral and cellular immune response against SARS-CoV-2.

Conclusions
B-cell epitopes on spike protein and T-cell epitopes in the nucleocapsid protein were predicted and analyzed in the current study. A total of 63 linear B-cell epitopes on spike protein were forecasted by ABCpred and BepiPred servers. Ten residues within spike protein (Gly 496 , Glu 498 , Pro 499 , Thr 500 , Leu 1141 , Gln 1142 , Pro 1143 , Glu 1144 , Leu 1145 , Asp 1146 ) were forecasted by DiscoTope 2.0 program as components of conformational B-cell epitopes. IEDB server was used for T-cell epitopes prediction, which gave rise to 81 and 64 peptides with binding capability to class-I and class-II molecule respectively. Four B-cell epitopes: Spike 315-324 (TSNFRVQPTE), Spike 648-663 (GCLIGAEHVNNSYECD), Spike 1064-1079 (HVTYVP AQEKNFTTAP), and Spike 333-338 (TNLCPF) were selected from the list based on their antigenicity score and surface accessibility. The T-cell epitopes Nucleocapsid 104-112 (LSPRWYFYY) and Nucleocapsid 66-75 (FPRG QGVPIN) on nucleocapsid protein could bind a wide spectrum of both HLA-1 and HLA-2 molecules. The final vaccine construct consists of CTL, Th, and B-cell epitopes that potentially protect individuals against SARS-CoV-2 by inducing both humoral and cellular immune response, which should be successively validated within both in vitro and in vivo models.