The entire sequences of the S protein of PEDV G2 genotype strains were obtained from the National Center for Biotechnology Information (NCBI) database. The antigenicity of all S proteins was predicted using the free online software VaxiJen v2.0.¹ Four representative PEDV strains with the best antigenicity from each G2 subgroup (G2a, G2b, G2c, and G2d) were chosen to identify antigen epitopes and design vaccine.

The sequence alignment of the S proteins from CH/SCZJ/2018, LW-L, CH/SDLQ/09/2020, and CH/SCST/04/2020 strains was performed using SnapGene tool,² and the mutant amino acid sites were then screened out. To determine the presence of a signal peptide region in the candidate antigen protein, the signal peptide of the S protein was predicted using SignalP-5.0 server³ (Almagro Armenteros et al., 2019).

To predict cytotoxic T lymphocyte epitopes (CTLs), the sequences of S proteins were submitted to the Immune Epitope Database (IEDB) server⁴ (Fleri et al., 2017). The 45 common Swine Leukocyte Antigen (SLA) class I molecules were utilized in the prediction by running the “IEDB-recommended” method for each allele. To estimate CTLs, 9-mer epitopes were screened using parameters such as TAP score > 1.0, IC50 < 500 nM, and proteasome score > 1.0. The dominant epitopes, which simultaneously appeared in three SLA-I alleles in each viral strain and had antigenicity values higher than the threshold (0.4) by using VaxiJen v2.0, were compared and incorporated for vaccine formulation. The common dominant epitopes of the four strains were used to create the final vaccine.

To predict helper T lymphocyte epitopes (HTLs), the sequences of S proteins were assessed using the online server NetMHCIIpan 4.0.⁵ Twenty-seven high-frequency human MHC II (HLA-II) alleles were used to predict epitope peptides. The epitopes with a length of 15 amino acid residues were utilized to estimate HTLs. The threshold for peptides with strong bindings was set to its default value. The dominant epitopes, which simultaneously appeared in three HLA-II alleles in each viral strain and had antigenicity values higher than the threshold (0.4) by using VaxiJen v2.0, were compared and incorporated for vaccine formulation. The common dominant epitopes of the four strains were used to create the final vaccine.

Linear B cell epitopes (LBEs) were predicted using the IEDB server with the method of Bepipred Linear Epitope Prediction 2.0 at the default threshold of 0.5. Additionally, the predicted epitopes were screened using VaxiJen v2.0. The common epitopes of the strains were used for vaccine construction.

The structures of S proteins were generated through modeling utilizing the SWISS-MODEL server.⁶ The reliability of the prediction results was evaluated by Global Model Quality Estimate (GMQE) (Yassim et al., 2022). Ramachandran plot and ProSA-web were used to assess the quality of the models (Waterhouse et al., 2018).

To predict conformational B cell epitopes (CBEs), the online software DiscoTope⁷ with a threshold value of −3.7 was used. The common regions of the predicted epitopes from the four strains were taken as the final conformational epitopes and visualized with PyMol.

The final subunit vaccine was designed by sequentially joining the resulting peptide sequences using suitable linkers. Porcine β-defensin-2 (PBD-2) (37 aa) was used as an adjuvant and fused to the N-terminal via the EAAAK linker. The CTLs, HTLs, and B cell epitopes were linked by AAY, GPGPG, and KK, respectively. The TAT sequence (11 aa) was attached to the C-terminal of the vaccine construct. The recombinant PEDV multi-epitope vaccine was named rPMEV in the study.

The secondary structure of rPMEV was predicted using SOPMA online analysis software⁸ (Deléage, 2017). The initial tertiary structure was predicted by Robetta server⁹ (Baek et al., 2021). To improve the structure quality and protein stability, the initial vaccine 3D model was refined by GalaxyRefine server.¹⁰ To validate the refined tertiary structure, a Ramachandran plot that could display favorable areas of backbone dihedra angels against amino acid residues in 3D structure was performed by the SWISS-MODEL workspace (see text footnote 6) (Waterhouse et al., 2018). In addition, ProSA-web was also employed in the final vaccine protein structure validation.

The online tool Protparam¹¹ was used to predict the physicochemical properties of rPMEV. To ensure the safety of the vaccine, the online software AllerTOP v2.0¹² was used to predict the allergenicity of rPMEV. Meanwhile, the VaxiJen v2.0 server was applied to evaluate the antigenicity of rPMEV.

To assess the interaction between the vaccine and the immune receptors, molecular docking of rPMEV and the TLR4 receptor (PDB ID: 4G8A) was performed by ClusPro server¹³ (Kozakov et al., 2017). Additionally, Ligplot and PyMol software were used to analyze the interaction interface residues.

