Immunoinformatics prediction of potential B-cell and T-cell epitopes as effective vaccine candidates for eliciting immunogenic responses against Epstein-Barr virus.

BACKGROUND: The ongoing search for viable treatment options to curtail Epstein Barr Virus (EBV) pathogenicity has necessitated a paradigmatic shift towards the design of peptide-based vaccines. Potential B-cell and T-cell epitopes were predicted for nine antigenic EBV proteins that mediate epithelial cell-attachment and spread, capsid self-assembly, DNA replication and processivity.
METHODS: Predictive algorithms incorporated in the Immune Epitope Database (IEDB) resources were used to determine potential B-cell epitopes based on their physicochemical attributes. These were combined with a string-kernel method and an antigenicity predictive AlgPred tool to enhance accuracy in the end-point selection of highly potential antigenic EBV B-cell epitopes. NetCTL 1.2 algorithms enabled the prediction of probable T-cell epitopes which were structurally modeled and subjected to blind peptide-protein docking with HLA-A*02:01. All-atom molecular dynamics (MD) simulation and Molecular Mechanics Generalized-Born Surface Area methods were used to investigate interaction dynamics and affinities of predicted T-cell peptide-protein complexes.
RESULTS: Computational predictions and sequence overlapping analysis yielded 18 linear (continuous) and discontinuous (conformational) subunit epitopes from the antigenic proteins with characteristic surface accessibility, flexibility and antigenicity, and predictive scores above the threshold value (1) set. A novel site was identified on HLA-A*02:01 with preferential affinity binding for modeled BMRF2, BXLF1 and BGLF4 T-cell epitopes. Interaction dynamics and energies were also computed in addition to crucial residues that mediated complex formation and stability.
CONCLUSION: This study implemented an integrative meta-analytical approach to model highly probable B-cell and T-cell epitopes as potential peptide-vaccine candidates for the treatment of EBV-related diseases.

At a glance of commentary

Scientific background on the subject

The design of peptide-based vaccines could engender the treatment of Epstein Barr Virus (EBV) infection since no cure exists till date. Hence, predicting highly probable B- and T-cell epitopes from a cohort of antigenic EBV proteins could be vital for the design and development of highly effective anti-EBV vaccines.

What this study adds to the field

Our thorough and combinatorial investigative approaches yielded B-cell and T-cell epitopes with favorable physico-chemical attributes that make them potential anti-EBV vaccine candidates. Interestingly, we were able to define a novel site on HLA-A*02:01 that could be further exploited for high-affinity T-cell epitope binding.

Epstein Barr virus (EBV) is a ubiquitous lymphotropic gamma-1 herpesvirus that belongs to the Herpesviridae family. EBV is also known as human herpesvirus 4 (HHV-4) which is widespread in humans infecting over 90% of the adult population globally. More specifically, in developing countries, EBV mediates an asymptomatic infection that usually occurs at the early infancy stage while there is a delay in acute infection until the early adulthood stage in developed countries with about 25% cases resulting into infectious mononucleosis (IM) with clinical features such as fever, swollen lymph nodes, severe fatigue, sore throat and immunologic dysfunction [1].

EBV is the principal cause of IM while it has also been implicated in several human cancers. About 125,000 new cases of IM are reported annually in the United States while about 200,000 new cases of malignancies associated with EBV occur yearly worldwide [2]. EBV has been associated with distinct forms of cancer such as nasopharyngeal carcinoma, stomach cancer, Burkitt's lymphoma, Hodgkin's lymphoma and some entities of B-/T-cell lymphomas [[3], [4], [5], [6], [7]]. Also, the virus has been linked with autoimmune diseases such as multiple sclerosis [8].

EBV has a life-long persistence in host and it is usually transmitted via the oral route. Upon close contact, viral particles are transferred from the saliva of infected individuals to the new host, entering the oral cavity after which the virus is amplified by a lytic (replicative) infection in the permissive cells of the oropharynx. Subsequently, the virus is shed in vast amounts into the throat followed by a simultaneous infection of the B cells, initiating a latent, persistent infection of the pool of B cells [[9], [10], [11]].

Amplification of the viral load follows the viral transmission to a native host and occurs via replicative (lytic) infection in a permissive cell type. This then has a life-long persistence in the host as an asymptomatic latent infection in a secondary cell type which is occasionally reactivated into the lytic cycle to release infectious virions that can be further transferred into a new host [9]. Regardless, the virus persists whether the infection is accompanied by a disease or not.

EBV is classified based on two types; 1 and 2, based on genetic variations in sequences for the EBV nuclear antigens (EBNA)-2 and EBNA-3. Comparatively, type 1 EBV strains exhibit higher prevalence and cell-immortalizing ability than the type 2 forms [12]. Although both types are detected across the world, type 1 is the most dominant in many populations [13].

The life cycle of EBV constitutes to widely distinct stages which include the latency and lytic stages. In the former, only a subset of genes is expressed which include six nuclear proteins, three latent membrane proteins, two EBV-encoded small non-coding RNAs and the BamHI A rightward transcripts (BARTs). The lytic stage typifies the viral replication stage (lytic reactivation) and involves the production of infectious viral particles. This constitutes about 80 proteins which are categorized into immediate-early, early and late proteins [8,[14], [15], [16]]. Lytic proteins are also classified based on their functionalities; capsid, tegument, envelope (glycoproteins - gp) among many others. These antigenic proteins mediate a wide range of function which includes B-cell transformation, viral attachment, membrane fusion, entry, replication and spread [8,[15], [16], [17], [18], [19]]. EBV antigens such as EBNA-1, -2, -3(A-C), -LP-, LMP1, LMP-2A, -2B and non-coding RNAs (EBERs and BARTs) are crucial for the efficient transformation of infected B-cells and subsequent lymphomagenesis [[20], [21], [22], [23]]. These sums the oncogenic mechanisms (Latency III) of EBV in cancer development. Moreover, B- and T-cell epitopic predictions for most of these antigens have been previously carried out [9,14,24,25].

There have been extensive studies on immune protection against EBV invasion, which involves both protective arms of the immune systems; innate and adaptive. A major component of the early innate immune response is the Natural killer (NK) cells, which produce interferon-gamma (IFN-γ) to delay or prevent the transformation of B-cells by EBV [26,27]. B-cells and T-cells, as components of the adaptive system also mediate strong immune responses since they are activated in the presence of a foreign pathogen. The presence of EBV antigens triggers potent B-cell mediated humoral immune responses. The induction of IgM antibodies against virus capsid antigen (VCA) with a subsequent switch to an IgG isotype characterizes acute primary EBV infection. This results in the formation of non-neutralizing IgG anti-VCA antibodies which exist for a lifetime. The major EBV glycoprotein (gp350) also facilitates the induction of neutralizing IgG antibodies, a process that only occurs after the infection has been resolved [28]. In addition to the aforementioned, other antibodies that target antigens that are non-neutralizing also do not appear until the primary infection is resolved. These non-neutralizing antigens include intracellularly located viral proteins [28,29]. Cytotoxic CD8 T cells also mediate strong adaptive immune responses triggered by EBV antigens. This response is directed towards the elimination of cells infected by the virus, which are recognized due to EBV peptide antigens bound by MHC-1 molecules located at the cellular surface; descriptive of the antigenic presentation pathway [30,31]. Previous studies have reported a dramatic pattern in IM patients whereby EBV CD8 T-cell epitopes derived from immediate-early or early antigens are recognized by about 50% CD8 T-cells [9,24]. Similarly, EBV peptide antigens bound to surface-expressing MHC-II molecules are recognized by CD4 T-cells but in a less-dramatic and focused manner after which the complex is committed into various cytokine-producing T-helper cells (Th) phenotypes that regulate immune responses. In most cases, antiviral responses from most CD4 T-cells that are EBV-specific are characterized by the production of IFN-γ and tumor necrosis factor-alpha (TNF-α). In contrast, only a small proportion of these CD4 T-cells produces IL-2 even though it is the typical and most expected Th1-mediated antiviral response [32].

