Article Information

• PDF (677 kb)
• Full Text with Large Figures
• Supporting Material
• Cited by in Scopus (29)

PubMed

• Articles by Maryna V. Krasovska
• Articles by Jiří Šponer

Related Articles

• A three-dimensional model of hepatitis delta virus ribozyme ...

  A three-dimensional model of hepatitis delta virus ribozyme based on biochemical and mutational analyses Current Biology, Volume 4, Issue 6, 1 June 1994, Pages 488-498 N.Kyle Tanner, Sophie Schaff, Gilbert Thill, Elisabeth Petit-Koskas, Anne-Marie Crain-Denoyelle and Eric Westhof Summary Our model supports the pseudoknotted structure and rationalizes several observations relating to the lengths of the various stems and the sequence requirements of the single-stranded regions. It also provides insight into the catalytic mechanism of the HDV ribozyme. We specifically propose that residues C75, U20 and C21 form the basis of the catalytic region and are close to the cleavable phosphate. Summary | Full Text | PDF (1699 kb)

• Structural Roles of Monovalent Cations in the HDV Ribozyme

  Structural Roles of Monovalent Cations in the HDV Ribozyme Structure, Volume 15, Issue 3, 1 March 2007, Pages 281-287 Ailong Ke, Fang Ding, Joseph D. Batchelor and Jennifer A. Doudna Summary The hepatitis delta virus (HDV) ribozyme catalyzes viral RNA self-cleavage through general acid-base chemistry in which an active-site cytidine and at least one metal ion are involved. Monovalent metal ions support slow catalysis and were proposed to substitute for structural, but not catalytic, divalent metal ions in the RNA. To investigate the role of monovalent cations in ribozyme structure and function, we determined the crystal structure of the precursor HDV ribozyme in the presence of thallium ions (Tl). Two Tl ions can occupy a previously observed divalent metal ion hexahydrate-binding site located near the scissile phosphate, but are easily competed away by cobalt hexammine, a magnesium hexahydrate mimic and potent reaction inhibitor. Intriguingly, a third Tl ion forms direct inner-sphere contacts with the ribose 2′-OH nucleophile and the -S scissile phosphate oxygen. We discuss possible structural and catalytic implications of monovalent cation binding for the HDV ribozyme mechanism. Summary | Full Text | PDF (1076 kb)

• Symposium 9: RNA in Action

  Symposium 9: RNA in Action Biophysical Journal, Volume 94, Issue , 1 February 2008, Pages 310-311 Full Text | PDF (72 kb)

• …more

Copyright © 2006 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 91, Issue 2, 626-638, 15 July 2006

doi:10.1529/biophysj.105.079368

Nucleic Acids

Cations and Hydration in Catalytic RNA: Molecular Dynamics of the Hepatitis Delta Virus Ribozyme

Maryna V. Krasovska^*, Jana Sefcikova^†, Kamila Réblová^*, ‡, Bohdan Schneider^§, Nils G. Walter^†, ,   and Jiří Šponer^*, ,  

^* Institute of Biophysics, Academy of Sciences of the Czech Republic, 61265 Brno, Czech Republic
^† Department of Chemistry, Single Molecule Analysis Group, The University of Michigan, Ann Arbor, Michigan 48109-1055
^‡ National Centre for Biomolecular Research, Faculty of Science, Masaryk University, 61137 Brno, Czech Republic
^§ Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, 166 10 Prague, Czech Republic

Address reprint requests to Jiří Šponer or Nils G. Walter.

The hepatitis delta virus (HDV) ribozyme is an RNA enzyme, found in the RNA genome of a human pathogen, the hepatitis delta virus. The ribozyme catalyzes a site-specific transesterification reaction generating 5′-hydroxyl and 2′,3′-cyclic phosphate termini, and plays an essential role in the formation of antigenomic and genomic strands during the viral replication of multimeric intermediates 1. The HDV ribozyme was the first RNA enzyme for which direct involvement of its own nucleobase C75 in catalysis was suggested 2,3.

Several crystal structures of the ribozyme are currently available, including wild-type and C75U mutant precursor and product forms 4,5,6. The global structure of the HDV ribozyme is characterized as a double-nested pseudoknot, stabilized by stacking of five helical stems and an extensive network of hydrogen bonds among helices, loops, and joiners (Fig. 1) 4,5,6. Cytosine C75 is located close to the cleavage site. The most significant conformational difference in the otherwise very similar tertiary structures of the precursor and product forms (root mean-square deviation (RMSD)=2Å) has been identified as a collapse of stem P1 and loop L3 toward the center of the ribozyme 6, resulting in a deeper placement of cytosine C75 into the catalytic pocket, close to L3 in the product form. The flexible L3 is bridged by a noncanonical G25/U20 basepair that shows an anti-to-syn flip of G25 between the precursor and product structures.

Display large version of this figure

Figure 1

Secondary structure of the simulated genomic HDV ribozyme with structural elements color-coded. The product form lacks A-2 and U-1. Open arrow, cleavage site; broken lines, quadruple and A-minor motifs of A77 and A78.

The crystallographic data suggest that C75 is well poised to act as the general base to deprotonate the nucleophilic 2′ OH at the cleavage site 6. The mechanistic data are contradictory, with the latest study suggesting that C75 may act as the general acid to protonate the leaving group, which requires C75 to be N3-protonated just before cleavage 7. However, since a chemical modification was introduced at the scissile phosphate group in the latter study to accelerate leaving group departure, the reaction pathway may be affected. Alternatively, it cannot be ruled out that the transition state structure of the catalytic pocket differs from the ground state structure observed in x-ray diffraction studies, which may lead to the contradictory picture obtained by structural and mechanistic studies. Finally, the HDV ribozyme may utilize distinct C75-based catalytic strategies depending on the exact RNA construct and metal ion conditions 8,9.

