Several mutations in the SpikeVOCs are related to increased binding affinity, increased infectivity, immune system escape and increased transmissibility of SARS-CoV-2 [1]. From amino acid sequences (Figure S1), we performed the multiple sequence alignment of SpikeWT, SpikeB.1.617.2, SpikeBA.2, SpikeBA.5, and SpikeBQ.1 proteins (Figure S2), which shows that SpikeB.1.617.2 presents mutations in NTD (T20R, G75R, D80Y, T95I, G142D), RBD (L452R, T478K), in the SD1 subdomain (T547K, D614G, P681R), the SD2 subdomain (D614G, H655Y), and in S2 subunit (D950N and V1264L). Conversely, Spikes of BA.2, BA.5 and BQ.1 share the following mutations when compared to SpikeWT: NTD (T20I, del25/27, A28S, G339D), RBD (S371F, S373P, S375F, T376A, D405N, R408S, K417N, S477N, T478K, E484A, Q498R, N501Y, Y505H), SD1 subdomain (D614G, H655Y, N679K, P681H, N764K, D796Y). In SpikeBA.2, residues H69, V70, N440, L452, F486 are the same as the residues in the wild type, but there is the substitution Q493R only found in this Omicron sub-variant. From the multiple alignment of amino acid sequences, we built a phylogenetic tree to investigate which proteins are more closely related and share similarities. The phylogenetic tree shows that the Delta variant is more similar to the WT than to the Omicron sub-variants (Figure S3). On the other hand, among the Omicron sub-lineages, BA.5 and BQ.1 share more similarities with each other than with BA.2. Based on the analysis of the aforementioned mutations, we can observe that there are few differences between BA.5 and BQ.1. SpikeBA.5 and SpikeBQ.1 have the L452R substitution and deletions of the residues H69 and V70. Comparing the sequences of SpikeBQ.1 and SpikeBA.5, we observed few differences such as a deletion at Y145 in BQ.1 and the substitutions K444T and N460K. These three mutations are only found in SpikeBQ.1. Nevertheless, these two mutations resulted in the dominance of BQ.1 over BA.5 [2].

In order to investigate how mutations are related to structural changes of different Spike of VOCs, we compared the backbone root-mean-square deviation (RMSD) of the protomers of SpikeB.1.617.2 (from Q14 to S1146), SpikeBA.2 (from Q14 to S1144), SpikeBA.5 (from Q14 to S1142) and SpikeBQ.1 (from Q14 to S1141) (Figure S4). The average RMSD values of SpikeB.1.617.2, SpikeBA.2, SpikeBA.5 and SpikeBQ.1 converged with ~ 5.0, ~ 5.1, ~ 5.0 and ~ 4.9 Å, respectively. Our previous study resulted in RMSD values of ~ 5 Å for SpikeBA.1 and ~ 4 Å for SpikeWT [21]. Therefore, SpikeB.1.617.2, SpikeBA.2, SpikeBA.5, and SpikeBQ.1 visit a greater number of different structural conformations, similar to SpikeBA.1, when compared with the wild type.

Next, we calculated the backbone root-mean-square fluctuation (RMSF) of the amino acid residues of SpikeVOCs (Fig. 1). We attributed a RMSF cut-off value of 3 Å and classified the residues with high conformational flexibility as those with an average RMSF greater than 3 Å. To facilitate the location of amino acid residues in the Spike (Fig. 1a–c), we have mapped the residue numbering based on the domain architecture of the SpikeWT protein, as shown in Fig. 1b. SpikeB.1.617.2 mutations are located in NTD (T20R, G75R, D80Y, T95I, G142D), RBD (L452R, T478K), SD1 (T547K, D614G, P681R), SD2 (D614G, H655Y), on HR1 (D950N) and HR2 (V1264L) (Figure S1 and S2). These mutations influence the structural flexibility of SpikeB.1.617.2, as shown by the high RMSF value in NTD (14–23, 71–74, 113, 145–155, 164–165, 167, 212–218, 248–257), in RBD (474–488), in SD2 subdomain (622–626), in proximity of the S1/S2 cleavage site (679–690), in proximity of the S2’ cleavage site (810–812), between the FP and HR1 domains (836–846), in HR1 domain (941–942), and in HR2 domain (1140–1146) (Fig. 1d). The residues with high RMSF values in Omicron sub-variants are: (1) BA.2, in NTD (14–20, 66–76, 142–150, 161–163, 210, 243–255), in RBD (468–486), in SD2 subdomain (620–623), in proximity of the S1/S2 site (675–685), in proximity of the S2’ site (806–810), between the FP and HR1 domains (842), in HR1 domain (938, 940–941), and in HR2 domain (1138–1144) (Fig. 1e); (2) BA.5, in NTD (14–21, 67–73, 139–146, 152, 175–180, 207–210, 240–253), in RBD (469–482), in proximity of the S1/S2 site (673–683), in proximity of the S2’ site (803–808), in FP (827), between FP and HR1 domains (830–841), and in HR2 domain (1134–1142) (Fig. 1f); (3) BQ.1, in NTD (14–21, 68–71, 140–145, 177, 206–208, 241–251), in RBD (470–479), in SD2 subdomain (616–620), in proximity of the S1/S2 site (673–680), in proximity of the S2’ site (803–807), and in HR2 domain (1135–1141) (Fig. 1g).

