3-O-Sulfation induces sequence-specific compact topologies in heparan sulfate that encode a dynamic sulfation code.

Heparan sulfate (HS) is arguably the most diverse linear biopolymer that is known to modulate hundreds of proteins. Whereas the configurational and conformational diversity of HS is well established in terms of varying sulfation patterns and iduronic acid (IdoA) puckers, a linear helical topology resembling a cylindrical rod is the only topology thought to be occupied by the biopolymer. We reasoned that 3-O-sulfation, a rare modification in natural HS, may induce novel topologies that contribute to selective recognition of proteins. In this work, we studied a library of 24 distinct HS hexasaccharides using molecular dynamics (MD). We discovered novel compact (C) topologies that are populated significantly by a unique group of 3-O-sulfated sequences containing IdoA residues. 3-O-sulfated sequences containing glucuronic acid (GlcA) residue and sequences devoid of 3-O-sulfate groups did not exhibit high levels of the C topology and primarily exhibited only the canonical linear (L) form. The C topology arises under dynamical conditions due to rotation around an IdoA → GlcN glycosidic linkage, especially in psi (Ψ) torsion. At an atomistic level, the L → C transformation is a multi-factorial phenomenon engineered to reduce like-charge repulsion, release one or more HS-bound water molecules, and organize a bi-dentate "IdoA-cation-IdoA" interaction. These forces also drive an L → C transformation in a 3-O-sulfated octasaccharide, which has shown evidence of the unique C topology in the co-crystallized state. The 3-O-sulfate-based generation of unique, sequence-specific, compact topologies indicate that natural HS encodes a dynamic sulfation code that could be exploited for selective recognition of target proteins.

Introduction

Heparin/heparan sulfate (Hp/HS), two members of the superfamily of glycosaminoglycans (GAGs), play important roles in numerous molecular and cellular processes. Classically, Hp/HS have been studied extensively as modulators of coagulation; however, growing evidence indicates that these biopolymers may possess major clinical value as anti-infectives and anti-cancer agents [1], [2], [3], [4]. Hp/HS have also been implicated as major regulators in diverse processes such as axonal growth [5], neurofibrillary tangles [6], and native immune response [7], [8].

The diverse biological activities of Hp/HS are mediated by different proteins including those that are secreted or membrane-bound, either within or outside cells. The recent release of the first draft of GAG interactome suggests more than 500 protein partners of these highly anionic polymers [9]. Among these, a number of proteins preferentially bind to 3-O-sulfate (3OS)-containing sequence(s) in Hp/HS including antithrombin [10], fibroblast growth factor receptor [11], glycoprotein D [12], cyclophilin B [13], neuropilin-1 [5], heparin cofactor II [14], tau glycoprotein [6] and spike glycoprotein [15]. Each of these sequences are distinct owing to the structural pattern flanking the 3OS group, which suggests the possibility of a “sulfation code”. It is expected that more proteins will be identified in the future as preferred targets with their own 3OS-based sulfation code.

The 3OS group is introduced in the final step of Hp/HS biosynthesis and is generally rare in most biological tissues [16]. For example, clinically used Hp has only one in three chains carrying the antithrombin-binding 3OS sequence, which approximates to one 3OS-containing residue in nearly 150 residues. On the other hand, human follicular fluid contains ∼6 % GlcNS3S residues [17]. Interestingly, nature has designed seven isoforms of 3-O-sulfotransferase (HS3ST-1, -2, -3A, -3B, 4, -5, and -6) to introduce just one 3OS group in an HS chain. In comparison, only one and three isoforms are known to date for other sulfotransferases (STs), which introduce O-sulfates at the 2- and 6-positions, respectively. Further, a high level of regulation is apparent in the biosynthesis of Hp/HS despite being a fully template-independent process. For example, the expressions of different STs are tissue- and time- dependent [18], [19], [20], [21]. Such a regulation appears to help generate a fairly consistent phenotype, as noted in the expression of STs in a specific tissue.

The seven isoforms of HS3ST display significant differences in substrate specificities. Whereas HS3ST-1 preferentially generates the antithrombin-binding HS sequence, HS3ST-3B generates binding motifs for herpes simplex virus-1 glycoprotein D [16]. The specificities of other isoforms have been presented in terms of generation of either antithrombin- or glycoprotein D- binding motifs. Yet, obvious differences in specificities among the seven isoforms have also been clearly noted [22], [23], [24], [25]. Alternatively, the HS3ST family of enzymes appears to encompass some overlap of specificity, while also exhibiting important differences.

The rarity of 3OS groups in nature coupled with an abundance of HS3ST isoforms that present an overlapping, and yet varying substrate specificity, make this functional group an enigmatic structural modification. A burning question is why has nature engineered 3-O-sulfation? Growing evidence points to its role in generating unique “sulfation codes” that are recognized by different proteins in a sequence-specific manner. We reasoned that 3-O-sulfation may induce unique conformational properties that endow an Hp/HS sequence with novel protein recognition properties. More specifically, we reasoned whether incorporation of a 3OS group in some Hp/HS sequences, but not all, generate novel conformational states?

