An Immunoinformatics Approach for Identifying and Designing Conserved Multi-Epitope Vaccines for Coronaviruses

1. Introduction

Coronaviruses (CoVs) were thought to only cause common cold in humans until 2002-2003 when severe acute respiratory syndrome-coronavirus (SARS-CoV) struck the world health system. The SARS-CoV outbreak caused 8,098 infections and 774 deaths globally (~10% mortality) [1]. Nearly 10 years after the SARS outbreak, another coronavirus outbreak took place in the Middle East and the etiological agent was identified as Middle East respiratory syndrome coronavirus (MERS-CoV). The outbreak then spread to South Korea with the first case reported in 2015 involving a 65-year-old man who had recently travelled to the Middle East [2]. Compared to SARS-CoV, MERS-CoV caused the highest mortality rate of ~35% among the three outbreak strains with a total of 2,458 infections and 848 reported deaths. In December 2019, a novel coronavirus strain namely SARS-CoV-2 caused a large-scale viral outbreak in WuHan, China. Since then, it has spread rampantly throughout the world and caused respiratory distress in humans. The disease associated with SARS-CoV-2 was coined COVID-19 and resulted in the breakdown of the healthcare system in many countries. The outbreak was eventually announced as a global pandemic by WHO on 11 March 2020 [3]. As of July 2024, COVID-19 has affected more than 775 million people and caused more than 7 million deaths worldwide [4].

As the coronaviral spike (S) glycoprotein is located outside the viral particle and mediates the viral entry into the host epithelial cells, it is undoubtedly the main target of neutralizing antibodies (NAbs) upon infection, therefore, making it the most important therapeutic target and in vaccine design. However, the emergence of new SARS-CoV-2 Omicron variants has rendered the vaccines ineffective, with ChAdOx1 nCoV-19 (Vaxzevria, AstraZeneca) conferring almost no protection from 20-24 weeks after the second dose of vaccination [5]. It is also notable that the emergence of new variants of concern (VOC), such as Omicron has raised attention globally as the new variants can escape the neutralizing antibodies and have increased transmissibility due to the presence of more than 30 mutations as compared to the parental strain, SARS-CoV-2-Wuhan-Hu-1 [6,7,8,9]. In view of the rising concerns regarding the increased infectivity of the new variants and controversies about the effectiveness of the vaccines, there is an urgent need for a pan-coronavirus vaccine that can induce the synthesis of neutralizing antibodies and is more comprehensive in conferring protection against the newly emerging variants as well as future coronavirus outbreaks While many groups predicted and identified evolutionarily conserved epitopes in-silico [10,11,12,13,14], and some of them were validated in-vitro and in-vivo [15,16,17,18,19,20,21], this study scrutinized the conserved epitopes further. Many predicted CTL, HTL, and LBL epitopes were found in the close vicinity in the S glycoprotein. Instead of studying them individually, they were aligned into single and relatively long peptide sequences. This novel strategy of having multi-epitopes is expected to stimulate a stronger and multi-faceted immune response against coronaviruses, addressing the limitations of the current vaccines against the emerging variants.

In this study, the evolutionarily conserved epitopes in both human and animal coronaviruses were identified using unique immunoinformatics approaches. After a stringent scrutiny and selection, 52 CTL epitopes, 11 HTL epitopes, and 68 linear B-lymphocyte (LBL) epitopes were identified from 30 coronavirus sequences derived from human CoVs (hCoVs) responsible for the common cold, SARS-CoV, MERS-CoV and SARS-CoV-2. Subsequently, the predicted epitopes were aligned and assembled into 2 final composite peptide sequences that were found to be evolutionarily conserved in SARS-CoVs, bats, and pangolin coronaviruses. These 2 assembled epitopes are not only found to be conserved in many of the coronavirus strains, but they also possess HTL, CTL, and B-cell antigen binding sites, and they match a diverse array of HLA class-I and -II supertypes prevalent in the human population, indicating their potential to activate both T and B cells effectively. Although the epitopes were only found to be mostly conserved in SARS-CoVs, bats, and pangolin coronaviruses, their distinctive compatibility with human HTL and CTL, and B cells renders them high potential in vaccine development. Altogether, the discoveries not only pave the way for the development of a pan-coronavirus multiepitope vaccine to combat existing and novel coronavirus strains but immunoinformatics are highly applicable in universal vaccine development, especially in identifying immunogenic conserved epitopes in target antigens.

