Sialoglycan binding triggers spike opening in a human coronavirus

Expression and purification of trimeric HKU1 spike ectodomains

The sequence of a HKU1-A spike protein (GenBank: ADN03339.1) coding for the ectodomain (residues 12–1266) was cloned into the pCG2 expression vector with an exogenous CD5 signal peptide. At the 3′ end, the coding sequence was ligated in frame with a GCN4 trimerization motif (IKRMKQIEDKIEEIESKQKKIENEIARIKKIK)^35,36, a thrombin cleavage site (LVPRGSLE), an 8-residue long Strep-Tag (WSHPQFEK) and a stop codon. The furin cleavage site at the S1/S2 junction was mutated from RRKRR to GGSGS to avert cleavage of the spike protein (Supplementary Fig. 17). The resulting construct was used for transient expression in HEK293T cells and purified as previously described³⁷. In brief, after incubation of the cells for 5 days, spike glycoprotein was purified from cleared cell culture supernatants by affinity chromatography using StrepTactin beads (IBA) and eluted in 20 mM Tris-HCl, pH 8.0, 150 mM NaCl, 1 mM EDTA, 2.5 mM d-biotin. The W89A mutant protein was produced as described previously¹⁷.

Sample preparation for cryo-EM

For the apo complex, 3 µl of 4.3 µM HKU1 spike trimer was applied to QuantiFoil R1.2/1.3 grids that had been glow-discharged for 30 s on a GloQube (Quorum) at 20 mW power. The sample was applied at 4 °C and 95% relative humidity inside a Vitrobot Mark IV (Thermo Scientific). The grids were then blotted for 7 s with +2 blot force and plunge-frozen in liquid ethane. For the holo complex and W89A negative control, 7 µl of 4.3 μM wild-type or mutant HKU1 spike trimer was combined with 3 µl of 1 mM sugar, resulting in a final spike protein concentration of 3 μM and sugar concentration of 300 μM. The samples were then incubated at room temperature for about 10 min before vitrification, which was carried out as described for the apo sample.

Cryo-EM data acquisition

The apo and holo HKU1 spike samples were imaged on a Thermo Scientific Krios G4 Cryo-TEM equipped with a K3 direct electron detector and a BioContinuum energy filter (Gatan) using EPU 2 acquisition software. The stage was pre-tilted to 30° to improve the orientation distribution of the particles. A total of 4,207 videos for apo spike and 4,065 videos for the holo spike were collected at a super-resolution pixel size of 0.415 Å per pixel, with 40 fractions per video and a total dose of 46 electrons per Ų. Defocus targets cycled from −1.5 to −2.5 μm.

The W89A mutant HKU1 spike incubated with disialoside was imaged on a Thermo Scientific Glacios cryo-TEM instrument equipped with a Falcon 4 direct electron detector using EPU 2 acquisition software. The stage was pre-tilted to 30° to improve the orientation distribution of the particles. A total of 896 videos were collected at 0.92 Å per pixel with 40 fractions per video and a total dose of 42 electrons per Ų. Defocus targets cycled from −1.5 to −2.5 μm. A summary of all data collection parameters is shown in Supplementary Table 2.

Single-particle image processing

For the apo complex, patch motion correction, using an output F-crop factor of 0.5, and patch CTF estimation were carried out in cryoSPARC live³⁸. Micrographs with a CTF estimated resolution of worse than 10 Å were discarded, leaving 4,202 images for further processing. The blob picker tool was then used to select 9,144,772 particles that were then extracted in a 100-pixel box (Fourier binned 4 × 4) and then exported to cryoSPARC for further processing. A single round of two-dimensional (2D) classification was carried out, after which 183,886 particles were retained. Ab initio reconstruction generated one well-defined reconstruction of the closed HKU1 spike protein. Particles belonging to this class were then re-extracted in a 300-pixel box. During extraction, particles were Fourier binned by a non-integer value, resulting in a final pixel size of 1.1067 Å. Subsequently, non-uniform refinement was carried out on the extracted particles with C₃ symmetry imposed³⁹, yielding a reconstruction with a global resolution of 3.3 Å. Subsequently, each particle from the C₃-symmetry-imposed reconstruction was assigned three orientations corresponding to its symmetry-related views using the symmetry expansion job. A soft mask encompassing one S1^A domain was made in UCSF Chimera⁴⁰, and used for local refinement of the expanded particles, yielding a map with a global resolution of 3.8 Å.