To evaluate the structural properties and interaction in rPMEV and TLR4 receptor complex, molecular dynamics (MD) simulation was performed using GROMACS (Groningen Machine for Chemical Simulations) (He et al., 2021). First, the AMBER99 force field was used to generate the topology for the protein-ligand complex in all of the MD simulations. Then the protein was solvated in a cubic box of TIP3P water (Grifoni et al., 2020), with at least 1.0 nm distancing the protein from the box edge. To balance the charged protein complex, ions were introduced using a genion tool (Shukla et al., 2018). Moreover, in order to prevent steric conflicts and improper geometry, the solvated electroneutral system was also relaxed through energy minimization. Then, the system equilibration was performed with 100 ps of NVT [substance (N), volume (V), and temperature (T)] equilibration and 100 ps of NPT [substance (N), pressure (P), and temperature (T)] equilibration without restraints. Finally, to analyze the trajectory, MD simulations of 100 ns were performed for the complex system.

To detect the immune response of rPMEV to the host, the C-ImmSim Server¹⁴ was used to perform the immune simulation (Rapin et al., 2010). The simulation step was set at 1050 (Bhatnager et al., 2021). The time steps were 1, 63, and 126 (one time step corresponds to 8 h). The other parameters were used as the default simulation parameters.

For cloning, codon optimization was performed by the online tool Java Codon Adaptation Tool (JCat)¹⁵ (Grote et al., 2005). Escherichia coli (Strain K12) was chosen to express the vaccine protein. Codon adaptation index (CAI) and the guanine and cytosine (GC) contents were evaluated. For raising the expression efficiency of the final vaccine protein, pET28a (+) was selected as the vector to express the vaccine protein. Before inserting the nucleotide sequence of the vaccine into the vector, the XhoI and EcoRI restriction enzymes, which are not contained in the vaccine sequence, were inserted at the N- and C-terminals of the sequence, respectively. Finally, the modified nucleotide sequence of the vaccine was cloned by SnapGene tool (see text footnote 2).

The GeneBank accession numbers of the four PEDV strains are MH061342.1 (CH/SCZJ/2018), MK392335 (LW/L), MZ161085 (CH/SDLQ/09/2020), and MZ161081 (CH/SCST/04/2020), respectively. The GeneBank accession numbers of S proteins derived from four PEDV strains are AZL47228 (CH/SCZJ/2018), QBM00061 (LW-L), UDL09544 (CH/SDLQ/09/2020), and UDL09540 (CH/SCST/04/2020). The results of the antigenicity of all S proteins are shown in Supplementary Table S1, and the amino acid sequences of these proteins are shown in Supplementary material.

To reveal the characteristics of the S proteins from these four PEDV strains, sequence alignment of S proteins was performed. The results, as shown in Figure 2, show that there are 114 distinct loci among the four strains. The S1 region contains 85 sites, 71 of which are in the N-terminal domain (NTD) and 6 in the C-terminal domain (CTD). The S2 region had 29 sites, 4 of which are in heptad repeat region 1 (HR1) and 2 in HR2. The conserved areas of S proteins are the non-mutant areas.

After prediction, signal peptide sections were found in S protein of all four PEDV strains, as shown in Supplementary Figure S1. The signal peptide region of CH/SCZJ/2018 is 1–18 (MKSLTYFWLFLPVLSTLS) (Supplementary Figure S1A), the signal peptide sequence of LW/L is (MKSLTYFWLFLPVLSTLS) (Supplementary Figure S1B), the signal peptide sequence of CH/SDLQ/09/2020 is (MKSLTYFWLFLPVLSTLS) (Supplementary Figure S1C), and the signal peptide sequence of CH/SCST/04/2020 is (MKSLTYFWLFLPVLSTFS) (Supplementary Figure S1D). All of the signal peptide sequences were removed for the later epitope prediction of these four PEDV strains.

The CTLs (9-mers) of 4 PEDV strains were obtained by the IEDB server, and each virus strain has 45 kinds of CTLs of SLA-I alleles. All CTLs were screened with a TAP score > 1.0, IC50 < 500 nM, and proteasome score > 1.0, as shown in Figure 3A. The dominant CTLs for each virus strain that appear in three SLA-I alleles simultaneously and have antigenicity values over the threshold (0.4) were chosen, and they were ranked from small to large according to the sequence start position as shown in Figure 3B. The dominant epitopes of each virus strain are listed in Supplementary Tables S2–S5. The common dominant epitopes of the four strains, as shown in Figure 3C, were compared and integrated for vaccine design.

