Next generation sequencing shows diversity of Omicron sub-lineages of SARS-COV2 circulating in Jeddah, Saudi Arabia

2.1 Sample collection

Procurement of samples for whole genome sequencing was conducted at Alborg Laboratories, Jeddah, Saudi Arabia between February and April 2022. The inclusion criterion was dependent on RT-PCR positive specimens for SARS-CoV-2 with a Ct value less than 30. A Sum of 96 specimens was chosen of which residual nasopharyngeal swab specimens were used for downstream viral RNA extraction. The Illumina CovidSeq sequencing kit was used for the sequencing. (Bhoyar et al., 2021).

2.2 RNA extraction and sequencing

The automated Abbott M1000 nucleic acid extraction instrument employing the QIAamp Nucleic Acid Extraction Kit (Qiagen) was used. The amplicon-based Illumina CovidSeq standard protocol was utilized. The Illumina MiSeq instrument was used to sequence libraries using the 600-cycle v3 MiSeq Reagent kit.

2.3 Genome assembly, phylogenetic clustering, classification of lineage and mutation analysis

Raw fastq paired-end reads were submitted for assembly using Exatype pipeline version 1.7.12 for NGS output SARS-CoV-2 pipeline v1.7.12 that applies quality control and mapping of read to reference genome of isolate Wuhan-Hu-1 NC_045512.2 ( https://www.sars-cov-2.exatype.com). All assembled reads were then saved in the GISAID reference database ( https://www.gisaid.org/), and corresponding accession numbers were provided in supplementary Table S1. Consensus sequences were then submitted to Nextclade web tool v2.5 ( https://clades.nextstrain.org/) for dynamic monitoring of SARS-CoV-2 evolution for phylogenetic analysis, clade, and lineage assignment (Aksamentov et al., 2021). Utilizing Nextclade v2.5 and comparing against the Wuhan-Hu-1 wild type, mutation analysis was performed. The phylogenetic tree for the sequences (n-95) was generated and visualized on Auspice v2.39.0.

2.4 3D modeling of Omicron Nucleocapsid protein

The Uniprot database (Uniprot ID: P0DTC9) was used to obtain the reference amino acid sequence. Omicron-specific mutations in the generated sequencing data were introduced to the sequence as follows amino acid substitution S413R, G204R, R203K, and P13L along with the deletions E31-, R32-, and S33-. The i-TASSER server was used for structural prediction of the Omicron Nucleocapsid (Yang and Zhang, 2015). I-TASSER is a protein structure prediction method that uses a combination of threading, fragment assembly, and refinement techniques to generate high-quality 3D models of proteins from their amino acid sequences. It first identifies structural templates from the Protein Data Bank (PDB) and then assembles fragments of known protein structures to construct a full-length model of the query protein. The assembled model is then refined using various techniques to improve its overall structure and energy. The predicted structure of the mutant with the best C-score, TM-score, and RMSD was aligned with the SARS-COV2 Nucleocapsid protein (PDB ID: 8f5d) using PyMol. Protein structural alignment using PyMOL involves loading protein structures (PDB files), selecting reference and mobile structures, performing sequence and structural alignment (align command), refining the alignment (align command), visualizing the results (color and transparency adjustments), analyzing alignment metrics (RMSD, TM-score, GDT-TS), and saving the aligned structures (PDB format) (“The PyMOL Molecular Graphics System, Version 2.0, Schrödinger, LLC.,” n.d).

2.5 Molecular docking study of the interaction of SARS-COV2 nucleocapsid protein with human RIG-I protein

The open-access server of Patchdock (Schneidman-Duhovny et al., 2005) was utilized to predict and assess the interaction between SARS-COV2 Nucleocapsid protein and human RIG-I protein. The electron microscopic structure of the SARS-COV2 Nucleocapsid protein was retrieved from the RSCB Protein Data Bank (PDB ID: 8f5d) along with the crystal structure of the human RIG-I protein (PDB ID: 2qfd). Two molecular docking studies were run in parallel; the human RIG-I protein with non-mutated Nucleocapsid protein, and RIG-I protein with the predicted structure of the mutated Nucleocapsid protein harboring Omicron-specific mutations. The best poses were determined based on weighted scores. PyMol was used to analyze docked structures to identify interacting residues of human RIG-I protein with SARS-COV2 Nucleocapsid protein, and the predicted Omicron Nucleocapsid protein respectively as well as the number of polar interactions.

2.6 Ethical approval.

This study was approved by the National Committee of Biomedical Ethics at King Abdul-Aziz City for Science and Technology, Saudi Arabia (HAP-02-J-099) on 13 December 2022.

3.1 3D modeling of Omicron Nucleocapsid protein

Predictive structural modeling of the Omicron-specific Nucleocapsid protein relied on incorporating high-frequency mutations in the sequencing output P13L (96.8 %), R203K (98.9 %), G204R (98.9 %), as well as S413R (70 %), and the deletions E31-, R32-, and S33- (100 %) were introduced to the SARS-COV2 Nucleocapsid protein amino acid sequence. The mutated sequence submitted to the i-TASSER server generated several structures, and the selected structure exhibited the best estimated TM-Score of 0.85 ± 0.08, C-Score of 0.98, and estimated RMSD of 4.8 ± 3.1 Å. The PyMol alignment showed that the two protein structures were not identical thus the Omicron alterations have impacted the folding of the Nucleocapsid protein mainly in the loop regions (Fig. 6) (see Fig. 7).

Close

Fig. 6

(A) EM structure of SARS-COV2 Nucleocapsid protein (PDB ID 8f5d). (B) Predicted structure of Omicron Nucleocapsid protein (modeled by i-TASSER server). (C) The PyMol alignment between the SARS-COV2 Nucleocapsid protein (red color) and the Omicron Nucleocapsid protein (modeled by I-TASSER server) (blue color) showed that the two structures were not identical. The non-aligned regions are in yellow and indicated by the arrows.

Export to PPT

Close

Fig. 7

Molecular docking model of Nucleocapsid protein with human RIG-I protein as predicted by Patchdock.

Export to PPT

3.2 Molecular docking study of interaction of SARS-COV2 nucleocapsid protein with human RIG-I protein

Patchdock was used to conduct comparative molecular docking to predict how both versions of the Nucleocapsid interacted with the human RIG-I protein that was experimentally demonstrated to interact with SARS-COV2 (Chen et al., 2020). The models exhibiting best-weightedhted scores were selected for further analysis (Table 1).

To further analyze the docked structures, PyMol was used to visualize both interacting models. For the SARS-COV2 Nucleocapsid - RIG-I interaction a total of 15 RIG-I residues were shown to interact with SARS-COV2 Nucleocapsid with a total of 26 polar interactions. On the other hand, Omicron Nucleocapsid - RIG-I interaction showed a total of 8 RIG-I residues interacting with Omicron Nucleocapsid with a total of 10 polar interactions (Table 2).

SARS-COV2 Nucleocapsid was found to interact with the two N-terminal domains (CARDS) (Ser 26, Asp 63, Leu 64, Gln 70, Asp 81, and Glu 136), the loop region (Arg 191, Ser 191, Ser 193, and Ser 194) and the central RNA helicase domain (Glu 366, Glu 367, Ser 413, Asp 415, and Ser 416). On the other hand, the Omicron Nucleocapsid protein was found to exhibit a lesser number of interacting RIG-I protein residues spread across the loop region (Ser 191, Arg 192, Lys 230, Gly 233, and Gly 235) between the N-terminal domains (CARDS) and the central RNA helicase domain (Thr 260, Ala 261, and Gln 300). Interestingly, the Omicron Nucleocapsid showed no interaction with the (CARDS).

4.1 Mutational features of Omicron lineages

The 21K Omicron clade was designated a variant of concern since it has 62 defining non-synonymous mutations including 36 mutations of the Spike (S) gene. The 69/70 (Spike H69- and V70-) deletion in this variant causes S gene target failure in which the S-assay within the TaqPath test yields a negative result (A. Li et al., 2022) hence dubbed the “stealth” Omicron since it complicates surveillance. While the 3 amino-acid insertion in Spike of 'EPE' at position 214 was identified as an insertion hotspot in the N-terminal domain region of the Spike Protein (Gerdol et al., 2022). However, the insertion was detected in 15 out 25 Omicron 21K clade samples (BA.1 = 1 and BA.1.1 = 14).

The 21L Omicron clade shares 38 amino-acid mutations with 21K including 21 shared mutations in the S gene. While 21L Omicron lacks the defining 69/70 deletion found in 21K, it carries six additional S gene mutations, e.g., Spike: T19I, V213G, S371F, T376A, D405N, and R408S. Additionally, a 9- nucleotide deletion was reported at position 21633–21641 leading to three deletions and a substitution Spike: L24-, P25-, P26-, and A27S as described (Jung et al., 2022).

Higher affinity binding of the spike protein RBD with human ACE2 receptor has been linked to increased transmission (Zahradník et al., 2021). Considering the RBD region of the Spike (319–541), there are 14 mutations in the sequenced samples from Jeddah: Spike: G339D, S371F, S373P, S375F, D405N, R408S, K417N, N440K, S477N, T478K, E484A, Q493R, Q498R, and N501Y. Spike G339D, K440N, K493Q, and R498Q were predicted to be unfavorable for Omicron RBD binding with human ACE2 in comparison to the wild-type virus, while the A484E substitution was predicted to be favorable for the binding of Omicron RBD with human ACE2 by computational molecular dynamics simulations of ACE2 interactions with Omicron variant spike protein (Rath et al., 2022). According to Zahradnk et al. (2002), the substitutions Q498R and N501Y together have a greater affinity for the human ACE2 receptor. Additionally, a cluster of mutations at the furin cleavage site with the alterations H655Y, N679K, and P681H was linked to increased Transmission (Gong et al., 2021). R203K and G204R mutations in the Nucleocapsid (N) gene were also associated with greater virus loads. (Jung et al., 2022).

Several experimental studies have assessed the immune evasion characteristics of Omicron variants and concluded that three clusters of Spike substitutions are responsible of Omicron variants and concluded that three Spike substitutions cluster were linked with increased resistance to vaccine-induced humoral immunity and neutralizing antibodies: cluster one (S371L, S373P, and S375F), cluster two (N440K, and G446S), and cluster three (Q493, G496, Q498, and Y505H), in addition to previously characterized Spike amino acid substitutions K417N, S477N, T478K, E484A, and N501Y (Cao et al., 2022), (Cameroni et al., 2022). One of the frequently mutated genomic regions is ORF1a in which deletion of 3 amino acids at ORF1a L3674-, ORF1a S3675-, and ORF1a G3676- (corresponding to an NSP6 deletion 105–107) is speculated to facilitate escape from the innate immune response by impairing viral degradation capability of infected cells (Benvenuto et al., 2020). Subsequent investigations further demonstrated interaction/modulation of interferon signaling pathways with several SARS-COV2 proteins at multiple levels, thus SARS-CoV-2 escapes interferon suppression through interacting with pattern-recognition receptors and modulating antiviral pathways by, which in turn could lead to severe disease pathogenesis (Q. Liu et al., 2022a). Certainly, the most important characteristic of the Omicron variants is their immune evasion capability in both previously recovered COVID-19 patients as well as those vaccinated.

Further discussion of potential phenotypic consequences for high prevalence mutations (more than 65 %) in the study population is in Table 3.

Based on GISAID reports, the Omicron (BA.1 & BA.1.1) were initially identified in Saudi Arabia in December 2021. We demonstrated the value of phylogenetic techniques in tracking the evolutionary changes of several lineages and sub-lineages of the Omicron variation in Jeddah, Saudi Arabia, through time in our genomic surveillance research. Our findings clearly demonstrated the change from 21K to 21L Omicron, with sub-lineages obviously co-circulating in the studied population. Whole genome sequencing for surveillance purposes can help identify mutations that may drive viral evolution and immune evasion. The binding affinity with the hACE2 receptor is still strong in spite of many alterations in the spike protein (Pascarella et al., 2022). Certainly, animal investigations of these lineages and sub-lineages are necessary to further evaluate this conjecture about phylogenetic and genomic analyses. Animal models have played a vital role in understanding the progression and variant diversity of COVID-19. Non-human primates, genetically modified mice, Syrian hamsters, and ferrets have been employed to study the virus's interaction with the host immune system, evaluate the efficacy of vaccines and therapeutics, and investigate the emergence and evolution of viral variants. Each model offers unique advantages and limitations, but collectively they have provided valuable insights into the pathogenesis of COVID-19. Animal models have allowed researchers to characterize the viral life cycle, identify key host factors involved in infection, investigate the immune response to SARS-CoV-2, evaluate the efficacy of vaccines and therapeutics, and study the emergence and evolution of viral variants. While animal models have been crucial for understanding COVID-19, it is important to recognize their limitations. No animal model perfectly replicates the human disease, and there can be significant differences in the immune response, disease progression, and susceptibility to different variants between humans and animals (Muñoz-Fontela et al., 2020).

As evident by the phylogenetic tree, over the course of our study, the Omicron lineage has evolved, and the successive variants and their sub-variants appear to become more transmissible and potentially immune-evasive. Therefore, continued monitoring programs are important for early warning of variants of concern that may have a deleterious impact on existing vaccines and population immunity. This is crucially significant in relation to Jeddah city, a major travel and religious tourism hub and a direct access point to Makkah. In addition to year-round visitors, during Hajj season almost 2 million people converge in Makkah from more than 185 countries in mass gatherings for four to five days for the annual pilgrimage which poses a critical challenge to preventing the spread and import/export of SARS-COV2 variants of concern (Badur et al., 2022).

4.2 Predicted impact of Omicron Nucleocapsid protein mutations on interaction with interferon signaling

The nucleocapsid (N) protein of SARS-CoV-2 is a versatile protein that plays a critical role in the viral life cycle. In addition to encapsidating the viral genome, the N protein interacts with various host cell proteins to manipulate the host immune response. One of the key ways in which the N protein evades the host immune response is by disrupting the interferon (IFN) signaling pathway (Bai et al., 2021).

IFN signaling is a crucial component of the innate immune response to viral infection. IFNs are a group of cytokines produced by host cells in response to viral infection. IFNs bind to IFN receptors on neighboring cells, triggering the expression of interferon-stimulated genes (ISGs). ISGs encode a variety of antiviral proteins, such as protein kinases, antiviral enzymes, and cell surface proteins that inhibit viral entry.

The N protein of SARS-CoV-2 can interfere with IFN signaling at multiple levels, including by disrupting the RIG-I pathway (Chen et al., 2020). RIG-I is a sensor protein that detects viral RNA and triggers the production of IFNs (Kawai and Akira, 2008). The N protein can bind to and inhibit RIG-I, preventing it from binding to viral RNA and becoming activated. The N protein can also interact with other proteins in the RIG-I pathway, such as MAVS and TBK1, to disrupt IFN signaling. One way in which the N protein suppresses RIG-I signaling is by binding to the DExD/H domain of RIG-I (Chen et al., 2020). The DExD/H domain is the ATPase domain of RIG-I, which is essential for RIG-I to bind to viral RNA and become activated. By binding to the DExD/H domain, the N protein prevents RIG-I from binding to viral RNA and becoming activated.

The N protein can also disrupt the interaction between RIG-I and MAVS (Chen et al., 2020). MAVS is an adaptor protein that is essential for RIG-I to signal to downstream kinases. By disrupting the interaction between RIG-I and MAVS, the N protein prevents RIG-I from signaling to downstream kinases and activating the IFN signaling cascade. In addition, the N protein can inhibit the phosphorylation of IRF3 by TBK1 (Chen et al., 2020). IRF3 is a transcription factor that is essential for IFN signaling. By inhibiting the phosphorylation of IRF3, the N protein prevents IRF3 from translocating to the nucleus and activating the transcription of IFN genes.

By disrupting the RIG-I pathway, the N protein of SARS-CoV-2 helps the virus evade the host immune response and replicate efficiently. The ability of the Nucleocapsid protein to suppress RIG-I-mediated interferon production is thought to be one of the key ways in which SARS-CoV-2 evades the host immune response. By suppressing IFN production, the Nucleocapsid protein allows the virus to replicate and spread more efficiently (Q. Liu et al., 2022b).

The human RIG-I protein is composed of three main regions: two caspase recruitment domains (CARDs) at the N-terminus, a central RNA helicase domain, and a C-terminal regulatory domain (CTD). The CARDs are protein domains that recruit and activate caspase proteases, which is an important step in the initiation of the innate immune response. The helicase domain unwinds double-stranded RNA (dsRNA), which is a critical step in the recognition of viral RNA by RIG-I. The CTD regulates the activity of RIG-I and contains phosphorylation sites that can be modified by kinases and phosphatases. This phosphorylation can either activate or repress RIG-I signaling (Kawai and Akira, 2008).

In this study, we sought to predict the impact of Omicron-specific mutations in the Nucleocapsid protein on its interaction with human RIG-I protein and to infer the potential impact on its ability to evade the immune response through suppression of IFN production.

Molecular docking showed that Omicron-specific Nucleocapsid protein exhibited reduced overall interaction with RIG-I with a total of 8 amino acid residues and 10 polar interactions in comparison to the SARS-COV2 Nucleocapsid protein which showed 15 interacting amino acid residues and a total of 26 polar interaction. Thus, it can be postulated from this model that Omicron-specific Nucleocapsid mutations P13L, R203K, G204R, as well as S413R, and the deletions E31-, R32-, and S33- have led to reduced impact on RIG-I interaction and subsequently INF production. However, this still has to be experimentally validated. While Omicron has exhibited overall increased immune evasion capability in comparison to previous SARS-COV2 variants; that can be attributed to the combination of the aforementioned mutations in the presence of Spike mutations N501Y/E484K that promote immune escape (Wu et al., 2021).

Computational methods have made significant progress in predicting protein–protein interactions (PPIs) and determining their binding energies, but they still face several limitations compared to experimental approaches. Experimental methods, such as yeast two-hybrid assays, co-immunoprecipitation, isothermal titration calorimetry, and X-ray crystallography, provide direct and accurate information about protein interactions and their binding energies. Computational methods, on the other hand, utilize various data sources and algorithms to predict PPIs and their binding energies. These methods include sequence-based methods, structure-based methods, and docking methods. Despite their advantages in terms of speed and cost-effectiveness, computational methods face limitations in accuracy, binding energy prediction, applicability, capturing dynamic interactions, and lack of experimental validation. Therefore, it is crucial to integrate computational methods with experimental approaches to gain a comprehensive understanding of protein interactions (Jessulat et al., 2011).