Molecular Dynamics Study of Twister Ribozyme: Role of Mg(2+) Ions and the Hydrogen-Bonding Network in the Active Site.

The recently discovered twister ribozyme is thought to utilize general acid-base catalysis in its self-cleavage mechanism, but the roles of nucleobases and metal ions in the mechanism are unclear. Herein, molecular dynamics simulations of the env22 twister ribozyme are performed to elucidate the structural and equilibrium dynamical properties, as well as to examine the role of Mg(2+) ions and possible candidates for the general base and acid in the self-cleavage mechanism. The active site region and the ends of the pseudoknots were found to be less mobile than other regions of the ribozyme, most likely providing structural stability and possibly facilitating catalysis. A purported catalytic Mg(2+) ion and the closest neighboring Mg(2+) ion remained chelated and relatively immobile throughout the microsecond trajectories, although removal of these Mg(2+) ions did not lead to any significant changes in the structure or equilibrium motions of the ribozyme on the microsecond time scale. In addition, a third metal ion, a Na(+) ion remained close to A1(O5'), the leaving group atom, during the majority of the microsecond trajectories, suggesting that it might stabilize the negative charge on A1(O5') during self-cleavage. The locations of these cations and their interactions with key nucleotides in the active site suggest that they may be catalytically relevant. The P1 stem is partially melted at its top and bottom in the crystal structure and further unwinds in the trajectories. The simulations also revealed an interconnected network comprised of hydrogen-bonding and π-stacking interactions that create a relatively rigid network around the self-cleavage site. The nucleotides involved in this network are among the highly conserved nucleotides in twister ribozymes, suggesting that this interaction network may be important to structure and function.

The primordial RNA world relied on the capabilities of RNA as a biocatalyst as well as a genetic coding informational molecule.[1−3] Such RNA enzymes exist currently in the form of ribozymes, which are noncoding RNA molecules that can perform chemical transformations at rates that are sometimes comparable to those of contemporary protein enzymes.[4,5] Ribozymes are pervasive in nature and participate in a variety of biological functions.[4] The small ribozymes, which catalyze self-cleavage at specific sites, encompass eight known natural classes: hammerhead,[6] hairpin,[7] hepatitis delta virus (HDV),[8] Varkud satellite,[9]glmS,[10] twister,[5] pistol,[11,12] and hatchet.[12,13] The site-specific phosphodiester bond cleavage in these ribozymes usually occurs via a mechanism in which the phosphodiester bond to be cleaved is aligned for in-line attack by O2′ of the preceding nucleotide, producing a 2′,3′-cyclic phosphate and a 5′-OH that corresponds to the new 5′-terminus of the substrate chain. Each class exhibits different geometric constraints, as well as varying modes of nucleophile activation, transition state stabilization, and leaving group protonation. All of these factors impact the efficiency of catalysis.[14−16]

Twister ribozyme is one of the newest classes of small ribozymes to be discovered and was identified by means of bioinformatics, followed by biochemical verification.[5] It is found in various eukaryotic and bacterial organisms. The original study introducing twister ribozymes identified nearly 2700 examples adopting a common secondary structure. Its size and structural complexity are comparable to those of riboswitches.[5] The resolved structures of all twister ribozymes exhibit two pseudoknots, T1 and T2 (Figure ), that are distinct in their residues and divalent cation locations at the self-cleavage site.[17−19] Additionally, four nonpseudoknot stems are predicted in the secondary structures[5] along with 10 highly conserved nucleotides. Mutation of these nucleotides, along with compensatory base changes of the stems and pseudoknots supported the relevance of the structural model.[5]

Figure 1

(A) Structure of the env22 twister ribozyme (PDB entry 4RGE). The U–1-A1 self-cleavage site is depicted as sticks with nitrogen and oxygen atoms colored blue and red, respectively, and the crystallographic Mg2+ ions are depicted as pink spheres. The active site Mg2+ ion, Mg101, is depicted as a cyan sphere, and a nearby Mg2+ ion, Mg109, is slightly visible in pink behind the ribozyme upward to the right of Mg101. Pseudoknots T1 and T2 are colored magenta and orange, respectively, and the P1 stem is colored light green. (B) Schematic depiction of the secondary structure of the same ribozyme. The self-cleavage site is marked with a black triangle. The pseudoknots and P1 stem are depicted using the same color scheme as in part A. Pseudoknots T1 and T2 are labeled, as well. The interactions of the residues denoted as P1 are drawn according to PDB entry 4RGE. (C) Schematic depiction of the site-specific phosphodiester bond cleavage in the env22 twister ribozyme at the phosphate linking nucleotides U–1 and A1.

As in many other ribozymes,[20−28] general acid–base catalysis is thought to play a significant role in the self-cleavage mechanism of twister ribozyme. Experiments on twister ribozyme indicate rate constants that strongly depend on the Mg2+ ion concentration and the pH of the medium, plateauing at ∼1 mM Mg2+ and at a pH of ∼6.5.[5] The observed rate constant, kobs, of the twister ribozyme construct used in the experiments was extrapolated from slow reacting conditions to ∼1000 min–1 under physiological conditions of 1 mM MgCl2 and pH 7.4. Four papers have reported five different crystal structures of twister ribozymes.[17−19,29] Sequences used in these crystal structures come from either Oryza sativa (two different sequences) or uncultivated organisms from the environment (two different sequences). Common to all of these structures is a highly conserved G that has been implicated as the general base for activating the nucleophilic attack by U–1(O2′). All of these structures also exhibit interactions of the nonbridging oxygens of the scissile phosphate by either hydrogen bonding to base amino groups or chelation to a metal ion. A table highlighting salient catalytic features of the available twister ribozyme crystal structures can be found in the Supporting Information (Table S1).

In this study, we perform simulations of the env22 twister ribozyme, which was generated by annealing two synthetic RNA strands consisting of 19 and 37 nucleotides, denoted the substrate and ribozyme strands, respectively.[19] Three crystal structures of environmental twister ribozymes have been published with the resolution improving from 3.88 Å of env [Protein Data Bank (PDB) entry 4QJH][17] to 2.89 Å of env22 (PDB entry 4RGE)[19] and then most recently to 2.64 Å of env22 (PDB entry 5DUN).[29] Our study focuses mainly on the 4RGE structure with a minor investigation of the more recently determined 5DUN structure, which shares the same sequence with 4RGE; the principal difference in composition between the two constructs is that 5DUN has a methoxy at U–1 to inhibit self-cleavage while 4RGE has a deoxy at U–1. (Note that 4QJH is a different sequence.) The 4RGE crystal structure, as well as the self-cleavage reaction, is depicted in Figure . The cleavage site features a Mg2+ ion directly coordinated to the nonbridging pro-SP scissile phosphate oxygen at a distance of 2.2 Å. To be consistent with the conventional self-cleavage site numbering, we have renumbered all of the nucleotides in the crystal structure so that the cleavage site U5-A6 becomes a U–1-A1 cleavage site, analogous to the conventional terminology used for other ribozymes. This revised numbering will be used throughout this paper (Table S1). Moreover, nucleotide G43, which has its N1 within hydrogen-bonding distance of the pro-RP nonbridging oxygen of the scissile phosphate, was shown to strongly contribute to catalysis through mutation experiments.[19] The Mg2+ ion directly coordinated to the pro-SP nonbridging oxygen atom was thought to stabilize the developing negative charge on the scissile phosphate, suggesting the possibility of its active participation in catalysis.[19]