For the holo complex, patch motion correction, using an output F-crop factor of 0.5, and patch CTF estimation were carried out in cryoSPARC live³⁸. Micrographs with a CTF estimated resolution of worse than 10 Å were discarded, leaving 4,045 images for further processing. The blob picker tool was then used to select 956,697 particles that were then extracted in a 100-pixel box (Fourier binned 4 × 4) and then exported to cryoSPARC for further processing. Four parallel rounds of 2D classification were carried out, using an initial classification uncertainty value of 1, 2, 4 or 6. Subsequently, the well-defined spike classes were selected from each 2D run and combined. Duplicate particles were then removed, after which 169,728 particles were retained. Ab initio reconstruction generated two classes corresponding to the closed and 3-up spike trimer. Particles from these two classes were used as the input for a second round of ab initio reconstruction that produced two classes corresponding to the 3-up and 1-up spike trimer, although the latter seem to be a convolution of 1-up and closed particles. These two volumes were then used as initial models for a round of heterogeneous refinement. To avoid missing spike particles that may have been removed during initial stringent selection of 2D classes, heterogeneous refinement was carried out on a larger particle stack of 895,888 particles, from which only carbon classes had been removed from the initial stack. Heterogeneous refinement produced two well-defined reconstructions of the 3-up and 1-up conformations. Particles corresponding to the 3-up class were subjected to a single round of 2D classification and the clearly defined spike protein classes were selected. These were then re-extracted in a 300-pixel box. During extraction, particles were Fourier binned by a non-integer value, resulting in a final pixel size of 1.1067 Å. Subsequently, non-uniform refinement was carried out on the extracted particles with C₃ symmetry imposed³⁹, yielding a reconstruction with a global resolution of 3.7 Å. As a result of the apparent heterogeneity in the 1-up sample, an additional round of heterogeneous refinement was carried out on the 895,888-particle stack, using higher-quality initial models, namely the fully refined 3-up map and the 1-up map obtained from the second round of ab initio reconstruction. Heterogeneous refinement produced well-defined reconstructions of the 3-up and 1-up conformations. Particles corresponding to both classes were individually subjected to a single round of 2D classification and the clearly defined spike classes were selected. These were then individually re-extracted in a 300-pixel box. During extraction, particles were Fourier binned by a non-integer value, resulting in a final pixel size of 1.1067 Å. Subsequently, non-uniform refinement was carried out on the extracted particles with C₃ or C₁ symmetry imposed, yielding reconstructions with global resolutions of 3.56 and 4.13 Å for the 3-up and 1-up conformations, respectively. After global refinement, a soft mask encompassing one S1^A domain of the 3-up sample was made in UCSF Chimera. Local refinement was then carried out on the 3-up particles, yielding a map with a global resolution of 4.19 Å. The particles belonging to the 1-up reconstruction were subjected to another round of heterogeneous refinement, which produced two clear reconstructions of the closed and 1-up spike protein. Non-uniform refinement was carried out on both sets of particles with C₃ or C₁ symmetry imposed, yielding reconstructions with global resolutions of 3.68 and 4.68 Å for the closed and 1-up conformations, respectively. For the closed spike protein, each particle from the C₃-symmetry-imposed reconstruction was assigned three orientations corresponding to its symmetry-related views using the symmetry expansion job. A soft mask encompassing one S1^A domain was made in UCSF Chimera⁴⁰, and the symmetry-expanded particles were subjected to masked 3D variability analysis⁴¹. Local refinement was then carried out on the particles belonging to the best resolved cluster, yielding a map with a global resolution of 4.13 Å.