Metal ions critically contribute to RNA folding and function 10. Active tertiary conformation and catalytic function are remarkably sensitive to the concentration and type of cation(s) present 11,12. HDV ribozyme activity has long been known as exclusively dependent on the presence of divalent cations 12,13,14. The majority of available crystal structures of the HDV ribozyme in its precursor form 6 reveals a consistent picture of cation binding with two typical binding sites, formed by L3 and the major groove of P4. In contrast, there are nine RNA-bound Mg²⁺ cations suggested by the crystal structure of the wild-type product form (Figure 2A) 15. For comparison, the number of folding-specific Mg²⁺ cations indicated by solution experiments for various medium-sized RNAs is only in the range of 0–4 16,17. As demonstrated recently, assessment of cations in x-ray structures is complex. Some observed cations are related to the crystallization process (packing, etc.) and their occurrence in the crystal is highly sensitive to experimental details 18,19. Occasionally, water molecules or even sulfate anions were suggested to be mislabeled as magnesium dications 19. Variability of Mg²⁺ binding patterns in otherwise similar structures could also reflect competition of divalent and monovalent cations for the same binding sites. Specifically, bound divalent cation clearly seen in one x-ray structure may be substituted by fluctuating divalent or monovalent cations that would evade detection by x-ray crystallography 20.

Display large version of this figure

Figure 2

(A) ESP contour maps of precursor (left) and product (right) ribozyme crystal structures, contoured at −25kT/e. Nucleotides are color-coded as in Fig. 1; crystal positions of Mg²⁺ are shown as green balls. Position of Mg²⁺ closest to the active site in product crystal is marked by an asterisk. (B) Na⁺ binding maps contoured at 10 σ in simulations PreC41+A-2 (left) and ProC41+A-2 (right). Na⁺ binding sites located at the major NESP sites have the same numbering as their associated NESP sites. Sites NA7 and NA11 (Table S2) are not seen at this contour level due to low occupancy.

A critically located hydrated magnesium ion is thought to be an important chemical participant in the reaction mechanism 2,3,9,21. The change in divalent metal-ion preference upon alteration of the linkage at the scissile bond provided the first biochemical evidence for coordination of a divalent metal ion in the active-site region of the cis-(self-)cleaving genomic ribozyme 22. Two equivalent reaction mechanisms have been proposed wherein a hydrated magnesium ion could either donate a proton to the 5′-oxygen leaving group acting as the general acid or deprotonate the 2′ hydroxyl group acting as the general base, playing a complementary role to C75 so that both general acid and base catalysis are utilized in a concerted fashion 2,3,7. Consistent with such models, anticooperative interactions between a protonated C75 and a magnesium cation have been demonstrated by pH and magnesium titration studies 3. Although cleavage activity of the genomic HDV ribozyme can also be observed in the absence of magnesium ions, it is greatly reduced even at molar concentrations of sodium cations in comparison to millimolar magnesium concentrations 2,3,8,23. Nakano et al. have detected structural and catalytic ions with 125-fold and 25-fold contributions, respectively, to cleavage rate enhancement 8. Additional studies have tentatively suggested that the structural site shows an inner sphere interaction with a preference for Mg²⁺ and that the catalytic site has outer sphere binding with little preference for a particular divalent ion 9.

In addition to cations, hydrating water is thought to be an integral part of nucleic acid structure 24,25,26,27,28,29,30,31,32,33,34. Despite limitations imposed by force field approximations and limited simulation timescales (sampling), molecular dynamics (MD) techniques are well suited to analyze hydration, especially regarding predictions of highly specific long-residency hydration sites 20,27,28,29,30,31,32,33,35. MD simulations are capable of providing qualitative insights into cation binding to nucleic acids, including detection of complex cation binding pockets that contain delocalized cations 20,29,32,36.

Recently, we analyzed ∼120ns of explicit solvent MD simulations of the HDV ribozyme 37. Precursor simulations with unprotonated C75 revealed weak dynamic binding of C75 in the catalytic pocket with spontaneous transient formation of a hydrogen (H) bond between U-1(O2′) and C75(N3). This Hbond would be required for C75 to act as the general base. Protonated C75H⁺ moved deeper toward loop L3 (resembling its product location) and established a firm Hbonding network. However, a C75H⁺(N3)-G1(O5′) Hbond, which would be expected if C75 acted as a general acid catalyst without substantial structural rearrangement, was not observed. The simulations confirmed that loop L3 is dynamical and may serve as a flexible structural element, possibly gated by the closing U20/G25 basepair, to facilitate a conformational switch induced by a protonated C75H⁺. In this study, we i), substantially extend our simulations to more than 200ns, and ii), provide a detailed analysis of hydration, cation binding, electrostatic potential, and backbone dynamics of the HDV ribozyme.

Starting geometries were based on crystal structures of the C75U mutant precursor HDV ribozyme (Protein Data Bank (PDB) codes 1SJ3 and 1VC7, resolution 2.30Å and 2.45Å) and the wild-type product (PDB code 1CX0, 2.30Å) (Table 1) 4,6,15. The remaining crystal structures were inspected to reveal the most typical crystal cation and water binding sites, if available. The structure of the catalytically inactive C75U mutant precursor was modified using InsightII 38 to carry C75 of the wild-type or a protonated C75H⁺. Note that the lower-resolution crystal structure of wild-type precursor in the absence of divalents (PDB 1VC5, 3.4Å resolution) has essentially the same structure as the C75U mutant (RMSD of 0.30Å) 6.

We tested different variants for initial positioning of Mg²⁺ ions (Figure 2A, Table 1). Most simulations were carried out with two Mg²⁺ ions initially placed as in the precursor x-ray structure. Inclusion of nine Mg²⁺ cations (as suggested by the wild-type product crystal) results in a Mg²⁺ concentration of 0.1M. Such simulations provide an upper limit for possible Mg²⁺-related structural effects seen in nanosecond-scale simulations.

