Interaction Between Heparan Sulfate Oligosaccharide and the Receptor Binding Domain of Wuhan and Omicron Variants of the SARS-CoV-2 Spike Protein

1. Introduction

Heparan sulphate (HS) and its structurally similar heparin are linear, highly sulphated polysaccharides, belonging to the glycosaminoglycan (GAG) family. their repeating unit consists of a disaccharide composed of uronic acid (β-D-glucuronic or α-L-iduronic) 1→4 linked to α-D-glucosamine (GlcN). This backbone undergoes various modifications [1], such as partial N-deacetylation/N-sulfation, O-sulfation, and epimerization, resulting in a heterogeneous population of polysaccharide linear chains with variable length, composition, and degree of sulfation. Precisely, the composition of HS and heparin can vary in terms of N-acetylation/N-sulfation, levels of 6-O- and more rarely 3-O-sulfation of glucosamines and 2-O-sulfation of uronic acid residues [2,3,4]. Heparin, differently from HS, is characterized by shorter chains, with higher sulfation degree and higher percentage of iduronic than glucuronic acid.

In the extracellular matrix (ECM), HS is organised in proteoglycans (HSPGs) or glycoproteins, which are ubiquitously expressed on the cell surface and within the ECM [2]. The polyanionic nature of HS, resulting from its high density of negative charges due to multiple sulfate and carboxylate groups, enables it to interact with a broad spectrum of proteins, including cytokines, chemokines [5], grow factors [6], microbial adhesion factors [7], and viral and/or bacterial recognition proteins [8], that are characterized by clusters of solvent exposed positively charged Lys and/or Arg side chains.

Recent attention focuses on the role of HS in viral pathogenesis, particularly in the context of emerging respiratory viruses. Several pathogens, including herpes simplex virus [9], dengue [10], and respiratory syncytial virus [11,12,13], exploit HS to initiate the contact with host cells. Similarly, several evidence suggest that SARS-CoV-2 coronavirus utilises HS as an auxiliary attachment factor [14]. This virus, characterised by its single-stranded RNA genome and a lipid capsid studded with spike (S) glycoproteins, employs its receptor binding domain within the S1 subunit (S1-RBD) to engage angiotensin-converting enzyme 2 (ACE2), which is located on the surface of the host cells. While ACE2 serves as the primary entry receptor, structural and biochemical studies have indicated that initial interactions between the spike protein and cellular HS facilitates the viral engagement with ACE2 [15,16], through the formation of a ternary complex S1-RBD – HS – ACE2 [14,17].

The emergence of SARS-CoV-2 lineages with enhanced transmissibility and altered receptor affinities requires understanding how variants of SARS-CoV-2 spike protein interact with the HS. The Omicron variant harbours multiple mutations in the S1-RBD, most of which are proximal to the ACE2 attachment site, which may modulate its interaction with both ACE2 and eventually with its co-receptors like HS [18]. Significant differences on both kinetic and thermodynamic aspects of the interaction between heparin and wild-type (WT), Delta and Omicron S1-RBD, were observed [19]. However, the structural and functional implications of these changes remain poorly defined.

A similar scenario is also observed concerning glycosylation, which is a crucial post-translational modification that shapes protein structure and function [20]. The S protein of SARS-CoV-2 is heavily glycosylated, forming a glycan shield that masks over 60% of its surface, helping the virus evade immune detection [21]. It has been shown that, in the S1-RBD region, the N331 and N343 residues are two highly conserved glycosylation sites across all the variants of SARS-CoV-2 RBDs [22]. While the N331 is located in the highly flexible region connecting the RBD to the NTD, the N343 glycan extends across the RBD surface, bridging the two helical regions flanking the central beta-sheet core and it is figured out to interfere with a significant portion of the domain that in S1-RBD binds HS [21,23].

In this study, the molecular interaction between HS hexasaccharide (hexa , Figure 1) and the S1-RBD of both the WT and Omicron SARS-CoV-2 variants was investigated. Parallelly, the interference between the glycosyl moiety at position N343 and the HS binding site of S1-RBDs is characterized. Multiple independent molecular dynamics (MD) simulations were applied to predict the glycan conformation in unbound and bound state, as well as the geometry of the complexes between hexa and Omicron and WT S1-RBD, that have been generated previously by docking [17]. Proton saturation transfer difference (¹H-STD) NMR spectroscopy allows the experimental validation of Omi-S1-RBD-hexa complexes and the corresponding proton epitope binding. The initial ¹H-STD build-up slope (STD₀) was determined and compared with the corresponding simulated values, using the RedMat application [24]. In summary, this study characterizes the structural and dynamic aspects of the interaction between HS oligosaccharide and Omicron and the WT S1-RBD, as well as the potential role of the N343 glycosyl moiety in perturbing this contact. These findings enhance our understanding of the molecular mechanism by which the SARS-CoV-2 virus take advantage of the cellular HS to approaches and begin the infection of the cell host, and is propaedeutic for the design of new antiviral strategies against the fast-evolving viruses such as SARS-CoV-2 coronavirus.