For the W89A mutant HKU1 spike incubated with disialoside, patch motion correction was carried out in MotionCor2 (ref. ⁴²), implemented through Relion version 3.1.1 (ref. ⁴³). The motion-corrected micrographs were then imported into cryoSPARC for patch CTF estimation and further processing steps³⁸. The blob picker tool was used to select 215,843 particles that were then extracted in a 100-pixel box (Fourier binned 4 × 4). A single round of 2D classification was carried out, after which 38,838 particles were retained. Ab initio reconstruction generated one well-defined reconstruction of the closed HKU1 spike protein. Particles belonging to this class were then re-extracted in a 300-pixel box. During extraction, particles were Fourier binned by a non-integer value, resulting in a final pixel size of 1.2267 Å. Subsequently, non-uniform refinement was then carried out on the extracted particles with C₃ symmetry imposed³⁹, yielding a reconstruction with a global resolution of 5.1 Å. Subsequently, each particle from the C₃-symmetry-imposed reconstruction was assigned three orientations corresponding to its symmetry-related views using the symmetry expansion job. A soft mask encompassing one S1^A domain was then made in UCSF Chimera⁴⁰, and used for local refinement of the expanded particles, yielding a map with a global resolution of 5.4 Å.

The ‘gold standard’ Fourier shell correlation (FSC) criterion (FSC = 0.143) was used for calculating all resolution estimates, and 3D-FSC plots were generated in cryoSPARC⁴⁴. To facilitate model building, globally refined maps were sharpened using DeepEMhancer (version 0.13)⁴⁵, as implemented in COSMIC2⁴⁶, or filtered by local resolution in cryoSPARC.

Modelling

Initially, a homology model for HKU1-A spike protein was generated by Phyre 2 (ref. ⁴⁷) with the embecovirus OC43 spike structure  (Protein Data Bank (PDB) 6NZK)²⁴ as template. The HKU1-A spike homology model was rigid body fitted into the apo-state cryo-EM map using the UCSF Chimera⁴⁰ tool Fit in map. The crystal structure of HKU1-A S1^B (PDB 5KWB; ref. ²⁴) was used to replace the equivalent S1^B domain in the homology model owing to clearly wrong homology modelling. Models were refined by carrying out iterative cycles of manual model building using Coot⁴⁸ and real-space refinement using Phenix⁴⁹. The Coot carbohydrate module⁵⁰ was used for building N-linked glycans, which were manually inspected and corrected. The apo state was modelled first, owing to its highest resolution. Subsequently, the closed holo, the holo 3-up and the holo 1-up were modelled in that order, using previous models as a starting point. For the initial holo (closed) S1^A model, Namdinator⁵¹ was used for flexible fitting in a locally refined and unsharpened map for the closed holo S1^A. Model validation was carried out using Molprobity and Privateer^52,53,54.

Elbow⁵⁵ was used to generate ligand restraints for the 9-O-acetylated terminal sialic acid based on the ‘MJJ’ ligand in the OC43 spike cryo-EM structure (PDB 6NZK)¹⁶, after which atom names were manually modified to be consistent with the earlier standard MJJ model and general sialic acid atom numbering, and the O2-attached methyl linker atoms of the original MJJ ligand were trimmed. As there is no standard MJJ–SIA α2,8 linkage defined in software packages used at present, we used molecular dynamics-based restraints (see below) to model this glycosidic linkage of the disialoside. The following restraints were used for the glycosidic linkage between the terminal 9-O-acetylated sialic acid (ligand code MJJ) and the penultimate sialic acid (ligand code SIA) based on the most common solution conformer: bond distance C2–O8 of 1.38 Å (σ of 0.01 Å); bond angles of 109.5° for O8–C2–O6 and for O8–C2–C3, and 114.5° for O8–C2–C1 (all σ of 2.0°); dihedral angles of 295.0° for C1–C2–O8–C8 and of 122° for C2–O8–C8–C7 (both σ of 5.0°).

MD simulations

Starting structures of the molecular systems were built on the basis of the cryo-EM structures of HKU1 (this work) using the graphical interface of YASARA⁵⁶. The N-glycans were attached to the protein on the basis of data from quantitative site-specific N-linked analysis of HKU1 spike protein³¹. Models of the complexes with α-Neu5,9Ac-(2-8)-α-Neu5Ac-OMe were built on the basis of the holo and apo versions of S1^A (residues 14–299). The ligand was positioned manually into the binding site guided by interactions found in PDB entry 6NZK (hCoV-OC43). The HKU1 N1 sequence was taken from GenBank entry NC_006577.2.