Divalent cations are poorly described by pair additive potentials and also sample insufficiently in simulations 29,39. Since Na⁺ are better parameterized and demonstrate rather satisfactory sampling 20,29,40, we have also performed MD simulations in the presence of just sodium ions to identify and compare cation binding sites between the precursor and product forms of the ribozyme. Na⁺-only simulations are justified since the C75U mutant ribozyme that was prepared and crystallized in the absence of divalent metal ions (2.7Å, PDB 1VBX) has the same structure as it does in the presence of divalents (RMSD=0.31Å). Our unpublished terbium(III) footprinting data on the genomic cis-acting ribozyme also support the notion that divalent ions are not required for HDV ribozyme folding. Finally, the same HDV ribozyme shows residual cleavage activity in molar concentration of Na⁺ in the absence of divalents 3,8,41. Anyway, the simulations are too short to reveal an unfolding caused by the lack of divalent cations 34. In the Supplementary Material, we provide justification for using minimal neutralizing Na⁺ concentrations.

All MD simulations (Table 1) were carried out using the AMBER 7.0 program package 42 with the parm99 Cornell et al. force field 43,44,45. The RNA was solvated in a rectangular box of TIP3P waters 46 extended to a distance of ≥10Å from any solute atom. The simulated system was neutralized by a minimal number of sodium cations 47 initially placed by the LeaP module at points of favorable electrostatic potential close to the RNA. This corresponds to an ion concentration of ∼0.2M. In several simulations, sodium cations were initially moved away from the solute (after the initial electrostatic placement) to prevent trapped cations. The ions then spontaneously locate to the binding sites starting from bulk. The Sander module of AMBER 7.0 was used for the equilibration and production simulations using standard protocols (see, e.g., Reblova et al. 29 and Razga et al. The particle mesh Ewald method 48 was applied with a heuristic pair list update, using a 2.0-Å nonbonded pair list buffer and a 9.0Å cutoff. The particle mesh Ewald charge grid dimensions are products of powers of 2, 3, and 5, resulting in a grid spacing of ∼1.0Å. The direct sum tolerance of 10⁻⁵ and the nonbonded cutoff of 9.0Å lead to an Ewald coefficient of 0.30Å⁻¹ (see Supplementary Material for further details).The production runs were carried out at 300K with constant-pressure boundary conditions using the Berendsen temperature coupling algorithm 49 with a time constant of 1.0ps.

The trajectories were analyzed using the Ptraj module of the AMBER 7.0 package and our own scripts and visualized by the programs PyMOL 50 and VMD 51. Long residency cation-binding and hydration sites were identified by means of calculation of cation and water density maps by a Fourier-averaging method 52. Individual solvent particles were traced up to the distance cutoff of 3.4Å for water molecules and 2.5Å for Na⁺ from ribozyme electronegative atoms. The positions of cations were taken in regular time intervals and then Fourier-transformed into pseudoelectron densities. The solvent density contour maps were visualized using the program Xfit 53. The electrostatic potentials of crystal structures and a series of simulated averaged structures were calculated using the program Delphi 54 by solving the nonlinear Poisson-Boltzmann equation for ionic strength 0.2M, and were visualized using InsightII 38. Each atom was placed in a medium with a dielectric constant of 2.0 in the solvent inaccessible surface-enclosed volume, which was obtained using a probe radius of 1.4Å, whereas solvent was treated as a continuum with a dielectric constant of 80. Our backbone analysis was based on systematic monitoring of the backbone torsion angles followed by comparison with known backbone conformational families 55.

Locations of the deepest negative electrostatic surface potentials (located within the −25kT/e contour and further referred to as NESP sites) are similar in crystal structures of the precursor and product forms of the ribozyme (Figure 2A) and do not change significantly in our MD simulations. The widest NESP site 1 (Figure 2A) is found at the pocket formed by the compact fold of loop L3, which encompasses the cleavage site and catalytic residue C75. The strong bend of the L3 backbone at the C21-U23 segment results in clustering of phosphates. NESP site 1 is associated with the global ESP minimum in the precursor (∼−63kT/e), whereas it is reduced to ∼−45kT/e in the product, becoming a local minimum. The NESP weakening is a consequence of the departure of the scissile phosphate (G1) and rearrangement of the active site after the cleavage reaction. At the same time, the collapse of P1 and L3 toward the center of the product and the shift of C75 deeper into the active site result in a shallower catalytic pocket. The presence of the −1 phosphate strengthens the NESP at the active site by only 2kT/e (test calculations not shown). NESP site 5 (at the major groove of the P1/P1.1 helices with peak near G1/U37 pair in the immediate vicinity of the active site) is also significantly amplified in precursor.

Two divalent cation-binding sites (marked as MG1 and MG2 in this article) observed in precursor crystal structures are located within NESP sites 1 and 2. On the contrary, only one of the nine Mg²⁺ cations seen in crystal structure of wild-type product locates within NESP site 5. Further, in this product structure, the Mg²⁺ closest to the catalytic pocket is 7Å away from the NESP site 1 (Figure 2A). The calculated ESP is thus consistent with cation binding in precursor crystal structures, but inconsistent with the distribution of Mg²⁺ cations in the wild-type product crystal. Binding of cations in simulations is substantially determined by ESP (see below).

Mg²⁺ possesses an octahedral coordination sphere 56,57. Hexahydrated Mg²⁺ interacts with RNA via nonspecific electrostatic interactions, whereas interactions including direct contact (inner shell) between RNA and Mg²⁺ require partial dehydration of the cation 58,59. The crystal structures, however, do not reveal positions of water molecules. Six functional groups (four phosphoryl oxygen and two uracil keto oxygen atoms) are located within the outer sphere coordination distance from the MG1 metal ion in crystals of C75U precursor 6. Further, a short contact between U75(O4) and Mg²⁺ suggests likely inner sphere binding.

In the precursor simulation PreC41+U75A-2Mg (cf. Table 1), the whole coordination sphere was in basic agreement with the crystal structure, including the inner shell Mg…U75(O4) contact (1.98Å). However, additional Mg…G1(O2P) (1.87Å) inner shell contact was formed (Fig. 3) in contrast to outer shell (water-mediated) coordination suggested by the crystal structure, accompanied by shift of the cation by ∼2.0Å from the initial position. The octahedral ligand shell is completed by four water molecules, bridging the cation with phosphates of U-1, C22, and U23. The first hydration shell waters do not exchange with the bulk solvent in agreement with the experimentally measured microsecond residency time 60 and computational studies 20,29,61.

Display large version of this figure

Figure 3

First hydration shell of Mg²⁺ bound at the catalytic pocket in simulations PreC41+U75A-2Mg (left) and PreC41+A-2Mg before the G1(O5′)…Mg²⁺ water bridge was broken (middle). (Right) First and second hydration shells of the catalytic Mg²⁺ in the simulation PreC41+A-2Mg (after the G1(O5′)…Mg²⁺ water bridge was lost). The residues are color-coded as shown in Fig. 1. Mg²⁺ is represented by a green ball; dark blue and cyan mesh represent first and second hydration shells of Mg²⁺, respectively. Functional groups of the catalytic pocket directly bound to Mg²⁺ are black balls, those bound to the first hydration shell waters are gray balls, and those bound to the second hydration shell waters are white balls. See also Fig. S1.

A metal ion is also identified in the active site of the low-resolution x-ray structure of the wild-type (C75) precursor ribozyme 6, ∼3.0Å away from the Mg²⁺ position in the C75U mutant precursor. This difference in the metal ion location has been attributed to the loss of a favorable contact between the keto oxygen of U75 and the metal ion 6. In all wild-type precursor simulations, C75 is shifted by 2Å within the active site compared to U75. This shift creates enough space for Mg²⁺ binding in the catalytic pocket while avoiding a direct contact with the N4-amino group of C75 (Fig. 3). Indeed, both wild-type simulations PreC41+Mg and PreC41+A-2Mg reveal stable binding of Mg²⁺ in the pocket. The average displacement of the Mg²⁺ from its initial position is 1.8Å and 1.3Å in the two simulations, respectively, and Mg²⁺ forms a direct 1.85Å contact with U23(O1P). The first hydration shell now consists of five water molecules. Four of these waters are equatorial and the fifth is positioned axially with respect to the U23(O1P) atom while being bound to U20(O2). Only the U20 nucleotide interacts with the first Mg²⁺ hydration shell through a nucleobase atom (O2), consistent with an experimentally suggested role of the U20 base in Mg²⁺ binding by the active site of the HDV ribozyme 62. The equatorial water molecules are involved in Hbonding with phosphates of G1, C22, and U23, and only one water molecule does not have contact with the ribozyme (Fig. 3). It should be noted that competitive inhibition of the wild-type HDV ribozyme by cobalt hexammine has been used to suggest that the catalytic metal ion binds into a region of low charge density and/or by outer-sphere contacts, in contrast to our simulation results 9.

The G1(O5′)…H-O…Mg bridge, critical for ribozyme cleavage in mechanistic proposals involving general acid catalysis by a hydrated Mg²⁺, was formed only at the beginning of the C75 wild-type simulations (Fig. 3, Fig. S1). This water bridge disappeared in simulations due to rearrangement of U-1/G1 backbone, described in Figure 2A of Krasovska 37. After this rearrangement, G1(O5′) becomes occasionally hydrated by water molecules from the Mg²⁺ second hydration shell. Currently, we cannot decide whether this local backbone rearrangement is correct or reflects a force field imbalance. The MD backbone geometry basically corresponds to the RNA backbone family 29 55, whereas the x-ray geometry does not match any established RNA backbone geometry class.

Twelve water molecules reside in the second shell of a bulk Mg²⁺ ion with average distances of 4.25Å between water oxygen and the magnesium 56,57. The second hydration shell of the Mg²⁺ bound at the catalytic pocket includes on average only nine water molecules, whereas the outer coordination sphere is completed by contacts with the ribozyme, filling the whole catalytic pocket with a dense network of Hbonds (Fig. 3). Binding of Mg²⁺ in the catalytic pocket results in a significant increase in the residence time of water molecules in the second hydration shell (by almost two orders of magnitude to up to >10ns, Table 2). The average water residency time in the second shell of the bulk Mg²⁺ is ∼15ps with the present force field, in line with literature data using a specialized force field 56, see also Auffinger et al. 61.

In summary, our MD simulations are in basic agreement with the crystal structure of C75U precursor. The simulations reveal smooth accommodation of the divalent ion in the catalytic pocket of the wild-type (C75) precursor, and a dense network of long-residency water molecules bridging the magnesium and the RNA.

Protonation of C75 would be required for it to play a role as the general acid during catalysis. Protonated C75H⁺ moves toward its product-like location in all precursor simulations and establishes a stable Hbonding network 37. When placing the Mg²⁺ into the active site, the cation is expelled within 1ns, consistent with the known competition between C75 protonation and Mg²⁺ binding at the active site 3. No magnesium cation is observed at the cleavage site in wild-type and C75U product crystal structures. In a product simulation with an Mg²⁺ cation initially placed at the pocket, the cation is expelled during the equilibration.

The very swift expulsions of the ions indicate that the ion binding is substantially destabilized. Further details can be found in the Supplementary Material .

The dynamic loop L3 regulates the negative electrostatic potential of the catalytic pocket 37. In addition, the residues U-1, G1, U20, C22, and U23 contribute specific coordination sites for catalytic cation binding. The U20/G25 basepair forms a rigidifying bridge across L3 and supports a compact fold of the otherwise flexible loop. The bifurcated cis-W.C./W.C. (Watson-Crick) configuration 63 of the U20/G25 basepair observed in the precursor crystal structures is preserved in all simulations with the catalytic Mg²⁺ bound in the pocket. It has two substates, one of them additionally mediated by Na⁺ cation (Figure 4A) for ∼60% of the basepair lifetime. In the absence of Mg²⁺, the U20/G25 basepair is immediately rearranged into a trans-W.C./H. (Hoogsteen) or cis-W.C./W.C. basepair 37 (Figure 4B). The shape of L3 is significantly affected by the type of U20/G25 basepair (Fig. S2).

Display large version of this figure

Figure 4

Interactions of the U20/G25 basepair with cations in precursor simulations. Basepair edge oriented inside the catalytic pocket is directed downward. Gray and black balls represent Mg²⁺ and Na⁺ cations, respectively. (A) Two substates seen in precursor simulations with Mg²⁺. The U20(O4)-G25(N7) distance trajectory in simulation PreC41+A-2 is shown below. (B) Rearrangement of the U20/G25 basepair in the absence of Mg²⁺; arrows show transitions observed in simulations. Substates 1–3 were observed in the following simulations: 1), PreC41+ and PreC41+A-2; 2), PreC41+C75+, PreC41+U75, PreC41+A-2 and PreC41+U75A-2; and 3), PreC41+Trc and PreC41+U75.