The HTLs of four PEDV strains were analyzed by online NetMHCIIpan 4.0 server. The strong binder peptides with a % Rank less than 1.0 were included in this study, as shown in Figure 4A. For each virus strain, the dominant HTLs that appear in three HLA-II alleles at the same time and have antigenicity values greater than the threshold (0.4) were chosen and ranked from small to large based on the sequence start position, as shown in Figure 4B. The dominant epitopes of each virus strain are listed in Supplementary Tables S6–S9. The common dominant epitopes of the four strains, as shown in Figure 4C, were compared and integrated for vaccine design.

The LBEs of four PEDV S proteins were predicted by IEDB software. The peptide sequence was identified as an epitope when the peptide threshold in the BepiPred-2.0 method is larger than 0.5 (Supplementary Figure S2). The dominant LBEs were identified based on peptide epitopes with antigenicity values larger than the threshold (0.4), as shown in Figure 5A. The dominant LBEs of four virus strains were compared, as indicated in Figure 5B, and the common dominant LBEs for the four strains were employed to produce the final vaccine (Supplementary Table S10).

Based on target-template alignment, the 3D structure of S proteins was modeled by SWISS-MODEL server. The GMQE scores of the four S protein models are around 0.70. All models were verified by ProSA-web, which can calculate the Z-Score of models. The Z-Score of CH/SCZJ/2018 model is −10.62, the LW/L model is −10.71, the CH/SDLQ/09/2020 model is −10.81, and the CH/SCST/04/2020 model is −10.71, as shown in Supplementary Figures S3A–D, respectively. The Ramachandran plot drawn by SWISS-MODEL showed that the model residues of CH/SCZJ/2018, LW/L, CH/SDLQ/09/2020, and CH/SCST/04/2020 in the Ramachandran favored region are 94.98%, 95.68%, 95.66%, and 95.79%, respectively, as shown in Supplementary Figures S3E–H. These findings suggested that the S protein models of the four PEDV strains are of high quality and accuracy.

To predict CBEs of each strain, the created models of S proteins were submitted to the online software DiscoTope. All the CBEs of each strain are represented in Supplementary Table S11. The positions of selected conformational epitope residues for each S protein are displayed in Figure 6.

The rPMEV construct was developed by joining the best candidate antigen epitopes. There are 9 CTLs, 11 HTLs, 6 LBEs, and 4 CBEs in rPMEV, as shown in Figure 7A. The CTLs were linked by the AAY linker, which helps the epitopes produce suitable sites for binding to the TAP transporter and enhance epitope presentation (Ayyagari et al., 2022). The HTLs were linked by the GPGPG linker, which can simulate HTL responses and conserve conformation-dependent immunogenicity of helpers as well as antibody epitopes (Ayyagari et al., 2022). The LBEs and CBEs were joined using the KK linker, which can increase immunogenicity (Li et al., 2015). Furthermore, at the N-terminus, porcine β-defensin-2 (PBD-2) (37 aa) was added to the vaccine construct to increase the immunogenicity of the multi-epitope vaccine using the EAAK linker, which can effectively separate proteins while reducing their contact (Kumar Pandey et al., 2018). To enable intercellular delivery, a TAT sequence was appended to the vaccine construct at the C-terminal. The number of vaccine residues in the final designed multivalent vaccine is 508 aa, as seen in Figure 7B.

The secondary structure of rPMEV was analyzed using SOPMA server. The results show that there is 27.56% alpha helix, 23.03% extended strand, 5.51% beta turn, and 43.90% random coil, as illustrated in Supplementary Figure S4.

Based on the sequence, the Robetta server was used to model the tertiary structure of rPMEV, and a total of five models were output. The quality and potential errors in the five models were verified by ProSA-web, as shown in Supplementary Figure S5. The Z-Score indicates overall model quality. The Z-Scores of models 1–5 are −7.16, −7.2, −5.4, −6.69, and −5.52, respectively. Model 1 was adopted as the initial model of rPMEV (Supplementary Figure S6A) since it has a lower Z-Score and is therefore of higher quality (Supplementary Figure S6B). The Ramachandran plot drawn by SWISS-MODEL workspace reveals that the initial model has 93.68% Ramachandran favored region, 1.58% Ramachandran outlier region, and 0.00% Rotamer outlier region, as shown in Supplementary Figure S6C. To optimize the tertiary of the vaccine, the initial model was refined by GalaxyRefine Server, and there were five optimized models output by the server. Similar to the initial model selection, all the optimized models were further verified by ProSA-web as shown in Supplementary Figure S7. The Z-Scores of the optimized models 1–5 are −7.39, −7.25, −7.29, −7.36, and −7.31, respectively. The optimized Model 1 was chosen as the final tertiary structure of rPMEV (Supplementary Figure S6D) since it has the lowest Z-Score and is of the highest quality (Supplementary Figure S6E). The Ramachandran plot drawn by SWISS-MODEL workspace shows that the final tertiary structure is 95.06% of residues in the Ramachandran favored region, 0.40% in the Ramachandran outlier region, and only 0.99% in the Rotamer region, as shown in Supplementary Figure S6F.

