Evolutionary Analysis and Prediction of Peptide Vaccine Candidates for Foot-and-Mouth-Disease Virus Types A and O in Bangladesh.

Foot-and-mouth disease (FMD), an endemic disease of cloven-hoofed animals, causes an annual economic loss of US$60-150 million in Bangladesh. There is no cross-protection among the foot-and-mouth disease virus (FMDV) serotypes and vaccination escape mutation may happen. Peptide vaccine is a safer alternative. The aim of this study is to predict and map the B and T cell epitopes of VP1 proteins of FMDV serotypes O and A that were circulating in Bangladesh from 2011 to 2013. Using evolutionary and computational approach (BCPred, BepiPred, DiscoTope, ElliPro, and ProPred-I, IEDB analysis for MHC-I prediction), a total of 11 B and T cell epitopes were predicted. Also, the three-dimensional (3D) structure of VP1 protein showed that the predicted five epitopes residing on N- and C-termini can be considered as good vaccine candidates, and epitopes on the G-H loop can serve as receptor recognition sites for vaccine design. The scores of predicted epitopes of one method were cross-checked with other one for potential epitope mining. Within the VP1 antigenic sites, significant evidence of positive selection was present indicating evolution of VP1 under high immune surveillance.

Introduction

Foot-and-mouth disease (FMD) is a highly contagious and devastating disease affecting cattle and other bi-ungulate species that causes significant financial losses.1 Foot-and-mouth disease virus (FMDV), the causative agent of FMD, belongs to the genus Aphthovirus of the family Picornaviridae. FMDV is classified into seven serotypes (A, O, C, Asia1, SAT1, SAT2, and SAT3) having several topotypes with a limited cross-neutralization. The virus possesses high mutation rate and antigenic variability.

FMDV is a small non-enveloped virus containing a pseudo T3 icosahedral capsid. The FMD virion is composed of a positive RNA genome of 8000 bases having a large single open reading frame (ORF) flanked by highly structured 5′ and 3′ untranslated regions (5′ UTR and 3′ UTR, respectively) and a protein capsid, which is assembled from 60 units of four structural proteins (VP1, VP2, VP3, and VP4) encoded by N-terminal half of the ORF.2,3 The predominantly exposed surface proteins are VP1, VP2, and VP3, whereas VP4 appears to be a completely internal structural protein.4

VP1 is the extensively studied protein owing to its significant roles in virus attachment, protective immunity, and serotype specificity. Generally, the most significant antigenic site of FMDV has centered on the sequence between amino acids 140 and 160 and that of the C terminus of VP1,5,6 therefore, the analysis of the entire or partial coding sequence of the VP1 capsid protein of FMDV merits attention.

Traditional FMDV vaccines are formulated with inactivated virus, which requires the production of tremendous amounts of the infectious agent. However, there are several disadvantages in the use of traditional FMD vaccines: (i) requirement of a cold chain to preserve vaccine stability; (ii) risk of contamination with incompletely inactivated viruses; and (iii) lack of a defined chemical content with regard to both the composition of the viral antigen and the presence of cellular contaminants causing occasional cases of anaphylactic shock in target animals.7 In addition, the presence of non-structural viral proteins in vaccine preparations often makes the distinction between infected and vaccinated animals difficult.8,9

Peptide vaccines containing discrete selected epitopes could be an alternative that avoids the manipulation of the infectious virus. Besides, peptide vaccines is safe in compare to whole viral inactivation vaccine and DNA vaccine where incomplete inactivation or reversion9 of the virus and integration of DNA10 into host cell genome might be occurred. These synthetic immunogens associated with different adjuvants have improved the duration of the response, as well as the delivery and presentation of the antigen, which allows the modulation of the immune pathways leading to protective responses.11 Peptide vaccines may also result in poor immunogens,12,13 and therefore it is essential to detect and evaluate the importance of the sequences of amino acid residues recognized by T and B lymphocytes in order to induce an effective immunity.14,15