In this work, we present molecular dynamics-based evidence that the presence of 3OS group induces a unique compact (‘C’) form distinct from the well-established linear, helical form known for decades. We studied a library of 24 Hp/HS hexasaccharides with varying distribution of sulfate groups to comprehensively understand their conformational states. We find that when a 3-O-sulfated GlcN residue is flanked by two iduronic acid (IdoA) residues, a unique ‘C’ shaped topology is significantly populated. The C topology arises from inversion of a preferred glycosidic torsion, which facilitates reduction in electrostatic crowding, release of one or more bound waters, and bi-dentate sequestration of a cation by iduronic acids. MD studies of a 3OS-containing octasaccharide also generate a compact form, similar to that observed for hexasaccharides. Interestingly, literature-reported crystal structure of the 3-O-sulfated octasaccharide shows the occupation of this compact topology in the bound state confirming the existence of distinct, sequence-dependent topologies, different from the cylindrical rod form known so far. Overall, the 3OS-based generation of unique, sequence-specific, compact topologies indicate that natural HS encodes a dynamic sulfation code that could be exploited for selective recognition of target proteins.

Methods

Molecular dynamics (MD) simulations

The library of 24 variably sulfated hexasaccharides was built using GLYCAM web-server () and files for performing MD simulations were downloaded from the server. The structure of each sequence was first confirmed, which included the order of residues, the bond lengths and angles, and initial torsional angles. For MD simulations, ring puckering of IdoA was restrained using weak torsional restraints and puckering occupancies were ascertained before and after simulations using the BFMP method [26], [27]. MD simulations were performed using AMBER18 with GLYCAM06J force field, essentially as described earlier [28], [29]. Each sequence was neutralized with an appropriate number of Na+ counterions and then solvated using TIP3P water with periodic boundary conditions such that a minimal distance of 12 Å was maintained between solute atoms and the nearest box edge. A two-step minimization was carried out for each hexasaccharide, wherein the solute was restrained with a harmonic potential of 100 kcal/(mol Å2) in the first step. Here, water molecules were “relaxed” for 500 iterations of steepest descent and 1500 cycles of conjugate gradient minimization. In the second step, conjugate gradient minimization was performed for 2500 iterations without any restraint. In the next step, the entire system (Hp/HS sequence in TIP3P water box) was equilibrated by implementing a Berendsen thermostat to achieve a constant temperature of 300 K first and then a pressure of 1 atm. Finally, the system was equilibrated without any restraints. All minimization and equilibration steps utilized a 2 fs timestep for a total of 1 ns. MD production runs were performed for 1 μs during which the coordinates were collected and stored every 10 ps. Covalently bonded hydrogen atoms were constrained using the SHAKE algorithm, as described in AMBER [28], [29].

Analysis of MD production runs

Torsional (Φ,Ψ) space, solvation, hydrogen bonding, end-to-end distance (EED), sodium cation (Na+) binding, and potential energy calculations were performed using the CPPTRAJ package with AmberTools14 [30]. IUPAC definitions for HS were used for glycosidic torsional angles, moving from the non-reducing end to reducing end: Φ = O5-C1-O1-C4 and Ψ = C1-O1-C4-C5 [31]. Torsional density maps were prepared using the ‘displot’ function in Seaborn for bivariate frequency distributions with the default parameters. Selection of Φ,Ψ data points corresponding to each minima was enabled by the findpeaks function in SciPy. The proportion of the secondary minimum at each linkage was calculated by generating a 2D histogram using the histogram2d function from the NumPy library and dividing by the total number of frames (100,000). EED was calculated as mentioned previously to approximate the length of the GAG, respectively [32]. Na+ binding was modeled by establishing distance criteria, where a Na+ was observed within 3.5 Å of C6 of Residues 2 and 4. Both distances were measured in parallel using the DISTANCE command in CPPTRAJ throughout the trajectory. Na+ binding was measured as a binary function, that is, 1 if present at the time the trajectory was captured (every 0.01 ns) and 0 if not. To evaluate changes in like-charge repulsions between carboxylates and sulfates, the coulombic potential energy term was isolated from the force field for the functional groups of interest using the energy command in CPPTRAJ. Default parameters were used in calculating electrostatics. Only terminal oxygens of residue III 3-O-sulfate, residue IV carboxylate, and residue V 6-O-sulfate were included in the calculations. Number of hydrogen bonds (H-bonds) per functional group were calculated by including only hydrogen bond acceptors and donors, as shown in Table S2. Water loss during conformational changes relied on changes in these H-bonds along the MD trajectory. An H-bond was considered formed if the distance between donor (D) and acceptor (A) atoms was less than 3.5 Å and the D-H∙∙∙∙A angle was 180 ± 60°.

MD simulations of the 3-O-sulfated octasaccharide

Octasaccharide GlcNS6S(1 → 4)IdoA2S(2SO) (1 → 4) GlcNS6S3S (1 → 4) IdoA2S(1C4) (1 → 4) GlcNS6S (1 → 4) IdoA2S (1C4) (1 → 4) GlcNS6S (1 → 4) IdoA2S(1C4) was built using the GLYCAM web server following the process for hexasaccharides as described above. Residue puckers were restricted to those observed in the co-crystal structure (PDB:3INA). Puckering occupancies were assessed using the BFMP method [26]. The minimization, equilibration, and production MD protocols for this octasaccharide were identifical to those as described for the hexasaccharide library. The analysis of MD production runs for EED, Φ/Ψ, water loss, Na+ binding analyses, etc. was performed in a manner identical to that of the hexasaccharides. A binary function allowed for the monitoring of three separate Na+ binding interactions, which included residue III to I, residue V to III, and residue V to I.

Results

Rationale behind the design of the hexasaccharide library