The U20/G25 basepair is different in the product crystal structure due to the anti-to-syn flip of G25. The anti-conformation of G25 observed in the precursor leads to a larger size of the catalytic pocket, required to accommodate the Mg²⁺ and U-1/A-2 residues. An actual anti-to-syn transition of G25 between the precursor and product configurations upon cleavage would require disruption of the U20/G25 basepair and at least a partial unfolding of L3. Indeed, disruption of the basepair accompanied by unfolding of L3 was observed in simulation PreC41+Trc, initiated using the precursor crystal structure but deleting U-1 (Table 1) as well as in the last 2ns of simulation PreC41+U75, both in the absence of Mg²⁺. Similarly, simulations of product show formation of two different types of U20/G25 basepairs and the effect of its occasional disruption on L3 unfolding 37. However, G25 always remains in the syn conformation in product simulations and anti in precursor simulations. The timescale of the simulations is unlikely to allow a spontaneous flipping, even in the simulation PreC41+Trc, attempting to mimic the precursor to product switch.

We suggest that the U20/G25 configuration observed in the precursor crystal structure requires binding of Mg²⁺. Furthermore, after cleavage and expulsion of Mg²⁺ from the cleavage site, L3 becomes more dynamic and its unfolding could allow for the anti-to-syn transition of G25 on a longer timescale.