Currently, no medicinal cure for EBV exists neither are there any prophylactic or therapeutic vaccines available. Numerous vaccine-development efforts to curtail EBV infection have been directed towards these antigenic proteins, particularly gp350 since it is the primary target of naturally occurring neutralizing antibodies and mostly abundant gp present on the virus and in infected cells [2]. Several clinical trials have been carried out for vaccine development, including a phase 2 trial for an EBV vaccine from gp350 protein, which reduced the rate of IM but not viral infection [2,[33], [34], [35]]. Suggestively, this limitation could be due to variations in vaccine attributes such as efficacy and duration of protection towards EBV transmission and disease progression. From previous reports, most of the potential vaccines have a short duration efficiency in which they may delay EBV infection but possibly cannot prevent the increases in IM cases over a long period [36].

Dosing is also an important factor to be considered in EBV vaccination. Taking a cue from a closely-related Varicella-zoster virus (VZV), the efficacy of its approved vaccine is based on dosing as studies have reported that about 85% efficacy is recorded with a single dose while a double dose increased the vaccine efficacy to 94% [[37], [38], [39]]. Relative to EBV infection, a logical explanation could be that current vaccines with shorter duration are efficient until they lapse after which a loss of protection is recorded in congruence with increased susceptibility to EBV-disease such as IM at a later age in patients [33,34,40,41]. Presumably, given the short duration of protection for current EBV (gp350) vaccines, the need for increased single or double doses may be considered in future clinical evaluations [34,42].

The ability of EBV to evade immune response is another factor that may account for the reduced efficacy of potential EBV vaccines. While most vaccines are focused on gp350, certain antigens have been identified for their crucial roles in immune evasion. These include BGLF5, BZLF1, BNLF2a, BILF1, EBNA1, EBNA2 among others [10,43]. Possible inclusion of these antigens to consist of a prophylactic vaccine may further improve efficacy if evaluated in future clinical trials [44].

Vaccination has served as a viable therapeutic option for the prevention and treatment of several microbial infections with a common goal of inducing adaptive immune responses by presenting antigenic fragments to the immune system [[45], [46], [47]].

Over the years, the use of peptide-based vaccines has been promising due to their safety and easy preparation [[48], [49], [50], [51]] and are designed based on different B- and T-cell epitopes. This is possible due to the antigenic nature of proteins secreted by the virus [15,52,53], hence the prediction of antibody epitopes using in silico methods have evolved as a crucial step in vaccine design [[53], [54], [55], [56], [57]].

Epitopic prediction has been made possible via the implementation of immunoinformatics methods which contributes strongly towards the development of fast, reliable and effective peptide-based vaccines. These involve the identification of linear (continuous) or non-linear (non-continuous) epitopes for B-cells and T-cells for in silico vaccine design with a desirable reduction in cost and resources compared to immunological experimentations which are tedious and costly.

These predictions are usually centered on the physicochemical properties of these epitopes and their tendencies to provoke desired immunological responses mediated by B- and T-cells [58]. Therefore, immunoinformatics prediction of potential B and T-cell epitopes for viral antigens remains a prominent way-forward for peptide vaccine development [59,60]. Moreover, tools for predicting T-cell epitopes help identify allele-specific peptides, thereby increasing the possibility of selecting and characterizing peptides that are potential vaccine candidates.

These approaches have been successfully implemented for vaccine prediction against multiple HPV, Ebola, HIV, Nipah, and Zika viruses among others [31,53,56,[61], [62], [63], [64], [65], [66], [67]].

Moreover, immunoinformatics methods have also been employed to predict antibody epitopes for some major latent and lytic antigens that play important roles in EBV virulence [14,53].

Complementarily, in this study, we carried out exhaustive B- and T-cell epitopic predictions of antigenic EBV proteins involved in capsid self-assembly, viral DNA synthesis/replication, and oral epithelial cell attachment/cell-to-cell spread as efficient sub-unit vaccine candidates (Table 1) (Fig. 1). This study is highly essential and could serve as an avenue to achieve the design of highly effective and dynamic prophylactic vaccines having a characteristic mix of multiple EBV antigens with diverse functions over the lytic stage of viral infection [68].

Table 1

Showing detailed information on selected EBV proteins investigated in this study.

EBV proteins	Functions	Expression stage	Accession number	Residue count	References
BMRF2	Viral entry, attachment to oral epithelial cells, manipulates cellular actin cell-to-cell spread	Lytic/structural	P03192	357	[[132], [133], [134]]
BDLF2	Interacts with BMRF2, rearranges cellular actin to promote viral cell-to-cell spread in epithelial cells	Lytic/structural	P03225	420	[18,135]
BORF1 (Triplex Protein 1; TRX1)	Capsid assembly (minor capsid protein-binding protein; mCPBP). Required for efficient nuclear transport of Triplex protein 2 (TRX2)	Early lytic/structural	P03187	364	[18]
BcLF1	EBV capsid self-assembly (Major capsid protein; MCP)	Early lytic/structural	P03226	1381	[18,136]
BDLF1 (Triplex Protein 2; TRX2)	minor capsid protein (mCP) required for capsid assembly in the nucleus	Early lytic/structural	P25214	301	[136]
BFRF3	Small capsomere-interacting protein (sCP) involved in capsid assembly	Early lytic/structural	P14348	176	[18]
BXLF1 (EBV Thymidine kinase; TK)	Involved in the activation and re-activation of viral DNA replication	Early lytic	Q3KSQ2	607	[123,137,138]
BGLF4	Serine/threonine protein kinase. Mediates the phosphorylation of host nucleoporins. Enhances nuclear export of EBV lytic proteins, facilitates virion production	Early lytic	P13288	429	[[138], [139], [140], [141], [142], [143], [144]]
BMRF1	DNA polymerase processivity factor. Plays essential roles in lytic DNA activation and replication. Major phosphoryl target for BGLF4. Interacts with other EBV DNA-binding proteins	Early lytic	P03191PDB ID: 2Z0L	404	[83,[145], [146], [147], [148]]

Fig. 1

Schematic representation of the EBV virion showing structural components and their respective locations. Also, possible locations of antigenic proteins investigated in this study are respectively shown to depict their functionalities.

Interestingly, these antigenic EBV proteins have not been previously examined using immunoinformatics method despite their crucial roles in EBV virulence.

Immunogenicity of the potential T-cell epitopes was also predicted, which entailed their binding strength to human MHC-1 molecule; HLA-A*02:01. This allele (HLA-A*02) is usually reported as the most frequent across various populations and locations, hence our focus on the *02 haplotype in this study. Likewise, MHC-I HLA-B*35, MHC-II DRB1*13, and HLA-DRB1*15 molecules have been pin-pointed for their frequencies in human populations, according to previous reports [[69], [70], [71], [72]].

The response of the immune system to a pathogen in different populations is largely influenced by host genetic polymorphisms. This is because, in the human genome, the HLA genes are the most polymorphic. This has posed a serious challenge in vaccine design and population coverage since each allele interacts with a certain group of peptides [[73], [74], [75]]. Most HLA-I alleles have been grouped based on their ability to bind similar peptides, with this classification covering over 80% of the HLA-A and HLA-B alleles. This classification is referred to as supertypes which consist of HLA molecules that bind similar peptide sets [31].

Both classes I and II MHC molecules present T-cell epitopes that are recognized by the CD8 and CD4 subsets of T-cells respectively, indicating there are CD8 and CD4 T-cell epitopes [30]. Upon recognition of their corresponding epitopes, CD8 T-cells are translated into cytotoxic T lymphocytes (CTL) while primed CD4 T-cells are changed into helper (Th) or regulatory (T-reg) T-cells [76]. MHC-I presented peptides generated from the intracellular region are commonly recognized by the CD8 CTLs while those that originate from the extracellular region and presented by the MHC-II molecules are mostly recognized by the CD4 helper T lymphocytes (HTLs) as earlier stated. According to previous reports, epitopes of CD8 cytotoxic CTLs are highly immunogenic and their presence in vaccine formulation has been linked to the induction of a robust immune response. Epitopes of the CD4 T cell subclass also play important augmentative roles when added relative to the duration and strength of the vaccine [[77], [78], [79]].

