Development of a novel multi-epitope mRNA vaccine candidate to combat HMPV virus

ABSTRACT

Human metapneumovirus (HMPV) is one of the main pathogens causing severe respiratory infections in children, as a common cause of immunodeficiency-related deaths in children and elderly individuals, the prevalence of HMPV has been showing an increasing trend during the last years. However, no vaccines or effective treatment plans are available currently. In this present, based on candidate proteins highly associated with viral virulence and has promising protective potential, we screened for immunodominant cytotoxic T cells, helper T cells, and Linear B-cell epitopes from the most promising candidate Fusion protein, together with G, SH, M, and M2. All epitopes were predicted to have strong antigenicity by Vaxijen and pose no potential toxicity, allergenicity, or hormonology to human proteins by Toxinpred, Allerpred, and Blast analysis, meanwhile, high conservancy is demanded to cover different subtypes. adjuvants β-defensin II and Pam2Cys was attached with EAAAK linkers to improve vaccine’s efficiency. Then, calculation of physicochemical properties proved the protein vaccine as a product can stably exist in the human body. Besides, we assessed the docking between the vaccine and immune receptors to evaluate its ability to stimulate immune responses, and the dynamic simulation further confirmed that the vaccine can tightly bind with immune receptors, which approved that the construction has the potential to induce strong humoral and cellular immune response. Finally, the vaccine was constructed into a multi-epitope mRNA vaccine, the immune simulations suggest that this is a vaccine candidate for controlling HMPV infection.

KEYWORDS:

• HMPV
• fusion protein
• epitope docking
• molecular dynamics
• multi-epitope vaccine
• mRNA vaccine

View responses to this article:

Multi-epitope mRNA vaccine candidate to combat HMPV virus: Comment

Introduction

Human metapneumovirus (HMPV), first identified in 2001 from 28 young children from the Netherlands, was recognized as a leading cause of Acute respiratory tract infections (RTIs) in children, the elderly and other population with immunodeficiency.¹ Usually, as a seasonal epidemic, HMPV threaten is controllable, its circulation lasts from winter to next year’s spring along with other respiratory viruses and generates a limited number of cases, the average detection rate of individuals under the age of 18 between 2016 and 2021 was approximately 4%.² However, the outbreak of HMPV exceeded previous years in 2023, the number of HMPV antigen-positive patients reported by hospitals and clinical laboratories from various states across the United States reached 19.6% in March (average 3.6% from 2008 to 2014 from same organization), resulting in overcrowding in pediatric and intensive care units.^3,4

HMPV belongs to the Pneumoviridae Metapneumovirus family. With a disease spectrum similar to respiratory syncytial virus (RSV), HMPV is responsible for about 15% acute respiratory infections among young children and about 5% in adults.^5,6 According to the World Health Organization’s report, HMPV is responsible for 10% of deaths in children under the age of 5 in 2017. While in 2018, it is estimated that the number of children under the age of 5 infected with HMPV worldwide reaches 14 million, among which 640,000 require hospitalization and 8,000 result in death.⁷ Typically, HMPV infection results in mild symptoms including fever, cough, nasal congestion and so on, but individuals with compromised immune function, such as premature infants, when infected, may develop bronchiolitis, pneumonia, and even severe pneumonia requiring admission to the ICU. In general, HMPV causes a significant medical burden every year, especially in middle- and low-income countries. Thus, it is highly necessary to have a targeted intervention strategy to protect infants and young children. Unfortunately, there is no solid treatment or vaccine for HMPV so far, the lack of vaccines is a significant factor leading to a large number of hospitalizations during the epidemic season, while the scarcity of treatment options has resulted in cases of death among hospitalized patients.

Supportive treatments such as immunoglobulins, glucocorticoids, and other symptomatic treatment methods are considered the main method against HMPV infection. In clinical practice, ribavirin and polyclonal intravenous immune globulin (IVIG) are utilized for severe infections in cases, though there is currently no study exploring their solid effectiveness in patients and further RCT studies are needed, these therapies are still indicated based on experiences.⁸ There are also some methods that have demonstrated antivirus effects in mice, including monoclonal antibodies, RNAi, NMSO3, and heparin, etc.⁹ However, none of these methods have been approved for licensed clinical use.

Developing a vaccine seems to be a feasible way to alleviate the pressure of viral infection in the absence of effective treatments. Human metapneumovirus (HMPV) has evolved into two genetic lineages with antigenic distinctiveness, A and B. Each lineage contains two subtypes, A1 and A2, B1 and B2, A2 can be further divided to A2a and A2b.¹⁰ Different subtypes have distinct genetic sequences, but still demonstrate considerable homology, as evidenced by protein sequence similarities between different subtypes exceeding 80%，by selecting the epitopes that are conserved in both lineages, subsequent vaccine candidates can provide protection against various strains of HMPV.¹¹ The HMPV virus carries a non-segmented negative-sense RNA as its genome, which has a total length of 13kb. There are eight genes on this RNA, collectively encoding nine proteins, among which nucleoprotein (N), phosphoprotein (P), and large protein (L) participates in synthetization of RNA. And matrix protein (M) is involved in viral assembly and budding, while M2 is involved in RNA transcription. The rest three proteins are anchored to the viral lipid envelope, F protein mediates membrane fusion while G protein mediates attachment. Less is known about small hydrophobic (SH) protein, studies suggest SH protein may be a viroporin forming ion channels, but it is believed to regulate host immune response.^12,13 After many attenuated vaccines failed due to insufficient attenuation or poor protective effects, subunit vaccine became the main focus of research. The F protein plays a crucial role in the fusion of the virus envelope and cell membrane, which means that targeting the F protein can terminate the invasion and subsequent replication of the virus. Moreover, the F protein is highly conserved (more than 90%) among different strains of HMPV, which make F protein the most important target in vaccine development.^14–16 Meanwhile, though other proteins are not considered as protective vaccine candidates before, they are also believed to be involved in immune progress, their epitopes have been identified to participate in immune regulation, thus improving immune responses against HMPV infection, the G and SH protein are also considered weaker protective vaccine provides neutralizing antibodies.^17,18 Considering the importance of F protein in immune response of HMPV, we set F protein as a main target for further epitope vaccine design, at the same time, given that G, SH, M, and M2 proteins are also highly related to viral replication, some immunodominant epitopes could be screened from their sequences to enhance the vaccine construction.

