SARS-CoV-2 Variants Show a Gradual Declining Pathogenicity and Pro-Inflammatory Cytokine Stimulation, an Increasing Antigenic and Anti-Inflammatory Cytokine Induction, and Rising Structural Protein Instability: A Minimal Number Genome-Based Approach.

Hyper-transmissibility with decreased disease severity is a typical characteristic of the SARS-CoV-2 Omicron variant. To understand this phenomenon, we used various bioinformatics approaches to analyze randomly selected genome sequences (one each) of the Gamma, Delta, and Omicron variants submitted to NCBI from December 15 to 31, 2021. We report that the pathogenicity of SARS-CoV-2 variants decreases in the order of Wuhan > Gamma > Delta > Omicron; however, the antigenic property follows the order of Omicron > Gamma > Wuhan > Delta. The Omicron spike RBD shows lower pathogenicity but higher antigenicity than other variants. The reported decreased disease severity by the Omicron variant may be due to its decreased pro-inflammatory and IL-6 stimulation and increased IFN-γ and IL-4 induction efficacy. The mutations in the N protein are probably associated with this decreased IL-6 induction and human DDX21-mediated increased IL-4 production for Omicron. Due to the mutations, the stability of S, M, N, and E proteins decreases in the order of Omicron > Gamma > Delta > Wuhan. Although a stronger spike RBD-hACE2 binding of Omicron increases its transmissibility, the lowest stability of its spike protein makes spike RBD-hACE2 interaction weak for systemic infection and for causing severe disease. Finally, the highest instability of the Omicron E protein may also be associated with decreased viral maturation and low viral load, leading to less severe disease and faster recovery. Our findings will contribute to the understanding of the dynamics of SARS-CoV-2 variants and the management of emerging variants. This minimal genome-based method may be used for other similar viruses avoiding robust analysis.

Introduction

According to a report in the USA, the cases of the Omicron (B.1.1.529) variant of SARS-CoV-2 have risen significantly compared to the Delta (B.1.617.2) or pre-Delta variants [1]. The infection rate is particularly high in the age group of 18–49 years, and Omicron is highly transmissible even among fully vaccinated adults, where it is 2.7 to 3.7 times as infectious as the Delta variant. However, the death rate is low in the case of Omicron compared to the other variants [1]. Similar findings of high transmission rate and decreased disease severity were also reported in other countries, including South Africa, where the Omicron variant was first reported [2, 3]. The reported death rate in South Africa from Wuhan, Delta, and Omicron are 19.7%, 29.1%, and 2.7%, respectively [4]. A similar low death rate was also seen in the USA and other countries [1]. However, the cause of the high transmission rate and decreased disease severity of the Omicron variant is still not fully understood.

Reports suggest that Omicron can escape from antibody neutralization and vaccine protection due to mutations in its spike (S) protein, responsible for high transmission, and the attenuated replication of Omicron is associated with decreased disease severity and death [5, 6]. In most cases, the spike (S) receptor-binding domain (RBD) and angiotensin-converting enzyme 2 (hACE2) interaction are the focus of studies targeted at understanding the increased transmissibility of the Omicron variant, where some researchers have shown that the Omicron RBD strongly binds to hACE2 [7-9]. However, some other reports suggest a low affinity of the Omicron RBD to hACE2 [10, 11]. Therefore, further research is needed to understand the strength of the Omicron RBD–hACE2 interaction and the stability of this complex.

To understand the features of Omicron that allow it to infect a wide range of age groups, we reviewed the literature to analyze why the other pre-Omicron variants of SARS-CoV-2 are not infecting or causing severe disease in young adults. From the immunological point of view, there are many possibilities for why COVID-19 is less severe in young adults [12]. Chitinase 3-like-1 protein (CHI3L1) stimulates the expression of hACE2 and viral spike protein priming proteases (SPP) in older adults and comorbid conditions. This increased expression and availability of hACE2 correlate with increased transmission of SARS-CoV-2 and severity of COVID-19 [13]. Another report indicates that high serum concentrations of interleukin-17A (IL-17A) and type II interferon (IFN-γ) in children and young adults provide innate immunity against different SARS-CoV-2 variants [14, 15]. It is also reported that children under 15 years of age show increased expression of type III interferon (IFN-λ1) in nasopharyngeal mucosa upon SARS-CoV-2 infection that prevents virus entry into the body [16].