Structural flexibility of the trimeric SpikeB.1.617.2 (Delta), SpikeBA.2, SpikeBA.5, and SpikeBQ.1. a S1 and S2 subunits of the Spike trimer (PDB ID 7W94 [69]). b SpikeWT amino acid sequence showing the S1 and S2 subunits and the cleavage sites S1/S2 and S2’ [1]. c RBD and NTD of the Spike. Backbone root-mean-square fluctuation (RMSF) calculated from the molecular dynamics (MD) simulations of the d SpikeB.1.617.2, e SpikeBA.2, f SpikeBA.5, and g SpikeBQ.1

In our previous study, SpikeBA.1 showed high RMSF values for residues in NTD (14–21, 70–75, 142–149, 175–179, 208–215, 243–255), in RBD (442–444 and 468–485), in SD2 subunit (599–603, 635 and 654), in proximity of the S1/S2 site (672–690), in FP (828–838), in HR1 domain (938–944), and in HR2 domain (1139–1144) [21, 70]. When compared with SpikeBA.1 [21], our results for SpikeBA.2, SpikeBA.5, and SpikeBQ.1 show that these proteins are less flexible in the S2 subunit, but more flexible in the S1 subunit (mainly in the NTD and RBD). In our previous study, we defined two regions of the Spike RBD interaction interface with human receptor Angiotensin-converting enzyme 2 (hACE2), E1 (residues K417, L455, F456, and 470–490) and E2 (residues 446–453 and 493–505), which have different physicochemical features that impact the interaction with hACE2 [11]. E1 makes mostly hydrophobic interactions, while E2 interacts mostly with hydrophilic interactions with hACE2 [11]. The Spikes of BA.2, BA.5 and BQ.1 share the same mutations in NTD (T20I, del25/27, A28S, G339D), in RBD (S371F, S373P, S375F, T376A, D405N, R408S, K417N, S477N, T478K, E484A, Q498R, N501Y, Y505H), and SD1 (D614G, H655Y, N679K, P681H, N764K, D796Y). It is known that substitutions in E1 region (K417N, S477N, T478K, E484K) and E2 region (N501Y) contribute to increase the binding affinity between RBD and hACE2 [1, 9, 31]. All results are in total agreement with several previous studies, which show high RMSF values in the residues 470–480 (in E1 region), interacting or not with hACE2 [9, 11, 21, 70]. Several studies have shown that there are high RMSF values in the E1 region of the RBD of the WT, Alpha, Beta, Gamma, Delta and Omicron variants between residues 460 and 490 [9, 11, 71,72,73]. We compare our results with the molecular dynamics study of Zhou et al. and we noted that the RMSF patterns of the WT glycosylated Spike protein and the Beta, Delta and Omicron variants [74, 75] were similar to those obtained in our study. Interestingly, even when interacting with hACE2, the glycosylated Spike protein showed the same pattern of RMSF values also when not interacting with this enzyme [70, 74,75,76,77]. Molecular dynamics studies of the Spike protein from other VOCs also showed the similar RMSF pattern, indicating that the structural data are, on average, convergent to the same result, regardless of whether or not it is glycosylated or even interacting or not with hACE2 [70, 74,75,76,77]. These results are also in line with our previous studies that show the similar pattern of RMSF values for subvariant BA.1 [21].

Overall, we did not observe significant RMSF differences in the S1 subunit of the SpikeBA.1 [21], SpikeB.1.617.2, SpikeBA.2, SpikeBA.5, and SpikeBQ.1. However, we noted a significant decrease of the RMSF values in the proximity of the S1/S2 and S2’ cleavage sites of the SpikeBA.5 and SpikeBQ.1 when compared with other Spike VOCs studied in this paper. These results suggest that, for both BA.5 and BQ.1, after hACE2 interaction and S2’ cleavage, the S2 subunit of Spike tends to decrease its structural changes to expose the fusion loop, when compared with previous variants, decreasing the fusion of virus envelope and host cell membrane. Indeed, the Omicron variant has lower fusogenicity than WT and the other previous variants, indicating that the Spike protein may be losing its ability to make important structural changes in the membrane fusion process [22, 78, 79]. Park and collaborators investigated the reverse mutation in specific positions of the BA.1 variant, in order to evaluate how these mutations affect the fusogenicity [79]. In this regard, eleven reverse mutations close to the S1/S2, S2’ and HR1 region cleavage sites were analyzed [79]. Reverse mutations near the S2’ site and in the HR1 region had a stronger effect on fusogenicity, where BA.2, BA.2.12.1 and BA.5 exhibited comparable levels of fusogenicity [79]. Overall, these results combined with our data show how mutations are significantly affecting the Spike structure and, consequently, the fusogenicity of Omicron sub-variants.