Each system was positioned in a periodic rectangular cuboid simulation box (10 Å buffer around the solute) and the AMBER14 force field was selected, which uses ff14SB (ref. ⁵⁷) and GLYCAM06j (ref. ⁵⁸) parameters (including mixed 1–4 scaling). YASARA offers several automated workflows (termed experiments) for system setup. The ‘neutralization experiment’ was used to adjust the protonation states of the amino acids (pH 7.4)⁵⁹ and to solvate the system in 0.9% NaCl solution (0.15 M). The ‘minimization experiment’ (short steepest descent minimization followed by simulated annealing minimization until convergence is reached) was used to remove conformational stress in the system. Simulations were carried out at 310 K using periodic boundary conditions and the particle mesh Ewald algorithm⁶⁰ to treat long-range electrostatic interactions. Temperature was rescaled using a tuned Berendsen thermostat⁶¹. The box size was rescaled dynamically to maintain a water density of 0.996 g ml⁻¹ (‘densostat’ method for pressure coupling)⁶². Position restraints were active during the equilibration phase for at least 3 ns (at the beginning, on all protein heavy atoms, and then only on backbone atoms). To prevent dissociation of the ligand, distance restraints were applied to maintain the critical H bonds for binding (between atoms K80:O and SIA_2:N5; K80:NZ and SIA_2:O1A;  T82:OG1 and SIA_2:O1B). Production simulations were carried out using YASARA with GPU acceleration in ‘fast mode’ (mixed multiple time-step algorithm reaching 5 fs)⁶² on ‘standard computing boxes’ equipped, for example, with one 12-core i9 CPU and NVIDIA GeForce GTX 1080 Ti. Harmonic position restraints (stretching force constant = 1 N m⁻¹) were applied to protein backbone atoms of residues 48–65 and 264–299 of the S1^A system to prevent system rotation in the cuboid box and to deal with the ‘artificially loose end’ at residue 299. The average root mean square deviation of the protein Cα atoms was monitored to check the overall stability of the simulation.

To visualize the glycan coverage of the closed spike protein, the fully glycosylated ectodomain system (590,814 atoms) was simulated with position restraints on backbone atoms of residues 1080–1110 for 250 ns with a performance of about 4 ns per day. Molecular systems based on S1^A alone were smaller (approximately 32,500–56,200 atoms, depending on the size of the N-glycans attached) and were sampled for an accumulated timescale of approximately 20 µs for the Caen1 sequence (apo + disialoside ligand, 5 µs, 6 simulations; holo, 15 µs, 27 simulations) and 52 µs for the N1 sequence (apo, 12 µs, 13 simulations; apo + disialoside ligand, 23 µs, 34 simulations; holo, 17 µs, 22 simulations) with individual simulations reaching up to 1.6 µs. The performance was about 100–200 ns per day. Distances shown in Fig. 5b were calculated from an example trajectory (Extended Data Fig. 10a) between the following atoms: L28:O and S86:OG; P33:CG and F94:CA; T30:O and W89:NE1; S36:N and F74:O; R34:O and -D76:N. Additionally, the solvated disialoside ligand was simulated without the protein using YASARA (general molecular dynamics parameters used as described above) in a cubic box with side length of 37 Å for 10 µs at 310 K using GLYCAM06j parameters. These simulation data were used to identify low-energy conformers of the disialoside ligand, which were used to support the modelling of the reducing-end Neu5Ac residue into the local low-resolution cryo-EM density.

Conformational Analysis Tools ( was used for analysis of trajectory data, general data processing and generation of scientific plots. VMD⁶³ was used to generate molecular graphics.

Analysis and visualization

Spike interface areas were calculated using PDBePISA⁶⁴. Surface colouring of HKU1-A spike protein according to sequence conservation was carried out using Consurf⁶⁵ and visualized in UCSF ChimeraX⁶⁶. The UCSF Chimera MatchMaker tool was used to obtain root mean square deviation values, using default settings. Domain rotations were calculated with CCP4 (ref. ⁶⁷) Superpose⁶⁸. Figures were generated using UCSF ChimeraX⁶⁶ and BioRender.com. Structural biology applications used in this project were compiled and configured by SBGrid⁶⁹.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.