2. Materials and Methods

2.1. Molecular Dynamics Simulations

2.1.1. System Preparation and MD Simulation Protocol

The initial geometries that were submitted to MD simulations were selected from docking results reported in our previous work [25], where the RBD sequence employed in the study was 333-526, therefore the glycosylation at the N331 site is not considered.

For the Wuhan S1-RDB system (referred to as WT-S1-RBD-hexa), we started from the single complex structure that was obtained previously [25] , while for the Omicron S1-RBD system, the two top-ranked, that represent the two principal binding modes, were selected. (referred to as Omi-S1-RBD-hexa-I and Omi-S1-RBD-hexa-II, respectively) (Figure S1) The Amber18 [26] software was used for the systems preparation, MD simulations and analysis. We applied the Amber (ff14SB) [27] and the GLYCAM-06j-1 [28] force fields for the protein and ligand, respectively. The glycosylated moiety at N343 was also described with GLYCAM-06j-1. The hexa oligosaccharide was described in its fully deprotonated form, since at pH=7.0, all the ionizable group are dissociated.

Each system was solvated with a TIP3P [29] cubic water box edge of 15 Å. Na⁺ ions were added to neutralise the system and then Na⁺Cl^- ions were added to mimic the ionic strength of a physiological environment ([NaCl] = 0.15 M). The cut-off for no-bond interaction was set to 8 Å; the Particle Mesh Ewald method was used to treat long-range electrostatic interactions. Each system was minimised using three consecutive minimisation steps, each consisting of 20 cycles of the steepest descent algorithm followed by 880 cycles of conjugate gradient minimisation. Harmonic restraints (k = 50 and k = 10 kcal/mol·Å²) were applied to solute atoms relative to the input geometry in these first two steps, while the final step was performed without restraints, allowing full minimisation.

The systems were then heated to 300 K in 100 ps in NVT conditions (Langevin thermostat [30] with collision frequency $\gamma =1p{s}^{-1}$) with a time step of 0.5 fs. The hydrogen atoms were fixed using the SHAKE [31] algorithm and a positional restraint of 10 kcal/mol Ų were set for the solute atoms. To equilibrate the system density, we first performed an equilibration step in NVT conditions (100 ps, dt = 0.5 fs) with restrained solute atoms (5 kcal/molŲ) followed by an equilibration run in NPT conditions (Berendsen Barostat [32], p = 1 atm, dt = 0.5 fs and pressure relaxation time ${\tau }_P$ = 2ps) and a final NVT step with no restraints (100 ps, dt = 0.5 fs). After equilibration at T = 300 K, MD simulations were carried out in triplicate. The production runs (200 ns, dt = 2 fs) were performed in NVT conditions, using the Langevin thermostat [30] to maintain a stable temperature of 300 K. All the replicas had the same coordinates, but different velocities, randomly selected from a Maxwell distribution at 300 K. Frame coordinates were saved every 20 ps for a total of 10.000 structures for each run. For the analysis, each of the three independent replicas was merged into one single meta-trajectory, including 60.000 sampled geometries.

The oligosaccharide hexa in the unbound state was characterised by three independent MD simulations, following the same protocol conditions previously described for the hexa in bound state with the S1-RBD. For one of the replicas, the initial conformation of hexa was characterised by residues: IdoA2S, GlcA and GlcNS6S in chair ¹C₄, ⁴C₁, and ⁴C₁, in that order. The initial glycosidic conformation was set as follows: ${\varphi }_1/{\psi }_1$ = 41°/14°, ${\varphi }_2/{\psi }_2$ =-40°/-30°, ${\varphi }_3/{\psi }_3$ = 60°/30°, ${\varphi }_4/{\psi }_4$ = -39°/-33°, ${\varphi }_5/{\psi }_5$ = 41°/14°, in order from non-reducing to reducing end as reported in our previous work [17]. The glycosidic dihedral angles are defined as ${\mathsf{\varphi }}_{\mathrm{i}}={\mathrm{H}1}_{\mathrm{n}}\to {\mathrm{C}1}_{\mathrm{n}}\to {\mathrm{O}4}_{\mathrm{n}}\to {\mathrm{C}4}_{\mathrm{n}+1}$ and ${\mathsf{\psi }}_{\mathrm{i}}={\mathrm{C}1}_{\mathrm{n}}\to {\mathrm{O}4}_{\mathrm{n}}\to {\mathrm{C}4}_{\mathrm{n}+1}\to {\mathrm{H}4}_{\mathrm{n}+1}$. The hexa initial geometries for the second and third MD simulation replicas were obtained by adding or subtracting 10°, respectively, to the initial ${\mathsf{\varphi }}_{\mathrm{i}}$ and ${\mathsf{\psi }}_{\mathrm{i}}$ values defined previously. The MD simulation trajectories of these three independent replicas were merged into one meta-trajectory that includes 60.000 sampled geometries of hexa for the analysis. Each meta-trajectory (hexa in unbound and bound states with RBD) was analysed using CPPTRAJ [33].

