ICMR-National Institute of Cancer Prevention and Research (ICMR-NICPR) NOIDA has a High Throughput Viral Diagnostic Laboratory (HTL) facility to test SARS-CoV-2 samples from Delhi and some UP regions of India. Nasopharyngeal/oropharyngeal samples (NPS/OPS) were collected from these regions in a viral transport medium (VTM) and transferred to ICMR-NICPR. 24 specimens were included in the study, 14 specimens represented vaccination and 10 specimens were non-vaccinated (TS1). To comprehend the exposure of various SARS-CoV-2 variants, a selection of these samples was put through whole-genome sequencing (WGS). The major inclusion criteria for WGS (RT-PCR) were positive SARS-CoV-2 screening and a real-time PCR cycle threshold value (Ct value) of ≤ 30. To avoid any kind of dissemination of SARS-CoV-2 to healthcare personnel, samples were initially thermally inactivated at 56°C for 30 minutes.

The viral genomic content of SARS-CoV-2 was isolated from NPS/OPS specimens using 200 μL of the VTM sample from GB Pure Coronavirus RNA Isolation kit (Genuine Biosystem). This extraction protocol is based on magnetic beads for easy and fast isolation of RNA in less than 30 minutes. Nucleic RNA was isolated as directed by the manufacturer and eluted in 40 μl of elution buffer.

The inactivated samples underwent RT-qPCR analysis using a COVID-19 RT-qPCR kit (GENES2ME VIRALDTECT-II) in order to determine the viral load, as previously reported (11). SARS-CoV-2 containing samples/positive samples having Ct values ≤ 30 for RdRp Gene, E Gene, and N Gene ( Supplementary Figure S1 ) were subsequently used for whole genome sequencing (WGS) of SARS-CoV-2 (11).

SuperScript™ VILO™ Master Mix (TFS) was used for cDNA synthesis from each isolate as follows: 5X VILO^™ Reaction Mix 2 μl; 10X SuperScript^™ Enzyme Mix 1 μl; and DNase-treated total RNA (10-12 ng) ≤ 7 µl mixed properly in a MicroAmp™ Optical 96 well reaction plate (0.2ml), sealed and centrifuged. The sealed plate was loaded into the thermal cycler for cDNA synthesis, setting the program as follows: 42°C for 45 minutes, 85°C for 5 min and hold on 10°C. The rest of the protocol was followed as per the manufacturer’s instructions (details mentioned in supplementary information). The sequencing reads were aligned with the NCBI SARS-CoV-2 Reference Genome using Torrent Suite v. 5.18.1. Several plugins were used to learn more about the genetic differences found in the data. These included Coverage Analysis (v1.3.0.2), Variant Caller (v5.16.0.0) with the default settings of “Generic-S5/S5XL (540)-Germ Line-Low Stringency,” and COVID19AnnotateSnpEff (v1.3.0.2), a plugin made for SARS-CoV-2 that can predict the effects of a base substitution. To confirm the accuracy and consistency of the nucleotide calls, Integrative Genomic Viewer_2.14.1 (IGV) software was employed to visualize the TVC (torrent variant caller) Bam files for each sample. The genome sequences of all the SARS-CoV-2 samples were shared in GISAID and NCBI ( Supplementary Table TS2 ).

The nucleotide sequences of all samples were analyzed using the BAM (Binary Alignment Map) file in IGV software (IGV_2.16.1) for whole genome analysis. All 24 sequences were aligned to the SARS-CoV-2 reference genome (NC 045512.2) using MEGA 11 software (12). An aligned file was downloaded in the nexus format to construct a phylogeny tree by using iTOL ( Figure 1 ) (13). Additionally, lineage assignment was carried out using Pangolin (14). The multiple sequence alignment file was visualized for the phylogenetic analysis using the iTOL software. The same nexus file was used to construct the phylogenetic tree using the maximum likelihood method and visualization was carried out by the iTOL. Default settings have been used in MEGA 11 and iTOL analysis. The tree was made with appropriate parameters for optimal visualization.