To understand whether the presence of a 3OS group endows a Hp/HS chain with unique conformational states, we studied a library of 24 Hp/HS hexasaccharides. Theoretically, four uronic acid (UA) and twelve glucosamine (GlcN) monomers are possible owing to variations at the 2-, 3-, and/or 6-positions,a which can generate 48 unique Hp/HS disaccharide building blocks, or 110,592 theoretically possible unique hexasaccharides of AcidNRE–AmineRE class (NRE = non-reducing end; RE = reducing end). An identical number of hexasaccharide sequences are possible for the AmineNRE–AcidRE class. This represents an astounding configurational diversity for such a small sequence. This diversity is further exacerbated by the different puckers possible for IdoA, of which 1C4 and 2SO forms are the most populated [31].

It is now well recognized that most protein-binding Hp/HS sequences are tetra- to octa- saccharide long [12], [14], [15], [33], [34], [35]. Hence, we chose hexameric sequences, which allow placing the 3OS group in the middle of the chain. A single 3OS group was positioned on the internal GlcN residue III, which had variations at 2- (acetyl or sulfate) and 6- (sulfate or unsubstituted) positions (Fig. 1). Three other variations were introduced. First, residues II and IV were either GlcA or IdoA with the former in 4C1 form, while the latter in either 1C4 or 2SO form. We chose to include IdoA because it has been shown to exhibit a nearly equal proportion of the two puckers [36], which could offer a better glimpse into interesting conformational dynamics. Second, residue V was either GlcNS6S or GlcNAc6S. Third, the 3-position of residue III could be either O-sulfated or —OH; the latter forming the six control sequences. The residues at the two termini, i.e., residues VI and I, were maintained invariant as GlcA and GlcNS6S, respectively. Thus, the library could be thought as six sub-libraries, each containing one non-3OS-containing control sequence and three 3OS-containing sequences with variations in the internal tetrasaccharide (Fig. 1). Of the 3OS-hexasaccharides, 12 contained IdoA (either 1C4 or 2SO), while 6 had GlcA (4C1). We reasoned that this library design afforded significant variation in configurational and conformational space. It is important to note that some of these 24 sequences are not known to exist in mammals as yet. The known substrate specificity of the C5 epimerase requires a GlcNS on the non-reducing side of the acceptor GlcA, which disfavors generation of the GlcNAc → IdoA disaccharide units [37]. However, these sequences are likely to be synthetically accessible and could serve as interesting chemical biology and/or drug discovery leads.

Fig. 1

The library of Hp/HS hexasacchrides. The library can be thought of in terms of 6 subsets derived from variations in the sulfation pattern and uronic acid composition (IdoA (1C4 or 2SO) or Glc(4C1)) of the four middle residues (Residues II through V). The terminal (reducing end (RE) and non-reducing end (NRE)) residues were maintained invariant as GlcNS6S and GlcA, respectively. Residue III was either 3-O-sulfated or unmodified at 3-position, which served as a control (Ctrl) for each subset. Symbols for each monosaccharide follow the Symbol Nomenclature for Glycans (SNFG) system, while presence of sulfate groups is represented by a red circle (see legend at the bottom). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

MD simulations indicate occupancy of a novel compact topology in a sequence-specific manner

Recently, a number of MD studies on GAGs in water have been performed using three major forcefields, including GLYCAM, CHARMM and GROMOS [29]. Considerable insights into glycosidic torsional angles, ring puckers, inter-atomic interactions, intermolecular interactions and global shapes of GAG chains have been derived from these studies [38]. In fact, recent studies on large number of GAG disaccharides have further confirmed the classic works by Mulloy and co-workers [31] that glycosidic torsions Φ and Ψ vary within a relatively small range irrespective of sulfation pattern [39], [40]. This gives rise to the canonical ‘helical’ topology of Hp/HS chains, which is used as the primary Hp/HS topology in docking [41], [42], [43], although pucker and glycosidic bond flexibility is typically allowed in these studies [44], [45]. Yet, Hp/HS chains are thought to be linear cylinders, which implies that “sulfation code” is wholly configurational with only a marginal room for IdoA conformations.

To quickly understand whether 3-O-sulfation impacts overall Hp/HS topology, we analyzed a parameter of global ‘shape’ called the end-to-end distance (EED, see definition in Fig. S1A) [32]. The helical form of the common Hp/HS hexasaccharide (i.e., (IdoA2S–GlcNS6S)3) exhibits an EED of 26.4 Å in its static state [46], which is similar to that derived from its MD simulations [32]. We calculated the EEDs at every 1 ns, i.e., for 1000 frames over 1 µs, for each of the 24 sequences. The results indicated that the majority of sequences displayed EEDs of ∼27 Å. This implies that the majority of sequences occupy the helical, linear topology (Fig. S1B). However, several sequences exhibited major reduction in EEDs. Specifically, major reductions were observed for sequences 2, 3, 12, and 14. Interestingly, each of these sequences contain IdoA residues at positions II and IV, either in the 1C4 (Seq. 2, 3 and 12) or 2SO form (Seq. 14).

To understand this better, we analyzed the EEDs of MD frames as a function of time. While the EEDs were essentially invariant (∼27 Å) for all control non-3OS sequences, the profiles for sequences 2, 3, 12 and 14 were significantly different (Figs. 2 and S2). Sequences 2, 3 and 12 displayed a transition to an alternate conformational form of ∼16 Å EED, while sequence 14 preferred a ∼13 Å form (Fig. 2A). Sequence 2 also displayed an additional more compact form with an EED of ∼8 Å. This suggested that this group of four sequences exhibit a high tendency to occupy non-linear, perhaps non-helical, topologies. A comprehensive survey of EED dynamics over the MD runs showed that sequences 1, 4–6 and 13 also sample non-linear conformational forms (Fig. S2); however, the extent of their compaction is not as high as sequences 2, 3, 12 and 14.