The vaccine protein is composed of 508 amino acids and 7,506 atoms, with the molecular formula C₂₄₁₂H₃₇₁₄N₆₃₀O₇₃₆S₁₄. The molecular weight of rPMEV is 54 KD, indicating that it is a suitable vaccine because proteins with molecular weights less than 110 KD can be efficiently purified (Barh et al., 2013). The theoretical pI of rPMEV is 8.59, and it includes 37 negatively charged residues (Asp and Glu) and 43 positively charged residues (Arg and Lys). The instability index (II) of rPMEV is 34.74, and the vaccine is classified as a stable protein since it is smaller than the instability index value of 40 (a value above 40 predicts that the protein is unstable) (Zhao et al., 2019). The aliphatic index of rPMEV is 71.50, which is calculated as the proportional volume occupied by aliphatic side chains. The grand average of hydropathicity (GRAVY) is −0.244. Since the GRAVY ranges from −2 to 2, negative values indicate that the protein is hydrophilic (Sha et al., 2020), so rPMEV belongs to a hydrophilic nature.

To avoid autoimmune reactions, an allergenicity analysis of rPMEV was performed using AllerTOP v2.0. Both rPMEV and the nearest protein, UniProtKB accession number O1449, are defined as non-allergens. The antigenicity of rPMEV predicted by VaxiJen v2.0 is 0.654, which is more than the threshold of 0.4. Therefore, allergenicity and antigenicity studies suggest that rPMEV is a safe protective antigen.

To evaluate the interaction affinity between rPMEV and TLR4, molecular docking was carried out using ClusPro Server. The interaction of rPMEV with TLR4 produces 30 rPMEV-TLR4 complexes, and all of the docking data are displayed in Supplementary Table S12. The conformation with the least interaction energy (−1046.1 kcal/mol) obtained by the docking of the rPMEV-TLR4 complex is shown in Figure 8A. The interaction interface residues were analyzed by PyMol in 3D and Ligplot in 2D format, as seen in Figures 8B,C, respectively. The results demonstrate that there are 13 hydrogen bonds at the docking interface engaged in the interaction of the complex subunits.

To estimate the stability of the rPMEV-TLR4 complex, MD simulation was performed using GROMACS. During the MD simulation, energy minimization, temperature, pressure, RMSD, and RMSF were analyzed. As shown in Figures 9A,B, the temperature and pressure plots indicate that the system is around 300 K and 1.7 atmospheres with 100 ps of the time interval, indicating that the system is stable and the MD operation is successful. The RMSD against time of the rPMEV-TLR4 complex, which represents the structural fluctuation of the overall structure of the protein complex, reveals a large fluctuation during the 0–8 ns simulation. After 8 ns, the RMSD value is kept around 0.5 nm, indicating the stability of the vaccine complex (Figure 9C). In addition, the RMSF of the rPMEV-TLR4 complex, which indicates the flexibility of the residue in the docking complex, was also performed. The majority of amino acid residues in the rPMEV-TLR4 complex are below 2.0 nm, and only a few residues have larger changes, indicating that the complex has stability and stiffness (Figure 9D).

The C-ImmSim server was used to study the immune stimulation of the final vaccine. The simulation results reveal that the immune response after immunization is mainly composed of B cell (Figure 10A) and T cell populations. After the three immune responses, significant antibody titers are produced against the antigen (Figure 10B). Additionally, after the secondary and tertiary immune responses, the active TH cell populations and the populations of total and memory TH cells are also increased, as shown in Figures 10C,D, respectively. Furthermore, as illustrated in Figure 10E, the population of cytotoxic T cells also increases steadily after each immunization. Importantly, after the three stages of immunization, IFN-g (interferon-gamma) shows a robust response among the cytokines and interleukins (Figure 10F).

To produce a suitable plasmid construct containing the vaccine construct sequence, codon optimization was undertaken (Supplementary Figure S8). The results of codon optimization of the vaccine construct show that the improved vaccine has a CAI value of 0.947 and a GC-content of 52.56%, indicating that the final vaccine has a high potential for good expression in E. coli since CAI >0.8 is regarded as a good value for high-level expression, and the optimal GC content ranges from 30% to 70% (Nezafat et al., 2016). For cloning the vaccine construct into the E. coli pET28a (+) vector, restriction sites XhoI and EcoRI were added to the N- and C- terminals of the codon sequence, and then the sequence of 1,524 bases was inserted into the pET28a (+) vector using SnapGene software, as shown in Supplementary Figure S9.