2.1.2. Filtering and Analysis of the MD Trajectories: Bound State Condition of S1-RBD and Hexa

We filtered the meta trajectories of each hexa-S1-RBD complex to select the frames in which the oligosaccharide hexa is in bound state with the corresponding S1-RBD. Precisely, the hexa is in the bound state with S1-RBD if and only if the centre of mass (CoM) of the two central monosaccharides (C) and (D) is within 8 Å from the CoM of the proximal amino acid. (MDAnalysis 8 Library tool v 2.8.0 Python v 3.10.16). The percentage of frames of each meta-trajectory, that agree with this bound state condition, was 73% for Omi-S1-RBD-hexa-I, 53% for Omi-S1-RBD-hexa-II and 39% for WT-S1-RBD-hexa. The MD simulation filtered trajectories were analysed using CPPTRAJ [33].

2.1.3. Model Validation by Reduced Matrix (RedMat)

The application RedMat compares experimental STD₀ measured upon interaction between hexa and Omicron S1-RBD with the corresponding theoretical values, that were simulated using selected geometries of the complexes between Omicron S1-RBD and hexa [24]. The following parameters were set for the simulation of STD₀: the spectrometer frequency at 600 MHz, the correlation time (${\tau }_c)$ of the complex at 25 ns, and a cutoff distance of 15 Å. The analysis was performed on the meta-trajectory that was previously filtered using the bound state condition (see paragraph 2.1.2). To select representative ligand binding modes, we carried out a ligand cluster analysis on those frames of the complex Omi-S1-RBD-hexa-I and Omi-S1-RBD-hexa-II, that exhibit an R-NOE factor below 0.3, after fitting on protein residues of the frames.

2.1.4. Cluster Analysis

The cluster analysis was applied to the meta-trajectories upon filtering with the bound state condition and R-NOE smaller than 0.3 in the complexes Omi-S1-RBD-hexa-I, and Omi-S1-RBD-hexa-II. Differently, the cluster analysis was applied on the meta-trajectory upon filtering with only the bound state condition in complex WT-S1-RBD-hexa, in fact, no STD₀ were available in this last case.

The cluster analysis was performed using CPPTRAJ [33] considering all the carbon, oxygen, sulfur, and nitrogen atoms for the hexa ligand. Clusters are generated using the hierarchical agglomerative approach, setting a cut-off distance of 10 Å between the clusters. The distance is calculated as the average distance between members of two clusters. A cluster is considered representative if the population percentage is not less than 10%.

2.2. NMR Spectroscopy: Sample Preparation and Experiment

The SARS-CoV-2 spike protein S1-RBD of the Omicron variant, expressed in HEK293 cells, was purchased from GenScript (Netherlands). The protein solution (1.12 mg/mL in phosphate buffer pH 7.4) was exchanged with 20 mM HEPES-d18 pH 7.2 with 200 mM NaCl (D2O) using VWR^® centrifugal filters (10 kDa membrane, 0.5 mL).

The NMR experiments were performed with a Bruker Avance III spectrometer operating at 600.13 MHz and equipped with a high-sensitivity 5 mm TCI cryoprobe. The ligand was dissolved in 0.2 mL of the purified protein, and the final solution was transferred into a 3 mm NMR tube. The concentrations of ligand and protein were 5.4 μM and 540 μM, respectively (protein/ligand ratio = 1:100). The spectra were acquired at 293 K. Based on our previous work [17], the number of scans was 2048, the dummy scans were 4, and the relaxation delay was set to 6 s. A spin-lock pulse of 10 ms was employed to remove the protein broad signal. The selected frequencies for acquiring the on- and off-resonance spectra were 580 and 20000 Hz, respectively. Five different experiments were performed, varying the saturation time (0.5, 1, 2, 3, 5 s). The ¹H-STD NMR spectra were obtained by phase cycling subtraction of the on-resonance and off-resonance data acquired in interleaved mode. The STD intensities at the initial slope were calculated for the anomeric proton signals, since they are sharp and well-resolved. A mono-exponential equation, as proposed by Mayer and James [34] (Equation (1)), was used to fit the experimental build-up curves.