Prediction of discontinuous B cell epitopes requires knowledge of the three-dimensional (3D) structure of antigen. In the induction of humoral or cell-mediated immunity using synthetic peptides, understanding of the nature of T and B cell epitopes can play an important role.16,17 The aim of the present study is to screen and identify the B cell and T cell epitopes on structural proteins of FMDV serotypes O and A using different bioinformatics and computational tools in order to gain further understanding of the antigenic structure of the proteins for application in the design of peptide vaccines.

In our recent investigation, we have predicted the antigenic sites for FMDV Asia1 serotype.18 In the present study, eight type A and eight type O VP1 nucleotide and protein sequences of circulatory FMDV in Bangladesh were used to predict the sequence variability, potential antigenic sites, and for selection analysis of types O and A.

Materials and Methods

Amino acid sequence alignment and variability index

A number of sequences (Supplementary Table 1) of serotypes O and A FMDV circulating in Bangladesh in 2012–13 (KC795948–KJ175182) and reference sequences (NC011450 and NC004004) from NCBI database were used to predict the amino acid sequence variability of VP1 region. EBI, EMBOSS Bioinformatics Tools (Transeq) were used to translate the nucleic acid sequences into respective amino acid sequences. VP1 amino acid sequences and reference sequences were aligned by ClustalX (1.81) and analyzed using SeqMan II (Lasergene 8.0; DNAStar Inc., WI, USA) for discerning amino acids sequence variability. Protein variability server (PVS) was used to calculate protein variability index using Wu–Kabat variability coefficient.19 The variability coefficient is computed using the following formula: variability = N*k/n, where, N is the number of sequences in the alignment, k is the number of different amino acids at a given position, and n is the times that the most common amino acid at that position is present.20

Selection pressure analysis of FMDV VP1 protein

For type O dataset, VP1 region nucleotide sequences of 11 local isolates of FMD virus serotype O along with 17 previously reported type O VP1 region sequences from 2009 in Bangladesh were used in this study (Table 1). Besides this, 29 other sequences of ME-SA topotype were selected, which represent related lineages of different outbreaks. For type A dataset, 13 VP1 region sequences of local isolates from 2012 outbreak along with 27 Indian isolates of recent outbreaks were collected. The selective pressure on FMDV sequences was inferred using different analyses (fixed effect likelihood [FEL], internal fixed effects likelihood [IFEL], mixed effect model evolution [MEME], branch site random effect likelihood [REL]) in Datamonkey webserver.21 Tamura–Nei model was selected for all of the analyses. Initially, FEL22 and IFEL were done to find positively selected sites (dN > dS) with a P-value cut-off of 0.3.23 But in many cases evolution can be diversifying and episodic in some codon positions in some branches, but not in the rest of the codon positions in other branches; so overall methods like FEL and IFEL cannot detect them. Therefore, MEME, a method to look for evidence of both diversifying and episodic selection at individual sites, was used with a significance level of 0.1.24 To detect diversifying selection, branch site REL was implemented using a p-value cut-off of 0.05.

Table 1

Results of positive selection analysis on FMDV VP1 region.

ANALYSIS METHOD	SEROTYPE	POSITIVE SELECTION RESULTS
		DATA-SET/SEQUENCE	CODON POSITION
FEL	Type O	Each data-set	45, 46, 56, 110
	Type A		4, 24, 43, 46, 48, 58, 96, 108, 134, 143, 168, 171, 194
IFEL	Type O		45, 46
	Type A		4, 48, 134, 143, 168
MEME	Type O	BAN_LA_Du-135_2013	14
		BAN_FA_Do-12_2012	17, 18, 19, 45
		BAN_FA_Do-11_2012	45
		BAN_FA_Do-12_2012
		BAN_PA_Kg-20_2012
		BAN_LA_Sa-137_2013	198
		BAN_TA_Dh-185_2013	19

Modeling and validation of VP1 protein 3D structure

For 3D structure analysis, SWISS-3D MODEL was used18,25 applying homology modeling method. PDB id: 1fod1 for type O and 4iv1A for type A were selected as templates for maximum sequence identity and E value. To evaluate the quality of 3D structures, Ramachandran Plot 2.0 to check stereochemical properties of the structure26 and PROSA for energy criteria27 were used. To refine target model, molecular dynamics simulation was done with YASARA force field minimization server.28

Epitope prediction and mapping of VP1 protein of FMDV