Cleavage rates above the background level have been measured for the genomic HDV ribozyme in molar concentrations of Na⁺ ions only 3,8,41. In all simulations lacking the Mg²⁺ ion at the active site, an exceptionally strong Na⁺ binding is seen at the catalytic pocket, which accommodates up to two Na⁺ cations simultaneously. A higher occupancy and slower exchange of cations between the pocket and bulk solvent were clearly observed in precursor simulations (Table 3).

The occupancy of the pocket by Na⁺ is sensitive to arrangement of the electronegative atoms pointing into the pocket. Thus, unfolding of L3 in simulation ProC41+C75+ results in reduced occupancy (Table 3). Furthermore, occupancy of monovalent ions is affected by orientation of the U20/G25 basepair. For example, a U20/G25 trans-W.C./H. basepair shows long residency (6–13ns) inner-shell Na⁺ binding in the pocket, whereas a cis-W.C./W.C. basepair seen in simulation PreC41+ is only involved in outer shell cation binding (Figure 4B). This fact explains a relatively low occupancy of the pocket in simulation PreC41+. Protonation of C75 is another factor that disfavors simultaneous binding of two Na⁺ ions (Table 3).

Although in most simulations there was essentially no exchange of ions between catalytic pocket and bulk solvent (Table 3), this observation does not reflect incidental trapping of ions at the beginning of the simulations 64. The outcome of the simulations does not depend on the initial location of the ions, as some simulations were initiated with all ions shifted away from the solute (Table 1). A detailed analysis of the ion binding in one representative precursor simulation is given in Table 4.

The first hydration shell of Na⁺ shows residency times comparable to the second hydration shell of Mg²⁺ (Table 2). Although maximal binding times of water molecules in the ligand shell of bulk Na⁺ cations do not exceed 200ps, they significantly increase when the waters bridge the ribozyme backbone with bound Na⁺ ions (Fig. 5).

Display large version of this figure

Figure 5

Two Na⁺ cations (cyan balls) and their hydration shell (blue mesh) bound at the catalytic pocket in simulation PreC41+A-2. Positions of cations in the figure are determined by points with highest pseudo-electron density of Na⁺. The residues are color-coded as in Fig. 1; functional groups of the catalytic pocket directly bound to the Na⁺ are shown as black balls and those interacting with the first hydration shell waters as gray balls.

In contrast to divalents, two Na⁺ freely migrate along the wide catalytic pocket, making multiple direct and water-mediated long residency contacts with the ribozyme (Table 4). The hydration shells of two Na⁺ cations fill the whole pocket; however, no long residency hydration of G1(O5′) is observed. It thus does not seem that a specific Na⁺-stabilized water molecule is poised to facilitate protonation of G1(O5′), at least within the approximations of the present force field.