2. Materials and Methods

2.1. Coronaviral S Gene Sequence Retrieval and Sequence Conservation Analysis

Forty-two coronaviral peptide sequences of the S gene were retrieved from the NCBI GenBank (https://www.ncbi.nlm.nih.gov/genbank/) as listed in Table 1, Table 2 and Table 3. A total of 24 sequences encompassing SARS-CoV-2 and its variants were retrieved from the Genbank (Table 1.) with the latest variant being XBB.1.5. Six sequences of the other coronaviruses causing SARS, MERS and common flu in humans are listed in Table 2. Table 3., on the other hand, tabulates twelve sequences of coronaviruses isolated in bats, pangolins and birds.

In order to identify the conserved regions in the coronaviral S glycoproteins, the amino acid sequence of the S glycoprotein of SARS-CoV-2 Wuhan-Hu-1 strain (Table 1.) was used as a reference sequence to perform Clustal Omega multiple sequence alignments in the EMBL-EBI (https://www.ebi.ac.uk/). The alignment was based on the Percentage Identity Threshold of 80% in the amino acid sequences using Jalview 2.11.2.6 (https://www.jalview.org/). The evolutionarily conserved regions of the S glycoproteins were identified and subjected to antigenicity screening, selection, and assembly.

2.2. The flow of Prediction of Conserved HTL, CTL and Linear B-lymphocyte (LBL) Epitopes of Coronaviral S glycoproteins

The prediction was performed separately for (i) CTL, (ii) HTL, and (iii) LBL epitopes by referring to their respective databases. The flow of prediction of conversed epitopes is depicted in Figure 1. The conserved epitopes were individually screened and identified, and their antigenicity and toxigenicity were predicted using VaxiJen 2.0 and ToxinPred, respectively.

2.2.1. Prediction of Conserved CTL Epitopes

The conserved CTL epitopes were identified using NetCTL-1.2 (https://services.healthtech.dtu.dk/services/NetCTL-1.2/). A total of 30 amino acid sequences of human coronaviral S glycoproteins were uploaded to NetCTL-1.2 by following the default criteria, which entailed 9 amino acids in length with a minimum threshold of 0.75. The available HLA class-I supertypes provided by NetCTL-1.2 included A1, A2, A24, A26, B7, B8, B27, B39, B44, B58, and B62. The redundant epitope sequences were filtered and subjected to the subsequent screening. The selected epitopes were then subjected to in-silico antigenicity screening using VaxiJen 2.0 (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html). The screening was performed using the default settings and “Virus” was selected as the target organism as part of the prediction criteria. Antigens labelled as “Probable Antigen” were sorted and selected. Subsequently, the selected antigens were screened for their immunogenicity using the IEDB Class I immunogenicity web-based prediction tool (http://tools.iedb.org/immunogenicity/). The non-immunogenic epitopes were indicated with negative scores and removed, the remaining epitopes were later examined using ToxinPred (https://webs.iiitd.edu.in/raghava/toxinpred/protein.php) to eliminate the probable toxic CTL epitopes. Following that, the epitopes were selected based on their high frequencies across different strains and supertypes. Table 4. summarizes the prediction of conserved CTL epitopes.

2.2.2. Prediction of Conserved HTL Epitopes

In the conserved HTL epitope prediction, the amino acid sequences of the coronaviruses were screened using IEDB MHC-II (http://tools.iedb.org/mhcii/) with the following prediction conditions: (i) Percentile rank: 25%, 15-mers amino acid; (ii) prediction method: Consensus 2.22; (iii) HLA supertypes: HLA-DR, HLA-DQ, and HLA-DP. Table 5. summarizes the prediction conditions for the conserved HTL epitopes. Epitopes with a percentile rank lower than 20.0 were eliminated as it indicates that those epitopes capture less than 50% immune response [22].

Similar to the conserved CTL epitope prediction, the antigenicity and toxigenicity of the epitopes were analysed using the same methods as shown in Table 4. In addition, the HTL epitopes were screened using IFNepitope (https://webs.iiitd.edu.in/raghava/ifnepitope/predict.php) for their abilities to induce interferon synthesis. The prediction criteria included “Motif and SVM hybrid” as the prediction approach and “IFN-gamma versus Non IFN-gamma” as the prediction model. The “NEGATIVE” HTL epitopes were removed.

2.2.3. Prediction of Conserved LBL Epitopes

ABCPred (https://webs.iiitd.edu.in/raghava/abcpred/ABC_submission.html) and SVMTriP (http://sysbio.unl.edu/SVMTriP/prediction.php) were used in predicting the conserved LBL epitopes with criteria such as 16-mers amino acids in length for both tools and the threshold of 0.51 and 0.50, respectively as shown in Table 6. All the predicted epitopes were selected for further antigenicity and toxigenicity prediction. The antigenicity and toxigenicity of the epitopes were screened as described in Table 4.

2.3. Alignment of the Predicted Conserved CTL, HTL and LBL Epitopes and Allergenicity Prediction

To identify the locations of the epitopes identified in Section 2.2, SARS-CoV-2-Wuhan-Hu-1 S glycoprotein was used as a reference sequence for multiple sequence alignment. The overlapped regions of the epitopes were aligned, they were then assembled into long amino acid sequences containing the conserved CTL, HTL, and LBL epitopes. Subsequently, the assembled peptide sequences were screened for allergenicity using AllergenFP v1.0 (https://www.ddg-pharmfac.net/AllergenFP/) to identify the probable allergens in the assembled sequences. The probable allergenic sequences, if any, were eliminated.

2.4. Structural Visualisation of Assembled Epitopes

The structural information of the close (PDB: 6VXX) and open states (PBD: 6VYB) of the S glycoprotein were retrieved from the RCSB Protein Data Bank (https://www.rcsb.org/). By using the retrieved information, the structures of the assembled conserved epitopes were then visualized using the UCSF ChimeraX (https://www.cgl.ucsf.edu/chimerax/).

3. Results

3.1. Conserved Regions in the S glycoproteins of Bat and Pangolin CoV, hCoVs, SARS-CoV-2, SARS-CoV, and MERS-CoV

The S glycoproteins of coronaviruses consist of a signal peptide, a receptor binding subunit (S1) and a fusion subunit (S2) with a length ranging from 1105-1351 amino acids [23,24,25,26]. With reference to the SARS-CoV-2-Wuhan-Hu-1, in the S1 subunit, there is an N-terminal domain (NTD, 14–305 residues) and receptor-binding domain (RBD, 319–541 residues), whereas, in the S2 subunit, there is a fusion peptide (FP, 788–806 residues), heptapeptide repeat sequence 1 (HR1) (912–984 residues), HR2 (1163–1213 residues), transmembrane domain (TM, 1213–1237 residues), and cytoplasmic domain (1237–1273 residues) [27]. As the coronavirus S glycoprotein is located outside of the viral particle and mediates the viral entry into the host epithelial cells, it is undoubtedly the main target of neutralizing antibodies (NAbs) upon infection, therefore, making it the most important therapeutic target and essential in vaccine design.

Given the importance of coronaviral S glycoprotein in vaccine development, this study aimed to screen and identify the evolutionarily conserved sequences in animal and human coronaviral S glycoproteins. The amino acid sequence of SARS-CoV-2-Wuhan-Hu-1 S glycoprotein served as the reference sequence to all retrieved coronavirus sequences. Upon screening and alignment, the results revealed that bat coronavirus strains such as Bat CoV RATG13, Bat CoV ZXC21, and Bat CoV YN02, and pangolin coronavirus strains such as Pangolin CoV GX-P2V, Pangolin CoV GX-P5E, Pangolin CoV GX-P5L, Pangolin CoV GX-P1E, Pangolin CoV GX-P4L, and Pangolin CoV MP789 showed some level of evolutionary divergence compared to that of SARS-CoV-2-Wuhan-Hu-1 with an identity threshold above 80% (Supplementary Figure S1). The results also explain why bats or pangolins are deduced as the most likely reservoirs of SARS-CoV-2. In contrast, three selected avian coronaviral S glycoproteins were distantly related to the reference sequence due to their relatively high evolutionary divergence.

The amino acid sequences of SARS-CoV-2 and hCoVs causing the common cold were also aligned with the reference sequence (Supplementary Figure S2). The S1 regions of hCoVs, i.e., H-CoV-HKU1–genotype B, CoV-OC43, CoV-NL63, and CoV-229E showed insignificant similarities to the reference sequence. Interestingly, their S2 regions were relatively conserved particularly at residues S₈₁₅ - S₈₇₄, S₈₉₇ - S₉₃₄, S₉₄₄ - S₁₀₆₉, and S₁₂₀₇ - S₁₂₁₈. The residues S₈₉₇ - S₁₀₆₉ corresponded to the HR1 and HR2 regions of the S2 subunit, whereas the S1207 - S1218 region was part of the HR2 and TM domain. The conservation of the HR1 and HR2 regions was documented previously and suggested as the targets for the development of fusion inhibitor agents [28,29]. Furthermore, the alignment of S glycoprotein sequences of SARS-CoV-2 and its variants, MERS-CoV, and SARS-CoV revealed that SARS-CoV, SARS-CoV-2 and its variants were highly similar to each other (Supplementary Figure S3). This finding suggests that the emergence of SARS-CoV and SARS-CoV-2 might be due to the recombination of viral genomes between bat coronaviruses in their natural reservoir (bats) or the intermediate host (pangolin), or both. There were no observable evolutionarily conserved regions in the S glycoprotein sequence of MERS-CoV in relative to that of SARS-CoV. Altogether, the alignment of coronaviral S glycoproteins with the reference sequence revealed a high evolutionary relationship between SARS-CoV (Urbani), bat CoVs, and pangolin CoVs. It is suggested that the emergence of the highly contagious and pandemic-causing SARS-CoV-2 is highly attributable to genome recombination or mutations of the coronavirus in animal hosts such as bats [30,31,32]. The high evolutionary relationship among coronaviruses sheds light into the development of universal vaccines using conserved epitopes.

3.2. Prediction and Screening of Conserved CTL Epitopes of S glycoprotein

The conserved CTL epitopes of S glycoprotein were first screened and predicted based on 30 coronaviral S glycoprotein sequences. The NetCTL-1.2 of DTU Health Tech provides high sensitivity and specificity among the publicly available bioinformatics tools [33,34]. This web-based bioinformatics tool utilizes a combination of predictive algorithms including proteasomal cleavage, TAP transport efficiency, and MHC class-I affinity to acquire highly probable CTL epitopes in a given sequence. Given the easy accessibility, the epitopes were selected based on the available human leucocyte antigen (HLA) class-I supertypes provided by the algorithms, such as A1, A2, A3, A24, A26, B7, B8, B27, B39, B44, B58, and B62 supertypes.

HLA class-I is known to be responsible for presenting processed antigens to T-cell receptors. Generally, there are three classical HLA-class-I encoding genes (HLA-A, HLA-B, and HLA-C) and all of them are extremely polymorphic. The number of identified HLA alleles has grown exponentially over the past decades and is likely to increase with time. To date, there are over 36,000 sequences of highly curated HLA alleles deposited into the IPD-IMGT/HLA Database (https://www.ebi.ac.uk/ipd/imgt/hla/). Undoubtedly, the vast number of HLA alleles makes the epitope prediction significantly complex and impractical. Thus, in the mid-1990s, an allele-specific classification called HLA Supertype, in which the first 9 HLA class-I supertypes were described [35] and 3 more HLA class-I supertypes were added later hence the 12 HLA class-I supertypes in the latest update [36]. dos Santos Francisco et al. (2015) investigated HLA class-I supertype frequencies among 55 human populations and found that HLA supertypes A2, A3, B7, B27, and B44 were evenly distributed and not specific to only certain populations [37]. Half of the populations showed frequencies at 14-29% for A2, 14-32% for A3, 18-31% for B7, and 21-32% for B44. In contrast, HLA supertypes A1, A24, B58, and B62 had greater frequency variations among the studied populations. It is also worth mentioning that the A24 supertype was found at higher frequencies (40% on average) in SEA, PAC, AUS, NEA and AME meanwhile the A1 supertype had an average frequency of 21% in Africa, Europe, and Southwest Asia [37].

The prediction of conserved CTL epitopes was based on the HLA class-I supertypes to cover as many human populations as possible. The initial screening yielded 1,048,575 potential CTL epitopes that matched the 12 HLA class-I supertypes. The large number of epitopes was then streamed down based on their antigenicity, immunogenicity, and toxigenicity. The elimination was performed using VaxiJen 2.0, IEDB MHC Class I immunogenicity and ToxinPred, respectively. The remaining 2,114 epitopes were subjected to another round of screening based on the frequency of appearance in 30 coronavirus strains and 12 HLA class-I supertypes. After stringent selection and removal of redundant epitopes, 12 epitopes (Table 7) were chosen for further analysis during the epitope’s alignment step.

3.3. Prediction and Screening of Conserved HTL Epitopes

The IEDB MHC-II prediction tool was used to predict and identify HTL epitopes because it provides a remarkable performance score owing to the embedded IEDB consensus 2.22 method [38,39]. Thirty coronavirus S glycoprotein sequences were mapped to 27 most widely distributed HLA class-II alleles as described by Greenbaum et al. (2011) [40]. A total of 108,767 HTL epitopes were identified and selected based on their being in the top 20% of the consensus percentile rank, corresponding to their abilities to capture 50% of the total immune response [41].

The epitopes were subjected to antigenicity, IFN-inducing, and toxigenicity predictions. The total number of remaining peptides was 1,377, which rendered difficulties in epitope selection. Consequently, peptides that matched 50% or greater of the HLA class-II alleles and coronavirus strains were chosen. There were 52 epitopes retained (Supplementary Table S1) and they were subjected to sequence alignment with that of SARS-CoV-2-Wuhan-Hu-1 to locate their positions. All of them were highly related to SARS-CoV-2 and its variants. Among them, 3 epitopes including HTL3, HTL6 and HTL26 were also related to SARS-CoV (Urbani).

3.4. Prediction and Screening of Conserved LBL Epitopes

The identified conserved LBL epitopes represented potential antigen candidates for stimulating the humoral immune response. Generally, B-cell epitopes are divided into (i) linear and continuous or (ii) conformational and non-continuous (Figure 2.). Although the vast majority of B-cell epitopes are conformational (approximately 90%) [42,43], the prediction of conformational B-cell epitopes is not as established as the LBL epitopes. Thus, the LBL epitope prediction has gained the most attention, especially in epitope-based vaccine development.

In contrast to the conserved HTL and CTL epitope predictions, screening and prediction of LBL are exclusive of HLA class-I and -II alleles. This is because B cells recognize antigens via B-cell receptors (BCR), known as membrane-bound immunoglobulins (Ig). Immunoglobulins consist of a constant fragment (Fc) region at the stalk and a variable (V) domain at the top. Given its functions in antigen binding, the V domain is responsible for the enormous theoretical diversity (10¹³⁻¹⁵) of the BCR repertoire [44,45]. Despite the high plasticity and diversity of the BCR repertoire, several lines of evidence demonstrated high frequencies of shared BCR clonotypes or elements in human BCR [46,47].

On this account, the conserved LBL epitopes were screened and identified using the ABCPred and SVMTriP prediction tools. These tools are considerably accurate in their predictions [48,49] and the results are analyzable. A total of 4,238 peptides with thresholds of 0.5 and greater were obtained after eliminating the duplicated sequences. The antigenicity and toxicity predictions were performed as described earlier to exclude non-antigenic and toxigenic epitopes, leaving only 621 peptide sequences for further analysis. Subsequently, the number of epitopes was narrowed to 68 by retaining the peptides similar to the sequences or found in 50% or greater of the coronavirus strains (Table 8.).

3.5. Alignment and Assembly of the Identified HTL, CTL T, and LBL Epitopes

In this study, a total of 131 epitopes (12 CTL epitopes, 52 HTL epitopes, and 68 LBL epitopes) were identified. Generally, the selection criteria included (i) within the S1 or S2 region, (ii) conserved regions, and (iii) matching most HLA class-I and -II supertypes. The identified epitopes were aligned to the SARS-CoV-2-Wuhan-Hu-1 S glycoprotein to identify their positions (Supplementary Figure S4). They were then assembled and combined into 2 peptide sequences encompassing 39 and 34 amino acid residues, respectively (Table 9.). Interestingly, the sequences of both assemblies, i.e., Epi1 and Epi2 corresponded to the epitopes located within the S1 region of the S glycoprotein. Epi1 was located at the N-terminal domain (NTD) (S_256-294) while Epi2 was found in the RBD (S_492-525). In addition, they were relatively conserved among the pangolin and bat coronaviral S glycoproteins (Table 3.) except for avian coronavirus strains. Notably, Epi1 was 66.7% (26/39) and 43.6% (17/39) similar to that of SARS-CoV and MERS-CoV, respectively (Supplementary Figure S3), suggesting that the conservation of Epi1 renders it a promising candidate of a broad-spectrum, cross-protective vaccine that potentially offers prophylactic protection against multiple coronavirus strains.

Next, the allergenicity of Epi1 and Epi2 was determined by using AllergenFP v1.0. The results showed that Epi1 was a potential allergen whereas Epi2 was a non-allergen. Epi1 had the highest Tanimoto similarity index to a major allergen Pru av 1 (UniProtKB/Swiss-Prot ID: O24248) that causes birch pollinosis and oral allergy in patients allergic to cherry. Epi1 shared 9 amino acids with that of Pru av 1 peptide albeit scattered throughout the latter’s peptide sequence (Supplementary Figure S5). Nonetheless, the overlapping residues in the Epi1 and Pru av 1 peptides do not coincide with the known IgE binding regions of the Pru av 1 peptide, i.e., the P-loop region (44LEGDGGPGT52) [50,51], suggesting that the allergenic potential of Epi1 is mostly negligible. Furthermore, due to its relatively smaller molecular size, the 3D structure and physiochemical properties of Epi1 are not definitive, thus it is inconclusive to serve as an allergen after modification for vaccine development.

The assembled peptides matched a wide range of HLA class-I and HLA class-II supertypes (Table 9.). In terms of the HLA class-I supertypes, the peptides matched 7 out of 12 (58.3%) HLA class-I supertypes mainly the supertypes A and B, which are globally distributed in all human populations [52]. Of the 27 available HLA class-II supertypes, Epi1 and Epi2 matched 23 (85.2%) and 20 (74.1%) HLA class-II alleles, respectively. In terms of the conserved LBL epitopes, Epi1 and Epi2 matched 3 and 2 LBL epitope sequences, respectively.

3.6. Identification of the Locations of the Conserved Epitopes in Coronaviral S glycoprotein

As mentioned previously, the Epi1 and Epi2 are located at the NTD (S_256-294) and RBD (S_492-525), respectively (Figure 3a). Figure 3. indicates the close-state and open-state positions of Epi1 and Epi2 on the S glycoprotein.

The position of epitopes on an antigen contributes to its antigenicity and immunogenicity. The findings showed that Epi1 (cyan) was slightly embedded within the NTD (grey) and, therefore, was relatively less exposed than Epi2 (orange) located within the RBD domain (S_319-541). This justifies the conservation of Epi1. Nonetheless, the conservation of epitopes is not solely determined by the exposure to immune cells or antibodies, it also depends on the functional importance of the epitopes. Mutations in conserved epitopes possibly disrupt key processes, such as viral attachment, entry, and immune evasion, thereby compromising the viral infectivity and replication in host cells [53,54,55]. In this light, conserved epitopes are important to ensure the structural and functional integrity of viral particles.

4. Discussion

The 2019 SARS-CoV-2 pandemic revealed the unpreparedness of global healthcare systems to effectively respond to such a crisis, eventually leading to the breakdown of healthcare systems. The evolution and natural selection of coronaviruses are believed to contribute to the emergence of various VOCs with exceptional abilities to escape vaccine- and infection-induced immunity. The COVID-19 prophylactic vaccines are mainly based on the whole S glycoprotein subunit of SARS-CoV-2-Wuhan-1 due to its high antigenicity and immunogenicity [56,57,58,59,60,61,62]. However, the effectiveness of the COVID-19 vaccines is becoming less pronounced following the existence of immune-evading variants due to the perpetual gene mutations [63,64,65,66]. In addition, the immune imprinting induced by immunization and previous infections also reduces the efficaciousness of the vaccines against newly emerged variants [67,68]. In this light, a conserved multi-epitope approach has been adopted to develop pre-emptive vaccines against highly mutable coronaviruses by targeting the critical functional viral antigens. This strategy not only induces broad and long-lasting immune responses, it also prevents the constant review of vaccine formulations due to viral mutations. To achieve this, epitope identification and characterization are entailed to generate epitope maps depicting their antibody specificities in-silico prior to rigorous in-vitro and in-vivo empirical investigations [13,14,18,19,20,21,69,70,71].

Phylogenetically related zoonotic coronaviruses including distantly related avian coronaviruses were included in this study to identify and analyze the conserved regions. A significant genetic divergence was observed in avian and human coronaviruses compared to SARS-CoV-2. This divergence is especially notable when comparing avian coronaviruses, which belong to the Gammacoronavirus genus with human coronaviruses and SARS-CoV-2, which belong to the Alphacoronavirus and Betacoronavirus genera, respectively. This observation is consistent with the phylogenetic data reported by Gilbert & Tengs (2021) [72]. It is noteworthy that none of the avian coronaviruses has been reported to infect humans to date. In this light, the prediction of conserved epitopes of coronaviruses prioritizes those and their close zoonotic counterparts causing diseases in humans (hCoVs, MERS-CoV, SARS-CoV, and SARS-CoV-2).

Initially, 12 HTL epitopes, 52 CTL epitopes, and 68 LBL epitopes were identified and the majority of them were conserved across the coronavirus strains. The avian coronaviruses, hCoVs, and MERS-CoV were distantly related to SARS-CoV, SARS-CoV-2, bat CoVs, and pangolin CoVs. The evolutionary convergence among those coronaviruses is likely due to the different natural and/or intermediate hosts [73]. In addition, it is noteworthy that the evolutionary convergence also results in the host-cell receptor variations as observed in hCoVs of which the surface receptors responsible for viral adsorption are mainly surface peptidases and sialic acid-rich glycan-based receptors [74].

Many of the identified epitope sequences overlapped with one another, therefore, they were aligned and assembled into single peptide sequences. Two peptide assemblies, Epi1 and Epi2, were produced, consisting of HTL, CTL, and LBL epitopes. Those sequences represented residues S₂₅₆-₂₉₄ (SGWTAGAAAYYVGYLQPRTFLLKYNENGTITDAVDCALD) and S₄₉₂-₅₂₅ (LQSYGFQPTNGVGYQPYRVVVLSFELLHAPATVC) of the S glycoprotein. Epi1 was located in the S1 region, particularly the NTD; Epi2, on the other hand, is found in the RBD region. Coincidentally, Epi1 was similar to the peptide sequence reported by Meyer et al. (2023) [17]. The study reported that the S₃₆₉-₂₇₇ epitope could induce a high magnitude of CTL immune response after in-vitro stimulation [17]. Furthermore, an in-vivo study confirmed that the Epi1 sequence (YYVGYLQPRTFLLKY) could induce a robust antigen-specific IFN-γ-producing CTL response [15] meanwhile an in-silico study identified the SGWTAGAAAYYV motif of Epi1 as the immunodominant site for T-cell and humoral responses [12]. Collectively, Epi1 is a robust candidate for the development of multi-epitope vaccines. Epi2 (S₄₉₂-₅₂₅), on the other hand, has not been reported in any previous studies but its PYRVVVLSF motif was hypothesized to induce adaptive immunity such as the production of neutralizing antibodies [11].

Human leukocyte antigen (HLA) alleles are among the most gene-dense and polymorphic regions in the human genome [75]. HLA molecules are responsible for antigen presentation to T-cell receptors (TcR) on CTL and HTL, and therefore, can readily affect the vaccine-induced immune response [34,76,77,78,79,80,81,82,83,84]. In this light, it is important to retain antigenic epitopes that can interact and bind to HLA class-I and class-II molecules in vaccine design and development. Epi1 and Epi2 fulfill the characteristics of immunogenic vaccine candidates given their multiple T-cell and B-cell epitopes. Among the matched HLA class-I supertypes, the A*02 supertype in both Epi1 and Epi2 is prevalently found in almost all human populations [52]. In regard to the HLA class-II supertypes, more than 20 HLA class-II supertypes were identified, together with the LBL epitopes, Epi1 and Epi2 are expected to trigger cellular and humoral immune responses, thereby providing more comprehensive protection against coronaviruses [43,85].

Given the conservation of Epi1 and Epi2, they hold promises as lead antigens in universal multi-epitope vaccine development, particularly in fighting the upcoming mutants. This helps address issues concerning the constant gene mutations and immune evasion seen in coronaviruses. Epi2 consists of most of the important residues required to form tight binding with ACE2 receptors [86,87,88]. It also encompasses well-known mutation sites found in the currently circulating Omicron variant, i.e., N501 and Y505. The N501Y mutation can lower neutralizing antibody binding in-vitro [89], while the Y505H mutation reduces viral protein stability, affects viral infectivity and promotes immune evasion [90,91]. Given the importance of the mutations, including them in a vaccine formulation is likely to add to the relevance of the multi-epitope vaccine with the circulating coronavirus variants hence greater immune protection.

Incorporating multiple epitopes in a vaccine formulation can offer broader and more durable protection against a wider range of viral variants. The multi-epitope sequences identified in this study shed light on the ongoing development and applications of coronavirus vaccines. To better model the binding affinity and stability between the epitopes with their counterparts such as TLR 3 and TLR4 receptors, docking analysis and dynamics simulation can be employed, thereby improving their efficacy in in vitro and in vivo studies [92]. Additionally, the epitopes also can be incorporated into vaccine formulations by conjugating them with virus-like particles (VLPs) or nanoparticle-based delivery systems to enhance their abilities in inducing humoral and cellular responses [93,94]. Furthermore, the epitopes can be developed into multivalent vaccines consisting of promising flu antigens such as nucleoprotein (NP) of influenza A virus (IAV), which assembles into virus-like particles (VLP) for vaccine delivery [23,95]. To strengthen the multivalency of the vaccine, highly conserved matrix 2 ectodomain protein (M2e) of IAV can be added into the vaccine formulation. The IAV M2e is known for conferring partial protection in animal models against IAV [96,97]. Collectively, the aforementioned prospective applications highlight the versatility and the potential of the multi-epitope peptides in curbing coronavirus infections.

This study revealed the evolutionarily conserved regions in SARS-CoV-2, SARS-CoV, and some animal coronaviruses while highlighting the genetic divergence in MERS-CoV and hCoVs that exhibited nearly no sequence similarity in the S1 subunit. This finding underscores a distinctive genetic difference among coronaviruses. Despite the genetic variations, the HR1 and HR2 regions of the S2 subunit displayed some degree of sequence similarities across the coronaviruses, indicating their potentials as the therapeutic targets.

5. Conclusions

In conclusion, the evolutionarily conserved epitopes are present among animal and human coronaviral S glycoproteins. Overall, 132 candidates representing HTL, CTL and LBL epitopes with relatively low evolutionary divergence were identified. They were screened and filtered into two final peptide assemblies: Epi1 is composed of 4 HLA class-I, 5 HLA class-II, and 3 LBL epitopes, meanwhile, Epi2 consists of 2 HLA class-I, 8 HLA class-II, and 2 LBL epitopes. These two peptide sequences, located within the S1 subunit, retain high population coverage and conservation properties, hence broad applicability and high effectiveness as lead universal vaccine candidates. Notably, Epi1 also contains immunodominant CTL epitopes, which adds to its potential as a vaccine candidate. Collectively, the conserved epitopes provide a robust foundation for vaccine development. They are expected to stimulate broad-spectrum immunity to mitigate the impact of infections of coronaviruses.