$\mathrm{S}\mathrm{T}\mathrm{D}\left({\mathrm{t}}_{\mathrm{s}\mathrm{a}\mathrm{t}}\right)=\mathrm{S}\mathrm{T}{\mathrm{D}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\left(1-{\mathrm{e}}^{\left(-k_{sat}\bullet {\mathrm{t}}_{\mathrm{s}\mathrm{a}\mathrm{t}}\right)}\right)$ (Eq. 1)

The resulting values were expressed as relative STD₀ percentages by normalising the intensities against the most intense signal.

3. Results and Discussion

The hexasaccharide (hexa) shown in Figure 1 was used as a molecular probe to characterise the interaction between HS and the tip of spike protein S1-RBD of the Wild type (or Wuhan, (PDB ID: 6M0J [35]) and the Omicron (DB ID: 7WBP [36]) variant of the SARS-CoV-2 (Figure S1).

In our previous work [25], we observed that the S1-RBD of WT, Delta, and Omicron variants are characterized by different affinity for immobilized heparin. The hypothesis raised is that mutations in the virus S1-RBD can modulate the association and dissociation process with HS, thereby affecting the interaction with the ACE2 receptor. Omicron S1-RBD has a stronger affinity for HS compared to the Wuhan strain, as predicted by other groups. [19] To characterize the interaction between hexa and Omicron S1-RBD, we applied a computational workflow combining MD simulations with ¹H-STD NMR interaction experiments. Analogously, the interaction between hexa and WT S1-RBD, was analyzed using MD simulation, since no experimental STD₀ were available for that system. Starting from the relevant binding modes generated previously by docking [25], we carried out multiple independent molecular dynamics (MD) simulations. Additionally, a full conformational analysis of hexa in unbound and in bound states with S1-RBD belonging to both WT and Omicron variants were done, to detect potential conformational changes induced upon interaction with these S1-RBDs variants.

3.11. H-STD NMR Interaction Experiment Between Hexa and Omi S1-RBD

¹H-STD NMR experiments were conducted to identify the protons and residues of hexa that interact with Omicron S1-RBD, therefore assessing their proximity to its surface. STD₀ percentages (STD₀ %), for non-overlapping ligand protons (i.e., anomeric protons) were calculated as described in the Materials and Methods section and subsequently mapped onto the ligand structure to depict the interacting epitope map (Figure 2 and Figure S2). Comparison of the STD NMR spectra with the corresponding ¹H spectra revealed that all anomeric protons of hexa exhibit significant STD effects, indicating that each sugar unit is in contact with the surface of S1-RBD (Figure 2). The highest STD₀ % effect was observed for H6C (100%), followed by H1D (96%), while H1F, H1E, H1C, H1B and H1A displayed lower STD enhancements, ranging from 10% to 36% (Table S1). H6E also shows an appreciable contribution to the interaction (Figure 2 and Figure S2b). However, the corresponding STD₀ was not determined due to spectral overlap (Figure S2b).

Notably, the H6C and H6E present approximately opposite orientation in the 3D structure of hexa (Figure S3), therefore, the detection of these two signals at comparable intensities (Figure 2 and S2b) indicates presumably that hexa interacts with the site I of S1-RBD (Figure S4) with opposite orientation of the glycan relative to the S1-RBD, corresponding to two different epitope binding (or binding modes), in which H6C or alternatively H6E faces to the surface of S1-RBD [17].

3.2. Conformation of Hexa in Bound State with Omicron and WT S1-RBD

The binding epitope of hexa with Omicron S1-RBD is characterized by strong ¹H-STD signals originating from the units C and D of the oligosaccharide in bound state with Omi-S1-RBD (Figure 2 and Table S1). This supports the idea that GlcNS6S(C) and GlcA(D) are, on average, closer to the protein surface in comparison to the terminal residues. Hexa, when bound to the Omi-S1-RBD or WT-S1-RBD almost preserves the conformational properties that were observed in its unbound state. Specifically, residues: A, C, E and D adopt chair ⁴C₁ conformation, while IdoA2S (residues B and F) assume the chair ¹C₄, even if a minor population of the skew-boat ²S₀ is present (Figure S5 and Table S2). The inter-glycosidic bond geometry of hexa remains essentially unaltered upon interaction with the Omicron and Wild-type S1-RBDs, preserving their unbound state conformation; this can be seen comparing the most populated ${\varphi }_i/{\psi }_i$ states in Figure 3a, with Figure 3b,c,d. (see also, Figure S6, and Table S3).

of hexa in the unbound state (a), RBD-Omi-hexa-I (b), RBD-Omi-hexa-II (c), WT-RBD-hexa (d). The colour gradient (yellow to blue) is proportional to the ${\varphi }_i/{\psi }_i$sampled state, identifying (qualitatively) the preferred conformations. The most populated ${\varphi }_i/{\psi }_i$ state is represented as a white star.

This result is consistent with the structural and conformational properties of hexa, in which the strong electrostatic repulsion between sulfated and carboxylate groups, along with the corresponding steric hindrance, are dominant and enforce the molecule to preserves its linear and stiff shape, that are typical of its unbound state, despite the interaction with different S1-RBDs. As a result, the inter-glycosidic bond geometries remain largely unaltered upon interaction with the surface of S1-RBD, figuring out a weak and low specificity of interaction between this oligosaccharide and these S1-RBDs. Interestingly, the structural rigidity of hexa preserves the position of its charged groups (sulphates and carboxylates) that remain oriented in approximately opposite direction. In summary the conformation of hexa defines and preserves two wider faces, both potentially available to contact the positively charged patches that are solvent exposed in the site I of S1-RBDs (Figure S3).

3.3. Characterisation of the Interaction Between Hexa and S1-RBD of Omicron and WT Variants

3.3.1. Analysis of the MD Simulation Meta-Trajectories

Starting from the most representative docking poses obtained in our previous works [25] (Figure S1), we run three independent MD simulations for each of the modelled complex, labelled as Omi-S1-RBD-hexa-I, Omi-S1-RBD-hexa-II, and WT-S1-RBD-hexa. In all these systems, hexa remains preferentially bound to the shallow channel defined by residues R346, N354, R355, K356, R357 and R466, which were previously identified as site I [17,37] (Figure S4). Interestingly, the docking does not give solutions in which hexa binds Omi-S1-RBD or WT-S1-RBD trough site II (K424, R454, R457, K458, K462, R466) and III (R403, R408, K417, and K444. According to docking results, the Omicron S1-RBD complexes show hexa bound within site I in two reverse orientations, allowing the two opposite surface of hexa (Figure S3) to interact with the surface of S1-RBD. These two alternative modes of interaction are identified by the GlcNS6S(A) oriented toward the ACE2 binding site and proximal to R346 (Figure S1a). Differently, when the opposite face of hexa approach the site I of S1-RBD, the residue A is in contact with R355 and R466 (Figure S1b). In the WT-S1-RBD-hexa, the oligosaccharide is docked with only one orientation, in which residue A is proximal to R346 and K356, as shown in Omi-RBD-hexa-II, with (Figure S1c).

Analyzing the MD simulation meta-trajectories of the hexa in bound state with Omicron and WT S1-RBD, the ligand exhibits a significant mobility along the channel corresponding to site I (Figure S7a-c) as detected by the RMSD, while the secondary structure of protein S1-RBD remains stable in the spanned simulation time (Figure S7d-f and Figure S8). In fact, in all simulated hexa S1-RBD complexes, the RMSD (calculated on heavy atoms of hexa relative to the input structure) shows a significant change, ranging between 10 to 30 Å (Figure S7a-c). Since the ligand conformation remains stable throughout the simulations (Figure 3), this variation in RMSD values (Figure S7a-c) reflects only on the rotational and translational movements that hexa experiences when is in contact with the concave and shallow surface of site I of both Omicron and WT S1-RBD (Figure S9). Therefore, despite the multiple electrostatic interactions between the positively charged amino acids and the negative charges of oligosaccharide, the ligand retains enough kinetic energy to explore a limited region surrounding site I on the surface of S1-RBD.

To identify the amino acids that are in persistent contact with the ligand, we applied a 6 Å cut-off filter between the hexa heavy atoms and key residues of S1-RBD. The frequency of each contact was then calculated as a percentage on the whole meta-trajectory and reported in Figure 4 as contact heatmap.

The color gradient ranges from lighter shades (orange to white), indicating less frequent contacts, to darker shades (dark red or black), representing more persistent (stronger) interactions, according to the color scale. If an atom of the oligosaccharide forms more than one contact with the same amino acid, only the contact with the highest percentage was reported. For each complex, the blue stars indicate a contact that was already present in the initial pose. Furthermore, these interactions have been classified as ionic or polar according to atoms or groups that are involved. Contacts between positively and negatively charged groups are classified as electrostatic interaction. Contact involving polar atoms with an average distance $\le$ 4.5 Å are predicted to be polar interactions and/or potential hydrogen bonds (Table S4).

The Omi-S1-RBD-hexa-I system, shows alteration of its initial binding epitope, indicating that the bound glycan explores the surroundings of site I. The amino acids mostly involved in this contact are R357 (67%), followed by Y396 (43%), K356 (32%), R466 (27%), N394 (25%), and R346 (20%). R357 makes the most frequent contacts with the hexa and compared to the initial state, forms additional electrostatic interactions with glucosamine C, glucuronic acid D, and iduronic acid F, with frequencies ranging from 66% to 33% (Table S4). K356 forms new electrostatic interactions with the central part of the oligosaccharide (residue C). Two residues, S359 and N360, that are proximal to site I (Figure S4) and known to be involved in GAG-SARS-CoV-2 recognition [38], engage the central C and D residues of the glycan. Furthermore, a new polar contact involving N394 and Y396 and glycan residues E and F, respectively was established, while the interaction between E and R466 was lost.

In the Omi-S1-RBD-hexa-II system, the contacts that were established upon docking are preserved, despite that, new interactions were also gained (Figure 4b). R346 persistently engages the central sugars C and D, showing contact percentages of 51% and 67%, respectively (Table S4c). These contacts are further stabilized by less populated polar interactions T345-GlcNS,6S(C) (49%) and GlcNS,6S(C)-A344 (38%). Furthermore, residues K356 and R357, which are part of the core of site I, frequently interact with IdoA2S (B) (40%) and GlcNS6S(A) (27%). Interestingly, the heatmaps in Figure 4a,b predict a different binding epitope for the hexa in bound state with Omicron S1-RBD, through its two opposite surfaces (Figure S3), suggesting these two surfaces presents a comparable capacity to interact with the tested S1-RBDs.

In the WT-RBD system, the hexa moves across the central area of site I, displaying greater mobility than that observed in Omicron S1-RBD and hexa. This increased mobility is evident in the persistent contacts in Figure 4c, where they differ from the initial binding epitope, indicated by blue stars. R346 preserves the contact with GlcA(D) (85%), IdoA2S(F) (60%), and GlcNS6S(E) (57%). Analogously, K356 and R357 form persistent electrostatic interaction with the central residue GlcNS6S(C), whose populations are 72% and 55%, respectively. Interestingly, K356 shows weaker contacts with terminal residues GlcNS,6S(E) and IdoA2S(B), while R357 shows contacts with IdoA2S(B) and GlcNS6S(A) (Table S4e-f). This reveals how the bound glycan spans the surface of the S1-RBD that surrounds site I in the simulated time window. Furthermore, N354, R355, K356, R357, N360, and R466, exhibit a range of contact frequencies between 57% and 18% with the terminal units IdoA2S(B) and GlcNS6S(A). Interestingly, the contacts between R346 and GlcA(D), or GlcNS6S(C); between K356 and IdoA2S(B), or GlcNS6S(C); between R355 and GlcNS6S(A), as well as the contact between R357 and GlcNS6S(A), were previously found in Omi-S1-RBD-hexa-II.

An important finding that emerges from the contact analysis and the visual inspection of the MD simulation trajectories is the mobility of HS oligosaccharide, here represented as hexa, when in bound state with the Omicron S1-RBD, and particularly with WT S1-RBD. In all these systems, the mobility of hexa promotes a repeated engagement-and-disengagement of ion-pair interactions, between the anionic functional group of the oligosaccharide, and the cationic side chains of Arg, and Lys that belong to site I and its surroundings. Specifically, the sulfo- groups of glucosamines A, C, and E, as well as the carboxyl groups of iduronic acid B and F, establish electrostatic interaction with the same key positive and/or polar side chain of R346, K356, R357 and N394 in Omi-RBD-hexa-I (Figure 4a). This suggests that hexa alternately contacts these residues moving through the shallow cavity of site I in the Omi-RBD-hexa-I system. Analogously, the same glycan that engages the surface of Omicron S1-RBD, but with opposite orientation in comparison to the previous system, allows its sulfo- and carboxyl groups to contact R355 and K356 in the Omi-RBD-hexa-II (Figure 4b), or K356 and R357 in the WT-RBD-hexa complex (Figure 4c). These results indicate a remarkable ligand mobility, with the hexa undergoing to a combination of rotation and translation movements on the positively charged surface, depicting a protein binding epitope that appears more extended than the size of the hexasaccharide probe.

3.4. The Mutation G339D in Omi-S1-RBD Reduce the Shielding Effect That the N343 Glycosyl Moiety Exerts on Site I

To elucidate the role that the N343 glycosyl moiety exerts on the interaction between hexa and S1-RBD, the corresponding oligosaccharide was the built and attached on the Omicron and WT S1-RBDs [25]. This N-glycan corresponds to the biantennary core-fucosylated octasaccharide FA2G2 with the structure [GlcNAc¹ (α1$\to$6) Fuc⁸] (β1$\to$4) GlcNAc² (β1$\to$4) Man³ [(β1$\to$3) Man⁴ (α1$\to$6) GlcNAc⁵] [(β1$\to$6) Man⁶(α1$\to$6) GlcNAc⁷], which has been identified as the most abundant glycan sequence across all SARS-CoV-2 variants [39]. The structure was obtained from the CHARMM-GUI COVID-19 protein archive [40], this corresponds to a full glycosylated all atom model of homotrimeric SARS-CoV-2 spike proteins in open conformation (PDB ID: 6VSB [41]). After selecting the FA2G2 glycan, it was connected to each S1-RBD by linking the reducing end of GlcNAc1 to the ND2 atom of the N343. Interestingly, the MD simulations show a reduced flexibility of the FA2G2 in Omicron S1-RBD in comparison to WT S1-RBD (compare Figure S11a,b with Figure S11c). In WT-S1-RBD, the octasaccharide equally populated two major conformations, defined by the ${\chi }_1$ dihedral angle around the N343 side chain, according the following consecutive atoms$C-C_{\alpha }-C_{\beta }-C_{\gamma }$. The two states allowed for this dihedral are the trans and gauche⁺ corresponding to approximate values: 160° and 70°, respectively (Figure 5a).

Then, FA2G2 has a dual orientation: when ${\chi }_1$ populates the trans state, the glycan occupies a pocket close to the 366-371 helix and the 372–375 loop without making any contacts with the residues of site I, as shows by the glycan with green tubes in Figure 5c, and supported by the heatmap in Figure S12. At values of ${\chi }_1$ around 70°, the glycan moves closer to the site I, making contacts with both hexa and the residues of site I, primarily K356, as shown by the glycan with cyan tubes in figure Figure 5c (see also Figure S12). In this conformation we observed a partial detachment of hexa from site I, thus suggesting a competitive effect for site I between FA2G2 and hexa (see movie S1a). On the other hand, for Omi-S1-RBD-hexa-I and hexa-II systems (Figure 5b,d, and Figure S13, respectively), the glycan sampled essentially the trans conformation ${\chi }_1$ (Figure 5b, movie S1b), and consequently, it remains proximal to 372-375 loop leaving the site I unperturbed, as shown by the green glycan in Figure 5d (see also Figure S12). Also, in a minorly populated state, where ${\chi }_1$ is around 75°, (populated less than 5%), the octasaccharide is solvent-exposed and loses contact with S1-RBD, as shown by the cyan glycan in Figure 5d (see also Figure 13). Notably, the reduced mobility observed for FA2G2 in the Omicron systems can be attributed to the G339D mutation. In the WT-S1-RBS, G339 is small and apolar, allowing free movements of the octasaccharide within 366-371 helix to site I. In contrast, substitution with aspartic acid in the Omicron RBD introduces a larger, negatively charged side chain that impairs movements toward site I. D339 is constantly in contact with the GlcNAc¹, GlcNAc² and Fuc⁸ units of FA2G2, impairing its movements. (Figure S14) Moreover, in the Omicron simulations, interactions between FA2G2 and residues in the helix-loop region 336-375 reduce the flexibility of this region, stabilizing the N-terminal part of the S1-RBD sequence, compared to the wild-type simulations, as can be observed comparing the RMSF of the corresponding residues in Figure 6.

In addition, S371L, S373P and S375F mutations also contribute to this stabilization. In fact, it has been demonstrated [42] that these mutated residues can establish a higher number of hydrogen bonds with helix 366-371, reinforcing the local structure (as confirmed by the decreased RSMF values in Figure 6.

3.5. Selection of the Complexes Omi-S1-RBD-Hexa Using the Experimental STD₀

The RedMat application [24] allows the identification of geometries of the simulated complex that agree with selected experimental STD₀. The application calculates a score (R-NOE factor) for each frame, according to eq. 2. This score is based on the standard deviation between each calculated STD₀ values$\left(ST{D}^{calc,k}\right)$ and corresponding experimental value ($ST{D}^{exp,k})$ considering selected k-signals. An R-NOE smaller than 0.3 normally indicates a good agreement between the experimental and theoretical proton binding epitopes.

$$RNOE=\sqrt{{\displaystyle \frac{\sum {\left(ST{D}^{exp,k}-ST{D}^{calc,k}\right)}^2}{\sum {\left(ST{D}^{exp,k}\right)}^2}}}$$

All the STD-NMR resonances of the central residues of hexa (from sugar E to B, Table S1) were selected for comparison with the corresponding theoretical signals. Terminal sugars A and F were excluded due to their high conformational flexibility (Figure 3). For both Omi-S1-RBD-hexa-I and Omi-S1-RBD-hexa-II systems, the comparison between experimental and theoretical STD binding epitopes shows variable agreement. Only 25% of the RBD-Omi-hexa-I and 18% of RBD-Omi-hexa-II present R-NOE smaller than 0.3 (Figure S15). To identify the dominant binding modes that has been selected by RedMat, a cluster analysis was performed on the MD simulated meta-trajectories of Omi-S1-RBD-hexa-I and Omi-S1-RBD-hexa-I, corresponding to the bound state condition (paragraph 2.1.2) and with R-NOE $\le$0.3 (Table S6 and Figure S16). For each cluster, the hexa geometry that possesses the lowest R-NOE was selected as the representative member of the cluster and reported in Figure 7.

The cluster analysis applied on Omi-S1-RBD-hexa-I system reveals two most populated geometries, with populations of 43% and 39%, respectively (Table S6). These correspond to two distinct binding modes. In the first cluster (43%), hexa is positioned in the central part of site I (Figure 7a and Figure S16a). GlcNS6S(C) forms a hydrogen bond between the backbone of R357 and the C6-OSO₃⁻ group, while the COO⁻ group of GlcA(D) forms a salt bridge with R357. In the second cluster (39%), hexa is in a region proximal to site I, forming interactions with residues S359 and N360 (Figure 7b and Figure S16a). In this case, the C6-OSO₃⁻ group of glucosamine C engages in a hydrogen bond with N334, and the N-sulfo group forms hydrogen bonds with T333 and N360. Additionally, the COO⁻ group of GlcA(D) forms a salt bridge with R357 and a hydrogen bond with S359.

In Omi-RBD-hexa-II, two representative clusters were identified with population of 73% and 10%. The most populated of them is reported in Figure 7c (see also Figure S15b and Table S6). The ligand, which position is consistent with that of the original docking pose, preserves residue IdoA2S(F) oriented toward the ACE region, and in hydrogen bond contact with the backbone of L441 [IdoA2S-3-O-H --- O=C-L441], and ionic interactions with K444. By contrast the N-sulfo group of sugar C and the carboxyl group of sugar D establish ionic interactions with R346, while the C6-OSO₃⁻ group of sugar C forms additional hydrogen bonds with NH of the backbone of T345.

In summary, the interaction between hexa and Omi-S1-RBD presents different binding modes, that agree with the observed STD₀ signals, in which the central sugars C and D are preferentially in contact with S1-RBD. These binding modes are reported in Figure 7a,b,c, and indicate a low level of specificity of the interaction, as was previously figured out in Parafioriti et al. [17]. In all these complexes characterized by R-NOE ≤ 0.3 the main interactions involve glucosamine C and glucuronic acid D, although contributions from glucosamine E and iduronic acid F are also observed. In contrast, units A and B are frequently solvent-exposed (Omi-RBD-hexa-II, Figure 7c), which accounts for their lower STD₀ values (Table S1).

The cluster analysis on the WT-S1-RBD-hexa system, that was applied on the meta-trajectory filtered only with the bound state condition (paragraph 2.1.2), reveals two representative clusters (Figure 8a,b and Table S6c). The most populated of them (48%) shows the ligand bound to site I (Figure 8a). In this binding mode, residues GlcNS6S(A) and IdoA2S(B) interact with R357, while GlcNS6S(C) and GlcA(D) form electrostatic contacts with K356. Additionally, GlcNS6S(E) and IdoA2S(F) establish interactions with R346 and R345. In the less populated cluster (20%), the terminal residues from IdoA2S(F) and GlcNS6S(E) are solvent-exposed, whereas residues from A to D are located in the upper part of site I (Figure 8b), establishing electrostatic interactions with R357 and R466, K356 and R466, and K356 and R346, respectively. Both clusters are consistent with our previous study [17]. Notably, the distance between all anomeric protons and the protein surface is less than 5 Å, indicating that all residues are likely contributing to the binding, and the saturation transfer from the protein to the ligand is figured out to be efficient, supporting (qualitatively) the STD intensities observed in our previous study [17].

4. Conclusions

Heparan sulfate (HS) plays an important role in SARS-CoV-2 infection by promoting the active (open) conformation of the spike protein and acting as a collector of viral particles at the external cell surface. More infectious variants, such as Omicron, appear to have undergone selection for enhanced interaction with both the negatively charged cell-surface HS and the ACE2 receptor by increasing the overall positive charge of the spike protein. However, the mutations observed in the Omicron S1-RBD are primarily localized within the ACE2 binding interface and do not significantly alter the heparan sulfate (HS) binding region.

We characterized the interaction of a hexasaccharide (considered the minimal interacting epitope of heparan sulfate) with the Omicron S1-RBD, using a combination of multiple independent MD simulations and STD-NMR. The oligosaccharide present stiff conformation both unbound and in bound state with Omicron, as well as, in the WT-S1-RBD. The hexa oligosaccharide explores the shallow cavity corresponding to the site I of these S1-RBDs, that is identified by the following adjacent polar or positively charged side chains: R346, N354, R355, K356, R357 and R466. This glycan undergoes partial detachment, re-binding, alternating rotation, and translation movements, probing the site I and the surrounding area without reaching the ACE2 binding region.

To experimentally support the 3D model of such low-affinity glycan-S1-RBD complexes, the RedMat application was used. This computational workflow identified three principal modes of interaction between hexa and the tested Omicron S1-RBD, indicating that the recognition event is characterized by a low level of specificity. As already observed for the Wuhan strain [17], the glycan binds site I using both its major surfaces, particularly involving the central glucosamine C and glucuronic acid D, while the terminal residues A, B, D and F display greater variability and solvent exposure.

In the WT-RBD-hexa complex, clustering of the bound frames revealed interactions extending across the entire ligand, with all sugar units positioned within 5 Å of the protein surface.

Interestingly, the two variants of S1-RBD present distinct conformational behaviour of the biantennary N-glycosyl FA2G2 at N343, that correlates with a different ratio of trans/gauche⁺ states of the dihedral angle ${\chi }_1$; this corresponds to a different capacity of this glycan to shield the site I of S1-RBD, and/or detaching the bound hexa. Precisely, in the wild-type system, the FA2G2 glycan exhibits higher mobility, that corresponds to a comparable population of trans and gauche⁺ states of ${\chi }_1$; in this condition FA2G2 shields the site I of S1-RBD and thereby interfering with the interaction between hexa and S1-RBD. Differently, the G339D mutation in the Omicron S1-RBD confines the FA2G2 between helix 366–371 and loop 372–375. In this condition ${\chi }_1$ populate essentially its trans state stabilising the S1-RBD structure, which is further reinforced by the point mutations S371L, S373P, and S375F.

We can summarize that the multivalent interaction observed in Omicron S1-RBD reflected by the presence of multiple binding modes and the reduced perturbative effect that the FA2G2 glycan moiety applies on the site I of Omicron S1-RBD, may enhance the accessibility of this latest S1-RBD to cellular heparan sulfate (HS), potentially supporting the increased affinity that Omicron S1-RBD shows in comparison to WT S1-RBD towards HS oligosaccharides.