On the other hand, the S and nucleocapsid (N) proteins of SARS-CoV-2 induce anti-inflammatory IFN-γ production in the host upon infection [17, 18]. Furthermore, the nonstructural protein-1 (NSP1) and N protein block the type I interferon (IFN-β) induction in the host and attenuate antiviral immune responses [19]. SARS-CoV-2 ORF8 activates the IL-17 signaling pathway and induces the secretion of inflammatory factors [20, 21]. ORF8 also down-regulates MHC-Ι through lysosomal degradation and is involved in SARS-CoV-2 immune evasion [22]. Furthermore, ORF8 attenuates the IFN-γ mediated antiviral responses in COVID-19 [23]. The membrane glycoprotein (M) of SARS-CoV-2 suppresses the expression and activity of IFN-λ1 through inhibition of the RIG-I/MDA-5 pathway [24]. Therefore, a wide range of host–pathogen protein–protein interactions modulates the host immune response in SARS-CoV-2 infection and disease severity.

Since SARS-CoV-2 is a rapidly mutating virus with varying degrees of transmission and disease severity abilities, the mutations in the Omicron variant are responsible for its high transmission and reduced disease severity. In this context, we aimed to analyze the structural proteins (S, M, N, and E) and pathogenicity-associated important accessory proteins (ORF3a, ORF7a, ORF7b, ORF8, and ORF10) [21, 25] of Omicron for a better understanding of the role these proteins in providing higher transmission and decreased disease severity/death of Omicron as compared to other SARS-CoV-2 pre-Omicron variants. Additionally, our second objective was to develop a bioinformatics approach using a minimal number of genomes, avoiding robust and complex analysis, to predict various biological features of considered variants, perform comparative studies, and correlate the findings with observed pathological phenotypes to address our main objective.

Methods

Variants, Genomes, and Mutation Mapping

We considered four SARS-CoV-2 variants in our analysis. Wuhan (wild-type SARS-CoV-2) (RefSeq: NC_045512.2) and three variants, whose sequences were submitted to the NCBI GenBank between December 15 and December 31, 2021, and were randomly selected. The variants and their genomes are: Gamma (P.1) (B.1.1.28.1) accession no: MZ477769.1 (submission date: December 31, 2021), Delta (B.1.617.2) accession no: OL966477.1 (submission date: December 22, 2021), and Omicron (B.1.1.529) accession no: OL901854.1 (submission date: December 17, 2021). The genome and proteome sequences were checked for their completeness. The structural and non-structural protein sequences were identified from these genomes for each variant. Multiple sequence alignment (MSA) was performed using Clustal Omega [26], and the mutations were identified based on the RefSeq Wuhan (wild-type) sequence. Jalview V.2 [27] was used to visualize the MSA results.

Prediction of Pathogenic and Antigenic Properties

We utilized the MP3 tool [28] with default parameters to predict the pathogenicity scores using the genome sequence of each of the selected variants. Hybrid results (SVM + HMM) were considered, and the SVM scores were used in the analysis. We calculated individually the pathogenic score for each protein of interest: ORFab, full-length S protein, spike-RBD (amino acid residues 331 to 524 of S protein) [28-30], M, N, E, ORF3a, ORF7a, ORF7b, ORF8, and ORF10. The cumulative and average scores of all these 11 proteins for each variant were also calculated. For overall antigenic property analysis, amino acid sequences of each of these proteins from each variant were used, and the VaxiJen v2.0 server [31] with the threshold of 0.4 was applied. The same calculation made for pathogenicity prediction was also applied for antigenic property assessment.

Prediction of Cytokine and Interleukin Producing Peptides

We predicted pro-inflammatory inducing peptides using the PIP-EL web server with its default parameters [32]. IFN-γ-producing peptides were predicted using the IFNepitope server [33], selecting scan, motif, SVM hybrid method, and IFN-γ versus non-IFN-γ option with other default parameters. IL4pred [34] was used to predict IL-4-inducing peptides. We selected protein mapping and hybrid (SVM + Merci motif) based prediction methods and kept other parameters as default. For prediction of IL-6-inducing peptides, we considered IL-6Pred [35], selecting protein scan and random forest (RF)–based prediction methods keeping other parameters as default. Finally, we used IL17eScan [36] to identify the IL-17-inducing peptides by selecting a protein scan module and a DPC-based model. For selecting the final epitope-related calculations, we used a threshold value (score) > 0.05 for IFN-γ, IL-4, and IL-17-inducing peptides, and for IL-6 production, we used the cut-off value as > 0.3. The number of positive epitopes generated and the cumulative or average SVM scores of the total number of positive epitopes at the set cut-off for full-length S, spike-RBD, M, N, E, ORF3a, ORF7a, ORF7b, ORF8, and ORF10 proteins were used for further calculations. The total number of positive epitopes and their cumulative or average SVM scores were used for the final results for each variant.

