Comprehensive Phylogenetic and Immunogenic Analysis of SARS-CoV-2 and Related Coronaviruses

1. Introduction

Coronaviruses (CoVs) are a group of viruses that mostly infect humans and wide range of animal species [1]. They can cause different kinds of diseases, ranging from mild to severe respiratory, hepatic, and neurologic conditions [2,3]. SARS-CoV-2 is the third coronavirus that can pass from animals to humans, after SARS-CoV and MERS-CoV [4], and the seventh coronavirus that can infect humans [5]. It is the cause of the COVID-19 outbreak that emerged in December 2019 [6]. Many human coronaviruses like HCoV-229E, HCoV-NL63, SARS-CoV, and MERS-CoV are believed to have originated in bats [7,8]. Additionally, CoV-OC43 and CoV-HKU1 are thought to have passed from rodents to humans [9,10]. However, the details of how these coronaviruses evolve in bats are still unknown, and the exact pathways through which they are transmitted from natural reservoirs to humans remain to be fully elucidated.

Coronaviruses have four main types of proteins: spike (S), envelope (E), membrane (M), and nucleocapsid (N) [11]. In SARS-CoV-2, the genome is protected by the N protein, which is closely associated with the M protein, often called the matrix protein. The outer envelope of the virus is made up of the E protein [11,12]. The N protein helps the virus replicate and has two parts: the N-terminal domain and the C-terminal domain, both of which can bind to RNA on their own [13,14,15]. The E and M proteins are crucial for assembling the virus, its pathogenesis, and budding off from cells [16,17]. Meanwhile, the S proteins on the virus’s surface form trimers and are key for recognizing and entering host cells [18,19].

Molecular and phylogenetic analysis could yield an important understanding of how viruses evolve and spread [20]. Though this method may not generally provide immediate clinical solutions, it offers promise for developing targeted therapies [21,22]. From the onset of the COVID 19 pandemic till date, plenty studies have been conducted to understand the phylogenetic relationships and immunogenic potential of different coronaviruses. In one comprehensive study, Oliviera and colleagues identified conserved B and T cell epitopes in the SARS-CoV-2 nucleocapsid protein [23]. However, the study was limited to nucleocapsid protein, overlooking other important viral proteins, such as the spike protein which plays crucial role in immune response and vaccine efficacy [24,25]. In addition, failure to account for all relevant coronaviruses could limit the discovery of many potential cross-reactive epitopes. Similarly, Chen and colleagues conducted a phylogenetic analysis to trace the evolutionary origin and unique niche among human coronaviruses—SARS-CoV, SARS-CoV-2, and MERS-CoV [26]. While their findings provided insights into evolutionary history, they mainly come from phylogenetic analyses, which might not capture all the factors affecting how the virus evolves and adapts to hosts. Another recent study by Lopes used bioinformatics to explore possible cross-reaction between SARS-CoV-2 and other mammalian coronaviruses by identifying similar protein regions [27]. However, like Oliviera's work, it failed to account for all relevant viral strains.

In summary, while recent studies have shed powerful light on the phylogenetic interconnections among coronaviruses [28,29,30,31], they often focus on only a subset of viral proteins or strains [23,24,26,27]. Few studies have comprehensively addressed both the phylogenetic and immunogenic aspects of coronaviruses considering all relevant strains and functional proteins [32,33,34,35]. Furthermore, the detailed mechanisms of the evolutionary changes and the structural adaptations owing to glycosylation and receptor binding in coronaviruses remain underexplored.

In this study we aimed to address these limitations by examining all relevant coronaviruses and their protein sequences to gain insights into their ancestral relatedness, cross-reactive epitopes, and potential immune escape mechanisms. Briefly, we leveraged web-based tools to determine whether different coronaviruses structural and antigenic variations could yield insights into their evolutionary relationships and ability to escape the immune system. This study could provide crucial insight to support the further development of vaccines, therapeutics, and diagnostic tools, as well as for preparing for future outbreaks.

2. Materials and Methods

2.1. Coronavirus Sequence Data Source