The more recently determined 5DUN crystal structure has a slightly higher resolution and contains all of the catalytic residues in 4RGE, as discussed below. Similar to the 4RGE structure, the 5DUN crystal structure also exhibited a Mg2+ ion directly coordinated to the pro-SP nonbridging oxygen atom of the scissile phosphate at a distance of 2.1 Å.[29] The pro-RP nonbridging oxygen was found to be within hydrogen-bonding distance of N1 of the conserved guanine close to the self-cleavage site,[29] as in the 4RGE crystal structure.[19] Nuclear magnetic resonance (NMR) spectroscopy was used to measure a pKa value of 5.1 for a 13C2-labeled adenosine at the +1 position;[29] N3 of this adenosine was also discovered to be crucial for self-cleavage.[29,30] The pKa of 5.1 ± 0.1 is shifted toward neutrality in comparison to the pKa of 3.7 ± 0.1 for the same adenosine in the substrate alone, which was interpreted as being possibly relevant to the general acid–base mechanism.[29] Additionally, an acceleration of the scissile phosphate Sp phosphorothioate cleavage, as well as a partial rescue for the RP phosphorothioate cleavage, was observed in the presence of Mn2+ and Cd2+, suggesting the direct involvement of divalent metals in catalysis. Hence, an RNA catalysis mechanism that employs both metal ion- and nucleobase-assisted acid–base catalysis was proposed.[29] In addition, a synthetic construct lacking the P1 stem,[29] which was proposed to be essential,[5] has been shown to be fully functional, implying that the P1 stem is not required for function.

In the context of these experimental findings, we conducted an extensive molecular dynamics (MD) study based on crystal structures of the env22 twister ribozyme[19,29] to elucidate structural and equilibrium dynamical aspects of this new ribozyme class. Moreover, the importance of the Mg2+ ions for this ribozyme was addressed in this study. For this purpose, we propagated three independent MD trajectories of the 4RGEenv22 twister ribozyme crystal structure, including all ten crystallographic Mg2+ ions. To investigate the role of Mg2+ ions, we also propagated trajectories in which one specific Mg2+ ion, namely either the purported catalytic Mg101 or the nearby Mg109, was removed, or all Mg2+ ions were removed. In all of these replacements, the overall charge conservation of the system was ensured by placing additional Na+ ions in the bulk solvent. In addition, to assess potential differences in the equilibrium dynamics and possibly gain insight into the general acid–base catalysis mechanism, we deprotonated O2′ of U–1 that attacks the scissile phosphate and propagated another trajectory. Finally, a shorter MD trajectory was propagated for a system based on the more recent 5DUN(29)env22 twister ribozyme crystal structure for comparison. The 5DUN structure was found to produce results similar to those obtained from the 4RGE structure. All of these MD trajectories were analyzed to investigate the role of the Mg2+ ions in twister ribozyme and to examine possible candidates for the general base and acid in the self-cleavage mechanism. During completion of this work, a computational study of the O. sativa twister ribozyme, which has a sequence different from that of the env22 twister ribozyme studied herein, was published;[31] a brief comparison of these two studies is presented in the Conclusions.

Methods

A set of MD trajectories was propagated for the env22 twister ribozyme, with the total sampling time amounting to 10 μs (Table ). The initial coordinates were obtained from two different crystal structures of the env22 twister ribozyme (PDB entries 4RGE and 5DUN).[19,29] The sequences of 56 nucleotides associated with these two structures are identical except for two missing nucleotides in the 5DUN structure that are far from the active site. As mentioned above, for the sake of being consistent with the conventional self-cleavage site numbering of nucleotides 1 to −1, we renumbered the 4RGE and 5DUN nucleotide sequences by subtracting 5 from the PDB residue numbers if the original residue number was >5 or by subtracting 6 from the PDB residue numbers if the original residue number was ≤5.

Table 1

Summary of the Employed Systems and the Amount of MD Sampling on Each System

trajectory name	basis PDB structurea	no. of Mg2+ ions	length of the production trajectory (ns)
full run 1	4RGE	10	2000
full run 2	4RGE	10	2000
full run 3	4RGE	10	2000
Na+ at Mg101	4RGE	9	1000
Na+ at Mg109	4RGE	9	1000
no Mg2+	4RGE	0	1000
U–1(O2′–)	4RGE	10	500
5DUN	5DUN	8	500

The 4RGE and 5DUN crystal structures can be found in refs (19) and (29), respectively.

The simulations utilized the ff12SB force field[32−37] implemented in the AMBER14 suite of programs[38] using the following MD protocol. First, we replaced the dU5 residue, which prevented self-cleavage during crystallization of the 4RGE structure, with U and added hydrogens to the crystal structure using the tleap functionality of AMBER14. The 10 Mg2+ ions and the water molecules resolved in this crystal structure were retained. We neutralized the system with the addition of Na+ ions and immersed this complex in an orthorhombic TIP3P water[39] box, ensuring that the sides are at least 10 Å from the solute and also adding 0.15 M NaCl buffer. We minimized the energy of this system using the steepest descent and conjugate gradient methods. This minimization was followed by gradually increasing the temperature to 300 K in a canonical ensemble (NVT) over 200 ps while maintaining weak harmonic restraints on the solute and then equilibrating it in the isothermal–isobaric ensemble (NPT) for 10 ns with no restraints. Analogous equilibration procedures were utilized for the slightly perturbed systems described below, as well as for the 5DUN crystal structure except here a methoxy U–1 was replaced with U. The MD trajectories were propagated using the cuda implementation of the pmemd program.[40]

We propagated eight independent trajectories for six systems that differed in the starting crystal structure, the ions present, or the protonation state of the nucleophile. Starting from an equilibrated system with initial coordinates obtained from the 4RGE crystal structure of the env22 twister ribozyme, we propagated three independent trajectories for 2 μs each. These three trajectories are denoted as “full run 1”, “full run 2”, and “full run 3” for the analysis (Table ). In addition to these trajectories, we propagated a 1 μs MD trajectory, in which the putative catalytic Mg2+ ion, Mg101, was replaced with a Na+ ion. This system is denoted as “Na+ at Mg101” for the analysis (Table ). Mg101 is directly coordinated to the pro-SP phosphate oxygen of the U–1-A1 site and pro-RP phosphate oxygen of the A1-A2 site. (For the sake of simplicity, we refer to Mg101 as the catalytic Mg2+, although its catalytic role has not been fully verified.) In another independent trajectory, Mg109 of the crystal structure was replaced with a Na+ ion because of crystallographic uncertainty about whether it was a monovalent or divalent cation.[19] Mg109 is coordinated to the pro-SP phosphate oxygen of the A1-A2 site and the pro-RP oxygen of the A2-U3 site and may influence catalysis. For our analysis, this system is denoted “Na+ at Mg109” (Table ). In the crystal structure, Mg109 is 6.2 Å from Mg101. Furthermore, another 1 μs MD trajectory was propagated after all crystallographic Mg2+ ions had been deleted and replaced with additional Na+ ions added to the solution to neutralize the system in the absence of the Mg2+ ions. This system is denoted as “no Mg2+” in the analysis (Table ). In addition to these trajectories, a separate 500 ns trajectory was propagated on a system after deprotonation of U–1(O2′) to investigate potentially more catalytically relevant precleavage conformations. This trajectory is designated as “U–1(O2′–)” for the analysis (Table ).