In summary, there is a strong Mg²⁺ or Na⁺ binding at the active site in the wild-type precursor, whereas only Na⁺ cations show stable binding at the active site in the product and C75H⁺ precursor simulations. Despite 100% occupancy of the precursor pocket by 1–2 Na⁺ cations, the binding is variable and the cations are delocalized (Table 4), which would make them difficult to detect in crystal structures.

Divalent cation binding to J1.1/4 and P4 was suggested by Pb²⁺ induced scission at the MG2 site 65 and was also revealed in crystal structures of precursor and C75U mutant product forms. Location of this cation close to C41/C41H⁺ indicates that it is a plausible candidate for the pH-sensitive structural cation associated with C41 8. Mg²⁺ is located within NESP site 2 in the crystals (Figure 2A) and makes no inner-sphere contact to the ribozyme. In our simulations, the cation remains in the pocket and typically interacts via inner shell binding with A42(O2P) and five water molecules, which do not exchange with the second shell. The axial water has a fixed position and makes no contact to the ribozyme, whereas the equatorial water molecules easily rotate, alternately bridging the cation with C41(O2′) and A43(O2P) (Fig. S3). The second hydration shell of this Mg²⁺ is considerably more dynamical than the second hydration shell of the Mg²⁺ at the active site pocket (Fig. S3). The binding times of water molecules in the second shell are up to 5.5ns. All water molecules with binding time above 1ns mediate the C41H⁺/A43 basepair. Although this basepair is stable and water mediated in Na⁺-only simulations, the water binding times are much shorter in the absence of bound Mg²⁺ (see below).

In the absence of Mg²⁺, a prominent Na⁺ binding site is located at the NESP site 2 (Figure 2B). Three to five Na⁺ ions exchange in the MG2 site in our simulations, leading to a total occupancy close to 100%. The cations bind alternately via inner and outer shell binding to C41H⁺(O2′), A42(O2P), A43(O2P), and G71(N7) (Table 5). The MG2 site may act as a structural site and together with the MG1 site could be the site of the strong competition between magnesium and sodium, affecting HDV ribozyme catalytic activity 9.

There are several other remarkable Na⁺ cation-binding sites in the HDV ribozyme (Table 5). Cation binding regions are similar in all precursor and product simulations and are primarily determined by ESP (Fig. 2, Fig. S4). The single-stranded junction J4/2 bearing the catalytic C75 is associated with two Na⁺ binding sites to phosphate. A rather compact inner-shell binding site NA3 bridges the trefoil turn with the C19/U20 backbone (Table 5, Fig. 6). It fits NESP site 3 and is in agreement with metal ion-induced cleavage experiments that have observed scission sites at positions G74, C75, and G76 induced by Mg²⁺66, Pb²⁺65, and Tb³⁺67,68 cations.

Display large version of this figure

Figure 6

Na⁺ binding sites associated with the J1.1/4 junction. (A) Interphosphate distances C32(P)-A78(P) (top) and C75(P)-C21(P) (bottom) in product and precursor simulations. (B) Na⁺ (gray balls) binding in sites NA3 and NA4 represented by cation density maxima in simulation PreC41+.

The cation binding site NA4 is located at the arching junction of P1 and P3 and fits NESP site 4. This inner-shell phosphate Na⁺ binding site is further extended by cation binding in the major groove of consecutive A77 and A78 (Table 5, Fig. 6). Occupancies of sites NA3 and NA4 are affected by the C32/A78 and C75/C21 interphosphate distances, which in turn are sensitive to arrangement of joiner J4/2. These sites have higher Na⁺ occupancy (Table 5, Table S1) in product simulations with deeply buried and firmly bound C75 and J4/2 placed closer to P1 and L3 (Fig. 6).

Cation binding in double-helical segments is described in the Supplementary Material .

We have carried out several simulations with nine Mg²⁺ cations, which are described in more detail in the Supplementary Material . Assessing all the data including Na⁺ simulations that provide sufficient sampling, we suggest that our simulations are inconsistent with the overall positioning of nine Mg²⁺ cations in the wild-type product crystal structure. In contrast, dynamics of cations seen in simulations is entirely consistent with the two Mg²⁺ binding sites in the precursor and C75U product crystal structures.

The available crystal structures of the HDV ribozyme provide essentially no information about hydration. The simulations reveal a variety of hydration sites, including common phosphate hydration sites 25 and specific hydration sites related to noncanonical basepairs and complex elements of tertiary structure (Fig. 7).

Display large version of this figure

Figure 7

(A) Long-residency water bridges in noncanonical basepairs. (*) Binding times of water molecules are significantly longer when they belong to the second hydration shell of Mg²⁺ (in parentheses, see text). (**) Two cis-W.C./W.C. basepairs show different water binding times, which are significantly longer in the U20/G25 basepair exposed into the catalytic pocket (in parentheses). (B) Closed (left) and open (right) substates of A-minor type II. Water density at contour level >10 σ (gray mesh); typical orientation of the water molecules, maximal (t_max) and averaged (t_ave) water binding times, are shown.

G1/U37, trans-W.C./H., or cis-W.C./W.C. U20/G25, G40/G74, and C41/A43 (Figure 7A) represent known types of water-mediated basepairs commonly observed in crystal structures 63,69,70. Water-mediated basepairs show water binding times 0.5–5.5ns 27,29,32. Water binding times of a given basepair, however, can significantly vary depending on its environment. Thus, vicinity of hydrated mono- and especially divalent cation results in a significant increase of water binding times as observed for C41/A43 and cis-W.C./W.C. U20/G25 basepairs (Figure 7A). Interestingly, the trans-W.C./H. U20/G25 basepair was observed in our simulations as either water- or Na⁺-mediated, with the cation or water molecule bridging U20(O2) and G25(O6) (Figure 4B). The cation mediated form is more stable 37.