Fig. 2

Profiles of end-to-end distances (EEDs). EED as a function of MD simulation time (A) and frequency of occurrence (B) are shown for control 1 and sequences 2, 3, 12 and 14. Corresponding profiles for control sequences 2 → 6 and sequences 4 → 11, 13, and 15 → 18 are shown in supplementary Figs. S2 and S3. EEDs were calculated every 1 ns for the 1 µs MD simulations (i.e., 1000 frames). For frequency distributions, the bin size was set to every 0.2 Å.

To further parse the differences, we calculated the frequencies of MD frames with similar EEDs (every 0.2 Å). Whereas each of the six control sequences displayed a primarily unimodal EED profile with maximal length of ∼27 Å, 3OS containing sequences displayed a non-unimodal distribution (Fig. S3). The most interesting among these were sequences 2, 3 and 12, which exhibited an unusual bimodal or trimodal distribution with EEDs of ∼27, ∼16 and ∼8 Å (Fig. 2B). Interestingly, sequence 14 was observed to be in its compact form for most of its MD simulation time with an EED of ∼13 Å. Other 3OS containing sequences, i.e., 1, 4–6 and 13, also displayed a bimodal distribution, albeit with a different distribution of EEDs (Fig. S3). In contrast, 3OS containing sequences 7 → 9 and 16 → 18, which have GlcA residues flanking 3OS residue, display a primarily unimodal distribution similar to their corresponding controls 3 and 6 (Fig. S3).

The above results indicate that when a 3OS group is present, a hexasaccharide samples multiple topologies depending on the structure of the sequence. Hexasaccharides with 3OS group and IdoA residues considerably populate non-linear, compact topologies. In fact, differing levels of compaction are observed for different types of sequences. Further, hexasaccharides with 3OS group and GlcA residues do not exhibit any tendency to populate one or more C form(s). Thus, the results indicate that the precise sequence of a 3-O-sulfated hexasaccharide determines the tendency and extent of compaction. This was the first evidence that some 3OS groups may induce an unusual compact form in Hp/HS sequences. Such unique, three-dimensional, structural motifs likely form the “sulfation code”.

3-O-sulfation selectively modifies a Psi torsional angle

To identify the site of origin for compaction, we studied torsions Φ and Ψ for each glycosidic linkage of the 24 sequences (Fig. S4). The IdoA/GlcA → GlcN and GlcN → IdoA/GlcA linkages of all control sequences C1 through C6 exhibited Φ/Ψ preferences that were identical to that observed by NMR of heparin dodecasaccaride having the common sequence (IdoA2S–GlcNS6S)6 [46]. This preference was also noted for 3OS sequences (see Table S1 and Fig. S4). However, sequences 2, 3, 12 and 14 displayed additional Φ/Ψ preferences, especially at the internal IdoA → GlcN linkage (i.e., IV → III residues). For example, sequence 2 displayed a second Φ/Ψ preference at −80.8°/66.4°, which represents a massive change in Ψ of nearly 165° from its linear torsional preference (-80.8°/-99.2°) (Fig. 3A). Likewise, sequences 3, 12 and 14 also displayed similar torsional preferences (Fig. 3A, Table S1). Extraction of geometries corresponding to these Φ/Ψ preferences from MD simulations shows that this major change in preference transforms the helical, linear form (EED ∼27 Å) to the compact form with EEDs of ∼16 Å in sequences 2, 3 and 12, and 13 Å in seq 14 (see models in Fig. 3A). Although the exact topologies of these compact forms are structurally unique, we refer to each of these as ‘C’ topologies, which are radically different from the linear ‘L’ form.

Fig. 3

Maps of glycosidic torsions pinpoint the origin of the compact (C) topology. (A) Ramachandran plot (Φ and Ψ ranges) for the IV → III glycosidic linkage of sequences 2, 3, 12 and 14 (see Fig. S4 of all 24 sequences). IV → III linkage corresponds to the middle IdoA → GlcN linkage of Hp/HS hexasaccharides. Blue regions refer to highly preferred Φ and Ψ. Sequences 2, 3, 12 and 14 exhibit an unusual Φ/Ψ of ∼-80°/∼66°, which induces a unique, compact (C) topology (shown on the right). Also shown is the linear, helical topology corresponding to the Φ/Ψ of ∼-80°/∼-105°. (B) Populations (in %) at the Φ/Ψ = ∼-80°/∼66° for all 24 hexasaccharides calculated using NumPy. Non-reducing terminus residue is labeled VI and reducing end residue is labeled I. Molecular models were generated in PyMOL. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

To gain some quantitative insights, we calculated the torsional occupancies using EED histograms, which corresponded to the occupancy of the compact form (Φ/Ψ = ∼-80°/∼66°) (Fig. 3B). The population of these Φ/Ψ torsions was found to be small, or non-existent, in the large majority of sequences, except for sequences 2, 3, 12 and 14. Whereas the C topology was populated nearly 30–70 % in sequences 2, 3, 12 and 14, it was 0–10 % in others. GlcA flanking hexasaccharides were essentially devoid of this conformer (Fig. 3B and S4). Thus, the results convey the role of HS structural pattern in populating the novel compact topology.

Atomistic origin of the L → C transformation

We sought to reveal the foundational factors that drive the topological switch from helical L shape to the C form. At the thermodynamic level, multiple enthalpic and/or entropic factor(s) may be responsible. We hypothesized that negative enthalpic factors, such as electrostatic repulsion due to sulfate group crowding, may initiate the process. This could then be followed by favorable entropic factors, such as release of Hp/HS-bound water molecules, which may synergize to drive the L → C topological transformation.