Finally, a 500 ns trajectory was propagated on a system with initial coordinates obtained from the more recent 5DUN crystal structure of the env22 twister ribozyme. This crystal structure has a slightly higher resolution of 2.64 Å,[29] compared to the resolution of 2.89 Å for the 4RGE structure.[19] In addition, the average collective B factor for the U–1 and A1 nucleotides in this crystal structure (39.1 Å2) is significantly lower than that for the 4RGE structure (94.5 Å2). Nucleotides C15 and A51, which are at the chain ends and are far from the active site, were omitted from the construct during experiments for the more recent 5DUN structure. Coordinates for these two residues were obtained by overlaying the neighboring residues in the 5DUN structure with those in the 4RGE structure in VMD[41] and extracting their coordinates from the superimposed 4RGE structure. These coordinates were appended to the 5DUN structure, followed by optimization of the newly added bonds, using Maestro.[42] The goal of the trajectory based on the 5DUN crystal structure was to obtain a qualitative comparison of the two crystal structures in terms of structural and equilibrium dynamical properties. A superimposition of their active sites is given in Figure S1, and the root-mean-square deviation (RMSD) between the heavy atoms of the two crystal structures is 1.0 Å for the entire ribozyme and 2.2 Å for only the active site. The greater RMSD for the active site is due to a slight difference between the conformations of the U–1 base; this difference between the U–1 base conformations decreases after MD equilibration of both structures (Figure S1), leading to an RMSD of just 0.8 Å for the active site.

The post-trajectory analysis for all of these trajectories included the calculation of RMSDs, root-mean-square fluctuations (RMSFs), key distances, radii of gyration, and cross-correlation matrices, as well as the investigation of interactions around the self-cleavage site. In addition, the positional stability of the catalytic Mg101 was assessed through radial distribution functions for Mg2+ around the self-cleavage site. Finally, the radial distribution functions of Mg2+ and Na+ ions around the A1 and A2 phosphate oxygens were examined to shed light on the probable identity of the cation designated as Mg109 in the 4RGE crystal structure and to characterize a third metal ion, a Na+ ion, that persists in this region.

Results and Discussion

Analysis of Structure and Equilibrium Motions

To compare the extent of structural variation within the nucleotide chains in the env22 twister ribozyme, we analyzed the RMSDs of the backbone P atoms (Figure ). The ranges of the RMSD fluctuations for the trajectories with all crystallographic Mg2+ ions, “full runs 1–3”, are slightly larger than those for the trajectories in which Na+ replaces either Mg101 or Mg109 (∼2 Å vs ∼1.5 Å), although the difference is fairly small and probably not indicative of significant differences in structural variation. The RMSD fluctuations in the system with no Mg2+ ions are initially moderate (∼1.5 Å) but then increase slightly at later times. Even in this case, however, they do not exceed what is observed for the “full runs 1–3” trajectories with all crystallographic Mg2+ ions, although the “no Mg2+” trajectory was not propagated as long as the “full runs 1–3” trajectories. Thus, the RMSD analysis indicates minimal differences in structural variation due to the Mg2+ ions.

Figure 2

Time evolution of the RMSDs for the independent MD trajectories. The black, red, and green lines represent the data from three independent trajectories containing all crystallographic Mg2+ ions. The blue line represents the data from a trajectory with Mg109 replaced by Na+, the purple line the data from a trajectory with Mg101 replaced by Na+, and the cyan line the data from a trajectory with no Mg2+ ions.

To investigate the mobilities of the individual nucleotides or groups of nucleotides, we calculated the RMSFs of the backbone P atoms (Figure A). Not surprisingly, we found the chain ends to be the most mobile regions across all systems. Notably, the P1 stem was observed to unwind in all of the MD trajectories by losing the hydrogen-bonding interactions of its U–4-A50 and U–3-A49 central base pairs after ∼150 ns, allowing the ends to move freely in solution but with the active site remaining intact. This melting of the U–4-A50 and U–3-A49 base pairs can be understood because the flanking U–5-A51 and U–2·G48 “base/wobble pairs”, as conventionally drawn in earlier twister secondary structures, are not present in the starting crystal structures: the terminal UA “base pair” is already melted in the crystal structure to allow U–5 to form a cis major groove triple with A45 in T1, and the upper UG “wobble pair” is melted to allow U–2 to form a trans major groove triple with A44 and for G48 to cross-strand stack on C26 of the T1 pseudoknot. Interestingly, melting of this UG wobble allows U–2 to extend up toward the active site, where it stacks between the putative general base G43 and U–5, and allows U–1 to extend into the active site to nucleophilically attack the scissile bond. Recent kinetic studies of a synthetic env22 twister ribozyme, which is the same ribozyme studied herein, in which the potential P1 stem-forming nucleotides were deleted, reported that the P1 stem of the twister ribozyme made essentially no contributions to self-cleavage,[29] consistent with the unwinding of P1 that we observe herein in silico. Apparently, P1 is not needed to form the active site, nor are the U–5- and U–2-containing triples. Note that studies on a twister ribozyme from Clostridium bolteae suggest that P1 is important to cleavage,[5] indicating that more experimental studies are needed.

Figure 3

(A) RMSFs of the P atoms for the 19-mer substrate and the 37-mer ribozyme strands of the env22 twister ribozyme. The break at residue 15 arises because it is the 5′-terminal nucleotide of the 37-mer. The black, red, and green lines represent the data from three independent trajectories containing all crystallographic Mg2+ ions. The blue line represents the data from a trajectory with Mg109 replaced by Na+, the purple line the data from a trajectory with Mg101 replaced by Na+, and the cyan line the data from a trajectory with no Mg2+ ions. (B) Depiction of the regions with low RMSFs. The cyan nucleotides are U–1 and A1, which constitute the self-cleavage site. The green region contains C9 and A10, while the blue portions feature nucleotides 25–47. Note that active site residues −1, 1, 25, 27, 29, 30, and 40–43 are associated with a lower RMSF.

The U–1-A1 self-cleavage site appears to be the least mobile region in the 19-mer substrate strand throughout all trajectories (see Figure ), and the mobility is also relatively low for nucleotides C9 and A10, which stack on pseudoknot T2 that connects G7 and C8 to C32 and G31, respectively. Nucleotides 26–47 are less mobile than the rest of the ribozyme strand and encompass the end of pseudoknot T2 mentioned above and both ends of the pseudoknot T1 connecting nucleotides 26–28 and 45–47. All of the relatively less mobile regions are located in the central core of the structure in all of the trajectories, as depicted in Figure B. In all of our trajectories, active site nucleotides −1, 1, 25, 27, 29, 30, and 40–43[19] exhibit a collective average RMSF of 0.9 Å that is lower than the overall average RMSF of 1.2 Å (Figure and Tables S2 and S3). This observation suggests that the active site of the ribozyme is organized to perform the reaction.