The A-minor motif is the most important recurring RNA tertiary interaction formed by adenines interacting with the minor groove edges of G=C basepairs 71. Dynamical water insertion into A-minor motifs was recently reported for ribosomal kink-turns 34. A-minor motif in the HDV ribozyme involves the consecutive stacked A78 and A77, interacting with two W.C. G=C pairs in helix P3 (Fig. 1). A78 optimally fits into the minor groove of the G29=C18 basepair as an A-minor type I motif 71. This interaction is entirely stable in the simulations. A77 stacks on C75 and forms an A-minor type II interaction. Its position is therefore related to the positioning of C75 in the active site. As a result, this A-minor type II interaction occurs in two substates in simulations (Figure 7B). The x-ray-like geometry is stable in simulations PreC41+U75, PreC41+U75A-2, and ProC41+C75+, whereas in all other simulations it is alternating with an open A-minor interaction mediated by a long-residency (up to 3ns) C19(O2)/C19(O2′)/A77(N3) water bridge. Thus, insertion of a water molecule stabilizes the open conformation of the A-minor motif and is likely related to subtle regulation of base 75 positioning in the active site. A77 is the only base that makes vertical (stacking) contact with the 75 base. The closed A-minor interaction observed in simulation PreC41+U75 leads to maximal mutual overlap of the aromatic rings of U75 and A77 (Supplementary Fig. S5a ). The closed A-minor interaction was also observed in simulation PreC41+C75+, but in this case A77 and protonated C75H⁺ bases were shifted against each other. The base-base stacking was supplemented by sugar-base stacking between the A77 five-membered ring and the C75+ ribose, which also is a quite efficient dispersion-controlled interaction similar to base stacking 72. Presence of the open A-minor interaction and canonical C75 (e.g., in simulations PreC41+Mg, PreC41+, and ProC41+) is characterized by reduced base-base overlap and improved sugar-base stacking (Supplementary Fig. S5b ).

When analyzing results of MD simulations, it is important to assess their quality. The simulations appear to provide good agreement with the available x-ray structures of the HDV ribozyme, as judged, for example, by RMSD, positions of bases, etc. 37. The backbone is inherently more difficult to describe by force fields compared to interactions involving the rigid nucleobases, since i), the backbone is anionic (and thus polarizable), and ii), the constant point atomic charges might not work equally well for distinct combinations of backbone torsion angles. Thus, the quality of the backbone description is becoming one of the main issues for nucleic acids simulations. Recent B-DNA simulations reported unexpected α-γ flips of the B-DNA backbone. Such switches lead to long-lived backbone substates with concomitant changes in B-DNA geometry 40,73,74. Major problems with DNA backbone topology were reported for G-DNA loops 75.

In the helical segments of our simulated structures, the backbone is rather rigid, with backbone torsion angles close to those of canonical A-RNA. There are, however, backbone flips mainly in nonterminal residues of the double-helical segments. These flips entail simultaneous transitions of the α (from 295 to 155°) and γ (from 55 to 180°) torsion angles, which are compensatory and do not affect arrangement of the helix. Switching occurs stochastically in all analyzed simulations and is reversible. The resulting geometry matches a modified A-type backbone conformation commonly observed in RNA crystal structures and described as family number 24 55. Fig. S6 in the Supplementary Material represents the population of phosphate switches observed in double-helical segments in two representative simulations (PreC41+ and ProC41+). Lifetimes of the individual flipped substates range from 2 to 10ns, and the number of flipped nucleotides happens to be smaller at the end of the simulations than in the middle of them. Fig. 8 shows a typical reversible backbone switch. Similar behavior was noticed also for multiple 25ns simulations of the 23S rRNA sarcin-ricin loop motif 76. Thus, for RNA, we do not observe any cumulation of the backbone flips, clearly contrasting behavior of B-DNA simulations 40,74.

Display large version of this figure

Figure 8

Characteristic values of backbone torsion angles α (black line) and γ (gray line) connected with the phosphate switch in the canonical helix P1 (residuum C33) in the simulation PreC41+.

There is a wide range of other backbone geometries in the nonhelical segments of the HDV ribozyme (Fig. 9). Many of them fall into typical backbone families observed in crystal structures of 23S and 5S ribosomal RNAs 55. By contrast, the residues U23-U27 of loop L3 do not match any established backbone family, which may reflect their dynamic disorder, their uniqueness, or it may be due to a limited crystallographic resolution. In fact, many nucleotide conformations as refined in x-ray crystal structures cannot be classified as one of the typical conformational classes 55. Our simulations reveal multiple transient and permanent backbone switches in the nonhelical regions (Fig. 9) which, however, do not distort the ribozyme. Most such switches represent transitions between some known RNA backbone families. In many other cases, the switch shifts the backbone from crystallographically observed conformation that does not belong to an established backbone family to a conformation that can be assigned to a specific RNA backbone family. Thus, our simulations reduce the number of nucleotides with unidentified backbone geometry compared to crystal structures. Such shifts were observed, for instance, in the L3 region (Fig. 9) and were especially apparent in simulation ProC41+C75+ with unfolded L3, where the unfolded L3 backbone became mostly canonical (not shown). The simulations preserve key noncanonical backbone arrangements of several residues. For instance, C41H⁺ retained its backbone topology that is specific for low-rise (A-platform like) dinucleotide steps 55.

Display large version of this figure

Figure 9

Secondary structure of the HDV ribozyme, with nucleotide backbone families 55 shown in the starting C75U precursor crystal (left) and PreC41+ simulated structure (right). Gray and black upper cases represent canonical and established noncanonical backbone families, respectively; black lower case letters represent an unidentified backbone angle combination. Temporary backbone flips seen in the simulation are marked by open circles, whereas permanent backbone switches are indicated by gray circles.