Immunogenic prediction is important for surface antigenic presentation relative to T-cell induced cytotoxic cell death [30,80]. Interactions between the peptide-MHC complexes were also modeled and investigated using peptide docking and molecular dynamics simulation algorithms to identify potential epitopic candidates for peptide/multi-peptide vaccine development. Current findings will contribute to the development of possible prophylactic vaccine candidates against EBV infections.

Materials and methods

Retrieval of EBV-target sequences

The respective sequences of the selected proteins were obtained from the Universal Protein Knowledgebase (Uniprot) database [81] using corresponding accession numbers presented in Table 1. The FASTA formats were retrieved and used for antigenic epitopic predictions and structural modeling. Structural modeling was performed using the Local Meta-Threading Protein Fold Recognition Algorithm [82], which applies to non-resolved structures while the structurally resolved protein; BMRF1, was retrieved from the Protein Data Bank (PDB) with entry 2Z0L [83]. Validation of the modeled structures was done using PROCHECK [84] and RAMPAGE [85].

B-cell epitopic prediction and validation

Linear (continuous) B-cell epitopes were predicted from the retrieved antigenic protein sequences while the 3D structures of the respective proteins were used to predict non-linear (discontinuous) epitopes.

Integrative methods were used for linear B-cell epitopic predictions based on the physicochemical attributes of the innate antigenic sequences. These attributes include surface exposure (accessibility), flexibility, antigenicity, and hydrophilicity. These were estimated from the linear antigenic sequence using Emini Surface Accessibility [86], Karplus & Schulz flexibility [87], Kolaskar & Tongaonkar Antigenicity [88] and Parker hydrophilicity [89] predictive tools on Immune Epitope Database (IEDB) analysis resources [90]. High thresholds range of 0.9–1.0 were selected for epitopic prediction.

The B-cell epitope prediction method (BCPRED) was also employed for linear epitope prediction and scoring [91]. The epitope lengths were varied from 12 to 22 to obtain epitopes within the high threshold range of 0.9–1.0.

Using the 3D structures of the respective proteins, linear and discontinuous B-cell epitopes were predicted with ElliPro which functions based on structural protrusion measured by residual protrusion indices (PI), an approximation of protein shape together with the final neighboring residual clustering [92].

Multiple prediction methods were employed for cross-validation and end-point selection of the most probable antigenic EBV B-cell epitopes based on their physicochemical attributes.

Allergenicity predictions and mapping of probable IgE epitopes from the designed peptides were performed using support machine vector, motif-based and BLAST-search algorithms integrated with AlgPred web server, with a reported 85% accuracy [93].

Prediction of potential cytotoxic T-cell epitopes and MHC-1 processing

The NetCTL 1.2 method [94] was used to predict probable CTL epitopes from the selected antigenic EBV proteins. Moreover, this method integrates the prediction of the three main steps of the MHC-1 antigenic presentation pathways which include proteasomal processing, transporter-associated with antigen processing (TAP) and MHC-1 binding. While TAP transport efficiency is predicted using the weight matrix, MHC-peptide binding and proteasomal C-terminal cleavage were predicted using an artificial neural network [95,96]. Default values of 0.15 and 0.05 were used for weights on C terminal cleavage and TAP transport efficiency while the threshold for epitope identification was set at 1.00 to enhance predictive accuracy.

Prediction of MHC-1 immunogenicity

This was essential to further identify epitopic from non-epitopic peptides from the predicted sequences. This explains their ability to be recognized by T cells. According to previous studies, amino acid residues with large aromatic side chains have a higher propensity to induce T-cell immunogenicity. Moreover, highly critical positions in a peptide for provoking immunogenic responses are positions P4-6 [97]. Immunogenic determination of the predicted epitopes was carried out using an integrated IEDB method [97], which determines potential immunogenicity of a peptide-MHC complex (pMHC) based on amino acid properties and location in the queried peptide sequence.

Construction of selected peptide library and blind peptide-protein docking

The PEP-FOLD algorithm was employed for the de-novo structural prediction of selected epitopes. This stochastic procedure entails a Hidden Markov Model-derived Structural Alphabet (SA) [98] which relies on the sOPEP coarse-grained forcefield for 3D sampling and generation, [99] and the Monte-Carlo procedure for final refinement. This is the final step wherein clusters are identified and sorted while resulting conformations are scored [100].