Possible correlated or anticorrelated motions in the env22 twister ribozyme were examined by calculation of cross-correlation matrices. The correlated and anticorrelated regions of this ribozyme under different Mg2+ conditions are illustrated in Figure S2 and are enumerated in the caption of this figure. The overall trends were consistent across all trajectories, with only a few minor variations in the degree of correlation or anticorrelation. Overall, the consistency of the RMSFs and cross-correlation matrices across the independent trajectories with and without Mg2+ ions indicates that the Mg2+ ions do not significantly impact the equilibrium motions of the overall ribozyme. However, as discussed below, specific Mg2+ ions could play a catalytic role and could be important for retaining the precise preorganization of the active site required to facilitate catalysis through electrostatic interactions.

The number of hydrogen bonds forming within the ribozyme was determined as a function of time and analyzed within a histogram, as depicted in Figure . The differences across the four trajectories are negligible, with typically 89–94 hydrogen-bonding interactions being observed. The number of hydrogen bonds fluctuates around 92 for all trajectories. This similarity indicates that the absence of Mg101 or the other Mg2+ ions does not cause drastic alterations in the overall fold of this ribozyme on the microsecond time scale of the trajectories. In addition, to investigate the compactness of the ribozyme, we calculated the radius of gyration over the course of each trajectory. The radius of gyration remained relatively constant at around 17.5 Å for all trajectories, as did the longest interatomic distance through the ribozyme of around 30.9 Å, also indicating that Mg2+ is not required to retain the overall fold on the microsecond time scale of the trajectories (Figure S3). For comparison, the radius of gyration is 17.4 Å and the longest interatomic distance is 29.5 Å for the 4RGE crystal structure. As mentioned above, Mg2+ ions could still be significant for effective catalysis.

Figure 4

(A) Time evolution and (B) histograms for the hydrogen bonds formed between the nucleotides. The black, red, and green lines represent the data from three independent trajectories containing all 10 crystallographic Mg2+ ions. The blue line represents the data from a trajectory with Mg109 replaced by Na+, the purple line the data from a trajectory with Mg101 replaced by Na+, and the cyan line the data from a trajectory with no Mg2+ ions.

Another MD trajectory of 500 ns was propagated with the initial coordinates obtained from the more recent 5DUN crystal structure of env22 twister ribozyme.[29] This trajectory exhibited structural and equilibrium dynamical properties similar to those of all of our other trajectories, which were based on the 4RGE crystal structure. The analyses associated with this trajectory are overlaid on those based on the 4RGE structure in Figures S3–S5.

Analysis of Metal Ions

We analyzed the putative catalytic Mg2+ ion (Mg101 in the crystal structure) through certain distance parameters for the systems in which it was present. This Mg2+ ion is coordinated to the nonbridging pro-SP oxygen of the scissile phosphate of the U–1-A1 site. This interaction distance was measured to be 2.2 Å in the 4RGE crystal structure and 2.1 Å in the 5DUN crystal structure. In the MD trajectories based on the 4RGE crystal structure, we observed an average value of ∼2.0 Å for this distance, ranging from 1.8 to 2.2 Å, in agreement with the value from the crystal structures. This Mg2+ ion is also coordinated to the nonbridging pro-RP oxygen of the A1-A2 site, with a distance of 2.2 Å in the 4RGE crystal structure and 2.1 Å in the 5DUN crystal structure. Our simulations based on the 4RGE crystal structure yielded an average distance of ∼2.0 Å for this distance, with a range of 1.8–2.2 Å, also in agreement with the crystal structures. The time evolution of these distances is shown in Figure S6. The interaction distances between the catalytic Mg2+ ion and these two oxygen atoms do not exhibit significant variation across the different trajectories that contain Mg101. These interaction distances are also similar for the trajectory based on the 5DUN crystal structure (data not shown). In addition to these two nonbridging oxygen atoms, Mg101 is coordinated to U3(O4) and three water molecules (Figure ). This dual phosphate coordination is reminiscent of the coordination of the metal ion in the active site of the HDV ribozyme.[43]

Figure 5

Depiction of the active site of the env22 twister ribozyme for a representative configuration from the “full run 1” trajectory. G43 (in ball and stick representation) is hypothesized to participate in the self-cleavage mechanism by stabilizing the negatively charged transition state through the hydrogen bonds it forms with the pro-RP oxygen of the scissile phosphate, shown as black dashed lines. Mg101 is observed to coordinate to the pro-SP oxygen of the scissile phosphate and thus is a candidate to actively contribute to the self-cleavage mechanism. It also coordinates to the pro-RP oxygen of A2 and O4 of U3, as well as three water molecules (cyan dashed lines). A1 displays a syn conformation at all times.

We further characterized the interactions of Mg101 with the ribozyme and compared these interactions to those of a Na+ ion placed at position 101 in the absence of Mg101. For this purpose, we calculated radial distribution functions and running coordination numbers for Mg2+ or Na+ relative to the pro-SP (OP1) oxygen of A1 and the pro-RP (OP2) oxygen of A2 (Figure ). The data from the “full runs 1–3” trajectories illustrate that Mg101 exhibits a strong chelated interaction with these two oxygen atoms and is relatively immobile. For comparison, two other Mg2+ ions distal to the active site were found to diffuse into the solvent over the microsecond trajectories (see below). The low mobility of the catalytic Mg2+ ion is confirmed by an average RMSF value of 0.6 Å for Mg101 over all trajectories based on the 4RGE structure. In the absence of Mg101, a Na+ ion also exhibits a strong chelated interaction with these two oxygen atoms and a similar degree of immobility (Figure ). The strong interaction of Na+ in the absence of Mg2+ at the Mg101 site suggests that this region of the active site generates a highly negatively charged pocket that attracts positively charged ions. Given its proximity to the cleavage site, this negatively charged pocket, as well as the associated positively charged ion, could potentially serve a catalytic function for self-cleavage.

Figure 6

(A) Radial distribution functions for Mg2+ and Na+ ions relative to the pro-SP (OP1) oxygen of A1 and the pro-RP (OP2) oxygen of A2 in the region of Mg101. The radial distribution functions for Mg2+ (black and red curves, “full run 1”) form pronounced peaks at ∼2.0 Å, with a second peak indicating that another Mg2+ is located at a distance of ∼4.2 Å for A2(OP2). The results for “full run 1” are representative of those for “full runs 1–3”. The radial distribution functions for Na+ (blue and cyan curves, “Na+ at Mg101”) form smaller peaks at ∼2.1 Å because of the higher density of Na+ in the simulation cell, given that the ion density occurs in the denominator of the standard definition of the radial distribution function. See the Supporting Information for the standard definitions of the radial distribution function and running coordination number. (B) Running coordination numbers for Mg2+ or Na+ relative to the pro-SP oxygen (OP1) of A1 and the pro-RP oxygen (OP2) of A2 in the region of Mg101. A shared Mg2+ or a shared Na+ is coordinated to both oxygen atoms at ∼2.0 Å for the “full run 1” or “Na+ at Mg101” trajectory, respectively. The coordination number for Na+ to A1(OP1) is slightly greater than unity at ∼2.0 Å, indicating that another Na+ samples this spherical shell in some configurations. The second step for the “full run 1” trajectory indicates that another Mg2+ is located at a distance of ∼4.2 Å from A2(OP2). Note that there is no practicing consensus in the RNA literature about the assignment of OP1. Thus, we used the terminology in ref (19), which is to assign OP1 as pro-SP.