In total, there has been little progress made in HMPV vaccine development, a multi-epitope mRNA vaccine is a proven strategy against bacterial and viral infection, an mRNA vaccine uses mRNA as a carrier to encode exogenous antibodies and introduce them into body cells. By utilizing the expression system to synthesize antigens, they induce cellular and humoral immune responses. Multi-epitope mRNA vaccines select immunodominant epitopes from different proteins as exogenous antibodies, thereby further enhancing the specificity and protective ability of the vaccine. With immunoinformatic methods, we designed a multi-epitope mRNA vaccine based on the F protein as well as some epitopes from G, M1, M2, and SH proteins to address the potential threat of an HMPV epidemic,

Method and materials

Reference sequence acquisition

The NCBI virus database stores the gene sequence and protein sequence of each virus, and the reference protein sequences of human metapneumovirus are obtained from this database. The G protein (YP_009513272.1), F protein (YP_009513268.1), M protein (YP_009513267.1), M2–1 protein (YP_009513269.1), M2–2 protein (YP_009513270.1), and SH protein (YP_009513271.1) were selected for further analysis.

T cell epitope prediction

The NetCTL1.2 is a server，comprehensively predicts CTL epitopes by means of peptide-MHCI binding, C-terminal cleavage of the proteasome, and TAP transport efficiency.¹⁹ Twelve supertypes, including A1, A2, A3, A24, A26, B7, B8, B27, B39, B44, B58, and B62 are selected, and the default threshold is 0.75. We used the MHC-II Binding Predictions tool from IEDB to predict the HTL epitope via the IEDB recommendation algorithm (NetMHCIIpan 4.1 EL).²⁰ Vaxijen 2.0 server was employed to predict the antigenicity of the epitopes according to their physical and chemical properties of amino acids.²¹ Epitopes with antigenicity scores >0.4 were used for further screening. Moreover, the AllerTOP v2.0 was utilized to predict the allergenicity of epitopes, and the ToxinPred servers were utilized to predict the toxicity of epitopes.^22,23

For CTL epitopes, their affinities to MHC I molecules were further predicted by the MHC-I Binding Predictions tool based on the NetMHCpan EL 4.1 algorithm. In general, CTL epitopes with rank percentiles <0.5% were considered to have strong binding to the MHC I molecules.

For HTL epitopes, epitopes with rank percentiles <2% were considered to have strong binding to the MHC II molecules. Cytokines induced by HTL epitopes are important in adaptive immunity. IFNepitope, IL4Pred and IL10Pred servers were employed to analyze the ability of HTL epitopes to induce IFN-γ, IL-4, and IL10, respectively.^24,25

Prediction of linear B cell epitopes

The ABCpred server predicted linear B cell epitopes with a threshold of 0.51 and an accuracy of 65.93%.²⁶ Then, the antigenicity, allergenicity, and toxicity of the LBL epitopes were predicted through the VaxiJen 2.0 server, the AllerTOP v2.0 server, and the ToxinPred server. LBL epitopes with good antigenicity, non-toxicity, and non-allergenicity were used for further analysis.

Conservative analysis of epitopes and homology analysis with human protome

The complete genomes of 392 human metapneumovirus strains isolated from human were obtained from the NCBI virus database (https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/), and the human metapneumovirus genome database was established via BLAST v2.120+. The tBLASTn was used to query the conservatism of epitopes, and more than 80% of the epitopes are used for vaccine construction. PeptideMatch server was employed to identify homology between epitopes and human proteins, thus minimizing the risk of self-antigens immune reactions.²⁷

Vaccine construct and calculation of physicochemical properties