To test the above conjecture, we analyzed the L topologies of sequences 2, 3, 12 and 14 and their corresponding controls C1, C4 and C5. Installation of a 3OS group appeared to increase electrostatic crowding due to close proximity of 3-O-sulfate, carboxylate, and 6-O-sulfate of residues III, IV and V, respectively (Fig. 4A). In contrast, these groups appeared ‘isolated’ in the C form. To gain better insight, we plotted the repulsive electrostatic force calculated by GLYCAM using a CPPTRAJ script between these anionic groups across the entire MD trajectory [30] (Fig. 4B). The plots show that rotation in Ψ at the IV → III linkage induces a significant drop in Coulombic energy of ≥20 kcal/mol for sequences 2, 3, 12 and 14. In fact, the timing of these drops correlated well with the observed topological transformation from the L to C shape. Further, such decreases in electrostatic energy were not noted for the corresponding control sequences 1, 4 and 5 as well as other 3OS sequences (Fig. S5). Thus, like-charge repulsion may be one of the factors initiating topological transition in the distinct group of 3OS sequences.

Fig. 4

Bad electrostatics contributes to the induction of the compact form. (A) Sequence 2 exemplifies this phenomenon. The 3-O-sulfate, 5-carboxylate, and 6-O-sulfate groups of residues III, IV and V, respectively (red spheres) in sequence 2 generate electrostatic crowding, which is relieved by rotation around Ψ of the IV → III glycosidic linkage. (B) Profiles of changes in electrostatic energy in the above dense anionic cluster as a function of the MD simulation time. The decrease in electrostatic energy coincides with the observed change in EED as a function of time (see Fig. 2). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Release of HS-bound water molecules favors compaction

The release of protein-bound water molecules during protein folding has been documented extensively [47], [48]. Solvent water molecules provide many roles in protein folding, one of which is to enhance entropic energy resulting in stabilization of the folded state. For GAGs, the role of such water molecules, if any, remains unknown. To assess whether water molecules are involved in L → C topological transformation, we first calculated the number of intermolecular hydrogen bonds (H-bonds) at the individual residue level as a function of MD simulation time [9], [30], [32], [40] (Fig. S6). Residue II of sequences 2, 3, 12 and 14 displayed a significant decrease in the number of intermolecular H-bonds that coincided with the topological transition. This phenomenon was absent in control sequence 1, 4 and 5 as well as in other 3OS sequences (Fig. S6). Interestingly in addition to residue II, residue IV of sequence 14 also showed a loss in number of intermolecular H-bonds.

To understand the atomistic origins of water loss, we repeated the analysis with focus on the Hp/HS partners involved in intermolecular H-bonds. Although several polar atoms/groups were involved, the carboxylates and O-sulfates of residues II, III and IV contributed most H-bonds, which remained consistent over the MD time period. However, a decrease of ≥2 intermolecular H-bonds was observed for sequences 2, 3, 12 and 14 upon formation of the compact form, specifically at IdoA carboxylates (Figs. 5 and S7). This loss in bound waters was not observed for their corresponding control sequences as well as for other 3OS sequences (Figs. 5 and S7). The results support our conjecture that the novel compact topology formation manifest from both enthalpic and entropic contributions. Interestingly, a careful review shows that the H-bonding profile of the 3OS groups of sequences 2, 3, and 12 increases by ∼1H-bonds. This is absent in other 3OS sequences, such as 9 and 18 (Fig. S8).

Sodium cation binding stabilizes the compact topology

Since the early 1990s, the binding of metal cations to GAGs has been the subject of many studies [49], [50], [51], [52]. A common observation has been that divalent cations, e.g., Ca2+, preferentially bind to carboxylate groups of GAGs. A more recent MD-based investigation of chondroitin sulfate disaccharides concluded that Ca2+ prefer to engage GlcA carboxylates rather than O-sulfate groups [53].

We hypothesized that additional favorable enthalpic contributions may arise from metal ion binding to contribute to the L → C transformation. To assess this possibility, we evaluated the MD trajectory of sequence 3 for the presence of metal ion, i.e., Na+, within the coordination sphere (≤3.5 Å) of each carboxylate and sulfate group. Interestingly, one Na+ ion was found to be engaged by the carboxylates of IdoA residues II and IV in a series of MD frames (Fig. 6A). This represents a bidentate IdoA–Na+–IdoA interaction, which arises only if the two carboxylates are brought in close proximity due to compaction. Thus, the corresponding control sequence 1 did not exhibit such a phenomenon.

Fig. 6

Stabilization of compact form arising from formation of bidentate IdoA-Na+-IdoA interaction. (A) The carboxylates of IdoA residues II and IV of sequence 3 generate a favorable site for Na+ resulting in stabilization of the compact (C) form. (B) Of the 24 Hp/HS sequences, a stable, consistent bidentate IdoA-Na+-IdoA interaction occurs only for sequences 2, 3, 12 and 14.

To assess whether this phenomenon is generally applicable, we developed a CPPTRAJ script for identifying the bidentate COO––Na+–COO– interaction for all 24 Hp/HS sequences. Fig. 6B profiles the results of this analysis. Although several MD trajectories did show the presence of such a bidentate interaction, only those corresponding to sequences 2, 3, 12, and 14 displayed a highly consistent phenomenon that lasted in a prolonged fashion. In fact, the timing of the formation of the bidentate IdoA–Na+–IdoA interaction correlated well with the L → C transition (see Fig. 5).