For the “full runs 1–3” trajectories, we also analyzed the distance between the catalytic Mg2+ ion (Mg101) and the closest other Mg2+ ion, which was crystallographically labeled Mg109 but was reported to also potentially be a monovalent ion (Figure S7).[19] According to the crystal structure, the distance between these two ions is 6.2 Å.[19] In our trajectories that contain both Mg2+ ions, Mg109 remains at an average distance of 5.6 Å from Mg101, fluctuating between 5.2 and 6.1 Å (Figure S7). Thus, both ions exhibit limited mobility, and the distance between them is long enough that ion–ion repulsion is not dominant. Furthermore, our analysis of the overall structure and motions identified no significant differences between the trajectories with a Mg2+ (“full runs 1–3”) or a Na+ (“Na+ at Mg109”) for the tentatively assigned Mg109. Thus, the identity of the cation at position 109 does not appear to significantly impact the structure or equilibrium motions of the env22 twister ribozyme, even though it is the closest crystal structure cation to the catalytic Mg2+. An overlay of a configuration from the “full run 1” trajectory with Mg2+ as this cation and a configuration from the “Na+ at Mg109” trajectory is depicted in Figure S8. Note that a Mg2+ ion at position 109 was resolved in the 5DUN crystal structure, supporting the assignment of this ion as Mg2+ rather than Na+.

Additionally, we compared the interactions and mobilities of Na+ and Mg2+ at the site tentatively assigned as Mg109 in the crystal structure. For this purpose, we calculated the radial distribution functions and running coordination numbers of Mg2+ or Na+ ions around the pro-SP (OP1) oxygen of A2 and the pro-RP (OP2) oxygen of U3, where the electron density for a monovalent or divalent cation was observed and assigned as Mg109 in the crystal structure (Figure S9). Clearly, the Mg2+ ion is much more strongly coordinated to the A2(OP1) and U3(OP2) oxygen atoms than is the Na+ ion. Only 0.054% of the saved configurations from the “Na+ at Mg109” trajectory exhibited a Na+ ion close to these oxygen atoms (within 3 Å), and the Na+ ions in this region were continually exchanging with other Na+ ions, indicating a high mobility for the Na+ ions. On the other hand, in all three 2 μs trajectories with a Mg2+ at position 109 (“full runs 1–3”), this Mg2+ ion always remained bound at the crystal structure location. Thus, the Mg2+ ion occupies the site bound to both oxygen atoms, exhibiting a chelated interaction, whereas the Na+ ion does not consistently occupy this site, exhibiting a diffuse interaction, as characterized previously for other ribozymes.[24,44] The lower mobility and higher occupancy of Mg2+ at this site are consistent with the assignment of Mg109 in the crystal structure.

As mentioned above, for all of the trajectories based on the 4RGE structure propagated in this study, we observed one (Mg104) or two (Mg104 and Mg110) of all 10 Mg2+ ions leave the crystal structure location and move into the bulk solvent. Both Mg2+ ions moved away from their original sites in the “full runs 1–3” trajectories, whereas only Mg104 moved away in the “Na+ at Mg109” and “Na+ at Mg101” trajectories. Because Mg110 moves away in the later portion of the “full runs 1–3” trajectories, the lack of movement observed for the “Na+ at Mg109” and “Na+ at Mg101” trajectories may be due to the shorter simulation time for those two trajectories. Mg104 is only weakly coordinated to the ribozyme, with five of the six coordinations to water molecules in the crystal structure, and over the course of the trajectory becomes coordinated to a sixth water molecule, allowing it to move away from the ribozyme. Subsequently, it occasionally moves closer to the ribozyme and then leaves again. Similarly, Mg110 is coordinated to four water molecules in the crystal structure and over the course of the trajectory becomes coordinated to a fifth and then a sixth water molecule, allowing it to move away from the ribozyme. All other crystal structure Mg2+ ions remained in their original locations throughout all trajectories. None of our analyses identified a change elsewhere in the ribozyme at or after the times when Mg104 or Mg110 moved away from the ribozyme, suggesting that these metal ions are not essential for the structural integrity or the equilibrium dynamics of the env22 twister ribozyme.

We observed that the relative positioning of A1 and A30 varies depending on the Mg2+ ion content, which in turn could affect the geometry of the cleavage site (Figure S10). All Mg2+-containing trajectories indicated that Mg107 coordinates to U25(O4), A30(OP2), three water molecules, and either another water in the “Na+ at Mg109” and “Na+ at Mg101” trajectories or U28(OP2) in the trajectories with all crystallographic Mg2+ ions (“full runs 1–3”). After all Mg2+ ions were removed, U25 shifted while maintaining its Watson–Crick base pairing to A30, leading to a π–π stacking conformation for A1 and A30. The conformational differences caused by the absence of the Mg2+ ions that are not in the active site might be important, as demonstrated by the geometrical shift around the active site caused by the removal of Mg107 when all Mg2+ ions were removed (Figure ). Thus, metal ions distal to the active site might also play key roles in facilitating rate acceleration.

Figure 7

Relative positioning of A1, A29, and A30 for a representative configuration from (A) the “full run 1” trajectory, in the presence of Mg2+ ions, illustrating a “T-shaped” stacking interaction between A29 and A1 and (B) the “no Mg2+” trajectory, in the absence of Mg2+ ions, illustrating a π–π stacking interaction between A1 and A30.

It is noteworthy that upon removal of all crystallographic Mg2+ ions, including Mg101, which is thought to be catalytic, the structural and equilibrium dynamical properties of the ribozyme did not exhibit significant differences. For the trajectories in which all Mg2+ ions are removed, however, the distance between the P atoms of A1 and A2 increased from 5.0 to 5.6 Å. This change is a relatively small increase in terms of the overall fold, although the shorter distance in the presence of Mg2+ could stabilize the structure somewhat, which might impact catalysis. As discussed further below, none of our trajectories adopted the alignment for the nucleophilic attack of U–1(O2′–), suggesting that this crystal structure conformation may be precatalytic. Thus, Mg101 may shift somewhat in the self-cleavage process.

Analysis of the Self-Cleavage Site Geometry

Interactions present at the self-cleavage site in the env22 twister ribozyme are depicted in Figure . In all of the MD trajectories discussed prior to this subsection, all bases were in their canonical protonation states. Recent work[31] published during the review of this paper proposes that the nucleobases equivalent to A1 and G43 serve as the general acid and base, respectively, in a related twister ribozyme. To examine these possible acid and base candidates, we propagated additional 200 ns MD trajectories for each of three additional systems with varying protonation states: (1) G43(N1) deprotonated and A1(N3) protonated, (2) G43(N1) in its canonical form and A1(N3) protonated, and (3) G43(N1) deprotonated and A1(N3) in its canonical form. These trajectories did not exhibit any major conformational changes and therefore do not provide definitive evidence regarding catalysis, although such changes could possibly occur on a longer time scale.

