Integration of Experimental and Computational Methods to Characterize Glycoproteins

Integration of experimental and computational methods to characterize glycoproteins Reyes Núñez-Franco,1 Francesca Peccati,1 Ángel Torres-Mozas,1 and Gonzalo JiménezOsés1,2,* 1Center for Cooperative Research in Biosciences (CIC bioGUNE), Basque Research and Technology Alliance (BRTA), Bizkaia Technology Park, Building 800, 48160 Derio, Spain 2Ikerbasque, Basque Foundation for Science, 48013 Bilbao, Spain Introduction Glycans encompass a wide range of fundamental biological functions, spanning form structural roles to energy storage to regulatory purposes. The rich variety of their three-dimensional structures −arising from the diversity of the carbohydrate building blocks and their linear or branched arrangements− makes them the most information dense biopolymers on Earth.1 Understanding the complex functions of glycoproteins requires knowledge of their interactions at an atomistic level (molecular recognition), underlining the importance of structural data that is often elusive to experimental techniques owing to the high flexibility of glycans. Computer simulations have been successfully filling this gap for over two decades, complementing and interpreting experimental results from a variety of structural techniques (NMR, X-ray crystallography and cryo-EM being prevalent), and providing insights into the function of glycoproteins and glycoenzymes through a variety of approaches encompassing quantum and molecular mechanics.1,2 The global challenge represented by the recent SARS-CoV-2 pandemics has catalyzed a deeper integration of experimental and computational efforts opening the door to an unprecedented use of simulations for the study of glycoproteins and showing their power in accurately describing how glycans can fine-tune large-scale conformational transitions and complex recognition phenomena.3 In this contribution, we emphasize this advancement presenting selected examples from the literature that cover a wide range of glycoprotein and glycoenzyme families, discussing how different types of computer simulations have been used to expand our knowledge on glycoprotein functions. Rather than providing a description of the technical aspects of computer simulations, which have been recently reviewed elsewhere,1 we focus on examples that show how the integration of simulation and experiment is key to provide an atomistic interpretation of experimental data. Molecular dynamics simulations of glycoproteins Unlike other biopolymers such as proteins and nucleic acids, carbohydrate polymers are not characterized by well-defined secondary and tertiary structures, but present complex geometry distributions better represented by conformational ensembles rather than by single conformations.1 Their intrinsic flexibility, which makes them essentially invisible to most structural biology techniques, has made molecular dynamics (MD) the family of computational methods of choice for shedding light on their structure and function as a complement to NMR, crystallographic and cryo-EM data.1,2 Molecular dynamics for biomolecular systems propagate the classical equations of motion and provide a picture of the microscopic time evolution of the system. In the context of glycoproteins, the conformational preferences of glycans are determined not only by their nature but also by a variety of interand intramolecular glycan–glycan and glycan–protein interactions, providing a highly complex scenario for computational modeling.1 Traditionally, NMR chemical shifts have been used as structural restraints for MD simulations of biopolymers, penalizing conformations that show inconsistencies with experimental data.4,5 Distance constraints derived mainly from NOESY NMR experiments were initially used in combination with conventional restrained molecular dynamics (rMD), and were later replaced by weighted time-averaged restraints (MD-tar) which were proven to reduce constraint energies and deviations from the target distances compared to standard rMD.6 However, chemical shift restraints obtained for flexible regions arise from a weighted average of the underlying conformational ensemble, in such a way that it might be impossible for an individual geometry to satisfy these average constraints.7,8 The last decade has seen a significant shift in perspective on the role of MD simulations in the study of glycoproteins. The advent of graphics processing units (GPUs) has changed the field of molecular dynamics of biological systems, extending both the size and time scales accessible to routine simulations (millions of atoms and microseconds).9 This, together with the development of specialized force-fields for carbohydrates (GLYCAM10 and CHARMM11 to cite a few) has made it possible to perform long unbiased simulations, providing insights on the effect of glycosylation on the stability, flexibility and conformational preferences, and function of glycoproteins. In several of the examples that will be presented in this chapter, the power of MD simulations lies in the full characterization of the low-energy conformational ensembles accessible to glycoproteins beyond the static picture provided by cryo-EM and Xray structures. Despite these significant advancements, there is still no unique approach to the simulation of these complex systems as different force fields12 and ensembles (canonical, microcanonical, grand canonical, etc.) can provide different results, and the lack of convergence of extensive simulations suggests that insufficient conformational sampling can bias the results, leaving an open challenge for the development of the next generation of MD approaches for studying glycoproteins.12 It is generally accepted that glycans increase protein stability enhancing their solubility and resistance to proteolysis. However, a complex relationship exists between the stabilizing effect and the polarity and position of the glycans; coarse-grained MD simulations have been successfully employed to describe how glycosylation fine-tunes protein stability and can exert a stabilizing or destabilizing effect according to the local environment.13–15 For instance, analysis of 63 SRC Homology 3 Domain variants revealed how the thermostabilizing effect of glycans depends on the number of attached chains and the position of the glycosylation sites, while it has little dependence on the size of the attached glycans.14 A different study by the same authors focused on the glycosylation of the MM1 protein (a miniaturized model of GcMAF, a serum factor that stimulates the phagocytic activity of macrophages) showed that glycosylation can have a globally destabilizing effect by inducing short-range glycan−protein interactions that break long-range protein−protein interactions.15 Atomistic models have been recently used by Jiménez-Barbero and co-workers to unravel NMR data on the glycoprofile of the IgE high-affinity receptor (FcεRIα) expressed in human HEK293 cells.16 Determining glycosylation patterns (composition and relative abundance of glycan epitopes) under physiological conditions is an open challenge owing to the intrinsic heterogeneity of glycans. Integration of NMR data with MD simulations revealed the structural and stabilizing roles of the glycans attached to residues Asn39 and Ans139. Analysis of the conformational ensembles suggests that glycosylated Asn39 bridges the two immunoglobulinlike domains through interactions at the sugar residues at the α(16) and α(13) branches, while the latter remains buried between the domains and might contribute to stabilize the protein structure.16 Similarly, MD simulations have been used by the same group to interpret NMR data in a systematic study of the interactions of the SARS-CoV-2 spike protein RBD (receptor binding domain) expressed in human HEK293F cells with a wide panel of human lectins (Gal-3, Gal- 7, Gal-8, DC-SIGN, MGL, Siglec-10 and Siglec-8), identifying novel glycoepitopes.17 As the RBD presents two glycosylation positions (Asn331 and Asn343), lectin recognition on either site leads to different geometries of the 1:1 lectin:RBD complexes, which were explored through MD simulations. Apart from the expected interactions at the canonical carbohydrate binding site, these simulations revealed a large dynamical network of glycan-receptor contacts matching NMR data and highlighting the complexity of the recognition phenomenon.17 Microsecond MD simulations have been used for extensive conformational analyses of several glycoproteins, highlighting how coupling between glycans and protein motions is able to finetune their conformational ensembles. For instance, it has been shown that glycosylation shifts the conformational ensemble of the yeast multidomain protein disulfide isomerase (yPDI) toward more open, extended, structures facilitating access to the catalytic site and highlighting the collective functional role of glycosylation.18 When comparing the conformational ensembles obtained from MD simulations of the glycosylated and non-glycosylated systems (Figure 1), highly populated (i.e. low energy) basins corresponding to open conformations are observed in the glycosylated system that are not observed in the ensemble of the nonglycosylated counterpart. These simulations are key to identify this functional role as the available crystal structures of yPDI and homologous proteins (represented as symbols in Figure 1) all fall in a narrow low-energy region and provide only limited information on the conformational ensemble accessible to yPDI. Figure 1. Representation of the open and closed conformations in the MD-derived conformational ensembles of of yPDI.18 (A) 2D projection of the free energy surface of glycosylated yPDI. Symbols show the positions of known crystal structures of yPDI and other homologous proteins within the conformational space spanned by MD simulations. The most populated energy minima with open structures (1 and 2) are shown for the glycosylated system. Reproduced with permission from the American Chemical Society. The trimeric HIV-1-envelope (Env) spike protein has been the object of extensive MD simulations to assess the collective effect of glycosylation on accessibility and conformational ensembles.19 Several virus proteins are highly glycosylated, and the biological role of the socalled glycan shield (the layer of surface glycans) is primarily to mask the surface of the protein from immune recognition. The Env spike protein presents a high glycan density (about half of its mass) and MD simulations revealed a synergy between protomer and glycan motions which coordinate to shield the protein surface.19 This dynamic information complements the static picture obtained from experimental structures, which have been mostly determined for the isolated glycoprotein by X-ray crystallography. These structures were solved with a single protomer in the asymmetric unit, which makes them perfectly symmetric. MD simulations incorporate thermal motion and environmental effects, thus breaking the symmetry and revealing the coupling between glycan and protein flexibility.19 Glycans have also been shown to locally fine-tune the conformational ensemble of the V3 loop of the gp120 envelope trimer of a HIV virion and consequently modulate its binding affinity to CD4.20 Glycans surrounding the V3 loop induce a pronounced narrowing of the conformational ensemble relative to the non-glycosylated loop in a composition-dependent manner, highlighting their functional role in shaping an intrinsically flexible region that escapes analysis with structural biology techniques. MD simulations were used to compare the conformational ensembles assumed by the V3 loop (Figure 2) in i) the non-glycosylated form, ii) a glycosylated form with Man-9 glycans at five nearby asparagine residues (Glycosylated5-glycans), and iii) a mono-glycosylated form with a single Man-9 glycan linked at Asn295 (Glycosylated295). The ensembles derived from MD simulations show markedly different conformational preferences in the three cases, with a sharper though different geometry distribution of the two glycosylated forms. Figure 2. Conformations spanned by the V3 loops according to the nature and number of glycans.20 A) Conformations spanned by the V3 loop in the three forms: non-glycosylated, Glycosylated5-glycans and Glycosylated295; B) conformations spanned by the V3 loop in the nonglycosylated form, C) conformations spanned by the V3 loop in the Glycosylated5-glycans form, and D) conformations spanned by the V3 loop conformations in the Glycosylated295 form. Reproduced with permission from Public Library of Science. The SARS-CoV-2 pandemics has been a turning point in the field, putting under the spotlight how extensive MD simulations can provide fundamental insight into the dynamics and function of the glycosylated spike protein. The timely release of complete and fully glycosylated models of the spike protein21 based on early cryo-EM and glycomic/glycoproteomic data22–24 was crucial for the rapid development of a rich body of literature based on MD simulations which provided precious insight into the shielding and functional roles of glycans,25–31 enabling intensive studies for identifying drug candidates, which will be discussed in the following section.3,32,33 (Figure 3). Figure 3. Milestones in experimental (grey) and MD simulation studies (blue) of the SARSCoV-2 spike protein. Indeed, an early landmark report by Amaro and co-workers on multiple microsecond, full-atom MD simulations of the spike protein revealed not only the role of extensive glycosylation in shielding the surface of the trimeric spike, but also the structural role of glycans at positions Asn165 and Asn234 in promoting transitions between the “closed” and “open” conformations of the RBD, key for the infection process (Figure 4).34 In the “closed” conformation, the RBD is inaccessible and shielded by the glycans, while in the “open” conformation, it is available for binding to the angiotensin-converting enzyme 2 (ACE2) receptor. Figure 4. Representation of the glycan shield in the region of the RBD interacting with ACE2.34 Panels A and D show the evolution of the accessible surface area of the RBD in the “open” and “closed” conformation as a function of the probe radius size obtained from averaging over multiple MD simulations. Panels B, C, E and F show the side and top view of an overlay of evenly spaced MD frames from simulations of the spike protein in the “open” and “closed” RBD conformations. Reproduced with permission from the American Chemical Society. Computational assesment of glycoenzyme activity CAZymes are active enzymes that display a key role in the biosynthesis, modification, binding, and catabolism of carbohydrates and are found in all living organisms. As of October 2022, the CAZy database (http://www.cazy.org) compiles the available knowledge on 173 glycoside hydrolases, 115 glycosyltransferases, 42 polysaccharide lyases, 20 carbohydrate esterases, 9 ligninolytic and 7 lytic polysaccharides mono-oxygenases families. Glycoside hydrolases and glycosyltransferases are the two largest CAZymes families, which have gained the attention of many researchers in the last years. They are responsible for the cleavage and formation of glycosidic bonds, respectively, which often involve the posttranslational modification of glycoproteins. Due to the high flexibility and diversity of both proteins and carbohydrates, modeling these enzymes is one of the biggest challenges for computational chemical biology. Protein flexibility, substrate conformations, enzyme conformational switches and electronic reorganization upon covalent bond cleavage/formation are key factors affecting the way these enzymes work.50 Mechanistically, glycosyltransferases (GT) and glycoside hydrolases (GH) are classified as “retaining” or “inverting” according to whether the stereochemistry at the donor’s anomeric bond is retained (α→α or β→β) or inverted (α→β or β→α), respectively (Figure 7). Figure 7. Canonical retaining and inverting mechanisms for glycoside hydrolases (R2 = H) and glycosyltransferases. Glycoside hydrolases catalyze the hydrolysis of glycosidic bonds, whose retaining profile can take place through either the classical Koshland double displacement, the substrate-assisted or the neighboring group participation mechanisms. Conversely, inverting glycoside hydrolases operate through a single displacement mechanism. These different mechanisms have been widely analyzed with the aid of computational techniques.50–57 On the other hand, glycosyltransferases catalyze the monosaccharide transfer from an activated sugar donor to an acceptor. Whereas their inverting profile is clearly established to follow a SN2-like mechanism, their retaining one has been extensively discussed.58 Among other techniques, hybrid quantum mechanics/molecular mechanics (QM/MM) calculations have been crucial for shedding light on this debate.59 QM/MM simulations allow treating a small part of the system involved in the formation and breaking of chemical bonds (i.e. the substrate(s) and the catalytic residues) with computationally expensive but very accurate QM methods (typically in the Density Functional Theory framework), whereas the rest of the system (typically the protein and its first solvation shell) is treated with faster but less accurate MM approaches, typically a general or specialized force field. An example of this methodology applied to a glycoside hydrolase is show in Figure 8. Figure 8. Representation of the enzymatic model used for the QM/MM study of the hydrolysis of 4-nitrophenyl-β-D-fucopyranoside catalyzed by Thermus thermophilus β-glycosidase.60 (Reproduced with permission from Frontiers Media SA). The inverting enzyme N-acetylglucosaminyltransferase I was the first glycosyltransferase enzyme studied by (QM/MM) methods.61 The authors calculated the generally accepted concerted SN2 mechanism in which residue Asp291 acts as a general base, with an activation energy for the proposed reaction mechanism estimated to be ~19 kcal mol-1 consistent with the experimentally observed activation barriers (15-25 kcal mol-1).62 QM/MM combined with metadynamics calculations is often the technique of choice for studying the formation and breaking of chemical bonds by fully solvated enzymes and including thermal motion.63 Metadynamics is an enhanced sampling technique for MD simulations that estimates the Free Energy Landscape (FEL) of a process (i.e. a chemical reaction or a conformational change) as a function of several selected degrees of freedom, also referred to as collective variables (CV) (Figure 9).64 This method accelerates the detection of rare events, helps the system escaping from local free-energy minima and allows characterization of transition states.65 Figure 9. Free energy surface for the formation of the covalent glucosyl–enzyme intermediate in 1,3-1,4-β-glucanase calculated with QM/MM metadynamics.55 The minimum energy path (MEP) is shown as a white dashed line. Contour lines are separated by 4 kcal mol-1. The green dots on the axes labels represent the bonds being formed/broken. R/R’: reactants; TS: transition state; P/P’: products. (Reproduced with permission from the American Chemical Society). From the methodological point of view, Rovira and co-workers have provided useful insights into different aspects that should be considered when running this type of simulations.66 The authors have shown the importance of equilibrating the starting structure (typically a crystallographic one) with classical MD prior to the QM/MM MD calculations, paying special attention to the sugar conformations at this initial step. They also showed how to select and control collective variables for the metadynamics simulations. Many studies describing the use of QM/MM-based calculations for characterizing both glycosyltransferase67–77 and glycoside hydrolases50–57 have been reported, including the exploration of different reaction mechanisms, both for retaining and inverting profiles, as well as the estimation of the reaction activation barriers which are generally in good agreement with experimental measurements. Examples of these investigations are shown in Table 1. Table 1. Examples of glycoside hydrolase (GH) and glycosyltransferase (GT) enzymes studied by QM/MM calculations. Family Enzyme Experimental activation barrier (kcal mol-1) QMM/MM calculated activation barrier (kcal mol-1) Preferred mechanism Ref. Retaining GT Mannosyl glycerato synthase 18.7 17 Front-face stepwise 76 Inverting GT O-fucosyl transferase 16.0 17–19 SN2 77 Retaining GH β-L-Arabino furanosidase 16.4 17 Double displacement 51 Inverting GH GH43 endoarabinnase 14.0 17 Single-displacement 53 In addition to deciphering reaction mechanisms, MD simulations have been widely employed for analyzing other important aspects of the function of these enzymes. As mentioned before, conformational switches are essential for understanding the way these enzymes work. HurtadoGuerrero and co-workers have extensively studied the glycosylation preferences of polypeptide N-acetylgalactosaminyltransferases (GalNAc-T). These enzymes control the transfer of a GalNAc moiety from UDP-GalNAc onto Ser/Thr residues of proteins. These enzymes are constituted by a catalytic and a lectin domain connected by a flexible linker. MD simulations have shed light on the dependence of the glycosylation preferences of those enzymes dictated by the rotation of the lectin domain, which is exerted by the flexible linker located between the lectin and catalytic domains.78,79 They also demonstrated for the first time that a subset of GalNAc-T isoforms showed a GalNAc binding site in the catalytic domain. They identified the molecular basis for GalNAc-T4’s longand short-range preferences. Conventional MD simulations supplemented with a combination of steered and Umbrella Sampling MD simulations were used to model the unbinding process of the glycopeptide from the catalytic domain of GalNAc-T4, as well as for studying the interdomain-loop flexibility and the occurring hydrogen bond and hydrophobic interactions.80 These simulations were carried out for the wild-type GalNAc-T4 and a truncated-loop variant; simulations showed a flexible openlike structure for the mutant opposite to the more static closed-like structure observed crystallographically for the wild-type (Figure 10). Figure 10. Atomic fluctuation analysis (Cα only) of wild-type (left) and a lectin flexible loop (LFL)-truncated mutant (right) in complex with a glycopeptide and UDP obtained from 0.5 μs MD simulations.80 Rigid enzyme regions are shown as thin purple and blue ribbons whereas more flexible residues are shown as red thick ribbons. Reproduced with permission from the American Chemical Society. Park and co-workers studied the catalytic activities of Salmonella SseK1 and SseK2 glycosyltransferases through a combination of molecular docking, conventional and Gaussianaccelerated MD (GaMD) simulations, NMR, crystallography, kinetic assays and isothermal titration calorimetry (ITC) experiments among other techniques.81 Their results suggested that the helix-loop-helix domain found in SseK1 and SseK2 determines substrate specificity and that the lid domain is the one regulating the opening of the active site. MD simulations have also been employed to decipher the high specificity of human 8oxaguanine-DNA glycosylase (hOGG1) towards the non-canonical nitrogen base 8oxaguanine, one of the most common DNA lesions.82 This work showed the importance of active-site residues not only for the catalytic reaction but also for driving a conformational modification of the nucleotide bound to the active site. Moreover, MD simulations performed at different pH values were used for evaluating the pHdependence of the enzymatic and non-enzymatic functionalities of an endo-β-D-glucuronidase (heparanase, HPSE).83 A combined study including MD simulations, molecular docking and protein angular dispersion analysis (PADω)84 on the apo and holo forms provided insight into the optimal pH values (5.0 and 6.0, respectively) driving the required structural changes for the enzymatic activity of HPSE. The authors demonstrated that at these pH values, the protonated forms of HPSE show optimal distances between the E225 and E343 catalytic residues to enhance catalytic activity. A flap region in HPSE was identified as a modulator of an opening/closing movement of the binding cavity; in turn, acidic pH favors the closing of the binding cavity residues thus promoting substrate binding. Molecular docking analysis also pointed HSPE at these conditions as the optimum conformations to accommodate HSPG. Another key aspect determining glycoenzymes’ activity is their thermostability. Thermostable enzymes are ideal for industrial applications due to their durability, tolerance to a wider pH range and resistance to organic co-solvents. MD simulations are among the most used methods to guide thermostability engineering of proteins since they allow the identification of highly flexible regions. These regions are identified as engineering targets for stability enhancement. Niu and co-workers employed the Spatial Compartmentalization of Mutational Hotspots (SCMH) method, which combines alignment of homologous sequences, spatial structure compartmentalization and MD simulations, for designing a thermostabilized variant of βglucanase from Bacillus terquilensis.85 Acevedo and co-workers also used MD simulations to identify flexible regions of a new cold-adapted xylanase glycoside hydrolase.86 Those regions are identified as the putative structural elements defining the thermostability of this enzyme. Mutagenic libraries were designed to stabilize the protein by introducing mutations in the previously identified flexible regions. They identified twelve mutant clones with improved T50 (the temperature at which the enzyme loses 50 % of its activity), the best mutant showing a T50 increase of 4.3 ºC. Chen and co-workers employed a similar technique for identifying key flexible positions suitable for proline substitutions with the aim of thermostabilizing αglucosidase from Xanthomonas campestris (XgtA).87 A quintuple mutant of XgtA wild-type (V167P/A177P/A220P/A230P/T345P), showed a 3-fold increase in the value of half-life at 45 ºC.i MD simulations are very useful to unveil the role of solvation in molecular recognition. For instance, the recognition mechanism of human thrombospondin 1 (HsTSR1) protein by Caenorhabditis elegans POFUT2 (CePOFUT2) glycosyltransferase was studied with microsecond MD simulations, which showed that the binding event is mediated by dynamically structured water molecules, alleviating in this way the need for multiple direct protein–protein interactions and allowing substrate promiscuity.88 The role of interstitial water was assessed by simulating also the HsTSR1-CePOFUT2 complex in absence of explicit solvent but maintaining only the crystallographic waters. The disruption in the complex communication between HsTSR1 and CePOFUT2 together with the higher stability of the complex when explicit water molecules are considered in the simulation confirmed the important role of specific water molecules in this recognition process. GH13 is the major glycoside hydrolase family among GHs. Many of the enzymes belonging to this family show, apart from their normal catalytic activity towards the hydrolysis of polysaccharides, transglycosylation activity to various degrees. While all the enzymes belonging to this group cleave the O-glycosidic bond of polysaccharide chains, water or other carbohydrates are subsequently transferred to an acceptor in enzymes showing hydrolytic or transferase activity, respectively. Arreola-Barroso and co-workers studied the amino acid residue composition defining the reaction specificity between hydrolysis and transglycosylation activities.89 This approach was based on the conservation of coevolving residues. The authors identified the composition of amino acid contacts that were frequently found in either hydrolases or transferase enzymes and underrepresented in the other group. Hydrolytic and transglycosidic contacts were defined as those that are more frequently found in hydrolases and transferases, respectively. Of note, some hydrolases showed only a few hydrolytic contacts, indicating that the fraction of transglycosidic contacts might play a more prominent role in defining the specific catalytic function. Using this methodology, the same authors identified mutagenic hotpots modulating the hydrolysis/transglycosylation (H/T) for two members of the GH13 family, namely Thermotoga maritima α-amylase (TmAmyA) and Thermotoga maritima glucanotransferase (TmGTase). Different mutants for both TmAmyA and TmGTase were constructed to modulate the H/T ratio, discovering a TmAmyA triplet mutant (K98P/D99A/H222Q) and a TmGTase variant (M279N) which produced a two-fold decrease and a five-fold increase in the H/T ratio, respectively. MD simulations were employed to explain the effect of these mutations on the activities in structural terms, finding out changes in local and global flexibility that can account for the modified T/H ratios. Conclusions The synergistic integration of experimental and computational approaches is one of the landmarks of modern glycosciences. Traditionally, the heterogeneity, structural and chemical properties and high flexibility of glycoproteins – with vast implications for host-pathogen interactions, immune response, cell adhesion, inflammation, etc. – has been largely neglected in structural biology studies (i.e. crystallography) due to both the use of recombinant (i.e. nonglycosylated) proteins and the “transparency” (i.e. lack of enough resolution) of glycans to those techniques. However, the development of increasingly powerful cell culturing, mass spectrometry, NMR, cryo-EM and computer modelling techniques is providing and exponential growth of the knowledge about the composition, structure, dynamics and function of glycosylated proteins beyond the amino acid sequence. The dramatic increase in computing power experienced in the last years – evidenced by the advent of exascale supercomputers and GPU-accelerated software – now allows performing unbiased molecular dynamics simulations in the microsecond and even sub-millisecond timescale, progressively approaching the experimental time regimes. The unprecedented efforts carried out by the chemical modelling community as a whole during the COVID-19 pandemic, mostly directed to elucidating the role of the multiple glycans attached to SARSCoV-2 spike protein, and simultaneously finding potential drug candidates using it as a therapeutic target, have revealed the paramount importance of considering the often neglected “glycan shield” of proteins. In parallel, quantum mechanics combined with enhanced sampling molecular dynamics are providing growing evidence about the operational modes of glycosylated enzymes, and enzymes cataylzing (de)glycosylation reactions; this is allowing to (40) Kapoor, K.; Chen, T.; Tajkhorshid, E. Posttranslational Modifications Optimize the Ability of SARS-CoV-2 Spike for Effective Interaction with Host Cell Receptors. Proc. Natl. Acad. Sci. U. S. A. 2022, 119, e2119761119. (41) Huang, Y.; Harris, B. S.; Minami, S. A.; Jung, S.; Shah, P. S.; Nandi, S.; McDonald, K. A.; Faller, R. SARS-CoV-2 Spike Binding to ACE2 Is Stronger and Longer Ranged Due to Glycan Interaction. Biophys. J. 2022, 121, 79–90. (42) Cao, Y.; Choi, Y. K.; Frank, M.; Woo, H.; Park, S. J.; Yeom, M. S.; Seok, C.; Im, W. Dynamic Interactions of Fully Glycosylated SARS-CoV-2 Spike Protein with Various Antibodies. J. Chem. Theory Comput. 2021, 17, 6559–6569. (43) Ramírez-Salinas, G. L.; Martínez-Archundia, M.; Correa-Basurto, J.; GarcíaMachorro, J. Repositioning of Ligands That Target the Spike Glycoprotein as Potential Drugs for Sars-Cov-2 in an in Silico Study. Molecules 2020, 25, 5615. (44) Maghsoudi, S.; Taghavi Shahraki, B.; Rameh, F.; Nazarabi, M.; Fatahi, Y.; Akhavan, O.; Rabiee, M.; Mostafavi, E.; Lima, E. C.; Saeb, M. R.; Rabiee, N. A Review on Computer‐aided Chemogenomics and Drug Repositioning for Rational COVID‐19 Drug Discovery. Chem. Biol. Drug Des. 2022, 100, 699–721. (45) Sahoo, B. M.; Ravi Kumar, B. V. V.; Sruti, J.; Mahapatra, M. K.; Banik, B. K.; Borah, P. Drug Repurposing Strategy (DRS): Emerging Approach to Identify Potential Therapeutics for Treatment of Novel Coronavirus Infection. Front. Mol. Biosci. 2021, 8, 628144. (46) Creutznacher, R.; Maass, T.; Veselkova, B.; Ssebyatika, G.; Krey, T.; Empting, M.; Tautz, N.; Frank, M.; Kölbel, K.; Uetrecht, C.; Peters, T. NMR Experiments Provide Insights into Ligand-Binding to the SARS-CoV-2 Spike Protein Receptor-Binding Domain. J. Am. Chem. Soc. 2022, 144, 13060–13065. (47) Alvarado, W.; Perez-Lemus, G. R.; Menéndez, C. A.; Byléhn, F.; de Pablo, J. J. Molecular Characterization of COVID-19 Therapeutics: Luteolin as an Allosteric Modulator of the Spike Protein of SARS-CoV-2. Mol. Syst. Des. Eng. 2022, 7, 58–66. (48) Zhao, P.; Praissman, J. L.; Grant, O. C.; Cai, Y.; Xiao, T.; Rosenbalm, K. E.; Aoki, K.; Kellman, B. P.; Bridger, R.; Barouch, D. H.; Brindley, M. A.; Lewis, N. E.; Tiemeyer, M.; Chen, B.; Woods, R. J.; Wells, L. Virus-Receptor Interactions of Glycosylated SARS-CoV-2 Spike and Human ACE2 Receptor. Cell Host Microbe 2020, 28, 586601.e6. (49) Clark, A. J.; Gindin, T.; Zhang, B.; Wang, L.; Abel, R.; Murret, C. S.; Xu, F.; Bao, A.; Lu, N. J.; Zhou, T.; Kwong, P. D.; Shapiro, L.; Honig, B.; Friesner, R. A. Free Energy Perturbation Calculation of Relative Binding Free Energy between Broadly Neutralizing Antibodies and the Gp120 Glycoprotein of HIV-1. J. Mol. Biol. 2017, 429, 930–947. (50) Coines, J.; Raich, L.; Rovira, C. Modeling Catalytic Reaction Mechanisms in Glycoside Hydrolases. Curr. Opin. Chem. Biol. 2019, 53, 183–191. (51) McGregor, N. G. S.; Coines, J.; Borlandelli, V.; Amaki, S.; Artola, M.; Nin-Hill, A.; Linzel, D.; Yamada, C.; Arakawa, T.; Ishiwata, A.; Ito, Y.; van der Marel, G. A.; Codée, J. D. C.; Fushinobu, S.; Overkleeft, H. S.; Rovira, C.; Davies, G. J. Cysteine Nucleophiles in Glycosidase Catalysis: Application of a Covalent β-lArabinofuranosidase Inhibitor. Angew. Chemie Int. Ed. 2021, 60, 5754–5758. (52) McGregor, N. G. S.; Artola, M.; Nin-Hill, A.; Linzel, D.; Haon, M.; Reijngoud, J.; Ram, A.; Rosso, M. N.; Van Der Marel, G. A.; Codeé, J. D. C.; Van Wezel, G. P.; Berrin, J. G.; Rovira, C.; Overkleeft, H. S.; Davies, G. J. Rational Design of Mechanism-Based Inhibitors and Activity-Based Probes for the Identification of Retaining α-l-Arabinofuranosidases. J. Am. Chem. Soc. 2020, 142, 4648–4662. (53) Meelua, W.; Wanjai, T.; Thinkumrob, N.; Oláh, J.; Mujika, J. I.; Ketudat-Cairns, J. R.; Hannongbua, S.; Jitonnom, J. Active Site Dynamics and Catalytic Mechanism in Arabinan Hydrolysis Catalyzed by GH43 Endo-Arabinanase from QM/MM Molecular Dynamics Simulation and Potential Energy Surface. J. Biomol. Struct. Dyn. 2022, 40, 7439–7449. (54) Jin, Y.; Petricevic, M.; John, A.; Raich, L.; Jenkins, H.; De Souza, L. P.; Cuskin, F.; Gilbert, H. J.; Rovira, C.; Goddard-Borger, E. D.; Williams, S. J.; Davies, G. J. A βMannanase with a Lysozyme-like Fold and a Novel Molecular Catalytic Mechanism. ACS Cent. Sci. 2016, 2, 896–903. (55) Biarnés, X.; Ardèvol, A.; Iglesias-Fernández, J.; Planas, A.; Rovira, C. Catalytic Itinerary in 1,3-1,4-β-Glucanase Unraveled by QM/MM Metadynamics. Charge Is Not yet Fully Developed at the Oxocarbenium Ion-like Transition State. J. Am. Chem. Soc. 2011, 133, 20301–20309. (56) Cuxart, I.; Coines, J.; Esquivias, O.; Faijes, M.; Planas, A.; Biarnés, X.; Rovira, C. Enzymatic Hydrolysis of Human Milk Oligosaccharides. The Molecular Mechanism of Bifidobacterium Bifidum LactoN-Biosidase. ACS Catal. 2022, 12, 4737–4743. (57) Morais, M. A. B.; Coines, J.; Domingues, M. N.; Pirolla, R. A. S.; Tonoli, C. C. C.; Santos, C. R.; Correa, J. B. L.; Gozzo, F. C.; Rovira, C.; Murakami, M. T. Two Distinct Catalytic Pathways for GH43 Xylanolytic Enzymes Unveiled by X-Ray and QM/MM Simulations. Nat. Commun. 2021, 12, 367. (58) Lairson, L. L.; Henrissat, B.; Davies, G. J.; Withers, S. G. Glycosyl Transferases: Structures, Functions, and Mechanisms. Annu. Rev. Biochem. 2008, 77, 521–555. (59) Gómez, H.; Mendoza, F.; Lluch, J. M.; Masgrau, L. QM/MM Studies Reveal How Substrate–Substrate and Enzyme–Substrate Interactions Modulate Retaining Glycosyltransferases Catalysis and Mechanism. Adv. Protein Chem. Struct. Biol. 2015, 100, 225–254. (60) Romero-Téllez, S.; Lluch, J. M.; González-Lafont, À.; Masgrau, L. Comparing Hydrolysis and Transglycosylation Reactions Catalyzed by Thermus Thermophilus βGlycosidase. A Combined MD and QM/MM Study. Front. Chem. 2019, 7, 200. (61) Kozmon, S.; Tvaroška, I. Catalytic Mechanism of Glycosyltransferases: Hybrid Quantum Mechanical/Molecular Mechanical Study of the Inverting NAcetylglucosaminyltransferase I. J. Am. Chem. Soc. 2006, 128, 16921–16927. (62) Seto, N. O. L.; Compston, C. A.; Evans, S. V.; Bundle, D. R.; Narang, S. A.; Palcic, M. M. Donor Substrate Specificity of Recombinant Human Blood Group A, B and Hybrid A/B Glycosyltransferases Expressed in Escherichia Coli. Eur. J. Biochem. 1999, 259, 770–775. (63) Ardèvol, A.; Rovira, C. Reaction Mechanisms in Carbohydrate-Active Enzymes: Glycoside Hydrolases and Glycosyltransferases. Insights from Ab Initio Quantum Mechanics/Molecular Mechanics Dynamic Simulations. J. Am. Chem. Soc. 2015, 137, 7528–7547. (64) Gund, P.; Halgren, T. A.; Smith, G. M. Chapter 27 Molecular Modeling as an Aid to Drug Design and Discovery. Annu. Rep. Med. Chem. 1987, 22, 269–279. (65) Mendoza, F.; Masgrau, L. Computational Modeling of Carbohydrate Processing Enzymes Reactions. Curr. Opin. Chem. Biol. 2021, 61, 203–213. (66) Raich, L.; Nin-Hill, A.; Ardèvol, A.; Rovira, C. Enzymatic Cleavage of Glycosidic Bonds: Strategies on How to Set Up and Control a QM/MM Metadynamics Simulation. Methods Enzymol. 2016, 577, 159–183. (67) Tvaroška, I. Atomistic Insight into the Catalytic Mechanism of Glycosyltransferases by Combined Quantum Mechanics/Molecular Mechanics (QM/MM) Methods. Carbohydr. Res. 2015, 403, 38–47. (68) Ardèvol, A.; Rovira, C. The Molecular Mechanism of Enzymatic Glycosyl Transfer with Retention of Configuration: Evidence for a Short-Lived Oxocarbenium-like Species. Angew. Chemie Int. Ed. 2011, 50, 10897–10901. (69) Teze, D.; Coines, J.; Fredslund, F.; Dubey, K. D.; Bidart, G. N.; Adams, P. D.; Dueber, J. E.; Svensson, B.; Rovira, C.; Welner, D. H. O-/ N-/ S-Specificity in Glycosyltransferase Catalysis: From Mechanistic Understanding to Engineering. ACS Catal. 2021, 11, 1810–1815. (70) Gómez, H.; Polyak, I.; Thiel, W.; Lluch, J. M.; Masgrau, L. Retaining Glycosyltransferase Mechanism Studied by QM/MM Methods: Lipopolysaccharyl-α1,4-Galactosyltransferase C Transfers α-Galactose via an Oxocarbenium Ion-like Transition State. J. Am. Chem. Soc. 2012, 134, 4743–4752. (71) Lira-Navarrete, E.; Iglesias-Fernández, J.; Zandberg, W. F.; Compañón, I.; Kong, Y.; Corzana, F.; Pinto, B. M.; Clausen, H.; Peregrina, J. M.; Vocadlo, D. J.; Rovira, C.; Hurtado-Guerrero, R. Substrate-Guided Front-Face Reaction Revealed by Combined Structural Snapshots and Metadynamics for the Polypeptide NAcetylgalactosaminyltransferase 2. Angew. Chem. Int. Ed. Engl. 2014, 53, 8206–8210. (72) Gómez, H.; Rojas, R.; Patel, D.; Tabak, L. A.; Lluch, J. M.; Masgrau, L. A Computational and Experimental Study of O-Glycosylation. Catalysis by Human UDP-GalNAc Polypeptide:GalNAc Transferase-T2. Org. Biomol. Chem. 2014, 12, 2645–2655. (73) Trnka, T.; Kozmon, S.; Tvaroška, I.; Koča, J. Stepwise Catalytic Mechanism via Short-Lived Intermediate Inferred from Combined QM/MM MERP and PES Calculations on Retaining Glycosyltransferase PpGalNAcT2. PLOS Comput. Biol. 2015, 11, e1004061. (74) Albesa-Jové, D.; Mendoza, F.; Rodrigo-Unzueta, A.; Gomollõn-Bel, F.; Cifuente, J. O.; Urresti, S.; Comino, N.; Gõmez, H.; Romero-García, J.; Lluch, J. M.; SanchoVaello, E.; Biarnés, X.; Planas, A.; Merino, P.; Masgrau, L.; Guerin, M. E. A Native Ternary Complex Trapped in a Crystal Reveals the Catalytic Mechanism of a Retaining Glycosyltransferase. Angew. Chemie Int. Ed. 2015, 54, 9898–9902. (75) Yan, L.; Liu, Y. The Retaining Mechanism of Xylose Transfer Catalyzed by Xyloside α-1,3-Xylosyltransferase (XXYLT1): A Quantum Mechanics/Molecular Mechanics Study. J. Chem. Inf. Model. 2020, 60, 1585–1594. (76) Ferreira, P.; Fernandes, P. A.; Ramos, M. J. The Catalytic Mechanism of the Retaining Glycosyltransferase Mannosylglycerate Synthase. Chem. – A Eur. J. 2021, 27, 13998– 14006. (77) Piniello, B.; Lira-Navarrete, E.; Takeuchi, H.; Takeuchi, M.; Haltiwanger, R. S.; Hurtado-Guerrero, R.; Rovira, C. Asparagine Tautomerization in Glycosyltransferase Catalysis. The Molecular Mechanism of Protein O-Fucosyltransferase 1. ACS Catal. 2021, 11, 9926–9932. (78) de las Rivas, M.; Lira-Navarrete, E.; Daniel, E. J. P.; Companõn, I.; Coelho, H.; Diniz, A.; Jiménez-Barbero, J.; Peregrina, J. M.; Clausen, H.; Corzana, F.; Marcelo, F.; Jiménez-Osés, G.; Gerken, T. A.; Hurtado-Guerrero, R. The Interdomain Flexible Linker of the Polypeptide GalNAc Transferases Dictates Their Long-Range Glycosylation Preferences. Nat. Commun. 2017, 8, 1–11. (79) Lira-Navarrete, E.; de las Rivas, M.; Compañón, I.; Pallarés, M. C.; Kong, Y.; IglesiasFernández, J.; Bernardes, G. J. L.; Peregrina, J. M.; Rovira, C.; Bernadó, P.; Bruscolini, P.; Clausen, H.; Lostao, A.; Corzana, F.; Hurtado-Guerrero, R. Dynamic Interplay between Catalytic and Lectin Domains of GalNAc-Transferases Modulates Protein O-Glycosylation. Nat. Commun. 2015 61 2015, 6, 1–12. (80) de las Rivas, M.; Paul Daniel, E. J.; Coelho, H.; Lira-Navarrete, E.; Raich, L.; Compañón, I.; Diniz, A.; Lagartera, L.; Jiménez-Barbero, J.; Clausen, H.; Rovira, C.; Marcelo, F.; Corzana, F.; Gerken, T. A.; Hurtado-Guerrero, R. Structural and Mechanistic Insights into the Catalytic-Domain-Mediated Short-Range Glycosylation Preferences of GalNAc-T4. ACS Cent. Sci. 2018, 4, 1274–1290. (81) Park, J. B.; Kim, Y. H.; Yoo, Y.; Kim, J.; Jun, S. H.; Cho, J. W.; El Qaidi, S.; Walpole, S.; Monaco, S.; García-García, A. A.; Wu, M.; Hays, M. P.; Hurtado-Guerrero, R.; Angulo, J.; Hardwidge, P. R.; Shin, J. S.; Cho, H. S. Structural Basis for Arginine Glycosylation of Host Substrates by Bacterial Effector Proteins. Nat. Commun. 2018 91 2018, 9, 1–15. (82) Tyugashev, T. E.; Vorobjev, Y. N.; Kuznetsova, A. A.; Lukina, M. V.; Kuznetsov, N. A.; Fedorova, O. S. Roles of Active-Site Amino Acid Residues in Specific Recognition of DNA Lesions by Human 8-Oxoguanine-DNA Glycosylase (OGG1). J. Phys. Chem. B 2019, 123, 4878–4887. (83) Nagarajan, H.; Vetrivel, U. Demystifying the PH Dependent Conformational Changes of Human Heparanase Pertaining to Structure–Function Relationships: An in Silico Approach. J. Comput. Aided. Mol. Des. 2018, 32, 821–840. (84) Caliandro, R.; Rossetti, G.; Carloni, P. Local Fluctuations and Conformational Transitions in Proteins. J. Chem. Theory Comput. 2012, 8, 4775–4785. (85) Niu, C.; Zhu, L.; Xu, X.; Li, Q. Rational Design of Thermostability in Bacterial 1,31,4-β-Glucanases through Spatial Compartmentalization of Mutational Hotspots. Appl. Microbiol. Biotechnol. 2017, 101, 1085–1097. (86) Acevedo, J. P.; Reetz, M. T.; Asenjo, J. A.; Parra, L. P. One-Step Combined Focused EpPCR and Saturation Mutagenesis for Thermostability Evolution of a New ColdActive Xylanase. Enzyme Microb. Technol. 2017, 100, 60–70. (87) Chen, L.; Jiang, K.; Zhou, Y.; Zhu, L.; Chen, X. Improving the Thermostability of αGlucosidase from Xanthomonas Campestris through Proline Substitutions Guided by Semi-Rational Design. Biotechnol. Bioprocess Eng. 2022, 27, 631–639. (88) Valero-González, J.; Leonhard-Melief, C.; Lira-Navarrete, E.; Jiménez-Osés, G.; Hernández-Ruiz, C.; Pallarés, M. C.; Yruela, I.; Vasudevan, D.; Lostao, A.; Corzana, F.; Takeuchi, H.; Haltiwanger, R. S.; Hurtado-Guerrero, R. A Proactive Role of Water Molecules in Acceptor Recognition by Protein O-Fucosyltransferase 2. Nat. Chem. Biol. 2016, 12, 240–246. (89) Arreola-Barroso, R. A.; Llopiz, A.; Olvera, L.; Saab-Rincón, G. Modulating Glycoside Hydrolase Activity between Hydrolysis and Transfer Reactions Using an Evolutionary Approach. Molecules 2021, 26, 6586.