Integrative Machine Learning and Immunoinformatics-Guided Design and In Silico Validation of a Multiepitope DNA Vaccine Against Avian Metapneumovirus

1. Introduction

Avian Metapneumovirus (aMPV) is one of the currently neglected respiratory viral poultry pathogens even in the US. The aMPV causes several clinical diseases in various species of birds including the swollen head syndrome in chickens and rhinotracheitis in turkeys [1]. Although the main tropism of the aMPV is main targeting the respiratory system, the virus could do further damage to other tissues like the reproductive system. The aMPV infection may enhance the development of other infections in birds such as colibacillosis, causing serious damage to infected chickens [2]. The aMPV is classified under the family Pneumoviridae and genus Metapneumovirus. The virus is subclassified into several subtypes [2]. Subtypes A and B found mainly in Europe, Subtype C in the United States, Canada, China, France, and South Korea, and Subtype D found only in France [1]. Subtype D has not been reported since its first discovery in 1985, while Australia remains the only continent considered aMPV and overall, subtype B remains the most prevalent subtype of the current groupings [3].

Due to the lack of monitoring of the virus on the genomic levels, the virus has the ability to continue to mutate and evolve, particularly within the migratory birds. This hypothesis is supported by the most recent outbreak of aMPV subtypes A and B reported in the United States in 2024, where infection was extremely high, with mortality rates reported as high as 80% [4]. At the current time, there are no therapeutic interventions for aMPV, in which preventative control strategies rely heavily on diagnostic techniques such as RT-PCR and vaccination measures as aMPV resembles other respiratory diseases such as Newcastle disease virus and avian influenza virus [5]. While vaccines for aMPV have been developed and used since the 1980s, continuous outbreaks point towards vaccine failures with many vaccines coming in the form of live vaccines that research indicates can revert to a virulent form [1]. This point towards the urgency of the design and development of some new generation’s vaccines against aMPV to cope-up with the dynamic changes of the virus on the genomic level and to avoid vaccine escape and failure of the currently used vaccines. Ideally, the new generation of vaccines should be prepared against the currently circulating strains and subtypes of the virus in certain area. Thus, continuous monitoring of viruses not only aMPV but also other viral diseases of human, birds and animals provide a basic foundation for the design of any new generation of vaccines [6,7,8,9,10,11,12,13].

The genome of aMPV is a negative sense RNA with 8 genes encoding various proteins: nucleoprotein (N), phosphoprotein (P), matrix protein (M), fusion protein (F), second matrix protein (M2), small hydrophobic protein (SH), attachment protein (G), and large polymerase protein (L). Of the proteins the F and G proteins are integral for virus attachment and entry with the G protein playing a role in mediating the attachment to the host cell and the F protein promoting pH-independent fusion between the cell membrane and virus envelope [14].

The gene encodes the G protein has an extensive history of genetic mutations like point mutations, insertions and deletions, leading to the varying size due to truncation of the G gene [15]. The protein stands out due to its unique sequence among other paramyxoviruses, where the hemagglutinin-neuraminidase (HN) attachment differs significantly [1]. The protein is heavily glycosylated with both N- and O-glycosylated sites suggesting a reduced antigenicity through shielding of the virus protein with host-specified sugars [14].

On the other hand, the F protein is highly conserved, less variable than the G protein when comparing between the various subtypes of aMPV [16]. The native aMPV-F protein is present in an inactive form called F0 form, cleaved into two disulfide-linked subunits, F1 and F2, mediated by host cell proteases such as Furin, critical for the viral fusion to the host cell membrane as well as for viral entry [1]. The F protein can bind to the cellular receptors and facilitate fusion with the virus and host cell membranes without the involvement of the G or SH proteins, also having another capability of mediating receptor binding as well [5]. There have been several approaches to vaccine work, including live, inactivated, recombinant and DNA vaccines with the need for a higher level of protection against aMPV, with each variation having its advantages and drawbacks. The integration of the genome sequencing data of the currently circulating strains of the target viruses with the newly developed machine learning, molecular docking and immune simulation paved the ways for the development of many vaccines against various viruses of humans/animals and birds [17,18,19,20,21]. In the current study we incorporated the genome data and the machine learning tools into the design and development of a vaccines for aMPV against this emerging viral pathogen [22,23,24,25].

The ability to rapidly screen datasets for identification of epitopes, predict model structures, and visualize ligand-receptor interactions using artificial intelligence shows great promise for vaccine development, changing the time and cost of development immensely [26,27,28,29]. The use of the DNA vaccines is favorably looked towards with the use of these types of vaccine constructs for its low cost development as well as its short synthesization time and stability [30]. When combined, genome sequencing and machine learning tools provide strong tools for the development of next generation vaccines with high levels of precision and accuracy [20,31,32]. Various vaccine designs of this nature pave the way for innovation and better vaccines for humans and animals.