Candidates for the general base abstracting the proton from U–1(O2′) to initiate nucleophilic attack on the scissile phosphate and for the general acid donating a proton to A1(O5′) have been proposed for all four twister sequences crystallized to date (Table S1). G43 has been proposed to be the general base for the env22 twister ribozyme. In the 4RGE crystal structure, G43 participates in a sheared cis-Hoogsteen sugar edge pair with A2, which positions G43 near the nucleophilic 2′-OH of U–1 in the crystal structure. The distance between G43(N1) and U–1(C2′) is 5.7 Å in this crystal structure. We analyzed this distance for the systems with all bases in their canonical protonation states. An initial geometry optimization of the ribozyme was performed after substitution of the hydroxyl at U–1(C2′) of the 4RGE crystal structure, as well as the addition of hydrogen atoms, water molecules, and buffer. Following this initial geometry optimization, the distance between G43(N1) and U–1(C2′) remained similar to that of the crystal structure at 5.5 Å, and the G43(N1)–U–1(O2′) distance was 5.1 Å. Following equilibration, however, these distances increased to more than 6.0 Å, and the average G43(N1)–U–1(C2′) and G43(N1)–U–1(O2′) distances over all trajectories based on the 4RGE crystal structure were 6.3 and 7.1 Å, respectively. This increase in distance occurred for all seven independent trajectories based on the 4RGE crystal structure. The G43(N1)–U–1(C2′) distance is 6.6 Å in the 5DUN crystal structure, and the average G43(N1)–U–1(C2′) and G43(N1)–U–1(O2′) distances in the “5DUN” trajectory are 6.5 and 7.3 Å, respectively. Thus, the MD trajectories based on the 4RGE crystal structure evolved toward the 5DUN crystal structure in this region. Note that the hydrogen on G43(N1) remained hydrogen bonded to the scissile phosphate pro-RP nonbridging oxygen (OP2), with a G43(N1)–OP2 distance of 2.6 Å in the 4RGE crystal structure and an average G43(N1)–OP2 distance of 2.9 Å over all trajectories (Figure ). A potential limitation of the current crystal structures of the twister ribozyme is that the U–1(O2′) is either deoxy or methoxy and therefore is not positioned correctly for catalysis. For example, G43(N2) may hydrogen bond to the pro-RP nonbridging oxygen while G43(N1) hydrogen bonds to U–1(O2’) during catalysis.

Figure 8

Evolution of active site nucleotides U–1, A1, and G43 in the MD simulations based on the 4RGE crystal structure. (A) Crystal structure geometry after substitution of the hydroxyl group at U–1(C2′) and a local optimization of its position. (B) Structure following geometry optimization of the entire ribozyme. (C) Configuration obtained after MD for 2 μs for the “full run 1” trajectory. The G43(N1)–U–1(O2′) distance is larger in part C than in parts A and B. A1(O5′) is hydrogen bonded to a Mg2+-bound water molecule in part B but then rotates to point toward A1(N3) in part C; the A1(O5′)–A1(N3) distance is smaller in part C (∼4.4 Å) than in part B (∼5.4 Å). The U–1(O2′)–A1(P)–A1(O5′) angle is more linear in part A than in parts B and C. The 5DUN crystal structure is similar to the configuration in part C.

The larger G43(N1)–U–1(O2′) distances observed in the MD trajectories are not conducive to G43(N1) acting as the general base. Given these large distances, we next examined systems with deprotonated G43(N1), protonated A1(N3), or both. For the two trajectories in which G43(N1) was deprotonated, the average G43(N1)–U–1(O2′) distance was 7.2 Å when A1(N3) was protonated and 8.1 Å when A1(N3) was in its canonical form. Thus, deprotonation of G43(N1) did not result in shorter distances that would be conducive to G43(N1) acting as the general base, although longer trajectories could be necessary to observe such large conformational changes. Other candidates for the base were considered, but none seem satisfactory (Figure ). A detailed discussion about the alternative candidates is provided in the Supporting Information.

Figure 9

Active site region for a representative configuration obtained from the “full run 1” trajectory. Possible general base candidates for the deprotonation of U–1(O2′), which is marked with a red label, and possible general acid candidates for the protonation of A1(O5′), which is also marked with a red label, are annotated. Mg101 is shown as a pink sphere.

The self-cleavage reaction would require in-line attack by U–1(O2′) if the reaction were to proceed with inversion of configuration at the reactive phosphorus, as in other small ribozymes.[45−47] After substitution of the hydroxyl at U–1(C2′) of the 4RGE crystal structure, as well as the addition of hydrogen atoms, water molecules, and buffer, only the U–1(2′OH) atoms were optimized, leading to a U–1(O2′)–A1(P)–A1(O5′) angle of 148° (Figure A). As reported previously[19] from modeling based on the crystal structure, this configuration is favorable for in-line attack by U–1(O2′). However, after subsequent optimization of the geometry of the entire ribozyme, but prior to molecular dynamics, this angle decreased to 104°, which is less conducive to in-line attack (Figure B). The average angle during production for all trajectories starting with the 4RGE crystal structure, including the “U–1(O2′–)” trajectory and the three additional trajectories with different protonation states for G43(N1) and A1(N3), was ∼80° (Figure C). Initial modeling starting from the 5DUN crystal structure resulted in a U–1(O2′)–A1(P)–A1(O5′) angle of 90°, with an average angle of 77° over the “5DUN” trajectory. Thus, the initial geometry optimization of the system, prior to MD, moved toward a configuration that appears to be less catalytically active than the crystal structure. To test the sensitivity of the results on the water force field, we also propagated a trajectory using the TIP4P[48] water model. We found that the positions of the bases in the vicinity of U–1 were qualitatively similar to the observations in the other trajectories, and the P1 stem exhibited a degree of melting, as illustrated in Figure S11.

With regard to the general acid, A1(N3) is a candidate for the general acid if it bears a hydrogen to donate at physiological pH. The syn conformation of this base, which is consistent across all of our trajectories and is present in all of the twister crystal structures to date,[17−19,29] offers this possibility as it places N3 in the general vicinity of the leaving group A1(O5′) (Figures and 9). In addition, this A1 is essential for catalysis, as revealed by the A1G mutant having no catalytic activity,[19] and loss of activity with N3C variants of A1.[30] Gaines and York recently proposed a pKa shift for A1(N3) caused by a constellation of phosphates,[31] similar to that proposed for C75 in the HDV ribozyme,[49] thereby allowing A1(N3) to serve as the general acid. However, measurements of the pKa of this A reveal a pKa of just 5.1 that could not be assigned specifically to the N1 or N3.[29] This pKa value is not optimal for general acid–base catalysis and is inconsistent with the pKa values of 7 from experimental studies, albeit on different constructs.[18,30] As discussed below, other possible candidates for the general acid were suggested by our MD simulations.