Fig. 5

Intermolecular H-bonds between solvent water and specific groups of HS residues. H-bonds formed over the 1 µs MD trajectory are shown for residues II (A), III (B) and IV (C) of control 1 and 3OS sequences 2 and 14 (see Fig. S7 for sequences 3 and 12). Note residues II and IV are IdoA, while residue III is 3-O-sulfated GlcN. Intermolecular H-bonds were calculated using CPPTRAJ.

It is important to note that the tendency to occupy the C form-related Φ/Ψ space at the IV → III glycosidic linkage is theoretically accessible for all Acid → Amine linkages (see Fig. S4, Table S1). This implies that longer oligosaccharide may present opportunities for the formation of additional bidentate IdoA–Na+–IdoA interactions, which may or may not compete with each other. Overall, the results support the hypothesis that the C form is likely to be stabilized by IdoA–Na+–IdoA interaction, just as loss of bound waters favor the L → C transition.

Kinetics of the L → C transformation suggests a concerted process

A key question arises whether the three different contributors to the L → C transformation, including electrostatic repulsion arising from multi-anion crowding, intermolecular H-bonding between IdoA carboxylates and water, and simultaneous engagement of Na+ to two IdoA residues, work independently or in a concerted manner? To assess this, we identified the MD time frames that display the maximal changes in the two descriptors of the L → C transformation, i.e., EED (Fig. 7A) and Ψ of IV → III linkage (Fig. 7B). Whereas sequence 2 exhibits start of topological transformation at about 444 ns, the start times for sequences 3, 12 and 14 were ∼566, ∼729 and ∼104 ns. Interestingly, these time frames for completion of L → C transformation for each of these sequences, as evidenced by EED and Ψ (IV → III), were very short and in the range of 2 to 5 ns. More importantly, the changes in intermolecular H-bonding (Fig. 7D) and bidentate Na+ complexation (Fig. 7E) occur over nearly the same time frame. In fact, the bidentate Na+ binding phenomenon appears to arise at the end of the transformation phase following which the C topology remains stable for an extended period of time. A closer review of changes in the electrostatic repulsion between 3OS of residue III, carboxylate of residue IV, and 6OS of residue V indicates a more variable and subtle process (Fig. 7C) that occurs over a longer time period and is initiated much before the rapid L → C glycosidic reorientation. This suggests that bad electrostatic crowding may be a contributor but not essential for L → C transformation. Overall, the kinetics of L → C transformation suggests a concerted process.

Fig. 7

The L-->C transition is concerted. Each descriptor of the L → C transition including EED (in Å), Ψ of IV → III linkage (in degrees), electrostatic energy between anions of residues III and IV (in kcal/mol), number of intermolecular H-bonds between carboxylates of IdoA and water, and number of Na+ bound simultaneously to IdoA residues II and V (bi-dentate interaction) as a function of MD simulation time is shown for sequences 2 (A), 3 (B), 12 (C) and 14 (D). Only the brief window of the L → C transition, wherein significant changes in each of these descriptors, is shown for clarity. Profiles shown are 10-point moving averages calculated for 100,000 frames over 1 µs simulations.

MD of an HS octasaccharide supports the principles enunciated from HS hexasaccharides

To assess whether the novel compact topology identified for 3OS HS hexasaccharides can be found in nature, we reviewed 3D structures of proteins in complex with oligosaccharides in the protein data bank (rcsb.org). A pentasaccharide containing one 3OS group has remained the major ligand studied to date. This includes fondaparinux in complex with adeno-associated virus-DJ [54], bacterial polysialyltransferase [55], platelet factor 4 [56], staphylococcus TarM [57], thrombospondin-1 [58] and antithrombin [34]. Unfortunately, fondaparinux has a GlcA-GlcN-IdoA motif that is different from the motifs studied in this work. Crystal structures of several longer oligosaccharides are reported; however, none of these contain the key 3OS group. Only one crystal structure of a heparin lyase in complex with an octasaccharide with 3OS group (PDB: 3INA) has been reported to date [59]. Although the structure of this 3INA octasaccharide (GlcNS6S-IdoA2S-GlcNS3S6S-IdoA2S-GlcNS6S-IdoA2S-GlcNS6S-IdoA2S) is different from the library of our hexasaccharides, it displays a characteristic compact topology shown by the four 3OS sequences in our MD simulations (Fig. 8A). Additionally, the compact topology was stabilized by a key bidentate IdoA-Cation-IdoA interaction in which an arginine residue served as a cation (Fig. 8B). This provides strong support for the existence of unique C topologies in nature.

Fig. 8

MD simulations of a 3-O-sulfated octasaccharide. A crystallized co-complex (PDB: 3INA) of 3OS octasaccharide in complex with a heparin lyase shows the characteristic compact topology observed for the four 3OS-containing HS hexasaccharides. (A) PyMOL-genereated molecular models of the 3INA octasaccharide from the co-crystal structure and MD-derived compact topologies of Seq 2, 3, 12 and 14. (B) Arg83 forms a key salt-bridging interaction with adjacent IdoA COO– groups, stabilizing the compact form of the 3INA octasaccharide. (C and D) MD studies on the 3INA octasaccharide showing EED as a function of time (C) and frequency distribution (D). EEDs (shown for every 1 ns) and frequency distributions (over 1 µs) were calculated for all 100,000 frames using a bin size of 0.16 Å. (E) The sequence of the 3INA octasaccharide and the Φ/Ψ space for each of its linkages derived from MD run (Res I is RE; VIII is NRE). Blue regions correspond to higher probability of occurrence, while red dots correspond to crystallographically observed Φ/Ψ values for the 3INA octasaccharide. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