We have carried out a set of explicit solvent molecular dynamics simulations to investigate the structural dynamics of the HDV ribozyme, with special emphasis on its interactions with ions and water. Careful analysis of the backbone behavior in the simulations suggests that the force field is performing reasonably well. Specifically, there are no irreversible backbone flips in the canonical regions, whereas key backbone segments that are critical for the complex 3D fold of the molecule are fairly stable. These observations obviously do not rule out the occurrence of force-field related problems on a longer timescale; in fact, it is unlikely that simple classical molecular mechanics would fully describe an anionic sugar-phosphate with its intricate electronic structure.

The compact tertiary structure of the HDV ribozyme is associated with multiple regions of deep negative electrostatic surface potential, which lead to several prominent and often uniquely structured cation-binding sites. When integrating all the available experimental and computational data, we suggest that the two most significant cation binding sites correspond to Mg²⁺ cations bound at the active site (MG1) and the major groove of P4 (MG2) as seen in the precursor crystal structures. The simulations reveal that the ability of the active site to accommodate the catalytic Mg²⁺ is determined by the ESP and the size of the pocket. Very low (negative) ESP and a larger catalytic pocket in the precursor structure are favorable for Mg²⁺ binding. In accord with the crystal structures, the divalent cation remains stably bound at the catalytic pocket in simulations of the C75U mutant and C75 wild-type precursors. It makes direct and water mediated contacts with the scissile phosphate, including a temporary G1(O5′)…Mg²⁺ water bridge as necessary for the general acid catalysis by the hydrated Mg²⁺ cation. Protonation of C75 in the precursor results in product-like deep binding of C75H⁺ in the active pocket and appears to be incompatible with binding of the catalytic Mg²⁺ cation at the active site. These findings provide additional support for a role of C75 as the general base during catalysis complemented by concerted action of a hydrated Mg²⁺ ion as general acid. The Mg²⁺ cation is not stable at the active site of the product form of ribozyme, in agreement with the crystallographic data. This observation is likely due to the significantly higher (less negative) ESP and shallower catalytic pocket of the product relative to the precursor.

A bound MG1 Mg²⁺ cation also stabilizes the structure of loop L3 (particularly the U20/G25 basepair), which shows rearrangements in the absence of Mg²⁺. We thus propose that the collapse of P3 and L3 toward the center of the ribozyme and the related flip of the noncanonical U20/G25 basepair as observed in crystal of product ribozyme might be a consequence of Mg²⁺ expulsion from the active site after the cleavage. We believe that the picture of the Mg²⁺ binding emerging from our simulations is qualitatively correct, though it is fair to admit that description of divalent cations as simple van der Waals probes with central point charge of +2 is a crude approximation. Thus, details of Mg²⁺ binding should be assumed to be imperfect. In addition, the simulation timescale is very limited for divalent cations.

Simulations carried out entirely in the absence of Mg²⁺ reveal that the MG1 and MG2 sites are occupied by monovalent cations with occupancy close to 100%. The catalytic pocket (MG1 site) typically accommodates two monovalent cations simultaneously. In contrast to Mg²⁺, bound monovalents are fluctuating in the pocket and thus are not likely to be detected by x-ray crystallography. The exchange of Na⁺ between the pocket and the bulk is slow in the simulations, so that no statistics could be computed. However, simulations utilizing different initial ion distributions show that the ions in the pocket are not incidentally trapped. When comparing with literature data, the catalytic pocket in the precursor HDV ribozyme is the most prominent Na⁺ binding site in RNA characterized by MD simulations so far 20,29,36. The simulations also suggest the presence of two significant cation-binding sites located near the J4/2 junction, as indicated by Pb²⁺ induced scission experiments 65. The description of monovalent ions is considerably more accurate compared with divalents and we suggest that the simulations are fairly sufficient to deal with highly occupied monovalent binding sites in RNA molecules.

Structural dynamics of the HDV ribozyme is associated with long residency hydration sites, including several water-mediated basepairs. A water-mediated interaction is also suggested as a substate for the fluctuating A-minor type II interaction involving A77, which is coupled with positioning of the catalytically critical nucleotide 75 in the catalytic pocket. The long-residency water molecules usually make contacts with electronegative RNA atoms and regions, and thus compete with cations for the same binding sites, as exemplified in the trans-Watson-Crick/Hoogsteen U20/G25 basepair. We suggest that selectivity of the binding site with respect to water molecules, monovalent, or divalent cations is affected by the geometry of the binding site. Monovalent cations may readily substitute for divalents. This is to be considered when analyzing experimental structural data, since a bound divalent cation can be rather easily replaced by fluctuating monovalent cations or by a water molecule, which are both less likely to be detected. Description of base stacking, Hbonding, and water binding can be considered as the most accurate part of the force field, since these interactions can be well captured by the point charge electrostatic model and the Lennard-Jones potential 77,78,79.

We evidenced another kind of hitherto unreported combined cation binding and hydration events, which are especially apparent in the catalytic pocket. Many of the long-residency water molecules present in the pocket participate in the Na⁺ first ligand shell while acting as structural bridges to the RNA. Thus, when the Na⁺ cations occupy the pocket, the pocket is filled by a complex network of water molecules structured by the ions and the solute. The binding times of water molecules in the first-shell of bound Na⁺ cations range up to ∼7ns, which is almost two orders of magnitude longer than in the case of bulk Na⁺ cations. Therefore, Na⁺ plays an important role in structuring of water molecules at the active site in analogy to Mg²⁺ and could possibly participate in catalysis in the absence of Mg²⁺. Interestingly, the same conclusion can be drawn regarding water molecules participating in the second shell of the bound Mg²⁺ ions. These water molecules also form stable water bridges between the hydrated cation and solute, and their binding time is substantially enhanced compared to the second shell water molecules of bulk Mg²⁺.