Analysis of the MD trajectories with all bases in their canonical protonation states indicated that the position of A1(O5′) evolved over time, suggesting the possibility of two or more thermally accessible configurations. Following the initial geometry optimization of the 4RGE crystal structure, A1(O5′) became hydrogen bonded to one of the water molecules coordinated to Mg101, with a distance of 2.8 Å between the oxygen atoms in the hydrogen bond (Figure B). This hydrogen-bonding distance remained similar throughout the entire equilibration process, including simulated annealing and MD equilibration for 10 ns. Moreover, the distance between A1(O5′) and A1(N3) remained at ∼5.4 Å during equilibration. After 40 ns of the “full run 1” trajectory, however, A1(O5′) rotated to point toward A1(N3), but with an intervening Na+ ion, and remained in this position for the rest of the 2 μs trajectory. The other trajectories also exhibited this conformational change in A1(O5′), and the average A1(O5′)–A1(N3) distance over all trajectories based on the 4RGE structure with the bases in their canonical protonation states was 4.4 Å. In the 5DUN crystal structure, A1(O5′) is already rotated toward A1(N3), but with a distance of 4.8 Å in the crystal structure and an average distance of 4.0 Å for the “5DUN” trajectory. Thus, similar to the movement of U–1(O2′) discussed above, A1(O5′) moves from the 4RGE structure toward the higher-resolution 5DUN structure during the MD trajectories. In this case, however, A1(O5′) appears to move toward the proposed A1(N3) base, albeit with a Na+ ion or water molecule intervening. Similar distances were observed for the trajectories with A1 and G43 in noncanonical protonation states. For the trajectories in which A1(N3) was protonated, this distance was 4.4 Å when G43(N1) was deprotonated and 4.6 Å when G43(N1) was in its canonical state, suggesting that protonation of A1(N3) does not decrease this distance. Thus, the average A1(O5′)–A1(N3) distance is more than ∼4.4 Å even when A1(N3) is protonated, which does not support the possible role of A1(N3) as the general acid, although it cannot be ruled out.

As mentioned above, upon visual inspection of the data from different trajectories, we noted that a significant fraction of the configurations from the trajectories of the systems in the canonical protonation states feature a Na+ ion between A1(N3) and A1(O5′) located directly above its ribose (Figure ). This Na+ ion coordinates to N3, O4′, and O5′ of A1 (Figure ), as well as to U–2(O2′) and two water molecules. This Na+ ion, which is present in approximately 7 μs of our total of 10 μs MD sampling with different Mg2+ ion compositions, could be facilitating the stabilization of the negative charge on A1(O5′) after self-cleavage. The running coordination number of the Na+ ions around A1(O5′) is depicted in Figure A. This scenario suggests that A1(N3) could coordinate a metal ion rather than act as a general acid, which would account for the recent report that N3 of A1 is critical for ribozyme activity.[30]

Figure 10

(A) Running coordination numbers for Na+ relative to A1(O5′). All trajectories except for “full run 3” exhibited a Na+ near this O5′. (B) Representative configuration from the “full run 1” trajectory depicting the Na+ ion as a purple sphere between A1(N3) and A1(O5′) located directly above the ribose of A1. This Na+ ion is coordinated to A1(N3), A1(O4′), A1(O5′), U−2(O2′), and two water molecules, although the coordination to U−2(O2′) is not shown here.

In particular, this Na+ ion is not present for the trajectories in the noncanonical protonation states. For the trajectory with A1(N3) protonated and G43(N1) deprotonated, A1(O5′) is hydrogen bonded to A1(N3) in some configurations (11% of all the saved configurations from a 200 ns trajectory), and both A1(O5′) and A1(N3) are hydrogen bonded to a single bridging water molecule or to two different water molecules in the rest of the configurations (Figure S12). For the trajectory with A1(N3) protonated and G43(N1) in its canonical state, A1(O5′) is hydrogen bonded to U–2(O2′), and A1(N3) is hydrogen bonded to a water molecule. The H-bonding interaction of the protonated A1(N3) with a water molecule suggests the possibility of water serving as the acid in the self-cleavage reaction. One possibility is that A1(O5′) samples both the configuration in which it is hydrogen bonded to a Mg2+-bound water molecule and the configuration in which it is rotated to point toward A1(N3) (Figure B,C). The general acid could be the Mg2+-bound water molecule when A1(O5′) is in the first configuration. The MD simulations cannot rule out any of these possibilities.

Molecular Interaction Network around the Self-Cleavage Site

The catalytic pocket for the env22 twister ribozyme has been defined to contain nucleotides −1, 1, 25, 27, 29, 30, and 40–43 (Figure ).[19] Among these nucleotides, the mutation of G43 was fully deleterious[19] and also significantly altered the pH dependence of the twister ribozyme,[18] underscoring its importance in catalysis, possibly because of its persistent hydrogen-bonding interaction with the pro-RP oxygen of the scissile phosphate (Figure ). Reduced catalytic activity was observed upon mutation of C41 and A42, suggesting that they also contribute to folding, structural stability, or catalysis. Moreover, as proposed previously,[19] we detected a water-mediated hydrogen-bonding network involving C41(N4) and A42(N6), along with U–1(O2′) (Figure ). This hydrogen-bonding network could provide structural stability in this vicinity. We observed a shared hydrogen-bonding water between C41(N4) and U–1(O2′) in 48.0% of all saved configurations from the “full runs 1–3” trajectories, whereas a common hydrogen-bonding water involving A42(N6) and U–1(O2′) was observed in only 16.0% of these configurations. This hydrogen-bonding network could explain the experimental result that C41U and A42G mutations reduce the catalytic activity, with a greater reduction for the C41U mutation.

Figure 11

Set of nucleotides lining the catalytic pocket for a representative configuration from the “full run 1” trajectory. The backbone ribose and phosphate of each nucleotide are omitted for the sake of simplicity. U–1 and A1 are displayed as thicker sticks. The water-mediated hydrogen bonding among C41(N4), A42(N6), and U–1(O2′) is shown with black dotted lines. The catalytic Mg101 is shown as a pink sphere.

G43 has been proposed to stabilize the transition state in the self-cleavage reaction through two hydrogen-bonding interactions with the nonbridging pro-RP oxygen on the scissile phosphate.[19] Our simulations illustrated that the pro-RP scissile phosphate oxygen (OP2) is involved in a set of bifurcated hydrogen bonds with G43 (Figure and Figure S13). One of these hydrogen bonds is formed with a hydrogen on G43(N2), while the other hydrogen bond is formed with the hydrogen on G43(N1). This persistent hydrogen bonding directly to the cleavage site supports the hypothesis that G43 might stabilize the transition state. Moreover, we observed G43 to establish a noncanonical sheared cis-Hoogsteen sugar edge[19] with A2 through two hydrogen bonds (Figure S13). The first hydrogen bond is formed between the other hydrogen on G43(N2) and A2(N7), while the second one is formed between G43(N3) and one of the hydrogens on A2(N6). An additional hydrogen bond is often formed between G43(O2′) and the other hydrogen on A2(N6), as well as between the hydrogen on G43(O2′) and A44(O5′). Moreover, A44 interacts with A2 through π–π stacking. Overall, these interactions combine to create a relatively rigid network around the self-cleavage site involving A1, G43, A2, and A44 (Figure S13 and Table S2).[19] As discussed above, however, G43 is a strong candidate for the general base in the self-cleavage reaction, most likely requiring a conformational change that brings G43(N1) closer to U–1(O2′).

The interaction network involving A1, G43, A2, and A44 also includes additional contributions from A29, C33, G24, U25, and A30 (Figure S14). Further details about these interactions are provided in the Supporting Information. All of these residues are among the highly conserved nucleotides in twister ribozymes. Thus, this highly interconnected interaction network around the self-cleavage site could be significant for the proper function of twister ribozymes. The only highly conserved nucleotide that does not participate in this interaction network around the self-cleavage site is G31. It forms a base triplet within pseudoknot T2, likely providing structural stabilization through the hydrogen-bonding interactions associated with this base triplet (Figure S15 and the Supporting Information for further details).