To assess whether the octasaccharide exhibits similar behavior in MD, we built the starting linear topology using GLYCAM-Web and subjected it to 1 µs MD simulation, as described for the hexasaccharide library. Visual inspection of the 100,000 MD frames of the octasaccharide showed that ∼100 ns into the simulation, the linear topology transformed into a compact form. This L → C transformation was accompanied with an EED decrease from ∼33 Å to ∼9 Å (Fig. 8C). Frequency distribution of the different conformers formed during the MD run indicated that the compact topology was the most abundant form (Fig. 8D). Closer inspection of the different glycosidic torsions revealed that the compact topology was generated by two Ψ reorientations at the III → II and V → IV linkages (Fig. 8E). Of these, the crystal structure presents only the latter Ψ reorientation. In fact, the V(IdoA2S) → IV(GlcNS6S) linkage exhibits a Ψ of 63.5°, which is remarkably close to the value presented by 3OS-containing hexasaccharide sequences at ∼63° (Fig. 8E and Table S1).

To further address the similarity of atomistic mechanisms behind the topological transformation in 3INA octasaccharide to 3OS hexasaccharides, we analyzed the two key components that help stabilize the C topology, i.e., loss of intermolecularly H-bonded water molecules and gain of bidentate IdoA-Na+-IdoA interaction. Similar to the 3OS-containing hesasaccharides, the results show that IdoA2S residues I and III of the 3INA octasaccharide lose about 2 to 3 water molecules as a results of the topological transformation around 100 ns (Fig. 9A). IdoA2S residues V and VII do not exhibit this phenomenon (Fig. S10), which correlates well with no noticeable changes in glycosidic torsions near the non-reducing terminus of the 3INA octasaccharide (not shown). Likewise, the MD profile shows a robust IdoA-Na+-IdoA bidentate interaction for the 3INA octasaccharide in which three IdoA carboxylates take part (Fig. 9B). Of these, residues I and III form more consistent interaction between them than residues I and V. Overall, atomistic foundations behind the MD profile of the octasaccharide are very similar to that of sequences 2, 3, 12 and 14, thereby supporting the mechanics and generality of the L → C transformation.

Fig. 9

MD simulations of a 3INA octasaccharide supports principles enunciated from the study of HS hexasaccharides. The 3INA octasaccharide displayed a loss of 2–3 water molecules (A) and gain of bidentate IdoA2S-Na+-IdoA2S interaction (B) as the L → C transformation sets in around 100 ns. Intermolecular H-bonds were calculated using CPPTRAJ for all frames of 1 µs MD trajectories. Shown here are residues I and III, while other residues are shown in Fig. S9.

Discussion

In this study, we discovered that a small group of 3OS-containing hexasaccharides significantly occupied one or more compact topologies (C forms) in the dynamical state, while also populating the linear, helical topology (H form) to some extent. Whereas the C form(s) were observed for sequences containing two IdoA residues, not all 3OS hexasaccharides with two IdoA residues sampled the C form to a significant extent. Also, none of the GlcA-containing 3OS hexasaccharides populated the C form. Likewise, none of the control sequences, which were devoid of a 3OS group, populated any of the C forms. These results were also found to hold parallelly for a naturally occurring octasaccharide containing one 3OS group. Our discovery challenges the field’s long-held understanding that irrespective of the sulfation pattern, HS exhibits a linear helical form with the possibility of small “kink(s)” around IdoA residues [27], [60], [61], [62], [63], [64], [65], [66], [67]. Our work shows for the first time that the conformational diversity of HS is much higher and detailed studies with configurationally diverse 3OS-containing oligosaccharides are critically needed to truly appreciate the static and dynamic nature of the HS sulfation code.

Within the class of 3OS sequences populating the compact form, certain sulfation pattern, e.g., sequences 3 and 12, may present one C topology, whereas others, e.g., sequences 2 and 14, may exhibit more (Fig. 2, Fig. 3). Alternatively, compact topologies are induced in a sequence-specific manner. This observation was also supported in MD studies of the 3INA octasaccharide, which contained one 3OS group and two 2-O-sulfated IdoA residues. Thus, 3INA octasaccharide and hexasaccharides 2, 3, 12 and 14 are structurally different. Yet, the 3INA octasaccharide exhibited a compact topology, which corresponded to the same change in Ψ as the hexasaccharides (Fig. 8). The observation that the 3INA octasaccharide populates a compact topology in the co-crystallized state with a protein lends experimental foundation to the concept of sequence-specific induction of compact topology.

One can think of natural HS as a repository of the most diverse “sulfation codes” that engineer selectivity of protein recognition. Yet, the sulfation code is not biosynthesized in the classic template-driven manner, as demonstrated for nucleic acids and proteins. To explain its origin, it has been suggested that the sulfation code could arise from the controlled spatiotemporal expression of the superfamily of N- or O-sulfotransferases [21], [68], [69]. Of these, the group of 3-O-sulfotransferases are most probably the most important. Not only is the 3OS group rare in HS [17], [70], [71] and has been implicated in mediating selective recognition of proteins [14], [15], [34], [6], nature utilizes seven different isoforms to engineer just one structural modification. Thus, the 3OS group has been thought of as the key driver of the sulfation code. Yet, how it engineers a unique sulfation code has remained elusive.