To predict continuous B cell epitopes on VP1 of FMDV, the amino acid sequences were analyzed using the DNAStar Protean system. The secondary structure was predicted using Garnier–Robson29 and Chou–Fasman methods.30 Surface properties of the VP1 protein, such as hydrophilicity, flexibility, accessibility, and antigenicity, were analyzed by Kyte–Doolittle,31 Karplus–Schulz,32 Emini,33 and Jameson–Wolf34 methods, respectively. Based on the results of these methods, the peptides with good hydrophilicity, high accessibility, and flexibility and strong antigenicity were selected.

The peptides located in α-spiral and β-sheet regions, which do not readily form epitope regions, were excluded. Besides, BCPred35 and BepiPred36 were used to evaluate linear B cell epitopes. Discontinuous or conformational epitopes were predicted by DiscoTope and ElliPro server with default settings.37

T cell epitopes were predicted using ProPred-I and IEDB analysis for MHC-I prediction. At present, about 60 full-length validated cattle MHC class-I cDNA sequences are available (http://www.ebi.ac.uk/ipd/mhc/bola), and based on that data, alleles were selected for MHC-I prediction. ProPred-I is a web-based tool that predicts binding peptides for MHC class-I alleles. IEDB analysis resource using Average Relative Binding (ARB) was utilized to predict IC50 values for peptides binding to specific MHC molecules.38 The peptides that have IC50 value less than 500 nm are considered as binders and those with IC50 value greater than 500 nm are considered as non-binders.17

Results

VP1 amino acid sequence variability of FMDV serotypes O and A in Bangladesh

The deduced amino acid sequences of both types O and A (Supplementary Table 1) were aligned using ClustalX software (Supplementary Fig. 1). In both serotypes, high variability was observed in the hypervariable region of the G–H loop between amino acid positions 130–160 (for type O) and 121–151 (type A), which contributes to major antigenic sites on the capsid coding region.39,40 In contrast, the RGDLXXL motif was found to be conserved in both serotype A and serotype O viruses and contained the RGDLGQL and RGDLQVL motif, respectively. Protein variability was calculated using Wu–Kabat variability coefficient in PVS (Fig. 1). Variability was mainly observed in the G–H loop exceeding a threshold value of 1.000 (Supplementary Table 2). Besides, a consensus sequence for both types was obtained from PVS analysis.

Figure 1

Wu–Kabat protein variability plot for VP1 type O (A) and type A (B) protein sequences. Variability is mostly observed in 130–160 position for type O and 121–151 position for type A.

Evidence of selection analysis on FMDV VP1 region

For type O dataset, FEL revealed four codon positions (Table 1) under positive selection. IFEL also corroborated evidence of positive selection on codons 45 and 46. For type A dataset, FEL revealed 13 sites under positive selection throughout whole VP1 and five of them were corroborated by IFEL. MEME revealed seven individual isolates of type O undergoing diversifying and episodic positive selection in different codon positions (Table 1). Isolates of 2012 outbreak in Bangladesh were consistently found with a 45Q* and 46N* substitution, which can be the consequence of positive selection under immunological surveillance. MEME revealed seven codons under diversifying and episodic positive selection for all type A local isolates; among them, codon positions 143V* and 168V* were positively selected and constantly deviated among all local isolates. All local isolates constantly deviate from other sequences in the dataset by 143V and 168R. As residue 143 is situated in the G–H loop, this change should have a major biochemical change in antigenic properties.

Construction and validation of FMDV VP1 protein 3D structure

To predict probable conformational epitopes and find distinct antigenic determining sites, knowledge of protein 3D structure is prerequisite. 3D structures for both types were constructed using SWISS MODEL. Templates (PDB id: 1fod1 for type O and PDB id: 4iv1A for type A) were taken considering better identity (88.72% – 1fod1 and 84.49% – 4ivA) and E-value (4.25e−97 for 1fod1 and 6.53e−84 for 4iv1A). The Qualitative Model Energy Analysis (QMEAN) value for type O 3D structure was 0.64 and was 0.58 for type A estimating the reliability of the model.41 Using PDB files of targets, Ramachandran Plots 2.0 for both types were developed (Fig. 2), and it was found that 98% of the residues in both models reside in the favored region and have good overall quality. The PROSA analysis of the model showed maximum residues to have negative interaction energy, and the overall Z-score is −3.71 for type O and −4.35 for type A model (Fig. 3).

Figure 2

Estimating 3D model quality by Ramachandran plot where most residues are in the favored region. Plot (A) shows type O VP1 model and plot (B) shows type A VP1 model.

Figure 3

PROSA Z-Scores for VP1 model of type O (−3.71) and A (−4.35), which reside in a zone of experimentally derived overall Z-Score. Plot (A) shows type O-VP1 model and plot (B) shows type A-VP1 model.

Using PyMol visualization tool, target models were superimposed on template models. Several loops (N-termini, BC, EF, FG, GH, C-termini) were labeled (Fig. 4). RMSD (Root-Mean-Square Deviation) for the whole model and loops were calculated using SuperPose version 1.018 to find the deviation among templates and target models (Table 2). Molecular dynamics simulations were done with YASARA force field minimization server for refining target models.28 For type O, the force field energy before minimization was −64360.1 kJ/mol and after minimization was −95945.5 kJ/mol. For type A, energy was −65015.2 kJ/mol, which became −89721.2 kJ/mol after minimization.

Figure 4

Protein 3D model of VP1 protein of FMDV serotype; type O (A) and type A (B) show different loops with presumptive positions.

Table 2

RMSD values for VP1 protein and different loops in type O and type A 3D model.

FMDV-TYPE O					FMDV- TYPE A
REGION	POSITION	RMSD VALUE (A°)			REGION	POSITION	RMSD VALUE (A°)
		ALPHA CARBON	BACKBONE	ALL			ALPHA CARBON	BACKBONE	ALL
VP1	All	0.10 (195 atoms)	0.12 (780 atoms)	0.27 (1511 atoms)	VP1	All	0.08 (168 atoms)	0.10 (672 atoms)	0.45 (1304 atoms)
N-termini	1–30	0.090 (30 atoms)	0.12 (120 atoms)	0.47 (248 atoms)	N-termini	1–24	0.0 (24 atoms)	0.11 (96 atoms)	0.45 (191 atoms)
BC Loop	45–55	0.09 (11 atoms)	0.11 (44 atoms)	0.13 (74 atoms)	BC Loop	43–45	0.05 (3 atoms)	0.05 (12 atoms)	0.10 (24 atoms)
EF loop	84–102	0.1 (19 atoms)	0.11 (76 atoms)	0.15 (147 atoms)	EF loop	70–76	0.07 (7 atoms)	0.08 (28 atoms)	0.10 (45 atoms)
FG loop	104–113	0.09 (10 atoms)	0.1 (40 atoms)	0.15 (79 atoms)	FG loop	109–116	0.01 (2 atoms)	0.181 (4 atoms)	0.846 (9 atoms)
GH loop	130–160	0.11 (31 atoms)	0.1 (124 atoms)	0.25 (237 atoms)	GH loop	121–151	0.066 (18 atoms)	0.088 (74 atoms)	0.081 (124 atoms)
C-termini	190–203	0.11 (6 atoms)	0.15 (24 atoms)	0.23 (43 atoms)	C-termini	185–189	0.039 (3 atoms)	0.061 (12 atoms)	0.07 (19 atoms)

Prediction and mapping of epitopes on VP1 protein

Variability, fragment mobility, and hydrophilicity are important features of antigenic epitopes. The existence of flexible regions like coil and turn regions provides powerful evidence for linear B cell epitope identification.42 In the present investigation, the secondary structure of the VP1 polyprotein was predicted using the methods of Garnier–Robson and Chou–Fasman. Prediction of continuous B cell epitopes was obtained by the methods of Karplus–Schulz (flexibility plot), Kyte–Doolittle (hydrophilicity plot), Emini (accessibility), and Jameson–Wolf (antigenicity) using DNAstar Protean system (Fig. 5). Based on the results obtained with these sequence-based methods, the potential seven B cell epitopes for type O and four for type A were predicted (Table 3). The results showed that the secondary structure of the VP1 polypeptide of FMDV consists of numerous regions of β sheet and turns.

Figure 5

Secondary structures, flexibility plot, hydrophilicity plot, surface probability plot and antigenicity index for VP1 polyprotein of FMDV type O (A) and type A (B).

Table 3

Continuous B-cell epitopes predicted by DNAStar Protean system.

FMDV TYPE		PREDICTION ITEMS
		HYDROPHILICITY	ACCESSIBILITY	FLEXIBILITY	ANTIGENICITY	GENERAL ANALYSIS RESULTS
Type O	Deduced Peptide positions	17–27, 35–39, 44–58, 60–64, 74–77, 98–104, 134–138, 158–161, 174–180, 191–206	1–15, 84–96, 182–188	5–16, 27–33, 67–70, 74–79, 83–89, 116–131, 142–159, 166–175, 181–189, 206–215	14–30, 32–39, 43–66, 91–114, 131–139, 141–147, 158–164, 169–185, 187–208	1–14, 28–31, 72–94, 118–128, 141–147, 170–172, 182–186
Type A		6–16, 78–81, 98–109, 111–123, 163–166, 174–177	77–88, 108–114, 173–180	5–8, 60–62, 67–72, 76–83, 85–87, 108–125, 172–188	6–12, 14–31, 35–60, 89–107, 123–138, 149–156, 160–175, 180–200	6–10, 76–89, 108–138, 172–180

Moreover, BCPREDS, BepiPred, ElliPro, and DiscoTope 2.0 were implemented to predict B cell epitopes. Both BCPREDS and BepiPred use sequence data to find continuous epitopes whereas DiscoTope and ElliPro use 3D data to find discontinuous epitopes. 24 linear epitopes were predicted by BCPREDS and BepiPred.43 But only six epitopes were selected as candidates that fully or mostly overlapped with ElliPro and DiscoTope prediction and were located in mostly conserved region (Table 4). For type O VP1, four epitopes were found to be promising with location in N-termini, E–F, G–H loop, and C-termini. For type A VP1, only two epitopes located in the E–F loop and the FG–GH loop sharing were promising (Fig. 6). To predict FMDV-specific CD8+ T cell responses in cattle, ProPred-I and IEDB T cell epitope prediction tools were used. Three epitopes for type O were chosen, which were located in N-termini, G–H loop, and C-termini. Two epitopes of type A were located in N-termini and the G–H loop (Table 5). These epitopes are considered as effective vaccine candidates because it was found that peptide sequences containing the G–H loop and C-termini of capsid protein VP1 in which a T cell epitope has been added were shown to provide immunity in model animals.12

Table 4

Prediction of VP1 B-cell epitopes with location and different programs scores.

	POSITION	LOCATION	SEQUENCE	BCPREDS SCORES	BepiPred SCORES (AVG.)	ElliPro SCORES	DiscoTope SCORES (AVG.)
FMDV type O	1–14	N-termini	ENYGGETQVQRRQH	0.738	0.509	0.851	1.788
	73–86	EF loop	WVPNGAPEAALDNTT	0.993	1.183	0.667	−7.74
	119–142	GH loop	CKYGEGAVTNVRGDLQVLAQKA	0.735	0.609	0.751	−3.65
	194–213	C-termini	VKQLLNFDLLKLAGDVESNPGO	0.779	1.151	0.858	−0.41
FMDV type A	71–96	EF loop	PEAALSNTGNPTAYNKAPFT	0.964	1.257	0.737	−9.07
	113–12 2	FG & GH loop	YSAASGRVRG	0.848	0.792	0.846	−3.70

Figure 6

Predicted B and T cell epitopes mapped onto the protein 3D structure using the PyMOL visualization tool. Green in (A) indicates VP1 protein structure, whereas ruby pink shows B-cell-specific epitope and wheat color shows T cell (CD8+) epitope. In (B), cyan shows VP1 protein of type A, lime-green shows B cell epitope and magenta shows T cell epitopes.

Table 5

Prediction of T-cell epitopes with location and different programs scores.

	POSITION	LOCATION	SEQUENCE	ProPred1 SCORES	IEDB (MHC-I(IC-50 UNIT))
FMDV type O	25–36	N-termini	VTPQNQINVL	5.70 (A20)	0.2 (Bola-N:00101)
	153–160	GH-loop	IKATRVTEL	6.21 (A20)	0.8 (Bola-N:01602)
	191–200	C-termini	VAPVKQLLNF	5.70 (A20)	1.1 (Bola-N:00101)
FMDV type A	29–37	N-termini	VKIGNVSPT	4.60 (A20)	0.3 (Bola-N:02301)
	123–135	GH-loop	DLGQLAARVAAQL	8.29 (A20)	0.6 (Bola-N:01001)

Discussion

FMDV infection of cloven-hoofed animals occurs throughout the world and causes huge economic losses. Extensive studies are focused on replacing conventional FMDV vaccines with peptide-based vaccines to avoid manipulation with viruses. Again, immunity to one serotype is not sufficient to provide immunity to other serotypes. Therefore, understanding of the amino acid sequences that can be recognized by T- and B-lymphocytes has attracted attention. This report conclusively presented preliminary information on the potential B cell and T cell epitopes concentrated on the G–H loop and the C-terminal region of VP1 proteins of FMDV types O and A isolated from Bangladesh in recent years.

In this investigation, multiple sequence alignment with NCBI Ref/seq sequences showed that 130–160 (for type O) and 121–151 (type A) positions on the G–H loop of VP1 between amino acid sequences contain significant variability. Protein variability index determined using the Wu–Kabat method also showed that the variability in those regions exceeded a threshold value of 1.00. Moreover, the conserved RGDXXL motif was also identified on the G–H loop from alignment and variability index. Consensus sequences from PVS for both types O and A were used for epitope prediction.

Sequence- and structure-based bioinformatics approaches were applied to predict continuous and discontinuous B cell and T cell epitopes. This report showed that epitopes 1–14 and 194–203 on N-termini and C-termini of type O VP1 (Tables 4 and 5), respectively, were previously identified by ELISA and Western blot42,44 and indicate auspicious vaccine candidates. Epitopes on the G–H loop of VP1 including the RGD motif may help in receptor recognition and can be considered as vaccine candidates. Besides this, epitopes on minor antigenic sites E–F and F–G loop have been considered in vaccine design for robust protection against FMDV.

T cell epitopes are important for understanding protective immunity mediated by CD8+ and CD4+ T lymphocytes. MHC-II-specific T cell responses may play an important role in protection against FMD, but no such prediction tool was found to identify MHC-II-specific T cell response rather than immunological assay. Alleles for MHC-I binding (BoLA) T cell epitope were selected from a previous study.45 Epitopes 153–160 of G–H loop, 194–200 of C-termini on VP1 of type O, and 29–37 of B–C loop and 123–135 of GH on VP1 of type A are considered as effective vaccine candidates as they showed protective immunity in a previous study.12 These results indicate that peptides on N-termini, G–H loop, and C-termini probably present the serotype-specific epitopes of FMDV of serotypes O and A.

Moreover significant evidence of positive selection was identified in both serotype O and A sequences on different positions. It has been previously identified that locations of codons under positive selection coincide closely with those of antigenic sites. However, these results suggest that arising antigenic variants benefit from such selective advantage in their interaction with the immune system, either during the course of an infection or in transmission to individuals with previous exposure to antigen. Thus, a comprehensive consideration of bioinformatics methodology and the molecular characteristics of the virus may improve the success rate of experiments aimed at mapping epitopes on VP1 protein.

Conclusion

The findings of this study therefore show that peptide epitopes for both B and T-lymphocytes can be defined precisely and mapped onto VP1 protein of FMDV types O and A collected from Bangladesh. These epitopes could be considered as candidates for designing peptide vaccines, which could provide specific protection against FMDV serotypes O and A found in Bangladesh.

Supplementary Table 1. FMDV serotype A and O VP1 sequences used in this study.

Supplementary Table 2. Protein variability data for FMDV serotype A and O VP1 protein calculated using Wu–Kabat variability coefficient.

Supplementary Figure 1. Alignment of deduced amino sequences for the VP1 proteins of FMDV serotype O (up) and serotype A (below) from Bangladesh local isolates. The RGDLXXL motifs are in the box. Red boxes show GH-loop where the inner small box indicates RGDLXXL motifs.