In order to identify key movements within the Spike trimers, a principal component analysis (PCA) was conducted by diagonalizing the covariance matrix of SpikeB.1.617.2, SpikeBA.2, SpikeBA.5 and SpikeBQ.1. The eigenvalues obtained from PCA allowed to calculate Schlitter entropy (ΔSSchlitter) of the SpikeB.1.617.2, SpikeBA.2, SpikeBA.5 and SpikeBQ.1 conformers, resulting in values of 74.8, 74.6, 75.3, and 74.3 kJ.mol−1.K−1, respectively. The previously reported entropy of trimeric SpikeBA.1 was 75.1 kJ.mol−1.K−1 [21]. Thus, the entropy values of the analyzed variants were similar but significantly higher than that of SpikeWT from a prior study (72.5 kJ.mol−1.K−1) [21]. This demonstrates that the Spike proteins of VOCs, in general, exhibit significant structural changes, which are quantified by ΔSSchlitter. Next, utilizing the first eigenvector (PC1, first principal component), it was observed that the major conformational changes in SpikeB.1.617.2, SpikeBA.2, SpikeBA.5 and SpikeBQ.1 primarily occurred in the NTD and RBD regions (Figs. 2 and 3). Cryo-EM studies have revealed that Omicron sub-variants exhibit reduced affinity for antibody interactions due to the mutations [6, 32,33,34, 38, 80]. Antibodies targeting the N-terminal domain (NTD) and receptor-binding domain (RBD) display limited binding to these sub-variants, particularly those specifically designed to recognize the RBD [5, 6, 21, 32, 81,82,83].

MD trajectories projected in eigenvector 1. The MD trajectory projected in eigenvector 1 (also called principal component 1, PC1) represents the main movements filtered from MD trajectory. The main movements were detected in the NTD and RBD. The panels represent structure superpositions of 100 snapshots projected in the eigenvector 1 (PC1) for a SpikeB.1.617.2 and b SpikeBA.2. The backbone conformational changes were quantified by calculating entropy using Schlitter method (ΔSSchlitter) at absolute temperature of 310 K. SpikeB.1.617.2, ΔSSchlitter = 74.8 kJ.mol−1.K−1; SpikeBA.2 = 74.6 kJ.mol−1.K−1. The small rectangle is located at the S1 subunit while the big rectangle is located at the S2 subunit

MD trajectories projected in eigenvector 1. The MD trajectory projected in eigenvector 1 (also called principal component 1, PC1) represents the main movements filtered from MD trajectory. The main movements were detected in the NTD and RBD. The panels represent structure superpositions of 100 snapshots projected in the eigenvector 1 (PC1) for a SpikeBA.5 and b SpikeBQ.1. The backbone conformational changes were quantified by calculating entropy using Schlitter method (ΔSSchlitter) at absolute temperature of 310 K. SpikeBA.5, ΔSSchlitter = 75.3 kJ.mol−1.K−1; SpikeBQ.1 = 74.3 kJ.mol−1.K−1. The small rectangle is located at the S1 subunit while the big rectangle is located at the S2 subunit

In general, the significant conformational changes in the RBD and NTD result in substantial decreases in antibody neutralization. One contributing factor is the reconfiguration of the RBD loop (residues 371–376), which disrupts the binding of neutralizing antibodies, ultimately hindering antibody neutralization [6]. The wide range of conformations in these regions is crucial as it may diminish recognition by neutralizing antibodies against SARS-CoV-2 [84]. Notably, high fluctuations in nearby loops 144–155 and adjacent loops 246–260, observed in different VOCs, contribute to decreased neutralizing antibody activity [85]. Structural alterations in the NTD and RBD significantly decrease binding affinity to neutralizing antibodies, facilitating immune system evasion. SARS-CoV-2 VOCs can evade neutralizing antibodies from convalescent and/or immunized individuals [1, 81, 82, 86]. Additionally, deletions 242–244 in Spike induce structural changes in the NTD, further reducing the efficiency of neutralizing antibodies [85]. Importantly, NTD and RBD contain crucial epitopes for neutralizing antibody recognition, and mutations in these domains lead to reduced recognition by vaccine-induced neutralizing antibodies [1, 5, 6, 40, 41, 87]. Numerous studies have demonstrated significant decreases in neutralization by antibodies against the Delta variant and Omicron sub-variants [1, 5, 6, 40, 41, 87]. Intriguingly, previous research has shown that reinfection risk with Omicron sub-variants is six times higher than other SARS-CoV-2 VOCs [88]. The mechanism of action of antibodies that recognize RBD is based on blocking the Spike/hACE2 interaction, preventing viral infection [89, 90]. However, although some antibodies interact with RBD, they do not prevent interaction with hACE2 but destabilize the homotrimer [91, 92]. Considering that the first principal component reveals the main movements during molecular dynamics simulation and that such movements are related to the structural changes of NTD and RBD, our results suggest a potential association between the loss of sensitivity to neutralizing antibodies and these structural changes. Our hypothesis is that the structural changes in these domains, which are the primary targets of neutralizing antibodies, could diminish the efficacy of vaccines against SARS-CoV-2. Therefore, the PCA results suggest that the pronounced conformational changes in the RBD and NTD of these variants are associated with diminished recognition by neutralizing antibodies.