The final vaccine consists of HTL epitopes, CTL epitopes, LBL epitopes, and adjuvants. AAY linkers, GPGPG linkers, and KK linkers were used to connect the CTL epitopes, HTL epitopes, and LBL epitopes, respectively. In order to enhance the immunogenicity of the construction, TLR4 agonist and TLR2 agonist were attached to the N-terminal and C-terminal of the vaccine, respectively. Further, A pan HLA DR-binding epitope (PADRE) that can bind to multiple MHC II molecules was incorporated into the construction to enhance the adaptive immunity induced by the vaccine. To figure out the physicochemical characteristics of the vaccine, the ProtParam tool was employed. The ProtParam tool can calculate the number of amino acids, theoretical PI, instability index, aliphatic index, half-life, and hydrophilicity values of the vaccine. Moreover, the VaxiJen 2.0 server and the ANTIGENpro server were utilized to evaluate the antigenicity of the vaccine. The AllerTOP v2.0 server and the ToxinPred 2.0 server were employed to evaluate the allergenicity and toxicity of the vaccine. To understand the probability of solubility of the vaccine when overexpressed in E. coli, the SOLpro server was used (http://scratch.proteomics.ics.uci.edu/).

Secondary structure, tertiary structure prediction, and refinement

The vaccine’s secondary structure was predicted with GOR4. The tertiary structure model of the vaccine was constructed by the Robetta server (http://robetta.bakerlab.org/). GalaxyWEB (http://sysbio.rnet.missouri.edu/3Drefine/) was then employed to refine the tertiary model. Three main indicators, namely GDT-HA, RMSD, and MolProbity were utilized to refine the tertiary model, which can reveal the model’s quality.²⁸ In addition, the quality of the model can be evaluated with Procheck, ERRAT, and Pro-sa Web servers. A model is considered to be of high quality if it has an ERRAT score above 85. The overall quality of the protein’s tertiary structure can be assessed using the Z-score, with a positive value indicating structural errors.^29,30

Molecular docking

Toll-like receptors (TLRs) play a crucial role in the recognition of different pathogen-related molecular patterns thereby contributing significantly to the innate immune response. The Crystal structures of TLR2 (PDBID:2Z81) and TLR4 (PDBID:2Z63) were obtained from the RCSB PDB Database. The ClusPro 2.0 server performed molecular docking between the vaccine and the TLRs, which used the PIPER16 program based on the on the Fast Fourier Transform (FFT) correlation approach to perform rigid-body docking. Further, the top rank result was submitted to the HADDOCK 2.4 server for refinement using default parameters. Meanwhile, we used the PDBsum server to analyze the interaction between vaccine and TLRs.³¹

Molecular dynamics simulation

Molecular dynamics (MD) simulation evaluated the stability of vaccine-receptor complexes. Then, a 100 ns molecular dynamics simulation was conducted by using Gromacs v2022.1 software.³² The Amber14SB_ armbsc1 force field and TIP3P water model were added to the whole system. Then, the anion (Cl-) was added to balance the charges of the whole system. Subsequently, energy minimization was performed based on steep algorithm until the maximum force ≤1000.0 kJ/mol/nm. NVT and NPT were performed 200ps and 1 ns, respectively. The modified Berendsen thermostat was employed to regulate the system’s temperature to 310k, while the Parinello-Rahman barostat was used to equilibrate the pressure to 1 atm. The hydrogen atom bonds were constrained by the LINCS algorithm and long-range electrostatic interactions were calculated by the Particle-Mesh Ewald summation scheme.

MM-GBSA calculation

The gmx_MMPBSA v1.56 tool calculated the binding energy between the vaccine and the receptors. In this study, we used the Molecular Mechanics Generalized Born Surface Area (MM-GBSA) methods to calculate the binding energy between the vaccine and the receptors via the last 10 ns trajectory.

Immune simulation and population coverage

Immune simulation can be used to evaluate the immune response of the human body after three injections of vaccine. C-IMMSIM server was employed to perform a progress of immune simulation, which identified epitopes via the PSSM methods.³³ Vaccination was injected on days 1, 28, and 56, respectively, and the simulation was performed for a duration of 350 d. The host HLA selection was set to HLA-A * 02:01, HLA-B * 07:02, and HLA-DRB1 * 01:01. The coverage of the vaccine in the world was computed with the Tepitool according to MHC corresponding of epitopes.³⁴

Construction of multi-epitope mRNA vaccine and secondary structure prediction

Jcat server (http://www.jcat.de/Start.jsp.) was utilized for reverse translation and optimization of the vaccine. To construct the conclusive multi-epitope mRNA vaccine, the DNA sequence of the vaccine is further modified, a Kozak sequence is attached to the 5’ end to improve RNA stability and translation efficiency,³⁵ signal peptide Tissue Plasminogen Activator (tPA) follows the Kozak sequence to facilitate the extracellular secretion of proteins and improve antigen presentation.³⁶ Meanwhile, an MHC I-targeting domain (MITD) sequence is attached to the 3’ end of the vaccine DNA sequence to guide the presentation of CTL epitopes,³⁷ then a TAA condon ends the translation. Besides, to stabilize the mRNA partial sequence of cytomegalovirus immediate-early gene Untranslated Region (UTR) was added as 5’ untranslated region (UTR), and partial sequence of the human growth hormone was attached as 3’ UTR.^38,39 The full DNA sequence was translated into mRNA using the Transcription Tool. The RNAfold server (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi) and mFold server (http://www.bioinfo.rpi.edu/applications/mfold) predicted the mRNA’s secondary structure.

Construction of vaccine vectors

XhoI and BamHI were selected from E. coli K12 for reverse translation and optimization of the vaccine sequence, gensmart tool was used for expression vector construction. XhoI and BamHI inserted the vaccine gene sequence into the pET28a (+) plasmid to construct expression vector⁴⁰(Figure 1).

Figure 1. Graphical abstract.

Results

Screening of epitopes

After being obtained from NCBI database, sequences of F, G, SH, M, M2–1, and M2–2 were put into the process of screening, Overall, all epitopes share some common screening principles, these epitopes are supposed to possess strong antigenicity (antigenicity >0.4), and proper immunogenicity, no toxicity, allergenicity, or homology to human proteins is allowed, besides, to provide protection against all subtypes, selected epitopes should be able to cover the majority of strains in 392 strains obtained from NCBI database. Then, epitopes that are deemed suitable and have high overall scores and rankings in the prediction are selected for further computation. Finally, 15 CTL epitopes (Table 1), 4 HTL epitopes (Table 2) and 9 LBL epitopes (Table 3) were selected. Most epitopes are from F and M proteins, none were from G proteins. In addition, the ability of inducing cytokines containing IFN-γ, Il-4, and Il-10 was considered when selecting HTL epitopes.

Table 1. Screen of CTL epitopes.

Protein	Start	Sequence	Antigenicity	Immunogenicity	Toxic	Allergic	Conservation	Score	MHC alle	Rank
F	30	TITEGYLSV	0.5971	0.03685	-	-	353/392	0.6504	HLA-A *02:06	0.13
								0.5281	HLA-A *68:02	0.17
								0.4682	HLA-A *02:03	0.27
								0.4316	HLA-A *02:01	0.33
	35	YLSVLRTGW	0.6288	0.07832	-	-	354/392	0.5733	HLA-B *57:01	0.44
								0.4501	HLA-B *58:01	0.37
								0.4396	HLA-A *32:01	0.14
	157	VLATAVREL	0.4093	0.2309	-	-	353/392	0.7337	HLA-A *02:03	0.08
								0.5376	HLA-A *02:01	0.24
	242	KLMLENRAM	0.6227	0.08979	-	-	356/392	0.3657	HLA-B *15:01	0.48
	243	LMLENRAMV	1.0043	0.06999	-	-	356/392	0.3160	HLA-A *02:03	0.48
	276	FGVIDTPCW	0.8004	0.1646	-	-	323/392	0.5519	HLA-B *58:01	0.29
								0.2783	HLA-B *53:01	0.31
	351	NINISTTNY	1.3653	0.03624	-	-	355/392	0.4902	HLA-A *30:02	0.18
								0.4718	HLA-B *15:01	0.33
								0.3946	HLA-A *26:01	0.18
								0.3250	HLA-A *01:01	0.34
								0.3192	HLA-B *35:01	0.4
	417	TVTIDNTVY	0.4045	0.21891	-	-	355/392	0.5565	HLA-B *35:01	0.2
								0.2844	HLA-A *30:02	0.49
								0.2627	HLA-A *01:01	0.43
								0.1722	HLA-A *26:01	0.43
	484	AEKGNTGFI	1.0937	0.09134	-	-	352/392	0.1890	HLA-B *44:02	0.5
	519	KPTGAPPEL	0.5202	0.1235	-	-	329/392	0.9611	HLA-B *07:02	0.03
								0.5289	HLA-B *35:01	0.22
								0.3580	HLA-B *53:01	0.23
M	12	IPYTAAVQV	0.7186	0.07995	-	-	350/392	0.9802	HLA-B *51:01	0.01
								0.4704	HLA-B *07:02	0.27
								0.3026	HLA-B *35:01	0.41
								0.2607	HLA-B *53:01	0.34
	53	KTLTITTLY	0.5848	0.22788	-	-	349/392	0.9669	HLA-B *57:01	0.04
								0.9355	HLA-A *30:02	0.01
								0.9354	HLA-B *58:01	0.04
								0.8030	HLA-A *11:01	0.08
								0.7408	HLA-A *32:01	0.03
								0.7183	HLA-A *01:01	0.09
								0.7029	HLA-A *03:01	0.16
								0.5264	HLA-B *15:01	0.27
								0.4308	HLA-A *30:01	0.21
								0.1982	HLA-A *26:01	0.38
	195	APYAGLIMI	0.4935	0.07045	-	-	350/392	0.7280	HLA-B *51:01	0.06
								0.5614	HLA-B *07:02	0.2
								0.1767	HLA-B *53:01	0.49
M2–1	22	KFNHNYWSW	0.9138	0.11095	-	-	349/392	0.7715	HLA-A *24:02	0.06
								0.7090	HLA-A *23:01	0.08
								0.5295	HLA-B *57:01	0.5
								0.3855	HLA-A *32:01	0.16
								0.3775	HLA-B *58:01	0.43
M2–2	57	YLENIEIIY	0.9347	0.43992	-	-	346/392	0.8639	HLA-A *01:01	0.05
								0.3091	HLA-A *30:02	0.42
								0.3080	HLA-B *35:01	0.41

Table 2. Screen of HTL epitopes.

Protein	start	Sequence	core peptide	Antigenicity	Immunogenicity	Toxicity	Allergic	IFN-γ inducer	IL-4 inducer	IL-10 inducer	Conservation	score	MHC alle
F	40	RTGWYTNV FTLEVGD	WYTNVFTLE	0.4378	0.19130959	-	-	+	+	+	355/392	0.2649	HLA-DPA1 *01:03/DPB1 *04:01
												0.1692	HLA-DPA1 *02:01/DPB1 *01:01
												0.2328	HLA-DPA1 *03:01/DPB1 *04:02
M	3	SYLVDTYQ GIPYTAA	VDTYQGIPY	0.4863	0.04596061	-	-	+	+	+	348/392	0.5837	HLA-DRB1 *15:01
	6	VDTYQGIP YTAAVQV	YQGIPYTAA	0.4467	0.17494668	-	-	+	+	+	348/392	0.2506	HLA-DPA1 *01:03/DPB1 *04:01
												0.2752	HLA-DPA1 × 03:01/DPB1 × 04:02
												0.1638	HLA-DPA1 × 02:01/DPB1 × 01:01
												0.4605	HLA-DPA1 × 01:03/DPB1 × 02:01
	184	EADQALTQ AKIAPYA	QALTQAKIA	0.8367	0.025362634	-	-	+	+	+	350/392	0.4803	HLA-DQA1 × 01:02/DQB1 × 06:02

Table 3. Screen of LBL epitopes.

Protein	Start	Sequence	Score	Antigenicity	Toxic	Allergic	Conservation
F	9	FSLLITPQHGLKESYL	0.74	0.5875	-	-	348/392
	39	LRTGWYTNVFTLEVGD	0.85	0.5023	-	-	355/392
	188	KMAVSFSQFNRRFLNV	0.74	0.7057	-	-	355/392
	206	QFSDNAGITPAISLDL	0.71	0.6216	-	-	352/392
	351	NINISTTNYPCKVSTG	0.84	0.9377	-	-	355/392
	433	EQHVIKGRPVSSSFDP	0.78	0.9249	-	-	332/392
M	134	KSVGKKTHDLIALCDF	0.73	0.6998	-	-	347/392
	182	SSEADQALTQAKIAPY	0.73	0.7844	-	-	350/392
	193	KIAPYAGLIMIMTMNN	0.84	0.6404	-	-	348/392

Vaccine construction and physicochemical properties calculate

To improve the effectiveness of the vaccine construct, β-defensin II and Pam2Cys were introduced to the vaccine sequence as adjuvants. Additionally, the vaccine also contains a sequence of PADRE to supplement the helper T cell (HTL) epitopes. The final vaccine consists of 2 adjuvants (Pam2Cys and β-defensin II), 15 CTL, 5 HTL (including PADRE), and 9 LBL epitopes. The vaccine’s final antigenicity is 0.6487, nontoxic and non-allergenic. The construction consists of 507 aas and has a 55.2 kDa molecular weight, its theoretical pI is predicted as 9.41. The half-life of the vaccine is estimated to be 30 hours in mammalian reticulocytes, in vitro, more than 20 hours in yeast, and more than 10 hours in E. coli. The protein’s instability index (II) is computed to be 22.50 which suggests the protein is stable. Furthermore, it’s 81.56 aliphatic index indicates that the vaccine has good thermal stability. The grand average of hydropathicity (GRAVY) for the protein was −0.051, indicating that it is hydrophilic as it is less than 0. The solubilization probability of the vaccine when overexpressed in E. coli was predicted to be 0.7936 according to SOLpro.

Secondary and tertiary structures prediction

The GOR4 server is utilized to predict the construction’s secondary structure, the sequence consists of 38.86% alpha helix, 18.15% extended strand, and 43% random coil (Figure 2). The vaccine’s tertiary structure was predicted using the RoseTTAfold, and Galaxyweb was then utilized to refine the structure (Figure 3a). The best model achieved a Molprobity score of 2.071, GDT-HA of 0.9763, and RMSD of 0.343. Ramachandran plots reveal that 375 residues (85.4%) lie in the most favored regions, with only 55 residues (12.5%) in the additional allowed regions, 4 residues (0.9%) in the generously allowed regions and 5 residues (1.1%) in the disallowed regions (Figure 3c). The Pro-sa web server analyzed the structure of the model, which achieved a Z-score of −7.6 (<0), indicating good quality (Figure 3d).

Figure 2. Prediction of secondary structure of vaccine.

(a) h for alpha helix, e for extended strand, c for random coil (b) percentage of second structure in vaccine sequence (c) blue for alpha helix, red for extended strand and purple for random coil.

Figure 3. Construction of multi-epitope vaccine.

(a) 3D structure of the vaccine (b) Protein sequence of the vaccine construction, β-defensin is lighted in red, PADRE sequence in yellow, HTL epitopes are lighted in blue, CTL epitopes in green and B epitopes are lighted in purple, Pam2Cys sequence is lighted in light blue, all the linkers are in black color (c) Ramachandran plot of the vaccine structure (d) Z-score of the vaccine structure

Molecular docking

The ClusPro 2.0 server was utilized to perform molecular docking between the vaccine and TLRs, while the preliminary molecular docking screening was run on HADDOCK. The vaccine-TLR2 complex has a center energy of −947.3 and a lowest energy of −1056.8, while the vaccine-TLR4 complex has a center energy of −819.5 and a lowest energy of −968.3. The HADDOCK refinement scored the vaccine-TLR2 complex to be −212.1 ± 2.3 and the vaccine-TLR4 complex to be −258.8 ± 2.1. The PRODIGY calculated the ΔG of vaccine-TLR2 to be −11.3 and Kd (M) at 25°C to be 5.4 × 10⁻⁹, which data are −16.9 and 4.1 × 10⁻¹³ separately for vaccine-TLR4 complex (Table 4). The interactions on the docking surface of the vaccine-TLRs complex are evaluated by PDBsum, results suggested plural salt bridges, hydrogen bonds as well as disulfide bonds are formed between proteins (Figures 4 and 5).

Figure 4. Docking analysis of vaccine-TLR2 complex, the docking result is on the left, and the analysis of the interactions inside the complex are on the right.

Figure 5. Docking analysis of vaccine-TLR4 complex, the docking result is on the left, and the analysis of the interactions inside the complex are on the right.

Table 4. Docking analysis of vaccine–TLR combination.

	TLR2	TLR4
	ClusPro
Center energy	−947.3	−819.5
Lowest energy	−1056.8	−968.3
	HADDOCK refinement
HADDOCK score	−212.1 ± 2.3	−258.8 ± 2.1
Cluster size	20	20
RMSD from the overall lowest-energy structure	0.6 ± 0.3	0.6 ± 0.3
Van der Waals energy	−99.8 ± 2.8	−105.7 ± 4.2
Electrostatic energy	−232.9 ± 8.5	−597.9 ± 26.7
Desolvation energy	−65.7 ± 1.0	−33.6 ± 1.6
Restraints violation energy	0.0 ± 0.0	0.0 ± 0.0
Buried Surface Area	2424.1 ± 64.4	3663.7 ± 26.2
ΔG (kcal mol^−1)	−11.3	−16.9
Kd (M) at 25°C	5.4 × 10^−9	4.1 × 10^−13

Molecular dynamics simulations

Gromacs was employed to possess molecular dynamics simulations to understand the binding affinity of the vaccine-TLRs complexes. The RMSD (root mean square deviation) served as an indicator of the structural fluctuations of the vaccine-TLRs complexes, The curves revealed both complexes tend to maintain a stable binding, and the RMSD of vaccine-TLR2 and vaccine-TLR4 remains stable after 10 ns in the 100 ns simulation, the average RMSD of vaccine-TLR2 is 1.284 while the average RMSD of vaccine-TLR4 is 1.382 (Figure 6a). The radius of gyration (Rg) was used to assess the degree of loosening in the complexes the curves of both vaccine-TLRs maintained relative equilibrium, the average Rg of vaccine-TLR2 is 3.801 while the average Rg of vaccine-TLR4 is 3.527, indicating both complexes gradually tighten in the 100 ns progress (Figure 6b). RMSF (root mean square fluctuation) analysis provided insights into the flexibility of the residues of vaccine and TLR proteins in complexes, the results demonstrated that residues around 400–500 of the vaccine-TLR2 exhibited relatively high flexibility while other residues kept a low level, the residues of vaccine-TLR4 complex always stay in a low level (Figure 6c,d). Meanwhile, the H-bonds curves showed both complexes tend to form more H-bonds during the interaction (Figure 6e). The gmx_MMPBSA v1.56 tool calculated the average ΔTOTAL binding energy of vaccine−TLR2 to be −91.73 and −126.59 kcal mol⁻¹ for vaccine-TLR4 complex (Table 5 and full results can be found in Table S1 and S2).

Figure 6. The results of molecular dynamics simulation.

(a) RMSD (root mean square deviation) plots of vaccine−TLRs, reflect the stability between the vaccine and TLR-2, TLR-4 receptor. (b) Radius of gyration (Rg) of vaccine-TLRs, suggesting the compactness of complexes (c, d) RMSF (root mean square fluctuation) of vaccine−TLRs, reflects the flexibility and fluctuation of the amino-acids residues in the side chain of docked complex Vaccine−TLR2(c), and TLR4(d). (e) H-bonds formed in complexes

Table 5. MMPBSA of vaccine-TLR2/TLR4 delta (complex - receptor - ligand).

Energy Component	Average(R2)	SD(R2)	Average(R4)	Average(R4)
ΔVDWAALS	−189.06	7.13	−229.19	9.61
ΔEEL	−522.57	49.76	−2696.96	67.77
ΔEGB	643.61	49.85	2831.71	69.32
ΔESURF	−23.72	0.92	−32.15	1.18
ΔGGAS	−711.63	50.56	−2926.14	69.31
ΔGSOLV	619.89	49.8	2799.56	68.84
ΔTOTAL	−91.73	8.69	−126.59	10.02

Immunization simulation

The C-IMMSIM tool accurately predicted immune responses in humans after three doses of vaccination, administered at 4-week intervals, the predicted results of the immune simulations generally align with the actual results of the vaccine injections. The second and third vaccinations induced significantly increased immune responses compared to the response after the first vaccination. The immune response exhibited elevated levels of IgM and combined IgG1+IgG2 antibodies, the peak titer of the antibody increases with the number of injections and gradually decreases after each injection, eventually reaching a stable presence in the body. (Figure 7a). The population of B cells and T cells also increased after each vaccination, and high levels of memory B cells and TH cells were generated simultaneously (Figure 7b–e). Meanwhile, the vaccination promotes the activation of the innate immune system, as evidenced by the activation and increase of macrophages (Figure 7h). Due to the introduction of HTL epitopes that have a strong stimulating effect on cytokines expression, the immunomodulatory factor, mainly IFN-g and IL-2, significantly increases after the vaccination (Figure 7i). These results indicate that the vaccine can elicit durable immune responses in both B cells and T cells, thereby demonstrating its effectiveness as a prophylactic measure.

Figure 7. Immune simulation after 3 doses of vaccination.

(a). Antibody levels induced by three doses of vaccine injection. (b) Population of B cells induced. (c) Levels of plasma B cells induced. (d) Population of helper T(TH) cells (e) Population of TH cells per state (f) Population of cytotoxic T(TC) cells induced (g) Population of Nature killer (NK) cells (h) Population of Macrophage (MA) cells (i) Levels of cytokines induced.

Population coverage

The distribution of MHC varies among different ethnic groups, and the population coverage of the constructed vaccine is calculated based on MHC alleles corresponding with epitopes in the sequence. The constructed vaccine has the capacity to provide coverage to 99.92% of the world’s population, and more than 99% in most regions (Figure 8). The graphical result is drawn based on pixelmap(https://pixelmap.amcharts.com/).

Figure 8. Population coverage of the vaccine.

Construction of multi-epitope mRNA vaccine and vaccine vector

Reverse translation and optimization were possessed on Jcat server, after improvement, a 2573 bp length cDNA was captured, which including full-length epitope vaccine as well as tPA and MITD sequence, its CAI value is 1.0 and its GC content is 49.98, both lies in an ideal range, suggesting good density and thermos-stability. After connecting with UTR sequence, the full-length DNA has 2624 bps. Then, RNAfold predicted its free energy of the thermodynamic ensemble (free energy of the secondary structure of vaccine mRNA) to be −796.96 kcal/mol, and the plasmid structure’s minimum free energy is −507.49 kcal/mol, graphical output can be found in supplement materials. The results predicted revealed the stability of the multi-epitope mRNA vaccine (Figures S1 and S2).

Finally, the DNA sequence of mRNA vaccine was inserted into Pet28a (+) plasmid between BamHI and XhoI with gensmart server to construct the expression vector (Figure 9).

Figure 9. In silico cloning of codon-optimized vaccine into E. coli K12’s expression system.

Discussion

HMPV has grown into a novel threat to public health these days, in the period from birth to 5 y old, almost 100% of the population experience infection once or more, after adulthood, people still face the threat of reinfection, especially those with poor immune response. Vaccination is considered the most effective method to combat epidemics, the final aim is to reduce RTIs in the target population including infants and the elderly. In earlier presents, Cseke et al. immunized cotton rats with soluble F protein, F DNA plasmid, or F protein after F DNA plasmid, all three groups were protected from HMPV infection, with the protection provided by the F protein being the strongest.⁴¹ F protein vaccines vectored on alphavirus replicon or other respiratory viruses have also shown varying degrees of protection in immunized animals.⁴² VLP-based vaccine is also a hopeful candidate, Cox et al. designed a VLP vaccine based on F and M proteins, similarly, Levy et al. developed a VLP vaccine based on F and G proteins, both vaccines provided protection for HMPV infected animal models.^43,44 Shaw’s team developed an mRNA vaccine coding F protein of HMPV and PIV3, which has entered Phase I clinical trials.⁴⁵ In the past decades, many attenuated live vaccines have been developed, most of which performed well in preclinical trials. rHMPV-Pa entered clinical trial but did not achieve the expected effects because of over attenuation.⁴⁶ Meanwhile, the knockout of G, SH, and M2–1, M2–2 proteins can reduce viral virulence, further confirming the role of proteins other than F in the immune process of HMPV infection.^47,48 The short discovery time and lack of sufficient attention have led to the absence of specific vaccines against HMPV to date, successful experiences in developing the Respiratory syncytial virus (RSV) vaccine indicated the important role of the F protein, however, led to the neglection of other proteins. Recently, doctor August’s team developed a mRNA vaccine mRNA-1653 based on F protein of both HMPV and Parainfluenza Virus Type 3. The safety and ability to induce antibodies of this vaccine have been verified in humans, further suggesting that mRNA is a safe and reliable platform, as well as the effectiveness of F protein as a target.⁴⁹ Thus, at this present, we managed to design a novel multi-epitope mRNA vaccine with immunoinformatic methods to face the threat of HMPV.

Advancements in immune informatics technologies, represented by epitope prediction, have brought new insights to vaccine development. Screening CTL, HTL, and B epitopes from target proteins and combining them with selected adjuvant sequences for design can enhance the safety and efficacy of peptide vaccines, methods such as dynamic simulations and immune simulations provide validation to vaccine design. Further animal experiments validate the protective effect of in silico designed vaccines on Streptococcus,⁵⁰ Shigella flexneri ,⁵¹ and other pathogens. To conclude results have proved the reliability of immunoinformatic methods.

Epitope selection is of paramount importance in the design of multi-epitope vaccines. In this report, our team predicted epitopes from mainly F proteins, as well as epitopes from G, SH, M1, M2–1, and M2–2. F protein is the most favorable target for vaccine development. F protein participates a vital role in recognition and fusion with host cells, and it has the strongest immunogenicity and proved protective response.^52,53 Meanwhile, most identified epitopes are located on F protein, the use of the F protein alone can induce the production of protective monoclonal antibodies, while currently no other proteins have been found to generate neutralizing antibodies.^54,55 However, proteins other than F are also involved in immune response, SH can inhibit the secretion of inflammatory cytokines, while G protein has the ability to decrease innate immune responses. Although viruses lacking G and SH proteins can still replicate, their replication efficiency will be greatly reduced.^56–58 M2–1 is a protein regulating RNA synthesis, it is highly associated with viral replication efficiency. And M2–2 helps HMPV escape from the immune system by inhibiting MAVS-dependent cellular response, knockdown of M2 proteins influences the immunogenicity of the virus.^47,48,59,60 Therefore, we mainly focus on the F protein epitope, and selected some predicted epitopes with better results from other proteins, and finally 9 LBL epitopes, 15 CTL epitopes, and 4 HTL epitopes were selected. Stimulating humoral immunity to generate a sufficient titer of antibodies is the focused point in vaccine design. However, simultaneously activating B cells and T cells can generate a stronger immune response, cellular immune response also responds to viral infection, and the ability to promote cytokine production of HTL epitopes mobilizes innate immunity to participate in viral defense. To provide protection for different serum types of HMPV, the epitopes are supposed to be conservative among 392 strains from 5 subtypes, due to high variability among different strains, there was no ideal epitopes found in G protein. The candidates were selected from prediction results by their composite score, antigenicity, and immunogenicity, also they are supposed to be not toxic or homologous with humanoid protein in case of potential harm. Besides, a PADRE sequence severs as a supplement to HTL epitopes, which is able to correspond with 15 MHC II alleles.

Appropriate adjuvants and linkers can enhance the efficacy of vaccines, GPGPG, AAY, and KK linkers are utilized to connect epitopes, these linkers are applied for better cleavage and enhanced immune response.^61,62 Two adjuvants are applied in the design, Pam2Cys is a short sequence agonizing TLR2/6 pathway and mediating a series of pro-inflammatory cytokine (including IL-1β, TNF-α, IL-12 and IL-17) releases.⁶³ Pam2Cys can enhance the humoral and cellular immune responses stimulated by mRNA vaccines without safety concerns, making it suitable as an adjuvant for mRNA vaccines.⁶⁴ While β-defensin II is a commonly used antiviral adjuvant that can enhance innate immunity and immune induction against viruses.⁶⁵ Thus, a sequence of Pam2Cys is attached to the C-terminal of the vaccine, meanwhile, β-defensin II is attached to the N-terminal both with EAAAK linkers.⁶⁶ Finally, the vaccine sequence contains 507 residues, its molecular weight is 55.2 kDa and its PI is predicted as 9.41, both are ideal for a vaccine construct.⁶⁷ With an instability index of 22.50 and an aliphatic of 81.56, the vaccine protein is considered stable and possesses high thermal stability, usually an instability index below 40 is considered stable, lower index perform better results, while higher aliphatic indicates higher thermal stability. High-level stability and thermal stability suggest that in-body synthesized vaccine proteins can stably exist, providing more opportunities for immune responses. With a hydrophilicity index of −0.051, the vaccine demonstrates its hydrophilic properties.⁶⁸ Meanwhile, the inspection of the constructed candidate revealed it has high antigenicity, and also immunogenicity without toxicity or allergenicity.

To ensure the vaccine can interact with immune receptors to induce immune response, molecular docking and dynamics simulation between the vaccine and TLRs are processed. Toll-like receptors (TLRs) are important immune receptors in viral inflammatory, the respiratory syncytial virus (RSV) virus, which is closely related to HMPV, can induce immune responses through TLR4 activation, TLR4-deficient mice are unable to clear RSV infection.^69–71 Additionally, as a molecule that recognizes a wide range of ligands, TLR2 also participates in the signal recognition of RSV, activating pro-inflammatory responses.^72,73 Thus, we simulated the docking between vaccine and TLR2/TLR4 to valid it has the ability to induce immune responses. Firstly, we predicted the vaccine’s tertiary structure and refined it, the model quality is acceptable through Ramachandran plotting and ERRAT quality analysis. In molecular docking, we observed that the vaccine binds tightly with TLR2 and TLR4, and plural bridges are formed in the interaction interface, these bridges include salt bridges and hydrogen bonds, disulfide bridges, etc., the more interactions formed, the more closely the molecular interaction is. Further MD simulation estimated the complex, during the 100 ns simulation, the vaccine and receptors tend to form stable interactions in simulated environments, RMSD (Root Mean Square Deviation) measures the difference between the simulated molecular structure and the experimental structure. The smaller the RMSD, the closer the simulation result is to the experimental result. both RMSD curves run flat at a low level, the average RMSD of vaccine-TLR2 is 1.284 while the average RMSD is 1.382 for vaccine-TLR4 complex. Rg (Radius of Gyration) describes the size and shape of the molecule. Smaller Rg values indicate a more compact molecular shape, while larger Rg values indicate a looser molecular shape. The Rg curves decreased smoothly, indicating the complexes tightened gradually until it reached a relative equilibrium state where it stabilized. Meanwhile, RMSF (Root Mean Square Fluctuation) measures the degree of variation in the positions of atoms within the molecule. Smaller RMSF values indicate more stable movements during the simulation process, and in the RMSF analysis of two complexes, only the vaccine-TLR2 complex has a slightly higher level of 400–500 residues of the vaccine, while the other residues maintain a lower level of RMSF, indicating that most of the residues are relatively stable in the simulation. Docking and analyzing of vaccine with TLRs indicates the candidate has an ideal affinity with receptors and the potential to stimulate strong immune responses, which is also validated by immune simulation. In immune simulation, the vaccine can elicit high titer antibody response, the total IgG and IgM levels exceed 65,000/ml after the third infection, which suggests that vaccination can generate a strong antibody response to neutralize the virus. Meanwhile, an increase in the population of B cells was observed after immunization, with multiple isotypes maintained at a higher level for an extended period of time following the final immunization. Similarly, vaccination leads to an increase in TH and TC cells, as well as elevated levels of cell factors such as IFN-g. Results align with the epitope selection and also suggests that the vaccine can induce strong humoral and cellular immune responses. Lastly, the protein sequence is converted into DNA and engineered to perform multi-epitope mRNA vaccine, and the full-length sequence was transferred to the pET28(a) plasmid as a precursor for synthesis in E. coli expression system.

There are still certain limitations in this study. Though many vaccines designed through epitope selection have achieved ideal protective effects equal to in silico predictions in animal experiments, still we cannot determine this candidate vaccine’s efficacy only in silico. Meanwhile, the experimentally validated epitopes can usually be obtained from epitope prediction tools, and they rank high in the overall evaluation. However, not all high-scoring peptides can be experimentally validated. Therefore, we have to appropriately increase the number of epitopes in the vaccine sequence, which will result in an increased total length, possibly affecting the stability of the final product and increasing the difficulty of acquisition.

Conclusion

In this present, we developed a novel multi-epitope mRNA vaccine candidate against HMPV with the assistant of immunoinformatic methods. Compared with previous researches, we choose multi-epitope mRNA vaccine as the formation, through incorporating immunodominant epitopes, we aim to achieve stronger targeting while minimizing undesirable effects. Meanwhile, a design with conserved epitopes can provide a broader range of protection among different subtypes, and a vaccine containing T cell and B cell epitopes is supposed to elicit both humoral and cellular immune responses. Besides, this design introduced epitopes from SH, M1, and M2 proteins to support F epitopes, by targeting other proteins with immune-modulatory functions, we can enhance the functionality of the vaccine. It should be pointed out that the final efficacy of the vaccine still needs to be determined through in vitro and in vivo experiments. But in summary, the design of this vaccine is a promising vaccine candidate that performs well in simulations. Given the limited research on HMPV vaccines, it can provide an effective strategy for the prevention of HMPV.

Author contributions

Conceptualization, F.Z, and S.M.; methodology, F.Z.; validation, S.M. and Y.X., and M.R.; formal analysis, P.Z, W.P, H.Y and Y.C; writing – original draft preparation, S.M, Y.X and C.T; writing – review and editing:J.C; supervision, P.P and J.C.; project administration, P.P and J.C; funding acquisition, P.P

Data available statement

The data that support the findings of this study are available from the corresponding author, Pinhua Pan, upon reasonable request.

Declaration of generative AI in scientific writing

No AI technology is used in scientific writing.

Submission declaration and verification

This article has not been published previously, and it is not under consideration for publication elsewhere, its publication is approved by all authors and tacitly or explicitly by the responsible authorities where the work was carried out.

Acknowledgments

Thanks to the support from HY.W, C.Z, and L.Z.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website at https://doi.org/10.1080/21645515.2023.2293300.

Additional information

Funding

This work is supported by the National Natural Science Foundation of China [No. 81770080]; Project Program of National Clinical Research Center for Geriatric Disorders [Xiang ya Hospital,Grant No. 2020LNJJ05]; The National Natural Science Foundation of China [No. 8210012334]; National Key R&D Program of China [No. 2016YFC1304204]; Key R&D Program of Hunan Province [No. 2022SK2038]; Project Program of central south university graduate education teaching reform [No. 2022JGB025]; The national key clinical specialist construction programs of China [Grant Number z047-02]; Research Project of Teaching Reform in Colleges and Universities in Hunan Province [2021jy139-2]; China Postdoctoral Science Foundation [2022M713520]; Project on the topic of emergency response to combat the COVID-19 epidemic in Hunan Province [2022SK2122]; Research Project on teaching reform of ordinary colleges and universities in Hunan province [2022JGYB037].