The present study aimed to design as multiepitope DNA-based vaccine, using top ranked B and T cell epitopes using two major proteins for aMPV (F and G). The designed DNA vaccine construct demonstrates a high antigenicity and strong binding affinity for chicken innate immune receptor (TLR7). In silico cloning as well as immune simulation were performed, showing great potential in stimulating an immune response through vaccination. The delivery of the DNA vaccines has hampered the success of this type of vaccines for decades however incorporation of the DNA vaccines with nanoparticles enhanced the delivery of the designed vaccines and proved to produce robust immune response in the vaccinated hosts [33,34,35,36].

The current data provides an in-silico blueprint in which, to validate the vaccine construct for application, further experimental studies will need to be performed, with in vivo studies and studies on the vaccination of chickens and other birds like turkeys and ducks which are commonly susceptible to aMPV.

2. Materials and Methods

2.1. Retrieval of Avian Metapneumovirus Protein Sequences and Multiple Sequence Alignment

The amino acid sequence of the Fusion (F) and Attachment (G) Proteins from various regions such as USA, South Korea, Brazil, China, and Netherlands in the last 15 years were retrieved from the National Center for Biotechnology (NCBI) GenBank database (https://www.ncbi.nlm.nih.gov/protein/). Sequences from subtypes A, B, and C were included. Any sequences with ambiguous bases, partial sequences, and duplicate isolates were excluded. Isolates of the same exact sequence were reduced to not overpower the sequence alignment. In total 53 sequences of the Fusion Protein and 68 sequences of the Attachment Protein were used to perform Multiple Sequence Alignment (MSA) using Geneious v2026.1.2 (www.geneious.com) and Clustal Omega v1.2.2 (https://www.clustal.org/people-search/) with default parameters [37]. Supplementary Table 1 provides all information about these sequences. The most conserved regions of the proteins were further considered for B cell and T cell (MHC class I and class II) epitope prediction. The Geneious software was then used to create a phylogenetic tree using PHYML, with 1000 bootstraps [38]. The phylogenetic trees of each protein can be seen in Supplementary Figure 1 and Figure 2.

2.2. Mapping of B-Cell Epitopes of the Fusion and Attachment Proteins in Avian Metapneumovirus

2.2.1. Prediction of Linear B Cell Epitopes

For the mapping of B cell epitopes using the generated consensus sequences for the two proteins (F and G), ABCpred was used (https://webs.iiitd.edu.in/raghava/abcpred/) [39]. The length of target epitopes was set to 16 and 18 mers with default parameters of 0.51 as the threshold. Epitopes were filtered based on matching the predicted linear epitopes from IEBD analysis resource server (https://tools.iedb.org/ellipro/) [40]. Further filtration of the epitopes were performed based on the antigenicity, allergenicity, toxicity. The antigenicity was assessed using the default threshold of 0.4 on the VaxiJen v2.0 server (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html) [41]. Allergenicity was assessed using the AllerTOP v.2.1 server (http://ddg-pharmfac.net/AllerTOP/) [42]. The toxicity was assessed using the default settings of the ToxinPred server (https://webs.iiitd.edu.in/raghava/toxinpred/multi_submit.php) [43].

2.2.2. Prediction of Discontinuous/Conformational B Cell Epitopes

Discontinuous/Conformational B cell epitopes were identified using IEDB’s prediction platform ElliPro (https://tools.iedb.org/ellipro/) [40]. The PDB files for the proteins were created using trRosetta (https://yanglab.qd.sdu.edu.cn/trRosetta/) [44]. The default settings were used with 0.5 as the minimum score and 6 angstrom for the maximum distance. Visualization of the epitopes was completed using Pymol version 2.5.4.

2.3. Mapping of T-Cell Epitopes of the Fusion and Attachment Proteins in Avian Metapneumovirus

Prediction of Cytotoxic T-Cell Epitopes (MHC Class I)

The cytotoxic T-cell epitopes were selected based on the chicken MHC I peptide-binding motifs in accordance with the BF24, BF212, BF215, and BF219 haplotypes or BF2∗02:01/BF2∗04:01/BF2∗21:01 as previously described [45,46]. Protein sequences were screened using the following patterns: x-(D or E)-x-x-(D or E)-x-x-E, x-x-x-x-(V or I)-x-x-V, x-x-x-x-(V or I)-x-x-x-V, (K or R)-R-x-x-x-x-x-Y, (K or R)-R-x-x-x-x-x-x-Y, x-R-x-x-x-x-x-(Y,P,L,F), x-R-x-x-x-x-x-x-(Y,P,L,F), and x-A/V/I/L/P/S/G-x-x-x-x-x-x-V/I/L. Of the epitopes matching the patterns, they are screened for the antigenicity, allergenicity, and toxicity as previously described [41,42,43].

Prediction of Helper T-Cell Epitopes (MHC Class II)

Due to the absence of chicken alleles in the IEDB database, the best human allele substitutes were used for MHC class II epitope prediction. Through the analysis of heatmap data in MHCcluster, the best substitute alleles for MHC class II are DRB1:1310, DRB1:1366, DRB1:1445, and DRB1:1482 [47]. The helper T-cell epitopes were identified using the IEDB server(https://tools.iedb.org/mhcii/) [48]. The binding predictor NetMHCII pan 4.1 BA (recommended binding predictor: 2023.09] was used, sporting the peptides in accordance to the inhibitory concentration (IC50) and percentile rank [49]. Peptides with an IC50<50nM and percentile rank<4 were chosen and further filtered via antigenicity, allergenicity, and toxicity as previously described [41,42,43]. MixMHC2pred (http://mixmhc2pred.gfellerlab.org/) using the MHC class II alleles Gaga_BLB1_002_01, Gaga_BLB1_012_01, Gaga_BLB2_002_01, Gaga_BLB2_012_01, Gaga_BLB2_012_02, was then used to predict the chicken alleles associated with the peptides [50].

2.4. Molecular Docking of Predicted T-Cell Epitopes to Chicken MHC-I and MHC-II Alleles

The chicken MHC allele sequence files (BF and BL) were retrieved from the Protein Data Bank (PDB). Molecular docking of the predicted MHC-I and MHC-II to the respective alleles was conducted using HDOCK (http://hdock.phys.hust.edu.cn/) [51]. Epitopes were filtered based upon their docking ability. The Pymol software was used to visualize the key residues, as well as information about the hydrogen bond donors and acceptors.

2.5. Assessment of Top Ranked Epitopes Affinity to Chicken Genome

The full amino acid sequence of the chicken (Gallus gallus) genome was retrieved from the National Center for Biotechnology (NCBI) GenBank database (https://www.ncbi.nlm.nih.gov/protein/) [52]. The sequence was then uploaded to Geneious v2026.1.2 (www.geneious.com). The top ranked epitopes for B and T cells were assessed for similarity to the chicken genome to prevent any unwanted autoimmune reaction to the designed vaccine. The predicted epitopes were mapped onto their respective protein structures using Pymol to assess surface exposure and estimate the accessibility of the epitopes to antibody recognition was done using solvent-accessible surface area (SASA) analysis. Top ranked epitopes were also assessed for conservation against the selected isolates for the project.

2.6. Assembly of Multi-Epitope Vaccine using Top Ranked Epitopes

The vaccine sequence was arranged starting with the Kozak sequence followed by a start codon. The primary area of the vaccine comprises the top ranked B and T cell epitopes based on the filtering values such as antigenicity, allergenicity, toxicity, molecular docking, and alignment to host genome. Linkers such as KK, GPGPG, and AAY were used to link the epitopes. Towards the end of the C-terminal end of the vaccine, chicken IL-18 is linked (Accession No. CAB96214) as an adjuvant as well as EAAK linker [53,54]. Prior to the stop codon, a 6 x His-tag (H) tag is added to the construct for expressions identification of the vaccine.

2.7. Analysis of Mulit-Epitope Vaccine Disulfide Bond Engineering

For the stability of the protein, it is essential to create disulfide bonds between the cysteine residues. The vaccine protein was submitted to Disulfide by Design v2.12 web server (https://cptweb.cpt.wayne.edu/DbD2/) for disulfide bridging [55]. The X³ value and the Cα-Cβ-Sγ angle were kept as the default parameters. Based upon previous studies, the X³ value must range between -87 to +97◦, as well as the energy score must not exceed 2.2 kcal/mol for the disulfide bridging. For bridging between possible residue pairs, mutations to cysteine residues occurred using the server’s “Create/View Mutant” option.

2.8. Assessment of the Physiochemical Properties and Stability of the Designed Multi-Epitope Pan-Avian Metapneumovirus Vaccine

The designed vaccine was assessed for its physiochemical properties using Protparam server (https://web.expasy.org/protparam/) [56]. The antigenicity was assessed using the default threshold of 0.4 on the VaxiJen v2.0 server (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html) [41]. Allergenicity of the vaccine was assessed using the AllerTOP v.2.1 server (http://ddg-pharmfac.net/AllerTOP/) [42]. The solubility of the vaccine construct was determined using the Innovagen Peptide solubility calculator (https://pepcalc.com/peptide-solubility-calculator.php). In terms of the stability and physical mobility of the atoms, the vaccine construct was analyzed using normal mode analysis (NMA) using the iMODS server tool (https://imods.iqf.csic.es/) [57]. Multiple descriptors such as eigenvalues, covariance, B-factors, and deformability and intrinsic motions of multiplex proteins were analyzed. Eigenvalues were used to assess motif stiffness and the biological targets’ efficacy measurement was used to predict deformity of the main chain.

2.9. Predictions of the Secondary and Tertiary Structures of Multi-Epitope Vaccine

The PDB file for the vaccine was created using trRosetta (https://yanglab.qd.sdu.edu.cn/trRosetta/) [44]. The secondary structure of the vaccine was predicted using PSIPRED (https://bioinf.cs.ucl.ac.uk/psipred/) [58]. The tertiary structure of the vaccine was predicted using AlphaFold Server (https://alphafoldserver.com/) [59]. The ProSA-web server (https://prosa.services.came.sbg.ac.at/prosa.php) was used to assess the quality of the model [60]. The stability of the construct was analyzed using Ramplot (ramplot.in) [61].

2.10. Molecular Docking of Designed Multi-Epitope Vaccine to Chicken Toll-Like Receptor (TLR)

Chicken toll-like receptor (TLR) was selected to be used for molecular docking with the vaccine construct. The chicken TLR7 (UniProt ID: A0A1L4FML6) was retrieved from Uniprot and the molecular docking was conducted using HDOCK (http://hdock.phys.hust.edu.cn/), showing the structures of the two molecules as well as the docking interactions [51]. Visualization of the docking was done by the Pymol server. A detailed analysis of the binding interfaces was characterized by the Tamarind Bio PDBsum tool (https://app.tamarind.bio/app) [62].

2.11. In Silico Immune Simulation of Designed Multi-Epitope Vaccine

To model the immune response to vaccination the C-ImmSim server (https://kraken.iac.rm.cnr.it/C-IMMSIM/index.php) was used [63]. The use of the server is to simulate the injection of the vaccine construct, modeling the three major immune organs (bone marrow, thymus, and spleen) and the influence of the construct on the B and T lymphocytes. The default parameters were used with 50 and 1000 simulation steps. It is proposed that the multiepitope vaccine construct is administered three times with 4-week intervals. Every step of the simulations is indicative of eight hours of real time, with periods equal to 1, 84, and 164 hours. Host cellular immune response as well as levels of cytokine expression induced by the vaccine construct are predicted by the server as well. The amino acid sequence of the vaccine was converted into a nucleotide sequence using VectorBuilder’s Codon Optimization tool (https://en.vectorbuilder.com/tool/codon-optimization.html). The sequence was inserted into the vector pcDNA3.1(+) with no overhangs accessed from the SnapGene website (https://www.snapgene.com/plasmids/image_consortium_plasmids/pcDNA3.1(%2B)) using the Insertion Cloning Action.

3. Results

3.1. Resuts of the analyssi of the currently circualting aMPV-F and G proteins strains

Protein sequences of the F and G proteins in aMPV were extracted to perform multiple sequence alignment in Geneious. Sequences feature various avian species such as chickens, turkeys, and ducks. The list of the sequences used, and sample information can be found in Figure 1. In the analysis of the multiple sequence alignment, similarity levels can be seen as 87.2% for F protein and 66.1% for G protein. This indicates a higher level of genetic similarity in the F protein sequence among avian hosts than the G protein sequence.

3.2. B Cell Epitope prediction and evaluation among the F and G proteins of the currently circaulting aMPV subtypes

Humoral immune response elicited by B cell epitopes, becoming a major factor in vaccine development. Linear/continuous B cell epitopes were predicted for aMPV proteins using the ABCpred computational tool. For the F protein, 54 epitopes were predicted and for the G protein, 46 epitopes were predicted. All the predicted epitopes were filtered for antigenicity, allergenicity, and toxicity as well as matching the predicted linear epitopes from the IEDB Ellipro tool. The remaining epitopes are shown in Table 1.

From the Ellipro server, the predicted conformational/discontinuous epitopes were also predicted using the 3D structure of the proteins, using the default settings of a minimum score of 0.5 and a minimum distance of 6Å. The respective epitope predictions for each protein can be found in Table 2. Visualization of the epitopes can be found in Figure 2 and Figure 3.

3.3. Cytotoxic T Lymphocyte (MHC Class I) and Helper T Lymphocyte (MHC Class II) Epitope Predictions and Assessments

Cytotoxic and helper T cells are important in the signaling and fighting infection on the cellular level. In Table 3, the predicted MHC class I epitopes based upon the published patterns are listed. All epitopes were screened for antigenicity, allergenicity and toxicity. For cytotoxic T lymphocytes, 27 epitopes were found for F protein, and 21 epitopes were found for G protein.

In Table 4, the predicted MHC class II epitopes are listed with the corresponding binding affinity (IC50 < 50 nM) for the aMPV proteins. Each of the epitopes were screened for antigenicity, allergenicity and toxicity. For the F protein, 3 helper T lymphocytes were found and for G protein, 1 epitope was found.

3.4. Assessment of Molecular Docking of MHC Class I and II Epitopesof the Currently Circualting aMPV Strains with Chicken Alleles

Using the Hdock server, molecular docking analysis was performed between MHC class I and II epitopes and chicken alleles (BF2 for MHC class I and Gaga_BLB1 and Gaga_BLB2 for MHC class II). On the server the chicken alleles were set as the receptors and the predicted epitopes from Table 3 and Table 4 were set as the ligands. In Supplementary Excel File S3 and S4, the binding affinity and confidence scores of the epitopes were listed. Visualization of the docking can be seen in Figure 4 and Figure 5, where H-bond donors can be seen highlighted in magenta, and H-bond acceptors can be seen highlighted in green.

3.5. Vaccine Construct Design

To design the final vaccine construct, the top ranked B cell epitopes and MHC class I and II T cell epitopes were chosen. The epitopes were selected based upon the antigenicity scores, non-allergenicty, non-toxicity, and binding affinity scores when applicable. All epitopes were screened against the chicken amino acid genome for any matches. Of the B cell epitopes, the highest percentage of matching was 43.75%, 75.00% for CTL, and 40% for HTL. This indicated an unlikely chance for the vaccine construct to cause autoimmune interactions. The epitopes were then also matched against the circulating strains of aMPV used in originally to see the individual sequences, subtypes, hosts, and geographical regions, seen in Supplemental Figure S5. In Figure 6, the the mapping of the selected B-cell, MHC class I, and MHC class II are shown in their respective docking locations as well as the solvent accessibility surface (SAS) regions.

Table 5. List of all selected aMPV-B and T cell epitopes incorporated in the final vaccine construct.

S.No	Protein	Start	Peptide	Antigenicity Score	Docking Score	Confidence Score
B Cell Epitopes
1	F	68	STELSLTQNALQELRTVS	0.7314
2		488	NIGAALAITILVVLSMSI	0.9328
3		500	VLSMSIIVGIAYYVVKKR	0.755
4		507	VGIAYYVVKKRKAKAKNG	0.7168
5	G	315	VNPREGGISPTAIQKN	1.627
6		339	DCTLSDPNEPPRICYQVG	0.8212
7		376	YGYACMATLYDTPFNC	0.6155
MHC Class I Epitopes
8	F	96	LSHRKKRFV	1.463	-235.5	0.8468
9		251	VRRKGFGIL	1.5903	-212.17	0.7762
10		252	RRKGFGIL	2.1777	-204.96	0.7501
11	G	48	SICISVEQV	1.4046	-186.24	0.6737
12		49	ICISVEQV	1.7104	-151.07	0.5053
13		50	CISVEQVKL	1.5918	-202.06	0.7391
MHC Class II Epitopes
14	F	238	SGQISLMLNNRAMVR	0.6209	-185	0.6682
15		240	QISLMLNNRAMVRRK	0.6937	-204.61	0.7488
16		241	ISLMLNNRAMVRRKG	0.607	-221.83	0.8079
17	G	20	QVLRRSKRILLGLVL	0.6683	-217.95	0.7956

Table 5. List of all selected aMPV-B and T cell epitopes incorporated in the final vaccine construct.

S.No	Protein	Start	Peptide	Antigenicity Score	Docking Score	Confidence Score
B Cell Epitopes
1	F	68	STELSLTQNALQELRTVS	0.7314
2		488	NIGAALAITILVVLSMSI	0.9328
3		500	VLSMSIIVGIAYYVVKKR	0.755
4		507	VGIAYYVVKKRKAKAKNG	0.7168
5	G	315	VNPREGGISPTAIQKN	1.627
6		339	DCTLSDPNEPPRICYQVG	0.8212
7		376	YGYACMATLYDTPFNC	0.6155
MHC Class I Epitopes
8	F	96	LSHRKKRFV	1.463	-235.5	0.8468
9		251	VRRKGFGIL	1.5903	-212.17	0.7762
10		252	RRKGFGIL	2.1777	-204.96	0.7501
11	G	48	SICISVEQV	1.4046	-186.24	0.6737
12		49	ICISVEQV	1.7104	-151.07	0.5053
13		50	CISVEQVKL	1.5918	-202.06	0.7391
MHC Class II Epitopes
14	F	238	SGQISLMLNNRAMVR	0.6209	-185	0.6682
15		240	QISLMLNNRAMVRRK	0.6937	-204.61	0.7488
16		241	ISLMLNNRAMVRRKG	0.607	-221.83	0.8079
17	G	20	QVLRRSKRILLGLVL	0.6683	-217.95	0.7956

The epitopes were linked using linkers such as KK, AAY, and GPGPC as well as the addition of the full-length chicken IL-18 gene (Accession No. CAB96214) as an adjuvant with a PEAK linker at the C terminal end of the construct. At the C-terminus for purification and identification of the vaccine expression, a 6 X His-tag (HHHHHH) is attached. For the construct as seen in Figure 7 the kozak sequence in orange starts, followed by the start codon in green while B cell epitopes are shown in red with a KK linker, MHC class I T cell epitopes are shown in purple with a AAY linker, and MHC class II T cell epitopes are shown in pink with a GPGPC linker. The MHC class II is linked to the IL-18 adjuvant in brown using HEYGAEALERAG while IL-18 is linked to the 6 x His-Tag in yellow using EAAAK followed by the stop codon in blue. In table 13 all the selected epitopes and their information are listed.

3.6. Disulphide Bond Engineering Results

In the analysis from Disulfide by Design v2.12 webserver showed one residue pair of GLY-ASP observed based upon the chi3 and B-factor energy parameters. The formation of the disulfide bonds was initiated with the introduction of cysteine residues, evaluating within the range of -87 to +97 chi3 values as well as an energy value <2.2. Figure 8 shows the steps to identify the residues for the formation of disulfide bonds to enhance the stability of the vaccine construct.

3.7. Results of the Physiochemical Properties of the aMPV Vaccine Construct

The weight of the vaccine construct is predicted to be 57939.52 Dalton with a theoretical isoelectric point of 9.40, indicative of an alkaline state of the vaccine. In total there are 51 negatively charged residues and 80 positively charged residues. The extinction coefficient measured at 280 nm in water was seen to be 32,540 with the assumption that all pairs of Cys residues form cystines. The aliphatic index was predicted to be 83.72, with a hydrophilic value of -0.200.

Using the iMODS server tool, the stability and physical movement of the vaccine construct was analyzed with normal mode analysis (NMA). In the figure the B factor in 9a is representative of the atomic fluctuations within the normal mode analysis compared to that of the experimental mode of vibration, with values <50 indicating stable proteins and value >50 indicating flexible proteins. In 9b, the variance plot shows the restricted motion through green and purple colors reflective of the cumulative variance. In Figure 9c the eigenvalue of the complex was seen to be 1.688728e-06. Figure 9d is the deformability plot, showing the capability of bending during motion while 9e is the covariance map showing the motions correlations with red indicating correlated movements, blue indicating anti-correlated movements, and white indicating no movement. Lastly, in Figure 9f, the elastic network shows the stable regions with dense springs and flexible region with sparse ones.

3.8. Results of the Secondary and Tertiary Structures of the aMPV Vaccine Construct

The secondary and tertiary structures of the multi-epitope vaccine were analyzed and modeled through PSIPRED and AlphaFold. Evaluation of the quality of the structures were done using ProSA and Ramplot. All of these can be seen in Figure 10.

3.9. Molecular Docking of the designed aMPV vaccine construct with the Chicken Toll Like Receptor (TLR7) Analysis

Molecular docking was performed between the aMPV construct and a toll-like receptor (TLR7) to analyze the potential immunogenic ability of the multiepitope vaccine. As previously discussed, TLR7 was selected from the ten toll-like receptors due to its role as a major immune sensor, utilizing various pathways to trigger the innate immune response. Using the Uniprot database, the sequence of TLR7 was retrieved and modeled using AlphaFold. The docking analysis was done using Hdock, where the results indicate a strong binding affinity between the vaccine construct and TLR7. Figure 11 shows the visualization of the construct binding as well as the docking interaction analysis of the amino acid residues and formation of hydrogen bonds and hydrophobic interactions as predicted by the PDBsum server in Tamarind Bio.

3.10. In Silico Cloning and Immune Simulation of the Designed aMPV Multi-Epitope DNA Vaccine

Shown in Figure 12 is the predicted immune response of the vaccine as seen through the B cell, T cell, and cytokine interactions. For cloning of the vaccine construct, the amino acid sequence was decoded using VectorBuilder for the expression within expression vector pcDNA3.1(+). The vector was cloned on SnapGene between the restriction enzyme cutting locations of BamHI and EcoRI, as seen in Figure 13.

4. Discussion

Avian metapneumovirus is a major emerging respiratory virus, with significant impact on the global poultry industry, with morbidity capable of affecting up to 10% of the population [1]. Many of the current vaccine for aMPV reduce clinical conditions and viral dissemination, but not all are equally effective in the prevention of infection and circulation, allowing for continued viral persistence of the virus [2]. This scenario triggers an urgent need for the development of effective vaccines that can cover currently circulating subtypes across the globe. As described in previously published articles, the use of machine learning has greatly improved, with multiepitope vaccines showing a long-term protection through induced trained immunity, but by also being safer and cost-effective [9,10,11,18,19].

In this study, the aim was to use machine learning immunoinformatic approach to develop a multiepitope based DNA vaccine, targeting two important proteins for avian metapneumovirus, F and G proteins. Using 53 isolates for the F protein and 68 isolates for the G protein, the vaccine construct targets the circulating isolates across the globe to provide maximum coverage. The vaccine construct was developed, showing promising physiochemical properties as well as high antigenicity and a tertiary structure suitable for further work and analysis. Using immune simulation, the designed aMPV vaccine construct was able to demonstrate a high potential for propagating humoral and cellular immune responses in silico. The data reflects as sequence conservation, seeming to have a potential advantage to be test in vivo. Our results showed that this vaccine construct potentially induce robust antibody production with class switching as well as a pertinent cytotoxic T-cell response, lasting help T-cell response, and immune activation seen by the cytokine profile.

The pipeline of the design and in-slico cloning of the multiepitope DNA vaccine are as follow : (1) retrieval of the amino acid sequences, (2) multiple sequence alignment generating a consensus sequence, (3) prediction of B cell (linear and discontinuous) and T cell epitopes (MHC class I and II), (4) molecular docking of T cell epitopes, (5) selection of top ranked epitopes and screening against the chicken amino acid genome, (6) construction of the multiepitope vaccine using linkers, (7) the use of IL-18 as adjuvant to the construct, (8) prediction of the physiochemical properties, (9) prediction of the secondary and tertiary structures, (10) molecular docking of the vaccine construct to chicken Toll-like receptor (TLR7) (11) immune simulation of the final vaccine construct and in silico cloning to assess the ability of the construct to elicit humoral and cellular immunity in chickens. For both B and T cell, predictions yield numerous epitopes per protein. Criteria filters were set to make sure the best epitopes were selected. For all the epitopes predicted, they were screened for allergenicity, antigenicity, and non-toxicity [64,65,66]. For B cell epitopes, the predicted epitopes were screened against the conformational epitopes predicted by ElliPro, in which overlap in residues indicates that the regions include properties from both linear and conformational B cell epitopes that is further carried to the multiepitope vaccine construct. When it comes to the MHC II epitopes, a percentile score of < 4 and IC50 < 50 nM was set as the criteria for filtering. Molecular docking was performed on all the MHC based epitopes, in which any that could not dock when using the HDock server were filtered.

After the filtering, top ranked peptides were selected based upon binding and high antigenic scores. Through the analysis of the conservation of the epitopes, it was seen the F protein had 53 isolates, there are 18 subtype A, 28 subtype B, 6 subtype C, and 1 mutant isolate, spanning over 6 avian hosts and five countries. For the G protein, with 68 isolates, there are 15 subtype A, 49 subtype B, 3 subtype C, and 1 mutant isolate, spanning over 5 avian hosts and 5 countries. It was then identified that two epitopes (RRKGFGIL and VRRKGFGIL) within the F protein have a 100% sequence identity across the 53 isolates. This shows the conservation that bridges the three subtypes and mutant, throughout all 6 avian species and 5 countries listed. It is important to note that epitope VGIAYYVVKKRKAKAKNG maintains a 94.44% sequence identity across 85. 71% of the subtype B isolates, based upon a single lysine to proline difference. On the other hand, the G protein epitopes show an exclusive specificity towards subtype B due to the high hypervariability of the G protein. Three epitopes ICISVEQV, SICISVEQV, and CISVEQVKL have a 100% sequence identity to subtype B isolates used, showing exceptional coverage. Of the targeted lineages, the used epitopes in the vaccine have crossed the host species barrier covering commercial poultry lines such as chickens and turkeys across international boundaries, based upon the computational analysis, to be further verified by in vivo testing. This conservation could theoretically minimize the risk of strain escape, in which further testing would be needed. The epitopes were constructed into the vaccine using linkers and adjuvants. For this vaccine, chicken IL-18 was selected as the adjuvant. It was selected based upon its importance as a cytokine in the innate and adaptive immune immunity, as well as it being used in various vaccine constructs against other viruses such as IBDV, NDV, and H5 influenza [53,54,67,68]. For preservation of conformational epitopes conformation, flexible linkers were used like GPGPG and EAAK. The 3D structure prediction was completed with AlphaFold, and a Ramachandran plot analysis along with other structural analysis were done to confirm proper folding and that the construct is solvent accessible.

A challenge presented in the study is that chicken alleles are not available in the IEDB server for MHC class I and II prediction. Alternative strategies were used to overcome this issue, using things such as using published binding motifs for MHC class I epitopes [45,46]. Specifically, MHC I has very little computational tools available in which using human surrogate models was not possible. For MHC class II, the human surrogate model was able to be used, leaving the default settings for the predictions. With this prediction, MixMHC2pred was able to be used to predict the corresponding chicken alleles (Gaga_BLB1_002_01, Gaga_BLB1_012_01, Gaga_BLB2_002_01, Gaga_BLB2_012_01, Gaga_BLB2_012_02) [50].

To predict the presences of cysteine residues containing a reactive thiol (-SH) group in the side chain that when incorporated into the protein would have major biochemical impact, a disulfide bond analysis was done. Disulfide bonds normally form due to the oxidation reaction between two cysteine residues, making the vaccine structure more stable. The spatial orientations of the two cysteines linked are determined by the conformational classification of the X³ dihedral angle (C-S-S-C), describing the rotation around the S-S bond [69]. At an angle approximately ±90◦, bonds show either the right- or left-handed conformations, and when approaching 180◦, it corresponds to a trans conformation, showing an arrangement of cysteines that is extended and linear. When the angle is close to 0◦, it indicated the cis conformation, which is energetically unstable due to its steric clashes and unfavorable geometry. As the resulting value was -80.08◦, the protein structure was making the right- or left-handed conformations with their corresponding energy values.

Normal mode analysis (NMA) was used to further validate the stability and flexibility of the vaccine construct. In the mobility plot, residues with higher peaks (0.5-2 arbitrary units) correspond to flexible loops or linker regions, while lower values (0-0.3 arbitrary units) correspond to stable structural regions. Bending during motion is shown in the deformability plot in the range of 0.05-0.15. The B-factor plot assesses atomic fluctuations where values below 50 indicate stable proteins and values above 50 indicate flexible proteins. Eigenvalues in the range of 1 × 10^-6 to 1 × 10^-4 are representative of stiffness of a protein, in which the vaccine protein at the higher end of this range shows a decent structural rigidity, and decent amount of energy would be needed for deformation. The variance plot shows individual (purple) and cumulative (green) variances, with lower individual variance values (<0.05) representing restricted motion of each mode, while the accumulation of modes was represented by a rise in the cumulative curve. Lastly, the stable regions seen with dense springs and flexible areas with sparse springs make up the elastic network.

In another study, a multiepitope vaccine candidate for infectious bronchitis virus in silico exhibited a strong binding to both TLR3 and TLR7 [70]. Following this, TLR7 was used for molecular docking. TLR3 was not used as it recognizes double stranded RNA while aMPV is single stranded. TLR7 is a primary endosomal receptor targeting single stranded RNA viruses like aMPV. The receptor is present on various cell types such as plasmacytoid dendritic cells, B cells, eosinophils, neutrophils, and monocytes/macrophages, making it a very crucial target for therapeutic research such as vaccine development. Previous studies have been published on similar multiepitope vaccines using TLR7 docking for viruses such as Newcastle Disease Virus and Avian Leukosis Virus, where string binding and likely robust interactions were seen in the proposed virus and TLR7. The interaction between the vaccine and TLR7 can lead to humoral and cellular immune response by stimulating the above-mentioned cells, promoting the production of cytokines, and accelerating the mobilization of effector cells. The docking simulation to TLR7 exhibited strong binding affinity between the receptor and the vaccine construct with a docking score of -264.26.

The last part of the study included an in silico immune simulation of the vaccine construct on the C-ImmSim server. The C-ImmSim is a very important server that provides predictive insights into the immunogenicity of the vaccine construct. The server also is a major limitation in this study as the server is optimized for human (mammalian) immune parameters. Distinct anatomical and physiological immune profiles such as the absence of encapsulated lymph nodes, unique splenic microenvironments, and highly specialized major histocompatibility complex profiles are seen in avian species. As the C-ImmSim server cannot replicate these features, the results are not a guarantee the vaccine can provide absolute protective immunity for avian species. These results are in-silico, pointing towards rather a proof of the concept, confirmation of a theoretical processing capacity, in which this can be treated as a baseline for needed in vivo validation in avian models.

The simulation shows that T cell populations such as cytotoxic T cells and helper T cells were activated, all important for humoral and cellular immunity. This simulation shows that the construct would be a good candidate for cloning as done using VectorBuilder and SnapGene. From this analysis the vaccine construct has the potential to elicit humoral and cellular immunity, confirmed through further studies with in vivo validation in avian models, which could change the role aMPV has, protecting chickens and other birds from infection. The early generations of the DNA vaccines have important disadvantages which is the low level of vaccine delivery to the target cells which impact the production of the immune response in the host after vaccination [71,72,73,74]. The new generations of the DNA vaccines have several modifications to overcome these problems including the encapsulation with some nanoparticles and the conjugation of some immune stimulates particularly the chicken IL-18 [20,54,75]. To improve the efficacy of this designed multiepitope DNA vaccine against aMPV, this vaccine construct should be encapsulated with some nanoparticle. This encapsulation will enhance the delivery of the vaccine into the target cells which will potentiate the immune response. This approach will overcome the low expression problem of the early generation of the DNA vaccines or those DNA vaccines that lack any nanoparticle [71,72,73,74]. Some recent studies suggesting incorporation of IL18 in the DNA vaccine construct against Newcastle diseases virus potentiate the humoral and cell mediated immunity suggesting its positive impacts as vaccine adjuvant [20,54,75]. There are many unique criteria of the aMPV vaccine construct designed in this study including (1) this vaccine construct potentially offers protection against aMPV based upon the observed conservation, (2) the vaccine construct will be encapsulated with nanoparticle which enhance the delivery of the vaccine into the target cells, (3) the vaccine include the chicken IL18 as an adjuvant which further enhance the immune response of the vaccinated host, (4) this vaccine could potentially protect multiple species of birds particularly chickens and turkeys which offer multispecies protection based upon conservation, to be validated by further testing, (5) the possibility of updating this vaccines to match the circuiting strains in case of any changes happened to the virus on the genomic level (6) the potential mass vaccination of this designed vaccine including in ovo administration. All these criteria make the designed new vaccine against aMPV an attractive vaccine for the poultry industry not only in the US but worldwide.

For future studies, the vaccine construct should be experimentally validated, homing in on studying the protein stability, safety and initial efficacy. Following that, further studies into in vivo evaluation can show the vaccine construct in a more complex system with a harder look at the safety and effectiveness of the immune response. Dosing can then be tested to see how the product does against various strains of aMPV, allowing for testing in animals such as chickens, turkeys, and ducks. In the event of more isolates, the construct can be updated using immunoinformatics to streamline all changes including, but not limited to, new epitopes and property testing.

5. Conclusions

A multiepitope DNA based vaccine for aMPV was successfully designed, utilizing 17 top ranked epitopes based upon, immunogenicity, non-allergenicity, and non-toxicity. The epitopes included B cell, MHC class I, and MHC class II from major surface proteins in the virus (F and G). High binding affinity was seen in the T cell epitopes to chicken alleles. The selected epitopes were linked together and had chicken IL-18 attached as an adjuvant. High binding affinity was also seen between the vaccine construct and Toll-like receptor 7. The vaccine was successfully cloned in-silico and when undergoing immune simulation, showed high immunogenic potential in both humoral and cell-mediated immunity. The vaccine construct is believed to provide possible protection various from circulating strains of aMPV. It is important to note the design of the vaccine is based upon computational prediction tools when assessing the high potential of the vaccine in terms of immunogenicity and structural behaviors. These characteristics are bound to vary in biological systems due to various factors such as the role in the cell, processing of the proteins, and modifications made post translation. Future work is proposed to validate the multiepitope vaccine construct in vitro, working towards animal models for full usage.