Testing the Effects of Nucleophile Deprotonation

Our analysis indicates that additional conformational changes would be required prior to self-cleavage. To determine whether deprotonating U–1(O2′) would induce such conformational changes, we propagated an additional MD trajectory for 500 ns, starting from the 4RGE structure with U–1(O2′) deprotonated, corresponding to the situation after the general base has removed this proton. To account for the additional charge of −1 caused by deprotonation, we added an extra Na+ ion to the bulk solvent to retain overall charge neutrality. Our analysis indicated that this trajectory exhibits behavior very similar to that of the trajectories in which U–1(O2′) is protonated (see “full runs 1–3”). The RMSD, RMSF, and radius of gyration analysis results are provided in Figures S3–S5. Analogous to the discussion above, the U–1(O2′)–A1(P)–A1(O5′) angle was 147° when only the 2′O– atom was optimized after the initial substitution in the crystal structure, but during the MD trajectories, this angle decreased to 104°, which is less conducive to in-line attack of the scissile phosphate. Most likely, additional conformational changes occur on a longer time scale but are inaccessible through microsecond MD trajectories.

Conclusions

This paper presents microsecond MD trajectories on the env22 twister ribozyme. The results illustrate that the overall structure and equilibrium motions of this ribozyme are not significantly influenced by the presence of Mg2+ ions on the microsecond time scale, although one crystallographic Mg2+ ion, namely Mg107, was found to influence the structure in the active site region. Moreover, the overall structure and equilibrium motions are similar for trajectories based on two related but different crystal structures of the env22 twister ribozyme[19,29] and for a trajectory that is propagated after deprotonation of the attacking O2′. In all cases, the active site region and the ends of the pseudoknots are less mobile than other regions of the ribozyme, most likely providing structural stability and possibly facilitating catalysis. The P1 stem is already partially melted at both its top and bottom in the crystal structures[19,29] and unwinds further in the MD trajectories, consistent with experiments indicating that this stem is not essential to function.[29] In addition, specific regions of the ribozyme exhibit consistent correlated or anticorrelated motions with other regions of the ribozyme in all trajectories. Furthermore, the simulations revealed a persistent interconnected network comprised of hydrogen-bonding and π-stacking interactions that create a relatively rigid network around the self-cleavage site. The nucleotides involved in this network, A1, A2, G24, U25, A29, A30, C33, G43, and A44, are among the highly conserved nucleotides in twister ribozymes, suggesting that this interaction network might be essential for the function of these ribozymes. The other highly conserved nucleotide, G31, forms a stacked base triplet with A38 and C8, leading to structural stabilization of the T2 pseudoknot and the P2 stem. These analyses suggest that this ribozyme could utilize a preorganized, relatively immobile active site combined with thermal equilibrium fluctuations of the overall structure to facilitate the self-cleavage reaction.

The interactions of specific Mg2+ ions were investigated to elucidate possible structural and catalytic roles.[19] The Mg2+ ion labeled Mg101, which is coordinated to the nonbridging pro-SP oxygen of the U–1-A1 site, the nonbridging pro-RP oxygen of the A1-A2 site, U3(O4), and three water molecules, has been proposed to play a catalytic role. Specifically, results of metal ion rescue experiments illustrating acceleration of SP phosphorothioate cleavage, as well as to a lesser degree RP phosphorothioate cleavage, in the presence of Mn2+ and Cd2+ support the catalytic importance of this Mg2+ ion.[29] The assignment of a nearby Mg2+ ion (Mg109), which is coordinated to the nonbridging pro-SP oxygen of A1-A2, the nonbridging pro-RP oxygen of A2-U3, and four water molecules, was uncertain in this crystal structure. Both of these Mg2+ ions, 101 and 109, remained chelated and relatively immobile during all trajectories in which they were present. Moreover, trajectories in which either of these Mg2+ ions was removed did not exhibit any significant changes in the structure or equilibrium motions of the ribozyme on the microsecond time scale. In contrast, the Mg2+ ion labeled Mg107, which coordinates to U25(O4), A30(OP2), three water molecules, and U28(OP2) or another water molecule, was found to significantly impact the active site structure. The absence of Mg107 results in structural changes that could inhibit the in-line attack of U–1(O2′) on the scissile phosphate bond, suggesting that this Mg2+ ion is both structurally and catalytically relevant.

These analyses of the Mg2+ ions, as well as the identification of a persistent Na+ ion located between A1(N3) and A1(O5′), provide catalytically relevant insights. Metal ion rescue experiments suggest an important interaction between Mg2+ and the nonbridging oxygens of the scissile phosphate. The identification of Mg101 coordinated to the pro-SP oxygen and the same other atoms in two different env22 twister ribozyme crystal structures, as well as the observation that it stays strongly chelated throughout a total of 10 μs of MD trajectories, is consistent with a catalytic role of this Mg2+. The Mg2+-bound water molecule could potentially play the role of the general acid. In addition, the presence of a Na+ ion close to A1(O5′) suggests that it could possibly stabilize the negative charge on A1(O5′) during self-cleavage and coordinate to key atoms in that region. The possibility that protonated A1(N3) could participate in the reaction[18,30] as a general acid is inconsistent with a shift of the pKa to only 5.1 in NMR experiments,[29] although different constructs were used in these two studies. Moreover, the average distance between A1(N3) and A1(O5′) is too large for such a protonation reaction in MD trajectories based on the 4RGE and 5DUN structures, particularly if a Na+ ion or water molecule intervenes, although the possibility cannot be ruled out. In terms of possible candidates for the general base, experimental studies strongly suggest that G43(N1) is the general base that deprotonates U–1(O2′),[18] and the 4RGE crystal structure suggests that G43(N1) is not too far from U–1(O2′), although the distance between G43(N1) and U–1(O2′) is larger in the MD trajectories and the 5DUN structure.

As this work was being completed, a computational study of the O. sativa twister ribozyme was published.[31] This ribozyme has a sequence different from that of the env22 ribozyme, and the crystal structure used as the basis for this study (PDB entry 4OJI) exhibited qualitative differences from the 4RGE and 5DUN structures studied herein. One of the major differences is that the 4OJI structure does not include a catalytic Mg2+ ion in the active site. As discussed above, the catalytic Mg101 remains strongly chelated during 10 μs of MD and could play structural, electrostatic, and possibly catalytic roles. In addition, the 4OJI structure is not in a conformation with U–1(O2′) aligned for in-line attack of the scissile phosphate, in contrast to the 4RGE structure. Another significant difference is that the P1 stem is partially melted in the 4RGE and 5DUN crystal structures and is found to unwind in the MD trajectories, whereas it is fully base paired in the 4OJI structure. Moreover, the melting of the UG wobble at the top of the P1 stem in the crystal structure allows conformational changes that facilitate the alignment of U–1 for nucleophilic attack on the scissile bond. Finally, the MD trajectories of the O. sativa ribozyme sampled configurations in which U–1 and G33 (G43 in env22) were stacked, whereas our trajectories of the env22 ribozyme did not sample this type of configuration over 10 μs, but rather G43 stacked on U–2. Additional experimental and computational studies will be necessary to further untangle the catalytic mechanism of the twister ribozyme.