Identification of Immune-Associated Human Proteins that Interact with ORF8, M, and N of SARS-CoV-2

The interacting human protein partners of ORF8, M, and N proteins of SARS-CoV-2 as described by Gordon et al. [37] and Enrichr-based gene set enrichment analysis (GSEA) [38] were used. To investigate the immune-modulating pathway for each of these SARS-CoV-2 proteins, we used the corresponding interacting human proteins in Enrichr. Under disease/drugs, LINCS L1000 ligand perturbations database-based results were used to provide cytokine, interleukin, and interferon-related regulation information [39]. Thirty, 15, and 47 human proteins that interact with M, N, and ORF8 of SARS-CoV-2, respectively, according to Gordon et al. [37], were used in this analysis. We identified the key integration of human proteins modulating cytokine, interleukin, and interferon-related pathways upon interacting with ORF8, M, and N of SARS-CoV-2. Furthermore, based on protein–protein docking using the HDOCK server [39], we predicted the mutational effects of ORF8, M, and N on binding activity and their possible input on regulating these cytokine-, interleukin-, and interferon-related pathways for each SARS-CoV-2 variant.

Protein 3D Structures Taken

Crystal structures of spike-RBD-hACE2 complexes for Wuhan (PDB: 6M0J), Gamma (P.1) (PDB: 7EKC), Delta (PDB: 7V8B), and Omicron (PDB: 7T9L) variants were retrieved from the PDB database (https://www.rcsb.org). Other crystal structures used are human DExD-box helicase 21 (DDX21) (PDB: 6L5N), human G3BP stress granule assembly factor 1 (G3BP1) (PDB: 4FCJ), human gamma-glutamyl hydrolase (GGH) (PDB: 1L9X), and SARS-CoV-2 ORF8 protein (PDB: 7F5F). For human stomatin (hSTOM) (UniProtKB: P27105), the AlphaFold-based structure AF-P27105-F1 available in UniProt was used. SARS-CoV-2 M and N protein sequences were identified from the specific genomes we have used. The 3D models of M and N were obtained by using the RaptorX server [40], followed by the GalaxyRefine server [41]. We also used the AlphaFold [42] and SWISS-MODEL servers [43] to model the mutant proteins. The SAVES v6.0 server (https://saves.mbi.ucla.edu) based PROCHECK tool [44] was used to analyze the stereochemical quality of the modeled protein structures, and the best models were selected based on the Ramachandran plot parameters (> 90% residues in the most favored region).

Protein–Protein Docking for Binding Strength Calculation

Static structures and web-based servers have been used in several studies to calculate the free binding energy [11, 45]. Additionally, there are several studies on spike-RBD–hACE2 interaction using molecular dynamics (MD) simulations for the SARS-CoV-2 variants [7, 46–49]. Therefore, we used a static approach, and to avoid redundancy, we did not perform MD simulation in our work.

We used the HDOCK server [39] for protein–protein docking. The UCSF Chimera X program [50] was used for 3D structure analysis and protein–protein 2D interaction maps were generated using LigPlot + v.2.2 [51]. We also used the commercial version of the Schrödinger Platform (https://www.schrodinger.com) for a second-line validation. For spike-RBD-hACE2 interaction analysis, for each variant, the crystal structure was cleaned for water molecules, ligands, ions, and crystallographic artifacts using Chimera X, and then the interactions/interacting residues between the spike-RBD and hACE2 were mapped using LigPlot + v.2.2/Schrödinger Platform. Next, we applied the HDOCK server and Schrödinger Platform and used the identified residues for protein–protein docking. The resultant best-docked complex was selected based on docking score and H-bonds to understand the binding strength.

We also used the HDOCK server and Schrödinger Platform to understand the interactions between (i) human STOM and four SARS-CoV-2 M protein variants, (ii) human G3BP1 and four SARS-CoV-2 N protein variants, (iii) human DDX21 and four SARS-CoV-2 N protein variants, and (iv) human STOM and four SARS-CoV-2 ORF8 protein variants. We performed blind docking for these interactions. The best-docked complex was selected based on docking score and H-bonds, and LigPlot + v.2.2 or Schrödinger Platform was used to identify the binding residues.

Analysis of the Effect of Mutations on Protein–Protein Interaction

Folding free energy (∆∆G)-based possible effect of SARS-CoV-2 mutations in corresponding protein stability was performed by DynaMut2 [52]. The 3D structures of the Wuhan (wild-type) variant of SARS-CoV-2 proteins was used to predict the stability of the mutant proteins of other variants.

Results

Omicron Showed Decreased Pathogenicity and Increased Antigenic Potential

Our analysis suggested a gradual decrease in pathogenicity and increased antigenicity as the new SARS-CoV-2 variants evolved. The overall pathogenic potential of Omicron was nearly 54% lower than that of the Wuhan variant, approximately 50% lower than the Gamma (P.1), and 20% of the Delta variant (Fig. 1A, C). In the context of only the structural proteins, we observed the same trend; i.e., the order of pathogenicity was Wuhan > Gamma > Delta > Omicron (Fig. 1B). The structural protein-based pathogenicity of Omicron was approximately 38% lower than that of Wuhan, 37% lower than Gamma, and 32% lower than the Delta variant (Fig. 1B, C) (Supplement Table S1). It is important to note that although the pathogenicity of spike RBD in Omicron was showing a negative result, it had very low pathogenicity. On the other hand, although the Omicron E protein showed a positive pathogenic value, it was less pathogenic compared to the E protein from other variants as per the MP3 hybrid prediction (Fig. 1A). However, these values did not affect the overall low pathogenicity of the Omicron variant (Fig. 2A) (Supplement Table S1).

Fig. 1

A Overall pathogenic properties of the four variants. B Overall pathogenic properties of structural proteins from the four variants. C The percentage of decreased pathogenicity of the four variants as compared to Omicron. D Antigenic properties of four variants.

Fig. 2

A Overall pathogenic, immunogenic, IFN, and IL induction abilities of four SARS-CoV-2 variants. B Pro-inflammatory epitope production scores of four variants. C IL-6-inducing epitope counts of four variants. D IL-17-inducing epitope counts of four variants. E Number of IFN-γ-inducing epitopes by four variants. F Number of IL-4-inducing epitopes by four variants.

We found a varied antigenic potential of individual proteins of these variants, and there was an overall slight increase in the antigenic potential in the following order: Omicron > Gamma > Wuhan > Delta (Figs. 1D and 2A). It is also important to note that the overall antigenic potential gradually increased in full-length S protein and spike RBD. We also observed a low antigenic score for Omicron ORF3a compared to other variants, and the antigenic scores of M, N, and ORF7a of Delta were low compared to other SARS-CoV-2 variants (Fig. 1D) (Supplement Table S1). Taken together, our analysis suggests a decreased pathogenicity and increased immunogenicity for the Omicron variant.

Omicron Showed Less Pro-Inflammatory and IL-6-Producing Epitopes as Compared to Other SARS-CoV-2 Variants

Next, we focused on the efficacy of pro-inflammatory cytokine and IL-6 production ability by Omicron as compared to the other SARS-CoV-2 variants. In our analysis, at the individual protein level, we observed that RBD, ORF7a, and ORF8 of Delta showed increased pro-inflammatory inducing peptide scores compared to the other SARS-CoV-2 variants. The E protein had the highest score compared to other variants, and the M protein showed a gradual decrease in the score in the order of Wuhan ≥ Gamma > Delta > Omicron (Fig. 2B) (Supplement Table S1). Studying the overall pro-inflammatory inducing peptide score for the four variants, we observed the following order: Delta > Wuhan > Gamma > Omicron (Fig. 2A, B) (Supplement Table S1). Therefore, the Omicron may produce less pro-inflammatory cytokines than the other variants.

We specifically focused on the IL-6-producing epitope counts and scores, and we observed the same trend in a very distinguished pattern. There is an overall gradual decrease in IL-6-inducing epitopes in the order: Wuhan ≥ Gamma > Delta > Omicron (Fig. 2A, C). At individual protein level analysis, this trend was observed mainly for the full-length S protein and to some extent in the case of the M protein, where the Wuhan, Gamma, and Delta may produce an equal number (n = 32) of IL-6-inducing epitopes; however, the Omicron M protein may produce a maximum of 27 epitopes (Fig. 2C) (Supplement Table S1). While the full-length S protein produces IL-6-inducing epitopes for all variants, it is interesting to note that the spike RBD of the Delta and Omicron variants may not produce IL-6-inducing epitopes. On the other hand, the N protein of Delta may generate IL-6-inducing epitopes but the N protein of other variants does not (Fig. 2C) (Supplement Table S1). Therefore, as per our results, the N protein of Delta is important for its disease severity through induction of cytokine storm, and in Omicron, a low IL-6 induction may be associated with reduced disease severity.

For IL-17-inducing epitope number and score, we did not see any differences across the variants. Only the N protein produced a total of seven IL-17-inducing epitopes with a score of 4.6 at cut-off > 0.5 for all variants; however, their positions vary due to mutations (Fig. 2A, D) (Supplement Table S1). Therefore, we presume that IL-17 may not have any specific role in disease susceptibility or severity in any of the SARS-CoV-2 variants we have used in this analysis.

Omicron Showed Increased IFN-γ and IL-4-Inducing Epitopes than Other Variants

To understand the anti-inflammatory cytokine production abilities of these four variants, we analyzed the IFN-γ- and IL-4-inducing epitopes of these four variants. The overall number and score of IFN-γ-producing epitopes by Delta were less than the other variants, and there were very small differences in these numbers and scores in the other three variants (Fig. 2A, E). At the individual protein level, although the full-length spike of Omicron generated more IFN-γ-producing epitopes, its RBD showed a lower number of IFN-γ positive epitopes as compared to the other variants. The ORF3a and ORF8 of Delta showed less IFN-γ-producing epitopes as compared to Omicron. However, the ORF7a of Delta generated more IFN-γ-producing epitopes than the other variants (Fig. 2E) (Supplement Table S1). As per our analysis, the order of overall IFN-γ induction ability was Omicron = Gamma > Wuhan > Delta. These results suggest a possible role of low IFN-γ-induction and disease severity by the Delta variant.

For IL-4, we observed an opposite trend of IL-6. There was a gradual increase of overall IL-4-producing epitopes and scores in the order: Omicron ≥ Delta > Gamma > Wuhan (Fig. 2A, F) (Supplement Table S1). S, N, and ORF3a were the three main contributing proteins to this difference. While the full-length S of Omicron may produce 58 IL-4-inducing epitopes, this number was 44, 46, and 52 for Wuhan, Gamma, and Delta, respectively (Supplement Table S1). While Wuhan, Gamma, and Omicron produced 50 epitopes, this number was 62 for Delta. For ORF3a, Delta may produce only seven IL-4-inducing epitopes, but the Wuhan, Gamma, and Omicron each produce twelve IL-4-inducing epitopes (Fig. 2A, F) (Supplement Table S1). From these data, it is also evident that IL-4 mediated anti-SARS-CoV-2 response was mainly regulated by S, N, and ORF3a, and the mutations in these three proteins were involved in immune evasion by the Delta variant and decreased severity in the case of the Omicron variant.

SARS-CoV-2 M, N, and ORF8 Protein Variants may Differentially Regulate Immune Response by Interacting with Human STOM, G3BP1, DDX21, and GGH

Next, we attempted to understand which SARS-CoV-2 proteins regulate interferon and interleukin production through host protein interactions and investigate the effects of mutations in these SARS-CoV-2 proteins in those interactions. From the literature, we found that the N protein of SARS-CoV-2 induces IFN-γ expression [17, 18], ORF8 regulates IL-17 signaling [20], IFN-γ mediates antiviral responses [23], and M protein inhibits the expression IFN-λ1 [24]. We modeled the mutant ORF8, M, and N proteins of the Omicron, Delta, and Gamma variants compared with the original Wuhan sequence.

For the M protein mutants, the LINCS L1000 ligand perturbations did not show any significant result (p > 0.05) related to the up- or down-regulation of IFN or ILs through interaction with any human proteins. However, it was found that the M protein may up-regulate IFN-γ, IFN-α, and IL-6 through interaction with human STOM (Stomatin), which regulated innate immunity (Supplement File 1: Table 1). In protein–protein docking analysis, we found the binding scores for hSTOM and M variants in this order: Delta ≥ Omicron > Wuhan = Gamma (Fig. 3A, B). Therefore, the Delta and Omicron variants showed higher IFN-γ and IL-6 production than the Wuhan and Gamma through hSTOM mediated pathway. However, since the M protein and STOM-mediated IFN-γ, IFN-α, and IL-6 production is not significantly (p > 0.05) enriched in Enrichr analysis, we should ignore this finding (Fig. 3A, B) (Supplement File 1: Table 1).

Fig. 3

The binding affinity between human STOM: M (SARS-CoV-2) variants, human G3BP1: N (SARS-CoV-2) variants, human DDX21: N (SARS-CoV-2) variants, and human GGH: ORF8 (SARS-CoV-2) variants. A Bar chart and B scatter chart.

On the other hand, the N protein of SARS-CoV-2 was found to interact with human G3BP stress granule assembly factor 1 (G3BP1) and human DExD-Box helicase 21 (DDX21) with a significant p-value (p < 0.05). The interaction of N with hG3BP1 and hDDX21 was associated with increased anti-inflammatory IL-4 production in human cells (Supplement File 1: Table 2). While we checked the binding of human G3BP1 with N protein variants, we observed that the binding strength of the Omicron N protein was comparatively lower than for the other variants. The degree of the binding strength followed the order: Gamma > Wuhan ≥ Delta > Omicron. Therefore, we presume that the N protein of the Omicron variant may not increase the IL-4 production through the hG3BP1 mediated pathway (Fig. 3A, B). In contrast, the Omicron N protein was found to bind more strongly to the human DDX21 than the N proteins from the other variants. Here, the binding strength order was: Omicron > Wuhan > Gamma > Delta (Fig. 3A, B). In Enrichr analysis, we found that ORF8 of SARS-CoV-2 interacted with the human GGH and down-regulated IL-17 and IFN-α (Supplement file 1: Table 3). Furthermore, the docking analysis suggested that mutations in ORF8 of the Delta variant slightly decreased the interaction between ORF8 and GGH. Therefore, the Delta variant may produce less IL-17 and IFN-α as compared to the other three variants (Fig. 3A, B).

Taken together, our analysis suggested that strong interaction of the Omicron N protein with human DDX21 may be associated with increased anti-inflammatory IL-4 production in cases of Omicron infection that reduced the disease severity.

Omicron Spike Showed Stronger Binding to hACE2, but its Stability is Low

To understand the increased transmutability of Omicron, we analyzed the binding affinity between spike-RBD and hACE2 as well as the stability of the various structural proteins of the four SARS-CoV-2 variants. The binding affinity with a higher negative score represents stronger binding) of the spike-RBD to hACE2 was highest for the Omicron variant compared to the other three SARS-CoV-2 variants. Both the HDOCK and Schrodinger-based PIPER pose scores showed the same trend. PIPER pose scores-based order of binding affinity was: Omicron > Delta > Wuhan > Gamma, and the HDOCK-based order was: Omicron > Gamma > Wuhan > Delta (Fig. 4A, B) (Supplement Table S2).

Fig. 4

The binding affinity of spike-RBD and hACE2 of four SARS-CoV-2 variants (A bar chart and B scatter chart). The stability (ΔΔG values) analysis of SARS-CoV-2 structural proteins (A bar chart and B scatter chart).

To understand the stability of SARS-CoV-2 proteins from the four variants, we first identified the mutations in four key structural proteins (S, E, N, and M) of the Gamma, Delta, and Omicron variants compared with the amino acid sequence of the original Wuhan (wild-type) variant (RefSeq). The MSA-based identified mutations are presented in Table 1 and Supplement File 2.

Table 1

Mutations in Structural Proteins of Gamma, Delta, and Omicron Compared to the Wild-Type (Wuhan) Variant of SARS-CoV-2.

Proteins	Gamma (P.1)	Delta	Omicron
Envelope protein (E)	NA	NA	T9I
Membrane glycoprotein (M)	NA	I82T	D3G, Q19E, A63T
Nucleocapsid protein (N)	P80R, R203K, G204R, T271I	D63G, S79I, R203M, G215C, D377Y	E31, R32, S33, P13L, R203K, G204R
Full spike (S)	L18F, P26S, D138Y, R190S, K417T, E484K, N501Y, D614G, H655Y, T1027I, V1176F	T19R, T95I, G142D, E156, F157, R158G, L452R, T478K, D614G, P681R, P812R, D950N	A67V, H69, V70, T95I, G142D, E156, F157, R158G, G339D, S371L, S373P, S375F, Q493R, G496S, Q498R, L452R, T478K, N501Y, Y505H, T547K, D614G, H655Y, N679K, P681R, N764K, D796Y, N856K, Q954H, N969K, L981F

Bold underlined residues are deleted residues

In mutation-based stability (ΔΔG values) analysis, all four structural proteins of Omicron were unstable compared to the Gamma variant and the wild-type Wuhan variant. The M protein in the Omicron variant was more unstable than in the other three SARS-CoV-2 variants (Fig. 4C, D). For the M protein, the Delta variant showed higher instability compared to the Omicron variant. However, the stability of the M protein of the Gamma variant was the same as for the Wuhan variant (Fig. 4C, D). The N protein of the Omicron variant showed nearly the same degree of instability as the Delta variant. However, the instability of the N protein from the Gamma variant was substantially lower compared to the Omicron and Delta variants (Fig. 4C, D). We observed a gradual decrease in stability in the following order: Omicron > Delta > Gamma (Fig. 4C, D). While we calculated the cumulative instability for all structural proteins, we observed that the Omicron proteins were less stable, and the mutations in the S and E proteins mainly contributed to this increased instability (Fig. 4C, D). For individual mutation-specific ΔΔG values, see Supplement Table S3.

Taken together, our results suggest that the Omicron S protein binds to the hACE2 more strongly as compared to other variants, and therefore, the transmission rate of the Omicron variant is higher. However, due to the higher instability of the S protein of the Omicron variant, the binding may not last long enough for systemic infection. Furthermore, the E protein of the Omicron variant is the most unstable, which may be associated with decreased viral particle production or maturation, resulting in a lower viral load, severity, and faster recovery in case of Omicron infections.

Discussion

It is well established that the severity of disease and death rate from the Omicron variant is comparatively lower, but the transmissibility is higher than seen for any other SARS-CoV-2 variant [1-4]. According to our analyses, we found a gradual decrease in the pathogenic properties of SARS-CoV-2 variants over time. The pathogenic properties decreased in the following order: Wuhan > Gamma > Delta > Omicron. Furthermore, in comparison to the Delta variant, the Omicron variant showed 20 and 32% lower pathogenicity at the genome and structural protein levels, respectively (Fig. 1A–C). The recent experimental findings by McMahan et al. [53] and Kawaoka et al. [54] also support our results where they have shown the pathogenicity of Omicron variants (BA.2 and BA.1) was lower as compared to early SARS-CoV-2 variants.

Recently, it was reported that the Omicron variant exhibits significant antigenic variation compared to other variants [55]. We also observed similar results in our sequence-based analysis. We found that the RBD, M, N, ORF3a, and ORF7a are the key proteins showing antigenic variation, whereas the Omicron RBD showed significantly higher antigenic properties than the other variants. The overall antigenicity found in our analysis was in the following order: Omicron > Gamma > Wuhan > Delta (Fig. 2A).

We skipped the B- and T-cell epitope analysis as our aim was not to identify peptides/epitopes based on vaccine developments or antibody evasion mechanisms here. Furthermore, several reports are available on B- and T-cell epitopes of SARS-CoV-2 variants. These previous reports support our gross findings. It has been shown that an important number of B- and T-cell epitopes in S, M, and N proteins remain unchanged among SARS-CoV-2 variants [56]. Furthermore, Omicron shows significantly conserved T-cell epitopes [57] and T-cell reactivity, however, due to host genomic diversity, the degree of reactivity varies [58]. Omicron S mutations occur in regions poorly targeted by CD4+ T cells but are more common in regions frequently targeted by CD8+ T cells [59]. Taken together with this previous B- and T-cell epitope-based findings and our overall antigenic property analysis, we propose that the Omicron variant could be the best possible attenuated vaccine candidate against COVID-19 in this current scenario.

It has been reported that the pro-inflammatory effect of the Omicron S protein is enhanced compared to other variants [60]. In our analysis, we found that the S protein may produce less pro-inflammatory and IL-6-inducing epitopes as compared to the Delta and other variants (Fig. 2B, C). Furthermore, the overall production of pro-inflammatory and IL-6-inducing epitopes are also lower in the Omicron variant than in other variants, and the order is Delta > Wuhan > Gamma > Omicron for pro-inflammatory epitopes and Wuhan ≥ Gamma > Delta > Omicron for IL-6-inducing epitopes. In addition to the S protein, the N protein plays an important role in the production of these peptides (Fig. 2B, C). On the other hand, we observed an increased anti-inflammatory IFN-γ and IL-4 induction ability of the Omicron variant compared to other variants in the following order: Omicron = Gamma > Wuhan > Delta and Omicron ≥ Delta > P.1Gamma > Wuhan, respectively (Fig. 2E, F). We also noticed that the S, N, and ORF3a proteins play an important role in IFN-γ and IL-4 induction differences in these variants. Additionally, we observed that the Omicron N protein could interact more strongly with human DDX21 than the equivalent protein of any other variant that may induce higher IL-4 production during Omicron infection (Fig. 3A, B). The human DDX21 interacts with the SARS-CoV-2 N protein [37] and induces the innate immune response in dengue virus infection [61]. Taken together, these results suggest that the milder disease severity associated with Omicron infections may be related to its increased ability to induce IFN-γ and IL-4 and reduced ability to induce pro-inflammatory cytokines and IL-6 as compared to other variants.

Spike-RBD binding to hACE2 is the key event in SARS-CoV-2 infection and transmission, and it is believed that a strong RBD-hACE2 interaction could be associated with hyper-transmissibility of the Omicron variant. In computational experiments, the binding affinity might be modified by employing either spike-RBD alone or the whole spike protein in a bound state with the hACE2 receptor. However, conflicting reports have shown both strong [7-9] and weak binding affinity [10, 11] of the Omicron RBD to hACE2 both in computational and in vitro studies.

In one of our previous reports, we have shown that the residue 493 in spike-RBD is important for hACE2 binding [46], and recently, through MD simulation (500 ns), the importance of this position has been confirmed by Socher et al. [47]. These authors have shown that the Omicron carrying the Q493R mutation in its spike-RBD exhibits higher flexibility in contact formation and stable binding with hACE2 as compared to the Wuhan or Delta spike-RBD [47]. Similarly, using MD simulation (500 ns), Lupala et al. 2022 have also reported the Omicron spike-RBD is more strongly bound to hACE2 as compared to the Wuhan variant [7]. In another study with 100 ns MD simulation, it was shown that pre-Omicron variants formed stable interaction between the spike-RBD and hACE2 in the following order of stability: Wuhan > Delta > Gamma [48]. This finding is also supported by a 100 ns MD simulation study by Cavani et al. (2022), where they showed a stronger binding affinity of the spike-RBD of the Gamma variant to hACE2 than for the Wuhan variant spike-RBD [49].

Consistent with these MD simulation-based observations, our analyses also showed that the Omicron spike-RBD had a stronger binding affinity to hACE2 than the spike-RBD of other SARS-CoV-2 variants, but that the Omicron S protein itself was less stable than the corresponding S protein of the other variants (Fig. 4C, D). Although the strong affinity of the Omicron spike-RBD to the hACE2 may generate high transmissibility, the interaction may not be sufficient for systemic infection due to poor stability of Omicron S protein leading to severe COVID-19. Attenuated replication of Omicron is associated with decreased disease severity and reduced death rates, and the E protein plays an essential role in this replication attenuation [6]. We observed that the structural S, N, M, and E proteins of the Omicron variant are less stable than the other variants, and importantly, the E protein of the Omicron variant is less stable (Fig. 4C, D). This unstable E protein may decrease the ability of replication or maturation of new viral Omicron particles, resulting in reduced viral load, disease severity, and faster recovery from Omicron infection.

Conclusions

In this study, we analyzed four genome sequences, one from each of the Gamma, Delta, Omicron, and the Wuhan (wild-type) RefSeq. We analyzed four structural and six accessory proteins from these genome sequences using various bioinformatics approaches to understand why the Omicron variant is more transmissible but causes less severe disease. Our analyses revealed many critical biological mechanisms in these aspects and showed that the SARS-CoV-2 variants lost their pathogenicity and inflammatory cytokine production ability and showed enhanced immunogenic and anti-inflammatory cytokine induction ability. Using these mechanisms, over time, through mutations in major structural and non-structural proteins, SARS-CoV-2 is trying to adapt to an attenuated co-existence with its human host. If this trend continues, emerging SARS-CoV-2 variants may show enhanced transmissibility but cause an even milder COVID-19 than the Omicron variant today. The bioinformatics strategy used in this analysis may also be useful for understanding the dynamics of other viruses and may help being prepared for any other future pandemics caused by viral or bacterial infections.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 1084 KB)

Supplementary file2 (DOCX 636 KB)

Supplementary file3 (XLSX 71 KB)

Supplementary file4 (XLSX 61 KB)

Supplementary file5 (XLSX 21 KB)