Corrected nucleotide sequences were translated to amino acid sequences using the Expasy Translation tool (https://web.expasy.org/translate/).

Mutation analysis was performed for all structural proteins of SARS-CoV-2. For this, these sequences were aligned using Multalign version 5.4.1 software. The default setting of Gap weight of 12, Gap length-weight 2, Consensus level high (90%) and low (50%) was selected. The result was generated in image format (http://multalin.toulouse.inra.fr/multalin/) (15).

PredictSNP tool was used to predict the pathogenicity of all mutations (https://loschmidt.chemi.muni.cz/predictsnp1/) (16). This server contains the prediction algorithms of several programs, such as MAPP, PolyPhen1 and PolyPhen2, SIFT, SNAP, and PANTHER these were utilized to achieve a consensus pathogenicity score. Thus, it provides a high degree of accuracy due to the consensus technique.

NetNGlyc-1.0 Server Output was used for the prediction of N glycosylation site (https://services.healthtech.dtu.dk/services/NetNGlyc-1.0/) (17) with a threshold potential of 0.5. The server uses Artificial Neural Network (ANN) and predicts the site based on the NXS/T sequence. The server also calculates the glycosylation ‘potential’ which is the average of nine ANN.

The selection of alleles was done from the available literature that was based on the association of alleles with COVID-19 confirmed cases across the world (18–20). For MHC class I, the following sixteen HLA alleles A02:01, A11:01, A24:02, B07:02, B08:01, B13:02, B18:01, B35:01, B40:06, B46:01, B51:01, B52:01, C01:02, C04:01, C6:02 and 07:01 were used. Cytotoxic T cell (CTL) epitope prediction, NetMHCPan 4.1 was used, which is based on an artificial neural network (ANN). Sixteen selected HLA class I molecules were selected for analysis.

B-cell epitope was predicted using ABCpredtool https://webs.iiitd.edu.in/raghava/abcpred/index.html server (21). This prediction tool is based on an artificial neural network. Default parameters of this prediction tool; window length/epitope length (16 amino acids) and threshold of 0.51 were used for prediction.

The complex structure of spike protein of SARS CoV-2 complex with human ACE 2 (hACE2) receptor was downloaded from the protein data bank (PDB ID: 6lzg) and used for mapping the mutations (22). Similarly, the complex structure of spike protein with antibodies S309 (PDB ID: 7xsw); Complex structure of gamma P.1 variant spike protein with Fab 4A8 (PDB ID: 8dls), and complex structure of omicron variant spike protein with Fab XGv282 (PDB ID: 7we7) were downloaded from the PDB database. Using the PDB-editor stand-alone tool, complex structures were changed to focus on the surface area of the spike protein with hACE2/antibody. The mutated residues were mapped on the complex structures using PyMOL.2(https://pymol.org/2/). For changing a specific residue and identifying the bonds in the structure “mutagenesis tool” of PyMOL was used.

Protein-protein docking was used to determine the effect of mutation on the binding affinity of RBD to ACE-2 receptor and antibody against RBD. The sequences of the RBD region were obtained from the sequencing data of individual samples having mutations in the RBD region. The structure was predicted by homology-based modeling using the I-TASSER (Iterative Threading ASSEmbly Refinement) server (https://zhanggroup.org/I-TASSER/about.html). Protein-protein docking was performed by docking the predicted structure of RBD with ACE-2 receptor structures (PDB ID - 1R42) using the HDOCK server (http://hdock.phys.hust.edu.cn/) and the structure of the complex, docking scores, and confidence scores were obtained.

The interaction of RBD regions (with mutation; already predicted by I-TASSER) with antibody 4A8 (PDB ID - 8DLS) was predicted by docking using the HDOCK server; and the structure of the complex, docking scores and confidence scores were obtained.

The genome sequences that have been compiled have been formally submitted to the Global Initiative on Sharing All Influenza Data (GISAID). The corresponding accession numbers, together with their respective sub-lineages, have been mentioned in the Supplementary Table (TS-2) . A phylogenetic tree based on the WGS of all 24 samples of SARS-CoV-2 is shown in Figure 1 . The genomic coverage for all samples included in the phylogenetic tree analysis was ≥98%. In addition to the 24 samples, we evaluated the sequences of spike protein of several globally circulating variants to gain a more in-depth knowledge of viral infections ( Supplementary Figure S2 ). The lengths and branches in the cladogram represent the evolutionary relatedness of the consensus sequence and samples. These sequences were grouped in 4 clades representing Omicron sublineages as BA.1.15, BA.2, BA.2.75.2 and BA.5. Maximum genomes (75%) belong to BA.2, followed by 8.3% of BA.2.75, 8.3% BA.2.75.2, 4.2% BA.1.15 and 4.2% BA.5. The genomes grouped in BA2.75.2 were more genetically related to XBB lineage.

All structural proteins, namely spike, nucleocapsid, membrane, and envelope proteins, had a total of 48, 9, 6, and 2 mutations, respectively. The frequency of each mutation within the samples is depicted in Figure S3 ( Supplementary Material ).

Two mutations were detected in the envelope protein, namely at locations T9I (all samples) and T11A (only four samples). The aforementioned mutations have been documented to exhibit a correlation with the ion-selectivity of envelope channels and a modification in pH sensitivity. Additionally, they have been associated with a decrease in cytokine production and cell death. The altered repercussions resulting from this mutation may perhaps account for the diminished efficiency in the release of the Omicron variant and the subsequent decrease in cellular damage (23).

The glycosylation pattern of membrane proteins can be influenced by mutations occurring at D3N (sample NICPR 294) (24) also found that the presence of D3G (sample NICPR 13) at the 3-8 position can be linked to the N-myristoylation site. The In-silico investigation of mutations Q19E and A63T exhibits uncertainty, as some analyses indicate that these mutations lead to the destabilization of the membrane protein’s structure, while other analyses imply that they contribute to its stabilization.

The nucleocapsid protein encompasses nine mutations, there is a mutation P13L and deletion of ERS at positions 31, 32, and 33 which are present at the N-terminal domain of the protein. This domain helps in modulating the RNA binding, and phase separation; thus, mutation at these sites affects the RNA binding. Mutations at positions 203 & 204, were observed in all the samples present in the linker region. At position 413, due to the change in codon, from AGT to CGT, alteration in amino acid, from Serine to Arginine, was seen in all the samples except one (i.e. NICPR13). This mutation is present in the C-terminal and it has no evidential effect on dimerization (25). F307L and D343G are unique mutations in samples NICPR293 and NICPR13 respectively, for which experimentation is required to establish their functional attributes.

The Omicron variant has a maximum number of mutations among the other Variants of concern. Spike being the largest structural protein, is highly mutated among all structural proteins. Similar to earlier observations, we also obtained mutations for BA.1, BA.2, and BA.5 ( Figure 2 ).

The genomes identified as NICPR1, NICPR4, NICPR5, NICPR6, and NICPR11 exhibit a shared profile of mutations or single nucleotide polymorphisms (SNPs) in spike. Their genetic lineage corresponds to the omicron variant BA.2, with the presence of additional SNPs (T19I, A27S, V213G, S375F, T376A, D405N, R498S, and Q493R) that distinguish them from the canonical lineage. Similarly, another set of samples denoted as NICPR2, NICPR3, NICPR7, NICPR8, NICPR9, NICPR10, NICPR12, NICPR14, NICPR15, NICPR16, NICPR17, and NICPR18 shared the identical set of SNPs ( Figure 2 ). Their genetic lineage aligns with the omicron variant BA.2, characterized by the inclusion of additional SNPs. However, noteworthy differences arise in the additional SNPs identified in these specimens, aligning them more closely with the subsequent XBB genetic lineage. Most of these mutations are related to the immune escape by reducing the neutralization activity of antibodies (26). In the later part of the study, we focused on the spike protein only.

The impact and nature of single nucleotide polymorphisms (SNP) were examined using PredictSNP for 48 amino acids ( Supplementary Table, TS-3 ), which revealed that most of the mutations (i.e.,41) are neutral in nature whereas seven are predicted to be deleterious ( Tables 1 , TS-3 ). Out of the total of 48 mutations, a subset of 13 mutations were not subjected to annotation and are presented in Table 1 . Out of the seven detrimental mutations discussed in TS-3, it is worth noting that five have been identified and documented for their involvement in immune evasion and increased susceptibility to infection. Two unique mutations (W152R and G339H) were observed that were not annotated according to PredictSNP. The impact of W152R mutation (NTD region) was studied in silico and found to interact with class III antibodies ( Figure 3A ). There is the formation of extra-polar interactions due to the presence of arginine residue ( Figure 3B ). We looked at the G339H and G339D mutation of the RBD region that falls in the group E epitopes how it affects the way ACE2 binds and how the immune system does not recognize it ( Figures 3C–E ).

Mutations at the RBM of the RBD-spike have been identified as Q493R, G496S, Q498R, and Y505H. Figure 4 also shows that the positions of Q493, G496, and Y505 have reverted. G496 reversion is also observed in the Omicron BA1 lineage. Q493 reversion mutation was observed among genotypes belonging to BA2.75 (NICPR20) and BA2.751(NICPR21). Q493R mutation in Omicron was observed as a compensatory mutation (gain of salt bridge between Omicron-spike and ACE-2) to K417N (loss of salt bridge between Omicron-spike and ACE-2) when compared with Delta-spike (3). The K417N is related to immune escape by decreasing antibody binding. The Q493 reversion would result in the loss of the salt bridge between Omicron-spike and ACE-2. Y505 reversion was only observed in BA2.75 (NICPR20), the Y505 position was also involved in ACE-2 binding through a Hydrogen bond.

The docking score of the RBD domain representative sample of different Omicron variant groups to the ACE-2 receptor has been mentioned in Table 2 . A lower docking score suggests a higher affinity for binding. The docking scores are computed using our knowledge-based iterative scoring functions, namely ITScorePP or ITScorePR. A lower docking score indicates a higher likelihood of a binding model. The confidence score determines the likelihood that the two molecules will bind: more than 0.7 indicates that it is highly likely; between 0.5 and 0.7 indicates that it is feasible; and less than 0.5 indicates that it is unlikely. In all the samples except sample 16, there is a decrease in docking score as compared to WT-SARS-CoV-2 and Delta, suggesting increasing affinity for binding to ACE-2 receptors. NICPR 19 and NICPR 21 have the highest binding affinity to the receptor.

Similarly, RBD interaction with antibodies of samples NICPR1, 13, 16, and 21 are nearly the same, NICPR19 and 293 have less, and NICPR294 has the lowest docking score, indicating the lowest, moderate, and highest binding to the antibodies.

SARS-CoV-2 and target cell glycosylation both have a significant influence on SARS-CoV-2 infection at various levels. Analyses of N-glycosylation in the spike glycoproteins revealed the presence of 22 sites for glycosylation ( Supplementary Table 3 TS 3 ). Except for four sites of N-glycosylation at position 17 (NLTT), 149 (NKSW), 165 (NCTF), and 343 (NATR), the rest were conserved in all the samples ( Table 3 ). N glycosylation site at position 17 (NLTT) was only present in one sample representing BA1.15. However, there was a loss of N165 (NCTF) glycan in the sample. This glycosylation site was also absent in the delta variant and only present in BA.1 and BA1.15. A reversal (absence) of the glycosylation at position 17 was observed in BA.2, BA2.75, and BA5 sub-lineages.

We have further checked the change in glycosylation extent ( Supplementary Table 4 ) for positions 149 (NKSW) and 343 (NATR). There are some mutations at nearby residues of the N-glycosylation site that change the potential of glycosylation, such as at position 149 NKSW to NKSR, the potential of glycosylation occurrence is increased to 0.6318 to 0.6946 (nearly 10%) as compared to the WT-SARS-CoV-2 sequence. At another position, 343 NATR to NATT the potential of glycosylation decreases from 0.5671 to 0.5497. For the same position, there is an increase in the glycosylation potential to nearly 10% for sample NICPR293 which has the same residue as the WT-SARS-CoV-2 sequence NATR with some mutations before this region (8 residues before). The interaction of these residues may be responsible for the increase in glycosylation potential.

In this part of the study, spike proteins of samples NICPR13, NICPR20, NICPR21, NICPR293, and NICPR294 were considered and compared with the reference sequence of SARS-CoV-2 (Wild Type or WT SARS-CoV-2 is used for representation). A large number of epitopes of different affinity were generated upon using the NetMHCPan 4.1 prediction tool, only strong binders to cytotoxic T cells were selected for further analysis. In a similar way for helper T cells, various epitopes in spike were compiled with higher affinity/Strong HLA binders to induce the immune response (shown in Supplementary Table 6 ). The total number of obtained epitopes for CTL, HTL, and B-cell epitopes from all selected samples are summarized in Table 4 and epitopes are compiled in Supplementary Tables 5 – 7 respectively.

The majority of anticipated epitopes are conserved in the Reference SARS-CoV-2 spike, and the allele-presenting characteristics of all samples are similar. Due to changes in amino acid sequences, some epitopes have changed in amino acid composition in comparison to WT, few epitopes remained only in the WT SARS-CoV-2, and this has also led to the formation of new epitopes. All the epitopes are mentioned in the Supplementary Table 5 . However, the new epitopes or altered epitopes within the selected samples are mentioned in Table 5 . These newly formed epitopes may play a role in generating immune escape in SARS-CoV-2-infected individuals.

HTL analysis of spike proteins in these samples’ results suggests that the majority of the epitopes are shared throughout the NICPR samples and exhibit conservation to the wild-type SARS-CoV-2 ( Supplementary Table 6 ). New epitopes or altered epitopes within the selected samples are mentioned in Table 6 . However, few epitopes are unique to BA.1.15 (sample 13), epitopes YSKHTPIIV and FSRLDKVEA. Two epitopes (FVIRGNEVS & IRGNEVSQI) were observed unique to samples 20, 21, 293, and 294. YHKNNKSRM epitope was unique to sample 294. The W152R mutation or the corresponding epitope seems a novel gain of helper T cell epitope (YHKNNKSRM) might compensate for the loss of helper T cell epitope due to the R346T/N (FNATRFASV) site ( Figures 5A–C ). The mutation R346T mutation loosens the binding of spike glycoprotein with class three neutralizing antibodies ( Figure 5B ) that is also similar to R346T ( Figure 5C ) The mutation also results in the loss of a potential epitope (FNATRFASV) as evident from the epitope predictions.

In addition to cellular immunity, humoral immunity mediates pathogen elimination in an antibody-dependent manner. To know the putative epitopes, the ABCPred server was used to further scan the SARS-CoV-2 spike protein for linear B-cell epitopes selecting the default parameter. Based on the above criteria, we obtained a total of 132 epitopes. Most of the B-cell epitopes are conserved (nearly 63%), while some are only present in WT-SARS-CoV-2 (Nearly 16%) ( Supplementary Table 7 ). Mutations in some epitopes have changed the antigenic score in the determinant (shown in Table 7 ). This could alter the epitope’s interaction with the antigenic determinant, which could aid in SARS-CoV-2 variants’ immunological escape from existing therapies. Unlike CTL and HTL-predicted epitopes, no new epitopes were discovered in our sample. Some overlapping epitopes are also present in WT-SARS-CoV-2 as well as our sequences, which can associate with some alterations in antigenic/immunogenic potential and may contribute to immune evasion ( Table 7 ). An epitope profiling study performed (27) on isolated epitopes from SARS-COV-2 infected individuals, has found overlapping regions with our predicted epitope findings, such as positions 19, 25, 257, and 547. Another study performed by found a conserved epitope 407VRQIAP412 (Sample NICPR-13) epitope and highly variable peptide/region among the SARS-CoVs 473YQAGSTP479 (Samples NICPR-13, 21, 293 & 20) (28).