The nucleotide and protein sequences of the four viral proteins were extracted in FASTA format from the Microbial Genomes and Protein Resources at the National Center for Biotechnology Information (NCBI) [https://www.ncbi.nlm.nih.gov]. Among the list are sequences for SARS-CoV-2, BatCoV BM48-3, HCoV-229E, HCoV-HKU1, Ty-BatCoV HKU4, Ty-BatCoV HKU33, HCoV-NL63, HCoV-OC43, Pangolin-CoV, BatCoV RaTG13, and Pi-BatCoV HKU5, MERS-CoV, and SARS-CoV-2 (refer to Supplementary Table 1 for details).

2.2. Antigenic Variation Prediction

To assess the variation among structural viral proteins from various coronaviruses, we performed pairwise sequence alignments using EMBOSS Needle software (version 6.6.0) available at https://www.ebi.ac.uk/Tools/psa/emboss_needle/ [36]. EMBOSS Needle aligns two input sequences globally using the Needleman-Wunsch algorithm. The alignment parameters included a gap open penalty of 10, a gap extend penalty of 0.5, and the BLOSUM62 protein alignment matrix. The software computes the optimal alignment along the entire length of the sequences, including gaps. We analyzed the percentage of gaps, similarity, and identity in the alignments, as well as the alignment length and score.

2.3. Antigenic Relationship and Structural Divergence

Amino acid sequences were aligned using ClustalW. We then used Molecular Evolutionary Genetics Analysis (MEGA) X software to construct phylogenetic trees for the spike glycoproteins, nucleocapsid, envelope, and membrane proteins [37]. The analysis employed the maximum composite likelihood method and the unweighted pair group method with arithmetic mean (UPGMA) [38,39]. This approach assessed the structural divergence of these proteins among the studied coronaviruses and evaluated their antigenic similarities.

2.4. Variation in Glycosylation Pattern of Spike Glycoproteins

Amino acid sequences of spike glycoproteins were analyzed using NetNGlyc 1.0 online software https://www.cbs.dtu.dk/services/NetNGlyc/ to predict glycosylation sites [40]. This analysis aimed to identify differences in the glycosylation patterns of spike glycoproteins, which may affect viral attachment to host cell surfaces.

2.5. Prediction of Cleavage Sites

The cleavage sites for 3CL proteinase as well as furin were predicated for the spike glycoproteins. We used NETCorona http://www.cbs.dtu.dk/services/NetCorona/ for the 3CL proteinase [41]. The NETCorona utilizes artificial neural networks to to predict cleavage sites based on known sites from seven coronavirus family members, reflecting the conservation of these sites within the family.

For furin cleavage sites, ProP 1.0 https://services.healthtech.dtu.dk/service.php?ProP-1.0 was utilized [42]. This tool predicts cleavage positions based on arginine and lysine propeptides using neural networks. Additionally, the SignalP 3.0 server, integrated into ProP 1.0, was used to predict the presence and location of signal peptide cleavage sites. Predictions were considered significant with a score greater than 0.5.

2.6. Epitope Prediction and Variation in Spike Glycoprotein

To identify and compare Cytotoxic T Lymphocyte (CTL) epitopes within the spike glycoproteins, we used the NetCTL 1.2 online tool http://www.cbs.dtu.dk/services/NetCTL/ (Larsen et al., 2007). This software predicts CTL epitopes based on three key factors: the efficiency of transporter associated with antigen processing (TAP), proteasomal cleavage patterns, and MHC class I affinity. These predictions reveal variations in epitope presentation across different spike glycoproteins, proteasomal cleavage, and MHC class I affinity [43,44].

3. Results

3.1. Phylogenetic Analysis and Pairwise Sequence Alignment of Coronavirus Proteins

The phylogenetic trees constructed from the structural viral proteins (spike, membrane, envelope and nucleocapsid) demonstrated the ancestral origin and distant evolutionary relationships of the newly emerged SARS-CoV-2 (Figure 1). It was found that SARS-CoV-2 was related to BatCoV RaTG13, Pangolin-CoV, Bat SARS CoV and BatCoV BM48-31. In all the trees, SARS-CoV-2 aligned with the same clade of BatCoV RaTG13, Pangolin-CoV, Bat SARS CoV and BaTCoV BM48-31. However, the phylogenetic analysis revealed that SARS-CoV-2, BatCoV RaTG13 and Pangolin-CoV originate from the same immediate common ancestor. Three bat-originated coronaviruses; HKU4, HKU5 and HKU33, and five human-originated coronaviruses; MERS, HKU1, OC43, 229E and NL63, were found to have divergent relationships with SARS-CoV-2. However, Ty-BatCoV HKU4 and Pi-BatCoV HKU5 showed ancestral relationships with MERS-CoV. Moreover, HKU33, of bat origin, always showed a phylogenetic relationship with the human-originated coronaviruses 229E and NL63. Three bat-originated coronaviruses (HKU4, HKU5, and HKU33) and five human-originated coronaviruses (MERS, HKU1, OC43, 229E, and NL63) were found to have divergent relationships with SARS-CoV-2, suggesting their potential for future zoonotic spillovers and emphasizing the need for continuous viral monitoring.

The results suggest that major proteins of SARS-CoV-2 were more aligned with BatCoV RaTG13 and Pangolin CoV structural proteins. Pairwise sequence alignment by EMBOSS Needle strengthened the phylogenetic relationship of SARS-CoV-2, BatCoV RaTG13 and Pangolin CoV. The identical and similarity patterns of SARS-CoV-2 proteins with their homologous proteins of other coronaviruses is indicated in Table 2. It has been revealed that SARS-CoV-2 proteins were highly similar and identical to BatCoV RaTG13 and Pangolin CoV than SARS-COV proteins. The spike, envelope, membrane and nucleocapsid proteins of SARS-CoV-2 showed approximately (89.8% and 97.3%), (100% and 100%), (97.7% and 99.1%) and (97.1% and 98.3%) nucleotide similarities with the respective protein molecules of Pangolin CoV and BatCoV RaTG13. Again, in the case of identity pattern, results showed that spike, envelope, membrane and nucleocapsid proteins of SARS-CoV-2 were (94.0% and 98.2%), (100% and 100%), (98.6% and 99.1%), and (98.1% and 98.6%) identical with the respective homologous proteins of Pangolin CoV and BatCoV RaTG13.

3.2. Antigenic Site Determination and Epitope Variation Analysis

The spike proteins of SARS-CoV-2, SARS-COV, BatCoV RaTG13, BatCoV BM48-31 and Pangolin CoV were used to determine the most antigenic sites by employing the NetCTL 1.2 epitope prediction tool. The NetCTL server, which gave a score well above the threshold value of 0.75, revealed the antigenic potential required to stimulate a protective response in host organisms (Flower et al., 2017). From the analysis, a total of 37 epitopes from the S proteins were found to be mostly antigenic in SARS-CoV-2, with almost 100% of peptides carrying more than the threshold value of the antigenic score of the NetCTL server (Table 3). Similarly, the other four coronaviruses (SARS-COV, BatCoV RaTG13, BatCoV BM48-31 and Pangolin CoV) also showed the presence of epitopic candidates in their spike proteins with values exceeding the antigenic threshold. Moreover, there were variation among the peptides present at the various epitopes in the spike proteins of the five coronaviruses including the recently emerged coronavirus SARS-CoV-2.

3.3. Comparative Analysis and Prediction of Spike Glycoprotein Glycosylation Sites

The glycosylation sites for the retrieved spike glycoproteins of the various coronaviruses were predicted (Table 4). A comparative analysis has demonstrated that 14 of SARS-COV-2 glycosylation sites can be found in other coronaviruses. Besides, three (3) novel glycosylation sites (NGTK, NXSN and NATR) were found in SARS-CoV-2. BatCoV RaTG13 had the most common sites (14) with SARS-COV-2, followed by Pangolin CoV (10), BatCoV BM48-31 (6), SARS-CoV (5), BatCoV HKU33 (5), BatCoV HKU5 and MERS-CoV (2), and CoV-229E, CoV-HKU1 and BatCoV HKU4 (1). None of the predicted glycosylation sites of CoV-OC43 and SARS-CoV-2 were common. These novel epitopes identified in SARS-CoV-2 can be critical for the development of peptide-based vaccines and should be prioritized for experimental validation.These glycosylation sites play crucial roles in viral infectivity and immune evasion, and their novel positions in SARS-CoV-2 could provide insights for antiviral drug design targeting glycosylation processes.

3.4. Predicted Cleavage Sites and Site Position of the Retrieved Spike Glycoproteins

The 3CL pro cleavage site was predicted for SARS-CoV-2, BatCoV BM48-3, BatCoV RaTG13, Pangolin CoV, SARS-COV, BatCoV HKU5, BatCoV HKU4, CoV-229E and CoV-NL63, whilst no 3CL pro cleavage site was predicted for MERS-CoV, BatCoV HKU33, CoV-HKU1 and CoV- OC43 (Table 5). By comparing the predicted cleavage site position of the other coronaviruses spike glycoproteins to SARS-COV-2 predicted cleavage sites, BatCoV RaTG13 had the nearest cleavage site position to SARS-CoV-2 cleavage site, followed by Pangolin CoV, BatCoV BM48-31 and SARS-CoV. This variation in cleavage sites can inform the development of specific inhibitors targeting these protease sites to hinder viral replication.

3.5. Predicted Furin Cleavage Site and Signal Peptide Cleavage Site

Predictive analysis reveal the presence of Furin cleavage site in SARS CoV 2, Pi-BatCoV HKU5, HCoV-NL63, HCoV-HKU1, HCoV-OC43, MERS-CoV but furin cleavage maturation in SARS CoV, BatCoV BM48-31, BatCoV RaTG13, Pangolin-CoV, HCoV-229E, Ty-BatCoV HKU4 ,Ty-BatCoV HKU33 was absent. This analysis confirms that among coronaviruses in the subgenus Sarbecovirus, Sars CoV 2, BatCoV RaTG13, Pangolin-CoV and Sars CoV, Sars Cov 2 is the only Sarbecovirusvirus with furin cleavage site SI/S2 (Figure 2A, B, C, D, E). The presence of this unique furin cleavage site in SARS-CoV-2 may account for its high infectivity and could be targeted in therapeutic strategies to mitigate viral entry into host cells.

Table 6. Summary of predicted furin cleavage site among coronaviruses.

Coronavirus	Furin Cleavage site position	sequence	Score	Signal peptide cleavage position
SARS-CoV-2	685	NSPRRAR|SV	0.620	13 and 14
BatCoV BM48-31	NONE		None	15 and 16
BatCoV RaTG13	NONE		None	13 and 14
Pangolin-CoV	NONE		None	14 and 15
Pi-BatCoV HKU5	745	TSSRVRR|AT	0.822	21 and 22
	884	TGERKYR|ST	0.507
HCoV-229E	NONE		None	21 AND 22
Ty-BatCoV HKU4	NONE		None	20 and 21
HCoV-NL63	863	LPQRNIR|SS	0.519	15 and 16
Ty-BatCoV HKU33	NONE		NONE	16 and17
MERS-CoV	751	LTPRSVR|SV	0.563	17 and 18
	887	TGSRSAR|SA	0.707
	1113	VKAQSKR|SG	0.512
HCoV-HKU1 AGW27881.1	759	SSSRRKR|RS	0.675	13 and 14
	758	SSRRKRR|SI	0.878
HCoV-OC43	757	SKNRRSR|GA	0.551	14 and 15
SARS CoV	NONE		NONE	15 and 16

Table 6. Summary of predicted furin cleavage site among coronaviruses.

Coronavirus	Furin Cleavage site position	sequence	Score	Signal peptide cleavage position
SARS-CoV-2	685	NSPRRAR|SV	0.620	13 and 14
BatCoV BM48-31	NONE		None	15 and 16
BatCoV RaTG13	NONE		None	13 and 14
Pangolin-CoV	NONE		None	14 and 15
Pi-BatCoV HKU5	745	TSSRVRR|AT	0.822	21 and 22
	884	TGERKYR|ST	0.507
HCoV-229E	NONE		None	21 AND 22
Ty-BatCoV HKU4	NONE		None	20 and 21
HCoV-NL63	863	LPQRNIR|SS	0.519	15 and 16
Ty-BatCoV HKU33	NONE		NONE	16 and17
MERS-CoV	751	LTPRSVR|SV	0.563	17 and 18
	887	TGSRSAR|SA	0.707
	1113	VKAQSKR|SG	0.512
HCoV-HKU1 AGW27881.1	759	SSSRRKR|RS	0.675	13 and 14
	758	SSRRKRR|SI	0.878
HCoV-OC43	757	SKNRRSR|GA	0.551	14 and 15
SARS CoV	NONE		NONE	15 and 16

4. Discussion

Human coronaviruses evolve from zoonotic transmission, and bats are commonly reported as reservoirs for such zoonotic viruses [45]. While COVID-19 is no longer a pandemic, understanding the genetic and structural characteristics of SARS-CoV-2 and related coronaviruses remains crucial for future preparedness. Our findings suggest that the unique features of SARS-CoV-2, including the novel glycosylation sites and distinct furin cleavage sites, may potentially contribute to the higher infectivity. This may ultimately provide valuable targets for developing broad-spectrum vaccines and antiviral therapies.

The output from the phylogenetics tree demonstrated that SARS-CoV-2, Pangolin-CoV and BatCoV RaTG13 share similar characteristics, and pairwise alignment indicated maximum similarity and identity between these variants. These findings are in accordance with previous studies that reported similar results [46,47,48,49]. Considering the diversity and similarities among coronaviruses, it is evident that frequent genetic recombination among various strains might lead to the emergence of new coronaviruses [50,51]. Like other coronaviruses, SARS-CoV-2 shares similarities with the six coronaviruses that infect humans but our phylogenetic studies suggest that SARS-CoV-2 is not a direct descendant of any of the six. However, it represents a divergence from coronaviruses that infect humans

The pairwise sequence analysis of SARS-CoV-2 spike glycoproteins with the other coronaviruses reveal a high similarity between RaTG13, Pangolin-CoV, Bat SARS CoV and BatCoV BM48-31, respectively. Although Pangolin-CoV and BatCoV RaTG13 revealed a 100% similarity in its envelope protein to SARS-CoV-2, a structural divergence was observed during the analysis of the nucleocapsid, membrane and spike proteins as BatCoV RaTG13 and Pangolin-CoV showed a similarity of 98.3% and 97.1% for nucleocapsid and 91.1% and 97.7% for membrane proteins. Accordingly, it is clear that the six human coronaviruses widely diverge from SARS-CoV-2 supporting the findings from the phylogenetic tree.

Antigenic epitope site determination analyses of coronavirus proteins were performed to determine the potential CTL epitopes that would interact efficiently with B lymphocytes to initiate an immune response against specific antigens [52]. Epitope sites were determined for the spike glycoproteins of SARS-CoV-2, BatCoV RaTG13, Pangolin-CoV, Bat SARS CoV and BatCoV BM48-31. SARS-CoV-2 was identified to have a total of 36 highly immunogenic CTL epitopes followed by BatCoV RaTG13 (38), Pangolin-CoV (41), Bat SARS CoV (40), and BatCoV BM48-31 (45) epitopes. The antigenic sites of SARS-CoV-2 were also compared with the other four coronaviruses spike glycoproteins. Higher levels of conservation among antigenic epitope sites of SARS-CoV-2, BatCoV RaTG13 and Pangolin-CoV were observed and this supports the claim that SARS-CoV-2, BatCoV RaTG13 and Pangolin-CoV evolve from a more recent ancestor [47,53]. Comparison between the antigen recognition site of SARS-CoV-2 and Bat SARS CoV spike glycoproteins showed a wide structural diversity among the two strains. Hence, it is evident that therapeutic agents used for treating severe acute respiratory syndrome diseases cannot be effective and efficient against COVID-19 [54,55]. Screening of predicted epitopes of SARS-CoV-2 spike glycoproteins with BatCoV RaTG13, Pangolin-CoV, Bat SARS CoV and BatCoV BM48-31 revealed six new CTL epitopes; XSNQVAVLY, NSFTRGVYY, NATRFASVY, STQDLFLPF, NSASFSTFK and ASFSTFKCY that may be responsible for the unique antigenic response of SARS-CoV-2. These novel CTL epitopes can be a target site for the development of peptide-based vaccine. The spike glycoprotein of SARS-CoV-2 is considered to be responsible for host specificity and is also the main receptor binding domain [56].

The SARS-CoV-2 spike glycoproteins were found to share fourteen out of its seventeen glycosylation sites with other coronaviruses, with majority of its glycosylation sites found in BatCoV RaTG13 (14 glycosylation sites) and Pangolin-CoV. Additionally, three novel glycosylation sites—NATR, NXSN and NGTK—were identified in the spike glycoprotein of SARS-CoV-2, suggesting that the virus might use a unique glycosylation mechanism to interact with its receptors. These glycosylation sites are potential targets for antiviral drugs that can interfere with viral attachment and entry into host cells. Analysis of the spike glycoproteins using the 3C-like protease tool revealed that the cleavage site for SARS-CoV-2 is 1000 amino sequence away from the start codons of the spike glycoproteins. 3C-like protease cleaves proteins to generate functional proteins such as single-stranded RNA-binding protein, RNA-dependent RNA polymerase (which is the main replicase of coronaviruses), helicase, exoribonuclease and endoribonuclease [57]. A comparative analysis of the SARS-CoV-2 cleavage site position revealed that BatCoV RaTG13 has the nearest cleavage site to SARS-CoV-2 at 998, whiles Pangolin-CoV has its site at 994. BatCoV BM48-31 and Bat SARS-CoV also have nearby cleavage sites at 988 and 972, respectively. In contrast, other human coronaviruses showed a wide variation in cleavage site position suggesting a major difference in their spike receptor binding protein and mechanism of replication from SARS-COV-2, Pangolin-CoV, BatCoV RaTG13, BatCoV BM48-31 and Bat SARS CoV. These variations might stem from gaps observed during pairwise alignment, as all these gaps appeared before the cleavage site position. Inhibiting the activity of 3C-like protease prevents replication since viral replication apparatus will not be produced [58,59]. Cleavage sites can serve as a target site for drug design [60]. Therefore, further research into designing 3C-like protease inhibitors against coronaviruses is warranted.

Furin is expressed at low level in the Golgi apparatus cells and has been shown to play a vital role in viral pathogenesis by cleaving polybasic or multi-basic sites such as those found in influenza virus subtypes H5, H7 and Mers CoV [61,62]. SARS-CoV-2 furin cleavage at the S1/S2 boundary primes the spike for an open conformation necessary for interaction to the ACE2 entrance receptor [63,64]. The presence of a unique furin cleavage site in SARS-CoV-2 at the S1/S2 boundary primes the spike for an open conformation necessary for interaction with the ACE2 receptor. This could explain the highly infectious nature of SARS-CoV-2 compared to phylogenetically closely related viruses of the Sarbecovirus subgenus. Targeting the furin cleavage site may be a promising therapeutic approach to block viral entry into host cells [62,65].

Taken together, future research should focus on validating the identified CTL epitopes and glycosylation sites through experimental studies, exploring their potential as targets for antiviral drugs and vaccines. Additionally, understanding the mechanisms behind the unique features of SARS-CoV-2, such as its furin cleavage site, can provide deeper insights into viral pathogenesis and inform the development of novel therapeutic strategies.

5. Conclusions

This study has revealed the genetic diversity among the structural proteins of SARS-CoV-2, BatCoV BM48-31, BatCoV RaTG13, Pangolin-CoV, Bat SARS CoV, Pi-BatCoV HKU5, HCoV-229E, Ty-BatCoV HKU4, HCoV-NL63, Ty-BatCoV HKU33, MERS-CoV, HCoV-HKU1, HCoV-OC43. The phylogenetic analysis of nucleocapsid, membrane, envelop and spike protein showed that SARS-CoV-2, BatCoV RaTG13 and Pangolin-CoV are closely related and evolved from an immediate common ancestor. Further characterization of spike glycoproteins using predicted epitopes, pairwise alignment, glycosylation site and cleavage site position proved that BatCoV RaTG13 is closely related to SARS-CoV-2 than the other coronaviruses. Also, SARS-CoV-2 showed six unique epitopes, three unique glycosylation sites and a unique furin cleavage site which can contribute to development of peptide-based vaccine as well as monitoring the consequences of glycosylation, enzyme cleavage and antigen binding variations during the process of SARS-CoV-2 infectivity.