A priori, one may theorize that an HS sulfation code is generated by either 1) the unique sequence of the oligosaccharide and manifested in the linear helical form and/or 2) conformational effects arising from occupation of different ring puckers, e.g., 1C4, 2SO, etc., in the linear helical form. Available co-crystal structural data support the first hypothesis wherein HS sequences exhibit linear, helical form with a complete turn every-four residues [60], [61], [62], [63], [64], [65], [66]. Many NMR-based studies in solution support the second hypothesis, wherein neighboring sulfation patterns impact ring puckering equilibria to engineer a sulfation code within the linear, helical form [27], [67]. Yet, none of these clearly identify the criticality of a 3OS group.

Interestingly, some clues on how a 3OS group could contribute to this phenomenon were already available in co-crystal structures of fondaparinux, the five-residue sequence with one 3OS group in the center of the chain. This pentasaccharide displays a “kink” in multiple co-crystal structures with antithrombin and platelet factor 4 [34], [56], which indicates a slight deviation from the linear, helical form. This departure from the linear form may not necessarily originate during protein binding because recent MD studies with 3OS-containing di- and penta- saccharides in the free state demonstrate the sampling of the kinked forms [39], [72]. With this work on longer oligosaccharides, it is clear that sequence-specific generation of compact topologies, i.e., dynamical sulfation codes, requires the presence of certain 3OS groups. It is important to note that 3OS group may not be the only contributor, especially when considering binding to proteins. For example, the ability of IdoA to switch ring puckers upon binding to antithrombin has been well-documented and likely to be an additional mechanism to explain selectivity of binding [73]. It is also important to not discard the possible contributions of IdoA2S, which are known to be important [74], [75], [76] and may or may not be similar to IdoA.

The compact topology-based sulfation code requires the formation of one or more bi-dentate IdoA-Cation-IdoA interaction, which stabilizes L → C transformation through a reorientation of Ψ at an appropriate IdoA → GlcN linkage (Fig. 6). Our observation that reorientation occurs at an IdoA → GlcN linkage, and not alternate GlcN → IdoA/GlcA linkages, agrees with earlier literature highlighting the flexibility of these linkages in HS [67], [77]. Interestingly, sequences 2, 3 and 12 carried IdoA in 1C4 form, while sequence 14 had IdoA in 2SO form (Fig. 1), which further conveys the importance of conformational bias in guiding overall topology.

At an atomistic level, our work resulted in delineating the underlying forces that direct the L → C transformation. These forces include 1) like-charge repulsion arising from sulfate group crowding (Fig. 4), 2) loss of HS-bound water molecules upon compaction (Fig. 5), and 3) favorable interactions of the bidentate IdoA-Cation-IdoA type following L → C transition (Fig. 6). Of these, electrostatic crowding while being clearly synchronous for sequences 2 and 3, was less synchronous for sequences 12 and 14 (Fig. 7). This suggests that the latter two forces, i.e., energetic gain from water loss and bidentate interactions, are likely to be more dominant contributors.

Compaction of proteins, which is associated with loss of bound water molecules, has been known to provide entropic benefits [47], [48]. In fact, it is known that protein folding relies on formation of a hydrophobic core that necessitates expulsion of water molecules. In fact, in a two-state protein folding model, stabilization of the ‘native’ state occurs following expulsion of “interstitial water and consolidation of the [hydrophobic] core” [78]. Our work demonstrates a rather similar phenomenon for HS sequences that undergo compaction. Expulsion of two or more water molecules occurs essentially simultaneously with bidentate IdoA-Cation-IdoA interaction, which could be thought of as consolidation of the core. To the best of our knowledge, this is the first report on the expulsion of bound water molecules during a topological transformation. In contrast, the cation binding property of HS has been well documented [49], [50], [51], [52]. Yet, the bidentate IdoA-Cation-IdoA “code” is novel.

A key question originates from this study on oligosaccharides whether larger HS chains exhibit compaction. Compaction or bending has been observed in HS and other glycosaminoglycans [79], [80], [81]. Yet, the extent of compaction is difficult to study because of polydispersity and heterogeneity. In fact, details on points of compaction or bending remain unknown. This work pinpoints atomistic origins and mediators of compactions, which could now be assessed at the polymeric level.

This work also conveys the importance of studying other members of the hexasaccharide library. Our library was designed to include IdoA and GlcA residues in a uniform manner. Natural and synthetic libraries could include a mix of IdoA, IdoA2S, GlcA and GlcA2S. It is very likely that these expand the concept of “dynamic sulfation code” including variations around the L → C transition.

Finally, our work does not attempt to offer one-to-one correlation between the substrate specificity of an HS3ST and a specific dynamic sulfation code. With the understanding gained so far, it is difficult to derive such as one-to-one correlation, even if it exists. Our work provides a reasonable foundation for supporting such a possibility and validating such as hypothesis will require much advanced work.

Overall, this study proposes a fundamental change in the concept of “sulfation code”. Rather than the sulfation code being a static configurational entity, it proposes a more dynamic conformational entity. Although we have studied 24 hexasaccharides and one octasaccharide, it is important to note that this is still a small number of HS sequences. Our work conveys the importance of understanding the topological terrain of 3OS oligosaccharides. Equally important would be to identify unique HS sequences that manifest the dynamic sulfation code and employ it in developing therapeutics or diagnostics against specific protein targets.

Author contributions

SGH performed analysis of MD experiments and contributed to paper writing; BN set up all MD experiments, reviewed and contributed to paper writing; URD supervised and guided research, revised and finalized the paper, and acquired funding. All authors have read and agreed to the publication of this work.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.