The best-ranked peptide models were then docked to HLA-A*02:01 (A2 supertype) which is a highly frequent human MHC-I allele [31] with a well-resolved X-ray crystal structure; 4UQ3 as retrieved from PDB [101]. Since the predicted epitopes could possess varying binding site affinities for the selected MHC-1 molecule (HLA-A*02:01), a blind peptide-protein docking approach was employed using the pep-ATTRACT algorithm (http://bioserv.rpbs.univ-paris-diderot.fr/services/pepATTRACT/) [102,103]. This incorporates a three-step procedure involving ATTRACT, iATTRACT and AMBER modules. The procedure is initiated by ATTRACT which executes rigid-body docking using the coarse-grained ATTRACT forcefield after which resulting models were ranked by their respective ATTRACT scores [104]. Afterward, the best 1000 models were firstly refined using iATTRACT [103], which is a flexible interface refinement method and then secondly by molecular dynamics simulation using AMBER14 [105]. Then, the final models were clustered and ranked by the average energy of top-ranked four models [106,107].

Taken together, this method performs a rigid-body global search of the surfaces of receptor proteins provided as inputs and identifies the most appropriate sites for peptide binding.

Moreover, the docked complexes with the highest global energy scores were analyzed and further simulated in explicit solvent using Molecular Dynamics (MD) simulation methods according to protocols employed in our previous studies [108,109].

MD simulation of peptide-protein complexes and interaction analyses

To further explore the interaction dynamics between the modeled complexes, we employed the GPU-PMEMD engine of AMBER18 MD package together with its interaction modules [110]. The complexes were parameterized using the FF14SB forcefield while the tleap module was used to neutralize the systems by the addition of Na+/Cl− counter ions. More so, the systems were explicitly solvated in a TIP3P water box of sizes 10 Å. Resulting complexes were subjected to an initial 2500 energy minimization steps with a restraint potential of 500 kcal/mol. This was followed by a second full minimization of 5000 steps without restraints. The systems were then heated gradually from 0 → 310K for 50ps to maintain a fixed volume and atomic number. The solutes were restrained with a potential harmonic restraint of 1ps and maintained at a collision frequency of 1ps. Long-range interactions were included using the SHAKE algorithm and Particle-Mesh Ewald (PME) method while a nonbonded interaction cut-off radius of 12 Å was employed. The systems were then equilibrated for 500ps time at constant temperature (310K) and pressure with fixed atoms, simulating the isobaric-isothermal NPT ensemble. The Berendsen barostat was used to maintain the pressure at 1 bar [111]. MD production run of 100ns was then initiated for the respective epitopic-HLA complexes while resulting trajectories were saved at every 2ps. The integrated CPPTRAJ and PTRAJ modules were used to analyze these trajectories at the end of the simulation.

Post–MD analytical methods employed in this study include root mean square deviation (RMSD) and root mean square fluctuation (RMSF). Cα RMSD measures time-based deviations and motions that occur across the protein backbone atoms relative to the starting structures. This provides a clue into the degree of stability exhibited by the complex from the resulting frames, particularly when bound by the respective T-cell epitopes. The mathematical expression for RMSD is given below:

From the equation, N depicts the total number of atoms in the system while Ri represents the vector position of the Ca atom of particle i. Alignment of the initial (o) and final conformations with the least square fitting was used to predict Ri.

Moreover, the equation for calculating RMSF is shown below; and accordingly, to obtain the standardized RMSF (sRMSF), the average RMSF is deducted from RMSFi (RMSF of the ith residue) and further divided by its standard deviation. The RMSF parameter is used to monitor the motion or fluctuations of individual residues that constitute the protein molecule, giving overall insights into the characteristic flexibility or rigidity of such molecule [112,113].

3D visualization was done on the Graphical User Interfaces (GUIs) of UCSF Chimera, (BIOVIA) Discovery studio and Schrodinger Maestro while resulting data-plots were done using the Origin analytical tool [[114], [115], [116], [117]].

MM/GBSA peptide-protein interaction analysis

The interaction-free energies between the predicted T-cell epitopes and HLA-A*02:01 were estimated using the Molecular Mechanics Generalized-Born Surface Area (MM/GBSA) method which also enabled the estimation of other energy components such as van der Waals, electrostatic, polar and non-polar solvation energies [108,[118], [119], [120]]. The equations below mathematically describe total binding energy (ΔG) for complex formation:

From the equation, G, G, and G represent the absolute free energies of the peptides, unbound protein and peptide-protein complex respectively. Decomposition of the total energy (ΔG) as shown in equations (2), (3), (4), (5) results in the solvation energy term (G) which sums up the polar (G) and non-polar (G) solvation terms; the gas-phase (E) energy contribution which also sums up the non-bonded [van der Waals (E) and electrostatics (E)] in addition to the bonded [internal energy (E)] energy terms, and the entropy term (TΔS), which was estimated using the quasi approximation method. The non-polar solvation term was determined by estimating the surface accessible surface area (SASA) using a water probe of 1.4 Å. Total binding energy (ΔG) for the peptide-protein complexes were computed in addition to other energy components. More so, to eliminate entropical effects, we selected frames at which the systems were structurally stable for MM/GBSA energy computations. 3D structures averaged over the 100ns simulation run were also used to investigate binding patterns and dynamics relative to complex formation.

Results

Sequence retrieval and 3D structural modeling

Accession numbers (AN) for the retrieved sequences are shown in Table 1, which also indicated the number of residues in the constituent proteins.

Also, the corresponding secondary structures of the proteins as modeled on LOMET webserver are shown in Fig. 2 including the crystal structure of BMRF1 obtained from PDB. Structural validation was carried out to verify the quality of the protein models used in this study. Ramachandran plots that estimate the percentage of residues in favored, allowed, expected and disallowed/outlier regions are shown in Supplementary Fig. S1.

Fig. 2

3D structures of EBV antigenic proteins. [A] BMRF2 [B] BDLF2 [C] BORF1 [D] BcLF1 [E] BDLF1 [F] BFRF3 [G] BXLF1 [H] BGLF4 [I] BMRF1.

Physicochemical characterization of linear B-cell epitopes

The selected epitopes were further analyzed for antigenic sequence characteristics using amino acid residue scale and Hidden Markov Models (HMMs) which incorporated integrative methods as earlier mentioned. The hydrophilicity of the predicted epitopes was accessed using the Parker hydrophilicity prediction method while graphical illustrations of corresponding residual hydrophilicity are presented in Supplementary Fig. S2. Regions with the highest score could correlate with high hydrophilicity. Threshold values of 1 were set to determine hydrophilic peptide sequences while those with the maximum and minimum scores were identified. Among the predicted epitopes, amino acid positions GEEDDDG54-60 for BXLF1 had the highest score while the minimum calculated score was −7.414 at positions FFIFLAL165-171 for BMRF2. Maximum and minimum hydrophilicity scores for other EBV antigenic proteins are presented in Table 2.

Table 2

Characterization of predicted EBV B-cell epitopes based on their innate physicochemical properties.

Physicochemical characterization
	Hydrophilicity		Surface flexibility		Surface accessibility		Antigenicity		Allergenicity
EBV proteins	Epitope	Score (Max/Min)	Epitope	Score (Max/Min)	Epitope	Score (Max/Min)	Epitope	Score (Max/Min)	AlgPred
BMRF2	CARGDHS197-203	4.571	LGSGSLA171-177	1.071	RYKLKK320-325	10.315	TVLVVCV225-231	1.3	Non-allergenic
	FFIFLAL165-171	−7.414	VVYMCCF336-342	0.876	VLVVCV226-231	0.068	AGARRNQ177-183	0.934
BDLF2	QDQNQNQ95-101	6.857	DQNQNQT96-102	1.133	QDQNQN95-100	4.805	VFIYYCC190-196	1.253	Non-allergenic
	LAFLALL198-204	−5.971	RMAAMAP146-152	0.899	CCYLAF195-200	0.070	NRNRTFN298-304	0.868
BORF1	GGSGGGS73-79	5.929	GGGSGGG72-78	1.142	RRRLQR9-14	6.535	VVHGPVV206-212	1.225	Non-allergenic
	LSLFILM129-135	−6.071	HCQIMNI286-292	0.892	GGVCFV245-250	0.078	GGGSGGG72-78	0.894
BcLF1	NGDHDHD785-791	6.700	QKSNSPR1197-1203	1.115	PRDRRE480-485	8.22	VPLVSLC637-643	1.251	Non-allergenic
	IFYYVFL811-817	−6.157	IFYYVFL811-817	0.882	ACVVSC1207-1212	0.052	GNMDEMA767-773	0.869
BDLF1	STHGTQN33-39	5.386	KLESNDV124-130	1.080	KNKYQP94-99	8.702	LKVVVSL3-9	1.227	Non-allergenic
	IVFPMVL132-138	−5.129	YFCIMVT208-214	0.883	CIMVTL210-215	0.106	RMRVMMT263-269	0.931
BFRF3	DTGSGGG155-161	6.357	GSGGGGQ157-163	1.147	DKRQRA75-80	4.673	SYLVFLT46-52	1.151	Non-allergenic
	RSYLVFL45-51	−3.214	ASAAAAV148-154	0.912	VAGAGA82-87	0.156	DTGSGGG155-161	0.898
BXLF1	GEEDDDG54-60	8.143	QKNGGKG204-210	1.126	KKRGRK447-452	6.122	VLLEVCL538-544	1.254	Non-allergenic
	VFPLMLL401-407	−6.086	LDCWILH385-391	0.880	VLLEVC538-543	0.070	GEEDDDG54-60	0.864
BGLF4	AEDGQDK127-133	6.757	TNSSSSG12-18	1.155	KKRFKE389-394	6.022	VVLLEVL344-350	1.25	Non-allergenic
	VVLLEVL344-350	−4.414	YGIMHLA307-313	0.884	AEVLAC349-354	0.117	DTMWGMG207-213	0.867
BMRF1	SEPEDKS308-314	6.629	RQQGSGG180-186	1.141	RQKQKH387-392	6.086	AHVVCSV237-243	1.249	Non-allergenic
	LSVILFN290-296	−3.686	CALMPYM143-149	0.899	GIIAVV271-276	0.052	TMEYDDK129-135	0.910

Surface flexibility prediction of the respective epitopes was predicted using the Karplus and Schulz flexibility method, which estimates atomic motions based on temperature or B-factor. While a high B-factor depicts a distorted structure, a low B-factor indicates a highly ordered or stable protein secondary structure [87]. Using threshold values set at 1, maximum predicted flexibility scores were 1.071 (LGSGSLA171-177), 1.133 (DQNQNQT96-102), 1.142 (GGGSGGG72-78), 1.115 (QKSNSPR1197-1203), 1.08 (KLESNDV124-130), 1.147 (GSGGGGQ157-163), 1.126 (QKNGGKG204-210), 1.155 (TNSSSSG12-18) and 1.141 (RQQGSGG180-186) respectively for BMRF2, BDLF2, BORF1, BcLF1, BDLF1, BFRF3, BGLF4 and BMRF1. Graphical illustrations coupled with the predicted minimum scores are shown in Supplementary Fig. S3. The Emini Surface Accessibility prediction method was used to determine the probable surface availability or exposure of the predicted epitopes on protein surfaces [86]. The threshold value was set at 1 which implies that epitopes with higher probability scores (>1) tend to localize on the surface of the respective antigenic proteins.

Among all predicted proteins, epitope RYKLKK320-325 from BMRF2 had the highest maximum score of 10.315. Moreover, predicted maximum scores for BDLF2, BORF1, BcLF1, BDLF1, BFRF3, BXLF1, BGLF4 and BMRF1 were 4.805 (QDQNQN95-100), 6.535 (RRRLQR9-14, RRLQRR10-15), 8.22 (PRDRRE480-485), 8.702 (KNKYQP94-99), 4.673 (DKRQRA75-80), 6.122 (KKRGRK447-452), 6.022 (KKRFKE389-394) and 6.086 (RQKQKH387-392) respectively. Predicted minimum scores are 0.068 (VLVVCV226-231) – BMRF2; 0.070 (CCYLAF195-200) – BDLF2; 0.078 (GGVCFV245-250) – BORF1; 0.052 (ACVVSC1207-1212) – BcLF1; 0.106 (CIMVTL210-215) – BDLF1; 0.156 (VAGAGA82-87) – BFRF3; 0.070 (VLLEVC538-543) – BXLF1; 0.117 (AEVLAC349-354) - BGLF4; and 0.052 (GIIAVV271-276) - BMRF1. Graphical illustrations showing the predicted surface accessibility of EBV are shown in Supplementary Fig. S4.

The ability of the epitopes to induce immune B-cell responses were also predicted using the Kolaskar & Tongaonkar antigenicity method. The epitopes with the maximum score were identified as those with greater probability to be immunogenic using a threshold of 1. As predicted, BMRF2 had the most antigenic epitope located at TVLVVCV225-231 with a score of 1.3 > 1 threshold set for predictions.

The allergenicity of the respective B-cell epitopes as predicted by AlgPred showed that all the proteins were non-allergenic as they do not match with any similar IgE epitopes. This is indicated in Table 2. Also, most antigenic epitopes for the other EBV proteins are presented in Table 2 together with their maximum and minimum scores. Graphical illustrations showing the degree of antigenicity per-residue are presented in Supplementary Fig. S5.

Prediction of linear B-cell epitopes

The most potentially antigenic B-cell epitopes selected were those with predictive values of 1 and are presented in Table 3. Using the BCPred server, a total of 16 epitopes were predicted for the nine antigenic EBV proteins, among which 14 had predictive scores of 1 and 13 were 20mer peptides. Two, with starting residue positions (62, 39) were predicted for BDLF2, three for BORF1 (positions 239, 71, 8); two for BcLF1 (positions 262, 475); two for BFRF3 (positions 157, 92); three for BXLF1 (positions 112, 2, 89); one for BGLF4 (position 10) and two for BMRF1 (positions 327, 349).

Table 3

Predicted linear B-cell epitopes of antigenic EBV proteins.

EBV proteins	Starting Position	Sequences	Epitope length	BCPred Score
BMRF2	90	SLRVYSTSTWVSAPCL	16	0.963
	195	IFCARGDHSVASLKET	16	0.939
BDLF2	62	PAAAAVPGPEPEPEPHPPMP	20	1
	39	VFYKAPAPRPREGRASTFHD	20	1
BORF1	239	GGAGGGGGVCFVGGLGLQPC	20	1
	71	AGGGSGGGSFSWSGWRPPVF	20	1
	8	DRRRLQRRIAGLLPPPARRL	20	1
BcLF1	262	KGVSTYTTAKGGEPVGGVFI	20	1
	475	NAAPAPRDRRETYSLQHRRP	20	0.998
BDLF1	189	ALPEVPGPLGLA	12	1
BFRF3	157	GSGGGGQPHDTAPRGARKKQ	20	1
	92	GGSSATPVQQAQAAASAGTG	20	1
BXLF1	112	GAQPPAPAHQKPTAPTPKPR	20	1
	2	AGFPGKEAGPPGGWRKCQED	20	1
	89	AVTSNTGNSPGSRHTSCPFT	20	1
BGLF4	10	SPTNSSSSGELSVSPEPPRE	20	1
BMRF1	327	RPRHTVSPSPSPPPPPRTPT	20	1
	349	SPARPETPSPAIPSHSSNTA	20	1

BcLF1 KGVSTYTTAKGGEPVGGVFI475-494 was selected for its characteristic surface accessibility with a BCPred predictive score of 0.998. B-cell epitopes with predictive scores of 0.96 and 0.94 were selected for BMRF2 being the highest scores obtained for the 16mer peptides at positions (90,195). Moreover, a 12mer peptide was predicted for BDLF1 with a score of 1 at position (189). 75% specificity was used for scoring indicative of its degree of accuracy [91]. These were cross-referenced with inherent or overlapping residues that constitute linear B-cell epitopic peptides predicted by the Emini Surface Accessibility and Kolaskar & Tongaonkar Antigenicity methods. Relatively, we checked for surface accessibility and antigenic properties possessed by the 12-, 16- and 20-mer B-cell epitopes predicted by BCPred for the respective EBV proteins. Our results are presented in Fig. 3.

Fig. 3

Sequential overlapping and analyses of predicted linear B-cell epitopes to identify surface accessible and antigenic attributes.

Prediction of B-cell conformational or discontinuous epitopes

Discontinuous (non-linear) B-cell epitopes were predicted from the 3D structure of the respective antigenic proteins using the Ellipro method which functions based on structural protrusion measured by the protrusion index (PI). Minimum scores were varied from 0.5 to 1.0 to succinctly explore potential epitopes on the protein conformation. Predicted non-linear epitopes and 3D positional representation are shown in Fig. 4. Predicted discontinuous epitopes correlate with some of the linear epitopes earlier determined (surface accessibility and antigenicity) as presented earlier in Table 2.

Fig. 4

Structural analysis of conformational or discontinuous B-cell epitopes (cyan). The locations of the respective epitopes (surface representation) are shown on the 3D structures of the selected EBV antigenic proteins. Corresponding amino acid sequences, as predicted, are also shown (cyan highlights).

Identification of potential T-cell epitopes

As earlier mentioned, the NetCTL method was used to predict epitopes for the selected antigenic proteins, a platform that incorporates the SNN, weight matrix and artificial neural network methods. Scoring of MHC-1 binding affinity was done against HLAs of the A2 supertype among which HLA-A*02:01 is reportedly the most frequent haplotype across various populations [[69], [70], [71], [72]].

Related antigenic processing stages were predicted using the proteasomal C-terminal cleavage and TAP scores, which altogether presents the possibility of identifying the most promising T-cell epitopes for the antigenic proteins. Peptides with nine residues (9mers) were selected since mostly presented peptides are of that length [121]. Possible MHC-1 ligands were identified (<- E) as indicated in Table 4 while the two highly ranked epitopes were selected and shown accordingly.

Table 4

Prediction of antigenic processing and presentation for potential T-cell epitopes of EBV proteins.

Antigenic proteins	Positions	Peptide sequences (T-cell epitopes)	MHC-1 binding affinity	Rescale binding affinity	Proteasomal C-terminal Cleavage	Transport affinity	Prediction score	MHC-1 binding
BMRF2	161–169	FMSPFFIFL	0.9130	1.361	0.9775	1.024	1.5588	< - E
	342–350	FVFTFCEYL	0.7902	1.1779	0.7815	1.103	1.3503	< - E
BDLF2	246–254	KVYTLIPAV	0.8028	1.1967	0.4223	0.74	1.297	< - E
	401–409	RLTDRSFPA	0.7341	1.0942	0.8579	−0.447	1.2006	< - E
BORF1	199–207	YLNKITTVV	0.8500	1.2671	0.973	0.230	1.4245	< - E
	157–165	LLYDSPATL	0.7634	1.1380	0.9742	1.306	1.3494	< - E
BcLF1	401–409	MQYSYYFPV	0.9188	1.3696	0.7904	0.756	1.526	< - E
	1327–1335	FLQEAFPAL	0.8944	1.3332	0.9758	0.901	1.5247	< - E
BDLF1	266–274	VMMTYLQSL	0.8772	1.3076	0.9626	1.335	1.335	< - E
	63–71	YLSDCTLAV	0.9317	1.3889	0.6197	0.451	1.5044	< - E
BFRF3	53–61	SQFCYEEYV	0.7511	1.1196	0.7174	0.5650	1.2554	< - E
	122–131	VAQSATPSV	0.4544	0.6733	0.8997	0.5510	0.8398	< - E
BXLF1	552–560	YLQFVLVDL	0.702	1.0464	0.977	1.076	1.2468	< - E
	548–556	RTLAYLQFV	0.788	1.1746	0.2347	0.565	1.2381	< - E
BGLF4	102–110	KLYDSVTEL	0.8824	1.3154	0.976	1.315	1.5276	< - E
	336–344	FLQFAAPKV	0.8245	1.2291	0.4544	0.155	1.305	< - E
BMRF1	144–152	ALMPYMPPA	0.9089	1.3549	0.6638	−0.3140	1.4388	< - E
	32–40	VLQVNLLSV	0.6918	1.0312	0.8818	0.3340	1.1802	< - E

Structural modeling of selected epitopes and HLA-A*02:01 docking analysis

The selected T-cell epitopes (peptides) were structurally modeled using the PEPFOLD server before docking to HLA-A*02:01. pepATTRACT blind peptide-protein docking method was used to identify the most appropriate binding sites on HLA-A*02:01 for the respective epitopes into which they were docked simultaneously. 51 clusters were obtained for each peptide-protein complex and the best-docked complexes were selected based on global energy scoring and shown in Table 5. Moreover, 3D models of EBV T-cell epitopes with the highest energy scores are shown in Fig. 5.

Table 5

Global energy estimations and MM/GBSA profiling of binding interactions between potential EBV T-cell epitopes and HLA-A*02:01. Estimations of associable energy components are also shown.

Antigenic proteins	Peptide-MHC-1 complexes	pepATTRACT	MM/GBSA estimations
		Global energy (kcal/mol)	ΔEvdw (kcal/mol)	ΔEele (kcal/mol)	ΔEGB (kcal/mol)	ΔESA (kcal/mol)	ΔGbind (kcal/mol)
BMRF2	FMSPFFIFL161-169	−18.44	−11.24	−58.83	61.37	−7.44	−10.18
BDLF2	KVYTLIPAV246-254	−14.50	−10.22	−12.95	17.20	−6.41	−7.25
BORF1	YLNKITTVV199-207	−15.08	−10.42	−11.55	15.85	−6.51	−7.43
BcLF1	FLQEAFPAL1327-1335	−20.42	−10.87	−9.78	14.44	−7.4	−7.69
BDLF1	VMMTYLQSL266-274	−14.37	−8.4	−15.03	17.02	−5.9	−7.59
BFRF3	SQFCYEEYV53-61	−13.89	−11.99	−4.78	11.98	−7.97	−6.38
BXLF1	RTLAYLQFV548-556	−14.79	−16.69	−58.64	64.04	−10.85	−13.47
BGLF4	FLQFAAPKV336-344	−17.40	−10.61	−67.05	69.94	−6.89	−9.1
BMRF1	ALMPYMPPA144-152	−18.54	−10.36	−14.07	16.84	−6.6	−8.88

: electrostatic energy; : van der Waals energy; Δ: total binding free energy; : polar desolvation energy; : non-polar solvation energy.

Fig. 5

3D structures of potential EBV T-cell epitopes for [A] BMRF2 [B] BDLF2 [C] BORF1 [D] BcLF1 [E] BDLF1 [F] BFRF3 [G] BXLF1 [H] BGLF4 [I] BMRF1.

As estimated, BMRF2 T-cell epitope FMSPFFIFL161-169 had the most favorable docking energy of −18.4 kcal/mol compared to FVFTFCEYL342-350.

Likewise, the most favorable peptide-HLA-A*02:01 energy scores based on site binding affinities include −14.5 kcal/mol, −15.1 kcal/mol, −20.4 kcal/mol, −14.4 kcal/mol, −14.8 kcal/mol, −17.4 kcal/mol, −13.89 kcal/mol and −18.5 kcal/mol respectively for BDLF2 (KVYTLIPAV246-254), BORF1 (YLNKITTVV199-208), BcLF1 (FLQEAFPAL1327-1335), BDLF1 (VMMTYLQSL266-274), BXLF1 (RTLAYLQFV548-556), BGLF4 (FLQFAAPKV336-344), BFRF3 (SQFCYEEYV53-61), and BMRF1 (ALMPYMPPA144-152).

Amongst all predicted T-cell epitopes, BcLF1 epitope BcLF1 (FLQEAFPAL1327-1335) had the highest docking score. Also, we observed that the binding of all predicted T-cell epitopes occurred majorly across two sites on HLA-A*02:01 designated Site 1 and 2, which are depicted in Fig. 6. The predicted Site 2 into which some of the epitopes were docked correlates with the binding site determined in HLA-A*02:01 crystal structure (PDB ID: 4UQ3). As predicted, T-cell epitopes for EBV proteins; BMRF2, BXLF1, and BGLF4 exhibited a higher affinity for Site 1 where they were docked while epitopes for BDLF2, BFRF3, BORF1, BcLF1, BDLF1, BFRF3, and BMRF1 had a higher binding propensity for Site 2 as predicted by pepATTRACT.

Fig. 6

Potential HLA-A*02:01 binding sites predicted for the respective EBV T-cell epitopes. [A] Site 1 and corresponding constituent residues [B] Site 2 and corresponding constituent residues. Also shown are the cohort of epitopes that bind with selective affinity to the respective sites as predicted by pepATTRACT blind-docking algorithm.

MD simulation and MM/GBSA-based interaction analyses

MD analyses were performed to obtain molecular insights into the binding and interaction dynamics of the predicted EBV T-cell epitopes and HLA-A*02:01. 100ns MD simulation runs were performed for the complexes and individual epitopes while resulting trajectories were analyzed.

This was performed using post-analytical metrics such as the root mean square deviation (RMSD) and root mean square fluctuation (RMSF) methods which were used to estimate trajectorial motions of constituent Cα atoms. Results showing estimated mean values for the epitopes and their corresponding MHC-1 complexes are presented in Table 6.

Table 6

Post–MD simulation analyses of structural motions and stability for the T-cell epitopes and peptide-protein complexes.

Molecular dynamic simulation analyses
EBV proteins	T-cell Epitope RMSD (Å)	T-cell Epitope-MHC complex
		Cα- RMSD (Å)	Cα- RMSF (Å)
BMRF2 [FMSPFFIFL161-169]	2.4 ± 0.4	4.2 ± 1.0	1.9 ± 0.8
BDLF2 [KVYTLIPAV246-254]	4.2 ± 1.0	6.8 ± 1.9	3.0 ± 1.2
BORF1 [YLNKITTVV199-207]	4.1 ± 1.0	4.1 ± 1.7	3.0 ± 1.3
BcLF1 [FLQEAFPAL1327-1335]	2.7 ± 0.5	6.4 ± 3.5	4.2 ± 1.8
BDLF1 [VMMTYLQSL266-274]	3.9 ± 0.6	5.5 ± 2.3	3.3 ± 2.0
BFRF3 [SQFCYEEYV53-61]	3.2 ± 0.3	4.3 ± 2.4	3.1 ± 1.1
BXLF1 [RTLAYLQFV548-556]	4.9 ± 1.0	2.8 ± 0.4	1.9 ± 0.6
BGLF4 [FLQFAAPKV336-344]	3.6 ± 1.5	9.0 ± 3.0	3.9 ± 1.7
BMRF1 [ALMPYMPPA144-152]	3.5 ± 0.5	5.6 ± 2.9	3.7 ± 1.4

Comparatively Site 1–bound BXLF1 T-cell epitope [RTLAYLQFV548-556] exhibited the highest RMSD with a mean value of 4.9 Å while the lowest RMSD of 2.4 Å was exhibited by BMRF2 [FMSPFFIFL161-169], indicative of its highly stable motion at the target Site 1 region of HLA-A*02:01. The trajectorial motions of the respective epitopes at the target Sites (1 and 2) are shown in Fig. 7. Moreover, it can be deduced from the plot that while their motions were highly unstable at the initial periods of the MD simulation, they seemingly attained stability towards the end of the simulation (∼50–80ns) which could mediate optimal binding and interactions at their target sites. Dynamical motions of the peptide-MHC complexes were also measured using the RMSD and RMSF metrics. Mean values are also presented in Table 6 while corresponding plots are shown in Fig. 7. BXLF1 [RTLAYLQFV548-556]-MHC complex was more stable while the highest structural instability was exhibited by the BGLF4 [FLQFAAPKV336-344]-MHC complex.

Fig. 7

Graphical plots depicting structural motions and residual fluctuations of potential EBV T-cell epitopes and protein complexes over 100ns MD simulation run [A] peptide-protein complex RMSD [B] peptide RMSD [C] peptide-protein complex RMSF [D] 3D representation of T-cell epitope binding on HLA-A*02:01 Sites 1 and 2.

Presumably, high RMSDs exhibited by the epitopes at the respective binding sites could correlate with systematic motions necessary to achieve optimal and high-affinity binding with target residues. Consequently, this could account for the complementary increase in structural motions and deviations observed for the peptide-MHC complexes until relative stability was attained.

Complementary interaction dynamics of the EBV T-cell epitopes were also investigated while binding affinities were measured using the MM/GBSA method. Findings revealed that predicted Site 1 HLA-A*02:01 binders; BMRF2 [FMSPFFIFL161-169], BXLF1 [RTLAYLQFV548-556] and BGLF4 [FLQFAAPKV336-344] had total binding energies () of −10.18, −13.85 and −9.1 kcal/mol, which are higher when compared to other EBV T-cell epitopes that are Site 2 binders as shown in Table 5.

Estimations of other energy components revealed that both van der Waals () and electrostatic energies () contributed to the while positive polar desolvation energy values () indicated that binding interactions were favorable in the hydrophobic regions. This is supported by the estimated negative non-polar energies () as shown.

Analysis of the intermolecular interactions between the T-cell epitopes and target MHC-1 molecule revealed a similar binding pattern across both sites as shown in Fig. 8.

Fig. 8

Analysis of intermolecular interactions mediated by potential EBV T-cell epitopes at their respective HLA-A*02:01 binding sites.

For Site 1 binders; amino acid residues Asp29, Asp30 and Gln115 (amongst other interacting residues) were crucial to the binding and stability of the epitopes; BMRF2 [FMSPFFIFL161-169], BXLF1 [RTLAYLQFV548-556] and BGLF4 [FLQFAAPKV336-344]. These include the formation of high-affinity hydrogen bonds with proximal residues of the T-cell epitopes. As shown, Site 1 residues; Asp29, Asp30 and Gln115 formed H bonds with Met2, Ser3 and Ile7 of BMRF2 [FMSPFFIFL161-169] respectively. Likewise, Asp30 and Gln115 recurrently interacted with Arg1 and Gln7 of BXLF1 [RTLAYLQFV548-556] respectively.

Similarly in BGLF4 [FLQFAAPKV336-344] complex, Asp29 formed two H bonds with Gln3 and Phe4 whilst also involved in H bond formation with Phe1. On the other hand, Asp30 existed in hydrophobic interactions with the hydrophobic side chain of Phe4.

Also, we identified important interactions that account for the stability of Site 2 binders at the target region (Fig. 8). A highly important residue that was uniformly involved in complementary H bonding with the respective EBV epitopes is Gln155. At the target Site 2, His70, Thr73, Lys246 and, Gln155 contributed to the binding of BcLF1 [FLQEAFPAL1327-1335] via H bond interactions while Gln155, Arg97, Lys66, Trp167 were involved in the binding of BDLF1 [VMMTYLQSL266-274]. Moreover, we observed that Lys146 and Arg65 played key roles in the binding of BDLF2 [KVYTLIPAV246-254] while BFRF3 [SQFCYEEYV53-61] binding was stabilized by interactions with His70, Gln155, Thr73, Arg97 and, Lys146. Gln155 and Lys66 were crucial for the binding of BMRF1 [ALMPYMPPA144-152] while constituent BORF1 [YLNKITTVV199-208] residues; Leu2, Ile5, Thr7 and, Val8 interacted with Site 2 residues; Thr163, Arg97, Thr73, and Gln155 respectively as shown in Fig. 8.

Discussion

Efforts to curtail the virulent activities of EBV have been ongoing requiring several therapeutic options ranging from the development of small molecule inhibitors to vaccination. While antagonists are directed towards crucial viral proteins [122,123], vaccines are developed to stimulate immunogenic B-cell or T-cell responses towards lytic and latent EBV proteins [2,15]. Although efforts to develop effective cures for EBV-mediated infections have been unsuccessful, a gp350 vaccine has made significant progress through clinical trials even though it does not curtail viral infections [2,[33], [34], [35]].

Among the various vaccine types, peptide or subunit vaccines have emerged as attractive candidates for treatment since there is a huge possibility of designing them based on essential information from viral pathogenic components. This is achievable via the implementation of immunoinformatics approaches which has been widely employed for B- and T-cells epitopic predictions from viral proteins [31,53,56,[61], [62], [63], [64], [65]]. The antigenic nature of viral proteins corresponds with their ability to elicit adaptive immunogenic responses by host B- and T-cells which could, in turn, lead to the death of infected cells. Molecular components (antigens) and not the whole pathogens are recognized by receptors on B- and T-cells; B-cell receptors (BCR) and T-cell receptors (TCR), a process that differs greatly for both cells. While BCRs recognize and bind surface-exposed antigens (epitopes) on viral proteins to secrete antibodies and mediate humoral adaptive immunity, TCRs recognize antigens presented on the surfaces of antigen-presenting cells (APCs) that are bound to major histocompatibility complex (MHC) molecules categorized into MHC-I and II. These are respectively recognized by CD8+ (cytotoxic) and CD4+ (helper) T-cells. Altogether, the pathway for antigenic presentation (at the extracellular surface) incorporates proteasomal degradation, TAP processing and MHC-1 [human leukocytes antigen – (HLA) in humans] binding [31,[124], [125], [126]]. Several crucial functional proteins (lytic and latent) constitute the EBV proteome and contribute to its virulence in host cells. Similar to EBV glycoproteins [8,53], proteins that regulate epithelial cell attachment, DNA synthesis/replication and, capsid assembly are highly important in the viral latent and lytic cycles. These proteins include BMRF2, BDLF2, BORF1, BcLF1, BFRF3, BXLF1, BGLF4, and BMRF1 (see Table 1). Therefore, immunoinformatics prediction of potential B-cell and T-cell (CD8+) epitopes for these proteins represent the highlight of this study.

Potential linear B-cell epitopes reported in this study were selected based on their characteristic surface-exposure/accessibility and antigenicity. These attributes are essential to facilitate their recognition and binding to paratopes on BCRs provoking humoral immunogenic responses. Therefore, linear (continuous) B-cell epitopes; BMRF2 [SLRVYSTSTWVSAPCL90-105], BDLF2 [PAAAAVPGPEPEPEPHPPMP62-81], BORF1 [DRRRLQRRIAGLLPPPARRL8-27], BcLF1 [KGVSTYTTAKGGEPVGGVFI262-281], BGLF1 [SPTNSSSGELSVSPEPPRE10-30] and [RPRHTVSPSPSPPPPPRTPT327-347] all possess high propensities for eliciting B-cell mediated immune responses. Since most B-cell epitopes are conformational or discontinuous, we also employed prediction algorithms for 3D conformational epitopes. As shown in Fig. 4, conformational epitopes for the respective EBV antigens are highly surface exposed. Antigenic and allergenic profiling was also carried out thereby presenting our predicted epitopes as potential candidates for the development of EBV B-cell vaccines. The presence of residues with hydrophilic, large and aromatic side chains could provide a clue into the characteristic surface-exposure and immunogenicity exhibited by the predicted epitopes [97].

T-cell epitope prediction was also carried out to identify the shortest peptides that constitute these antigens and can stimulate CD8 T-cells activation [56,127]. Herein, binding to a highly frequent human MHC I molecule (HLA-A*02:01), TAP processing and C-terminal cleavage were considered which resulted in the selection of T-cell epitopes; BMRF2 [FMSPFFIFL161-169], BDLF2 [KVYTLIPAV246-254], BORF1 [YLNKITTVV199-207], BcLF1 [FLQEAFPAL1327-1335], BDLF1 [VMMTYLQSL266-274], BFRF3 [SQFCYEEYV53-61], BXLF1 [RTLAYLQFV548-556], BGLF4 [FLQFAAPKV336-344] and BMRF1 [ALMPYMPPA144-152]. With a 75% threshold, higher predictive scores obtained could indicate the immunogenic propensities of these T-cell epitopes. Information on the peptide-blinding clefts of MHC molecules is also essential since they could be vital for the structure-based design of high-affinity and selective epitopes for peptide-vaccine design [56,128]. Without restraining epitopic binding to a particular site on HLA-A*02:01, two sites were identified using a blind peptide-protein docking algorithm, one of which corresponds to an experimentally defined site (PDB entry 4UQ3). These sites were identified based on a global search of the protein surface after which the peptides (conformation) were docked to the corresponding sites [102,103]. Investigating the binding dynamics, affinity, and stability of the T-cell epitopes at their target HLA-A*02:01 sites using MD and MM/GBSA methods further yielded interesting insights, particularly for the identification of crucial site residues that mediate high-affinity peptide-protein interactions.

MD simulation analyses revealed the dynamical motions of HLA-A*02:01 when bound by the predicted T-cell epitopes from the respective EBV antigens. According to previous reports, high structural motions are necessary for MHC-1 molecules and correlates with their ability to translocate the processed peptide (proteasomal C-terminal cleavage and TAP) from the endoplasmic reticulum (ER) lumen to the Golgi and ultimately to the extracellular space [30,128,129]. As shown in Fig. 7, the complexes demonstrated high RMSDs indicative of high structural motions which became more stable towards the end of the simulation period. These motion patterns were also reflected by the RMSF analysis which estimates residual fluctuations. While notably high motions were observed in HLA-A*02:01 bound by predicted T-cell epitopes from BcLF1 [FLQEAFPAL1327-1335], BGLF4 [FLQFAAPKV336-344], BMRF1 [ALMPYMPPA144-152] and BFRF3 [SQFCYEEYV53-61], low structural motions were observed in the BXLF1 complex. However, the T-cell epitope from BXLF1 [RTLAYLQFV548-556] exhibited the highest RMSD which could indicate its degree of motion at its target site (Site 1). Recurring inter-molecular interactions between specific HLA-A*02:01 site residues and the predicted T-cell epitopes further explained binding patterns relative to high-affinity binding and stabilization. Identification of a novel site (Site 1) for epitope binding on HLA-A*02:01 could also aid in the structure-based design of high-affinity immunogenic binders for the target MHC-I molecule. This is because the most favorable ΔGs estimated were for T-cell epitopes bound to Site 1; BMRF2 [SLRVYSTSTWVSAPCL90-105] (−10.18 kcal/mol), BGLF4 [SPTNSSSGELSVSPEPPRE10-30] (−9.1 kcal/mol) and BXLF1 [RTLAYLQFV548-556] (−13.471 kcal/mol) compared to those bound to Site 2.

As proposed in this study, while BMRF2 [SLRVYSTSTWVSAPCL90-105], BDLF2 [PAAAAVPGPEPEPEPHPPMP62-81], BORF1 [DRRRLQRRIAGLLPPPARRL8-27], BcLF1 [KGVSTYTTAKGGEPVGGVFI262-281], BGLF4 [SPTNSSSGELSVSPEPPRE10-30] and [RPRHTVSPSPSPPPPPRTPT327-347] depict potential epitopes that could elicit immunogenic B-cell responses based on physicochemical attributes, BMRF2 [FMSPFFIFL161-169], BDLF2 [KVYTLIPAV246-254], BORF1 [YLNKITTVV199-207], BcLF1 [FLQEAFPAL1327-1335], BDLF1 [VMMTYLQSL266-274], BFRF3 [SQFCYEEYV53-61], BXLF1 [RTLAYLQFV548-556], BGLF4 [FLQFAAPKV336-344] and BMRF1 [ALMPYMPPA144-152] have high propensities of provoking T-cell mediated cytotoxic immune responses towards infected cells.

An interesting futuristic approach is the possible derivation of highly selective antagonists from the predicted epitopes, as reported in a study by Abdul et al. [59]. This could present avenues to reverse-map highly interactive regions of the potential EBV epitopes and generate pharmacophore models that can be used to screen virtually for possible lead candidates. These derived compounds would presumably be able to achieve selective targeting and inhibition of crucial pathogenic machinery of EBV.

Also, more precise and accurate methods can be incorporated to distinguish epitopes from non-epitopes from the large pool of predicted epitopes. This is an effective approach as demonstrated in a previous study by Khanna and Rana where they implemented a machine learning model for determining epitopes of Mycobacterium tuberculosis [130].

Immunoinformatics predictions performed in this study covered multiple EBV antigens thereby presenting a possibility of developing highly effective anti-EBV vaccines consisting of a mixture of antigens from these pools of B- and T-cell epitopes [68].

Furthermore, while we propose the ability of these predicted epitopes to provoke efficacious immune responses, it is still necessary to validate their propensity using in vitro and in vivo techniques in the future. These could include binding assays for in vitro determination of peptide binding and/or ELISPOT or CFSE dilution assays for estimating T-cell recognition [131].

Conclusion

The recent paradigm shift towards peptide or subunit vaccine development has been highly resourceful, particularly with regards to provoking adaptive or cell-mediated immune responses. These are designed based on B-cell and T-cell epitopes derived as components of antigenic proteins secreted by viruses during their infectious life cycles. Although no significant progress has been made to date, vaccination remains a promising treatment option to curtail the virulence of Epstein Barr virus (EBV) implicated in infectious mononucleosis (IM) and several human cancers. In this study, we implement a computational meta-analysis paradigm that integrates immunoinformatics, blind peptide-protein docking, and molecular dynamics to predict potential B-cell and T-cell epitopes for antigenic EBV proteins that mediate epithelial cell–cell attachment/spread, capsid self-assembly, DNA synthesis, replication and, processivity. Predictive BCPred and Ellipro algorithms were used to rank linear and conformational B-cell epitopes based on their physicochemical properties such as surface-exposure or accessibility, flexibility, hydrophilicity and, antigenicity. More so, based on the Major Histocompatibility Complex (MHC) class 1 antigenic presentation pathway, potential T-cell epitopes were predicted using the integrative weight matrix and Artificial Neural Network (ANN) algorithm of the NetCTL 1.2 webserver. Peptide modeling and blind docking were used to identify probable binding sites on a highly frequent human MHC-1 molecule; HLA-A*02:01 while varying binding affinities for the peptide-MHC complexes were obtained. A novel high-affinity binding site (Site 1) was identified in addition to a structurally resolved site (Site 2). Dynamical studies revealed that the binding of the T-cell epitopes induced a high degree of structural motions in HLA-A*02:01 with favorable binding energy (ΔG) that correlates with their stabilities at their targeted sites. Our results suggest potential B- and T-cell epitopes that could serve as effective candidates for peptide-vaccine development towards the treatment of EBV-related diseases.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.