Sodium and Magnesium Ion Location at the Backbone and at the Nucleobase of RNA: Ab Initio Molecular Dynamics in Water Solution.

The interactions between Na+ or Mg2+ ions with different parts of single-stranded RNA molecules, namely, the oxygen atoms from the phosphate groups or the guanine base, in water solution have been studied using first-principles molecular dynamics. Sodium ions were found to be much more mobile than Mg2+ ions and readily underwent transitions between a state directly bonded to RNA oxygen atoms and a completely solvated state. The inner solvation shell of Na+ ions fluctuated stochastically at a femtosecond timescale coordinating on average 5 oxygen atoms for bonded Na+ ions and 5.5 oxygen atoms for solvated Na+ ions. In contrast, the inner solvation shell of Mg2+ ions was stable in both RNA-bonded and completely solvated states. In both cases, Mg2+ ions coordinated 6 oxygen atoms from the inner solvation shell. Consistent with their stable solvation shells, Mg2+ ions were more effective than Na+ ions in stabilizing the RNA backbone conformation. The exclusion zones between the first and second solvation shells, solvation shell widths, and angles for binding to carbonyl oxygen of guanine for solvated Na+ or Mg2+ ions exhibited a number of quantitative differences when compared with RNA crystallographic data. The presented results support the distinct capacity of Mg2+ ions to support the RNA structure not only in the crystal phase but also in the dynamic water environment both on the side of the phosphate moiety and on the side of the nucleobase.

Introduction

Ribonucleic acid (RNA) molecules support essential physiological activities in living systems, including control of gene expression,[1−3] encoding genetic information as a blueprint for protein synthesis,[4,5] transferring amino acids into ribosomes,[6] carrying out protein translation in ribosome complexes,[7,8] catalysis of biochemical reactions,[9,10] and regulation of cellular responses to environmental stimuli.[11−13] This wide range of RNA functions is accomplished by the capacity of RNA to fold into specific three-dimensional conformations. Because the high negative charge of the sugar-phosphate backbone works against RNA folding into a compact structure, positively charged metal ions promote folding by reducing the electrostatic repulsion between RNA phosphate groups.[14] The stabilizing effect of metal ions on the RNA structure, however, varies between different ion types. In particular, Mg2+ ions at millimolar concentrations enhance dramatically the stability of RNA tertiary structures that are otherwise only marginally stable in the presence of high monovalent cation concentrations.[15,16] Magnesium ions also enhance catalytic activity of ribozymes[17] and mediate catalytic processes in ribozyme active sites.[18,19] The fact that RNA molecules possess specific Mg2+ binding motifs[20] further implies that the role of Mg2+ ions in RNA folding is not limited to simple compensation of electrostatic repulsions, but it is driven by very specific contextual requirements.[21,22]

The interactions between metal ions and RNA have been previously studied employing a number of experimental techniques. The RNA binding sites of monovalent and divalent metal ions in solution and the crystal phase have been explored by nuclear magnetic resonance[23−25] and vibrational spectroscopy methods including Fourier-transform infrared and Raman spectroscopy.[26] The locations of the ions in the RNA crystal phase could be determined with improved spatial resolution by X-ray crystallography;[27,28] however, the structure of RNA molecules in the crystalline phase may differ from that in solution.[29] A serious problem is also the distinguishing of Na+ ions from the water molecules.[30]

Computer simulation of RNA dynamics is commonly performed with the use of classical molecular dynamics (MD)[31−33] and Monte Carlo algorithms.[34−36] The accuracy of these classical methods for MD simulations, however, depends critically on the appropriate choice of force field parameters.[37] Other limitations in resolving ion dynamics in classical simulations is the difficulty of including explicit water molecules. The very slow exchange of water in the tight Mg2+ hexahydrate complex is a particular problem.[38−41]

To complement the available experimental and theoretical methods, first-principles quantum chemical simulations can be employed for detailed characterization of the dynamical interactions between metal ions and RNA. For example, the interactions of Na+ and Mg2+ ions with the phosphate groups have already been studied with Born–Oppenheimer MD, employing density functional theory (DFT) [at the generalized gradient approximation (GGA) level] to calculate the forces between atoms.[29] The structural model that was used consisted of the sugar-phosphate backbone of single-stranded RNA (ssRNA), without taking into account the nucleobases.

In the present work, we complemented the model, studied earlier,[29] with guanine and cytosine bases in order to make it much more realistic. In this way, we have the whole molecular structure of ssRNA. The ab initio molecular dynamic simulations are focused on the interaction of Na+ and Mg2+ ions with electronegative centers of ssRNA, namely, the oxygen centers from the sugar-phosphate backbone and the oxygen centers of the nucleobases. We explicitly model all interactions at the ab initio level and present real-time quantum dynamics using water as a liquid phase and solvent. This provides important insights into the possible function of X-ray crystallographic structures (such as those deposited in the Protein Data Bank), for which, as a rule, no water phase is presented.[42]

Noteworthily, the present study is primarily focused on the local interactions of the cations with ssRNA fragments, including phosphate groups and parts of the nucleobases. Because the dimensionality of the quantum wave function in ab initio methods grows exponentially with the total number of atoms, we are restricted to work with simple dinucleotide ssRNA fragments. This precludes direct simulation of secondary and tertiary RNA structure transitions, which occur at the microsecond to millisecond timescale.[43] The employment of ab initio methods, however, grants us the ability to study the local interactions at the femtosecond to picosecond timescale with high fidelity,[44] without the use of the predetermined force fields of the classical MD, and allows for direct extrapolation of our conclusions based on these precise simulations to physiological water-based solvents in living organisms.

Methods

Ab initio Born–Oppenheimer MD simulations of ssRNA with Na+ or Mg2+ countercations in water solution are performed using the freely available CP2K/Quickstep package.[45] Kohn–Sham density functional theory[46,47] with GGA is applied based on the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional.[48] The PBE was chosen as a tradeoff between the lower-precision LDA functional and the hybrid ones, which are too computationally expensive for the current task. Also, the studied interactions in the present article, between the charged Na+ and Mg2+ ions on one hand and the phosphate groups, H2O (O) atoms, and guanine on the other hand, are expected to be electrostatic in nature, which are properly described by the PBE functional. For all atoms, the basis set DZVP-MOLOPT-SR-Goedecker-Teter-Hutter (GTH) is employed, which is suitable for molecular studies in the condensed phase.[49] For reducing the computational cost, Gaussian and plane-wave expansion sets are used for expanding the electronic wave functions.[50,51] Only the valence electrons are explicitly considered. Their interaction with the remaining ions is described using the pseudopotentials of GTH.[52,53] The charge density cutoff of the finest grid level is equal to 400 Ry. The number of used multigrids is 5.

All simulations are carried out in the NVT ensemble with a timestep of 1 fs. The temperature was set to 320 K using canonical sampling through a velocity rescaling thermostat.[54,55] For each simulation, statistics is obtained after the initial 1 ps, which is required for thermal equilibration of the molecular system (Figure S11). The studied systems consist of a dinucleotide ssRNA fragments containing either two cytosine bases (CC) or one guanine and one cytosine (GC) and the metal countercations, two Na+ or one Mg2+, necessary to neutralize the negative charges of the phosphate groups. The dinucleotide abbreviations used throughout the article follow the standard FASTA format.[56] The whole system is hydrated in a periodic box with dimensions of 12.645 Å × 17.3 Å × 17.3 Å. For simulations with CC 2Na+, GC 2Na+, CC Mg2+, or GC Mg2+, the box contains 106, 104, 110, or 107 water molecules, respectively, which assures that the density of the water phase is equal to 0.99 g/cm3. Periodic boundary conditions for all axes are used. The dinucleotide skeleton is oriented along the shortest axis (12.645 Å). The only spatial structure that the dinucleotide can form in this case is a straight ssRNA chain. It should be noted that such a periodic system has certain limitations concerning its conformational plasticity – the RNA chain cannot form secondary and tertiary structures, for example, loops, and it also cannot form double-stranded stems. Thus, we simulate an “infinite” straight chain of single-stranded RNA (ssRNA). Also, the simulation time of 100 ps is relatively short to capture all possible conformations of the dinucleotide. We emphasize that when the definition (ssRNA) is addressed in this work, concerning our DFT simulations, the dinucleotide model that we used should be considered. Despite the limitations, as can be seen further, the model is able to describe the primary modes of local interactions between the solvated counterions and RNA. After the initial placement of the metal ions, there are no structural constraints imposed on any of the atoms. Visualization and analysis of the MD are performed with visual MD version 1.9.3.[57] Radial distribution functions (RDFs) between elements are computed with the default bin width of 0.1 Å. Coordination numbers in the inner solvation shells were computed using custom scripts with a cutoff at 3.2 Å or 2.8 Å for Na+ or Mg2+ ions, respectively. All simulations are started from geometry-optimized structures. After 1 ps, thermal equilibration of the systems is observed (a thermostat time constant of 50 fs or 50 steps). The temperature of the system remained stable (±30 K) during the whole simulations. After the 20th ps, the thermostat was practically switched off (a time constant of 1 ps or 1000). Also, after the 20th ps, no drift of the potential energy of the systems was observed.

Four major conformations of the RNA backbone at the phosphate moiety are distinguished:[58]

The absolute value of the angle between the adjacent C–O bond and the plane formed by the phosphodiester bond was considered >20° for g and ≤20° for t. Scattergrams were computed with the use of custom scripts in Microsoft Excel and Wolfram Mathematica.

Six ssRNA simulations with metal ions are performed. The initial location of cations and the duration time of the simulations are as follows:

CC 2Na+ simulation for 102.68 ps with two bound sodium ions, directly coordinated at phosphate group P1 or P2 in a CC ssRNA fragment.

CC 2Na+ simulation for 103.90 ps with two solvated sodium ions, coordinated through their inner water shells at phosphate group P1 or P2 in a CC ssRNA fragment.

GC 2Na+ simulation for 103.23 ps with two bound sodium ions, directly coordinated at the oxygen atom of guanine or the hydroxyl group attached to the pentose ring in the 2′ position of guanosine in a GC ssRNA fragment.

CC Mg2+ simulation for 101.68 ps with one bound magnesium ion, directly coordinated at phosphate group P2 in a CC ssRNA fragment.

CC Mg2+ simulation for 105.45 ps with one solvated magnesium ion, initially coordinated through its inner water shell at phosphate group P2 in a CC ssRNA fragment.

GC Mg2+ simulation for 101.68 ps with one bound magnesium ion, directly coordinated at the oxygen atom of guanine in a GC ssRNA fragment.

Results and Discussion

Mobility of Na+ Ions

The interaction of sodium ions with RNA was studied in three simulations, each of which with a time duration of over 100 ps (see the selected snapshots in Figure ). The interaction of sodium ions with the phosphate groups from the RNA backbone was studied in a model consisting of an ssRNA fragment containing two cytosine bases (CC) in water solution, whereas the interaction with the guanine base was studied in a model with an ssRNA fragment containing one guanine and one cytosine base (GC).

Figure 1

Selected snapshots from the sodium-containing ssRNA simulations: (a) two bound Na+ ions directly coordinated at phosphate groups P1 and P2 in a CC ssRNA fragment; (b) two solvated Na+ ions interacting through their inner hydration shells with phosphate groups P1 and P2 in a CC ssRNA fragment; (c) two bound Na+ ions directly coordinated at O5 from the guanine base or O6 from the hydroxyl group attached to the pentose ring in the 2′ position of guanosine in a GC ssRNA fragment.

In the CC 2Na+ simulation with two bound sodium ions, directly coordinated at phosphate group P1 or P2 in a CC ssRNA fragment, both sodium ions stayed near their initial positions (Figure a). One of the sodium ions, denoted as Na+1, remained directly bonded to an oxygen atom, O1 or O2, from phosphate group P1 for the whole period of the simulation (1–102 ps). The other sodium ion, Na+2, interacted directly with oxygen atom O3 from phosphate group P2 for most of the simulation duration. However, for a short period of time (82–90 ps), it was solvated in water and interacted with both oxygen atoms O3 and O4 through water molecules from its inner hydration shell. The RDF for the Na+1-O1 distance during the time intervals 1–75 and 97–102 ps and for the Na+1-O2 distance during the interval 75–97 ps exhibited a maximum at 2.35 Å (Figure S1a–d), which indicates a direct contact between Na+1 and an oxygen atom from phosphate group P1. The sodium ion Na+2 was also directly bound to oxygen atom O3 from P2 during the intervals 1–82 and 90–102 ps, with the corresponding RDF peaks at 2.55 and 2.75 Å (Figure S1e,f,h). For a short time period 82–90 ps, however, Na+2 was solvated and interacted with oxygen atoms O3 and O4 through a water molecule as evidenced by the corresponding RDF peaks at 3.65 and 4.25 Å (Figures a and S1g). Distances between the phosphate groups (P atoms) and the sodium ions throughout the simulation are presented in Figure a. It should be noted that the distance to the P atom is roughly 1 Å longer than the distance to the nearest O atom of the PO4– group. Experimental results obtained with X-ray crystallography indicate Na+-O(P) distances in the interval 2.65–2.75 Å when the ion is directly bonded to the phosphate group.[59,60] Comparison of the ab initio results is made with already performed classical MD simulations employing CHARMM27 and AMBER (ver. 2.0).[29] Water molecules were described by TIP3P potential. The optimal distance Na+–O(P) in the classical MD simulation with AMBER is 2.35 Å, the same as the distance in the ab initio MD. The location of the cation in the case of classical MD (AMBER) differs from that obtained with the ab initio method. Calculations with AMBER suggest bidentate configuration (Na+ bonded to two oxygen atoms of a phosphate group), while DFT results suggest monodentate configuration. The simulation with CHARMM also suggests monodentate configuration. Both classical methods suggest Na+–O(P) shorter by 0.2 Å by the DFT method.

Figure 3

Dynamics of the distances from the Na+ and Mg2+ ions to particular phosphorus and oxygen atoms of the systems: (a) Na+1-P1 and Na+2-P2 from the CC 2Na+ simulation with two bound sodium ions; (b) Na+1-P1 and Na+2-P2 from the CC 2Na+ simulation with two solvated sodium ions; distances to the closest phosphate groups are presented; (c) Na+1-O5 and Na+2-O6 from the GC 2Na+ simulation with two bound sodium ions; (d) Mg2+-P1 from the CC Mg2+ simulation with one bound magnesium ion; (e) Mg2+-P1 from the CC Mg2+ simulation with one solvated magnesium ion; (f) Mg2+-O5 from the GC Mg2+ simulation with one bound magnesium ion.

In the system containing CC and 2Na+, the two sodium ions remain completely solvated during the whole simulation, Figure b. The sodium ion Na+1 interacted with oxygen atom O1 from phosphate group P1 during the interval 1–20 ps (Figure b) through a water molecule as indicated by the RDF peak at 5.05 Å (Figure S2a,b). In the rest of the simulation, 20–103 ps, Na+1 remained fully solvated with RDF peaks at 6.95 and 8.05 Å from O1 and O2, respectively (Figure S2c). In the last time interval (20–103 ps), the distance between the Na+ ion and an oxygen atom of the phosphate group shorter than 6.5 Å indicates the interaction through two water molecules, and a higher distance indicates the interaction through three water molecules. The other sodium ion, Na+2, interacted with O3 from P2 through a water molecule for most of the simulation duration, with the corresponding RDF peaks at 4.75 or 4.65 Å, respectively, for the time intervals 1–76 and 82–103 ps (Figures b and S2d,e,g). Na+2 interacted with the phosphate group through two water molecules for a short time in the interval 76–82 ps (Figures b and S2f).

In the GC 2Na+ simulation with two bound sodium ions, one of the Na+ ions interacted directly with oxygen atom O5 from the guanine base, while the second Na+ ion was coordinated to O6 from the hydroxyl group attached to the pentose ring in the 2′ position of guanosine in a GC ssRNA fragment (Figures c and 3c). The sodium ion Na+1 was directly bound to O5 from guanine, with an RDF peak at 2.55 Å for the time intervals 1–31 and 41–103 ps, respectively (Figure S3a,b,d). Thermal fluctuations led to intermittent detachment of Na+1 from O5 in the interval 31–41 ps as evidenced by two RDF peaks at 3.15 and 2.75 Å (Figure S3c). The other sodium atom Na+2 was stably attached to oxygen atom O6 in the pentose ring of guanosine with an RDF peak at 2.45 Å for most of the simulation time, except for a short period, 31–33 ps, during which thermal noise created two RDF peaks at 2.85 and 3.15 Å (Figure S3e–h).

Mobility of Mg2+ Ions

We modeled the interaction of magnesium ions with phosphate groups from the RNA backbone or with the guanine base in the same molecular systems consisting of CC or GC ssRNA fragments in water (Figure ). The simulations were performed for at least 100 ps.

Figure 2

Selected snapshots from the magnesium-containing ssRNA simulations: (a) one bound Mg2+ ion directly coordinated at phosphate group P2 in a CC ssRNA fragment; (b) one solvated Mg2+ ion interacting through its inner hydration shell with phosphate groups P1 and P2 in a CC ssRNA fragment; (c) one bound Mg2+ ion directly coordinated at O5 from the guanine base in a GC ssRNA fragment.

In the CC Mg2+ simulation with a bound magnesium ion, Mg2+, directly coordinated at phosphate group P2 in a CC ssRNA fragment, the magnesium ion stayed firmly attached to the phosphate group (oxygen atom O3) (Figures a and 3d) with an RDF peak at 2.15 Å for the whole period of the simulation 1–101 ps (Figure S4a,b). The inner hydration shell of Mg2+ invariably consisted of 5 water molecules in addition to the O3 center. One of the water molecules, denoted as W1, located on the opposite side of O3 along the Mg2+-O3 axis was exchanged with another water molecule W2 at the end of the time interval 28–29 ps (Figure a). The RDF for the Mg2+-W1O or Mg2+-W2O distance measured between Mg2+ and the water oxygen atoms during the time interval 1–28 ps exhibited a maximum at 2.15 or 4.35 Å, respectively (Figure S4c,d). In the rest of the simulation duration, 29–101 ps, the RDF peaks for Mg2+-W1O or Mg2+-W2O were reversed, 4.35 or 2.15 Å, which indicated the exchange of W1 with W2 (Figure S4e,f). The observed exchange of water molecules from the inner hydration shell of Mg2+ is consistent with dynamic instability of water molecules located opposite to the oxygen atom from the phosphate group.

In the CC Mg2+ simulation with a magnesium ion solvated by 6 water molecules, the Mg2+ ion interacted through a water molecule from its inner hydration shell with one of the bridging oxygen atoms from phosphate group P1 or P2 (Figures b, 3e, and S5). Initially, the Mg2+(H2O)6 complex interacted with 3′ bridging oxygen atom O8 from phosphate group P2 during the interval 1–7 ps (Figure b) through a water molecule, W4, that formed a hydrogen bond as indicated by the W4H1-O8 RDF peak at 2.55 Å (Figure S5e,f). Then, in the time interval 7–42 ps, the whole complex moved away from the phosphate groups (Figure b) with RDF peaks at 3.45, 3.25, 4.25, and 3.75 Å for W3H1-O7, W3H2-O7, W4H1-O8, and W4H2 O8 distances, respectively (Figure S5c,g), indicating the interaction of the Mg2+ with the phosphate group through two water molecules. For the rest of the simulation, 42–105 ps, the complex approached the RNA backbone again and interacted with the 5′ bridging O7 atom from the P1 moiety through a water molecule, W3, that formed a hydrogen bond as indicated by the W3H2-O7 RDF peak at 2.85 Å (Figure S5a,d).

In the GC Mg2+ simulation with a bound magnesium ion to the guanine base, Mg2+ remained directly attached to oxygen atom O5 of guanine for the whole simulation duration 1–101 ps with an RDF peak at 2.15 Å (Figures c, 3f, and S6a,b). The inner hydration shell of Mg2+ retained the same 5 water molecules in addition to the O5 atom and allowed for the intermittent formation of hydrogen bonds between the oxygen atom of water molecule W5 and the hydrogen atoms from the amine group of the cytosine base stacked below guanine during the time intervals 1–7 and 22–101 ps (Figures c and S6c,d,f). The observed Mg2+-mediated interaction between vertically stacked guanine and cytosine bases is consistent with the essential role of Mg2+ ions in facilitating high structural complexity and folding arrangements that allow RNA molecules to perform diverse life-supporting cellular functions.[22]

In one example, the experimentally determined distance between the directly bonded Mg2+ ion and the O6 oxygen atom of guanine has been experimentally determined by X-ray crystallography.[61] The distance between (Mg2+-O6G) for the Mg2+ 201 ion from PDB 2A43 is 2.37 Å. Classical MD simulations with CHARMM27 and AMBER (ver. 2.0)[29] for magnesium ions, interacting with RNA moieties, show similar positions of the cations with respect to the phosphate groups. For the directly bonded Mg2+ ion, monodentate coordination to one of the O atoms of the PO4– group is observed. The most commonly observed Mg2+-O distance is found to be 0.1–0.3 Å shorter than that for the DFT simulations (this work).

Solvation Shell Structures of Na+ and Mg2+ Ions

In their solvated complexes, sodium ions have a poorly defined radius of about 3.2 Å (Figure a) for their inner solvation shell due to frequent exchange of water molecules.[29] This radius cutoff at 3.2 Å corresponds to a local minimum between the first and second solvation shells visualized through histograms of oxygen atoms surrounding the Na+ ion (Figure S7). The number of coordinated oxygen atoms in the inner solvation shell of Na+ ions varied from 3 to 8 (Figure S8), which on average was higher for solvated Na+ ions compared with Na+ ions coordinated to oxygen centers of RNA (Table S1). In the simulation with two Na+ ions bound to phosphate groups from the CC ssRNA backbone, the average coordination number was 5.0 ± 0.7, whereas for the simulation with two solvated Na+ ions, the average coordination number was somewhat higher, 5.5 ± 0.8. Interestingly, in the GC ssRNA simulation, the Na+1 ion bound to the guanine base exhibited an even lower average coordination number of 4.7 ± 0.5 (toward water oxygen atoms and the carbonyl group of guanine), whereas the Na+2 ion bound to the hydroxyl group of the sugar backbone had an average coordination number of 5.4 ± 0.6, which was closer to the coordination number of solvated ions. The histograms for these simulations are shown in Figure S7, and the variations of the coordination number with simulations time are presented in Figure S8.

Figure 4

Histograms of the O atoms in the solvation shells of the Na+ ion (a) and Mg2+ ion (b) averaged over all performed CC and GC ssRNA simulations. The illustrated widths of the exclusion zones (gray) and the second coordination shells (pink) were determined in an extensive crystallographic PDB survey by Leonarski et al.(66) The exclusion zone in the crystal phase for Na+ is 2.8–3.6 Å, solvated: 3.1–3.3 Å (not depleted). The second coordination shell for Na+ in the crystal phase is 3.6–4.9 Å; solvated: 3.3–5.2 Å. The exclusion zone in the crystal phase for Mg2+ is 2.2–3.8 Å, solvated: 2.7–3.4 Å. The second coordination shell for Mg2+ in the crystal phase is 3.8–4.6 Å, solvated: 3.4–4.9 Å.

We can compare our theoretical results with experimental data from the literature. For example, fully hydrated Na+ ions in complexes with RNA motifs are studied with X-ray crystallography.[62−64] It is observed that sodium ions coordinate 5 water molecules in the first solvation shell (Na+8—PDB 434D and Na+19—3ND4).[62,63] A structure with 6 water molecules was also reported (Na+102—PDB 466D).[64] Generally, for both bonded and completely solvated Na+ ions, the most frequently observed (via X-ray crystallography) solvation numbers are 4, 5, and 6.[60] X-ray crystallographic results, however, may differ from the biological systems as they represent frozen in time system and ion solvation is altered. Classical molecular mechanics (MM) simulations give different coordination numbers for the sodium ion, dependent on the force field.[65] For solvated Na+ ions, the coordination number can range from 7.1 (CHARMM27) to 5.3 (AMBER99). These variations are caused by the different parameters in the used force fields. With ab initio methods, however, such large variations are avoided as the properties of the different systems are calculated from the wave function. In this regard, for the solvated Na+ ion, we found a coordination number of 5.5 (this work), while an earlier study with the BLYP functional and the combined BLYP/MM method gave average values of 5 and 5.2, respectively,[65] which is closer to the AMBER99 force field.

Solvated magnesium ions have a well-defined radius of 2.8 Å (Figure b) for their inner solvation shell.[29] This radius cutoff at 2.8 Å includes a part of the completely depleted region of the exclusion zone between the first and second solvation shells of Mg2+ ions (Figure S9). In the solvated state, the Mg2+ ion exhibited an average coordination number of 6 oxygen atoms throughout the whole time period of the simulation (Figure S10). In the CC ssRNA simulation, the Mg2+ ion bound to a phosphate group was able to exchange a single water molecule from its inner solvation shell (Figure a). During the exchange, the Mg2+ ion transiently (less than 1% of the simulation time) exhibited a coordination number of 5 oxygen atoms. The exchange of the water molecule is obviously governed by local fluctuations of the water environment. Our simulation is relatively short (100 ps) for this result to be used to determine the lifetime of a water molecule in the Mg2+ solvation shell. Other studies indicate this time to be in the millisecond range.[67] In the GC ssRNA simulation, the Mg2+ ion bound to the guanine carbonyl oxygen also transiently exhibited for a very short time a coordination number of 5 oxygen atoms but did not exchange any water molecules (Figure c). X-ray crystallographic data and MM simulations with either CHARMM or AMBER support a solvation number equal to 6 for Mg2+ in both bonded and hydrated states.[60] The ab initio dynamics of Na+ and Mg2+ ions in solution supports qualitatively the results reported in an extensive crystallographic PDB survey by Leonarski et al.(66) but also reveals some quantitative differences. In the crystal phase, the reported exclusion zones between the first and second solvation shells were 2.8–3.6 Å for Na+ and 2.2–3.8 Å for Mg2+. Furthermore, the reported widths of the second solvation shells were 3.6–4.9 Å for Na+ and 3.8–4.6 Å for Mg2+ (Figure ). In comparison, the histograms obtained from the quantum dynamics in solution revealed that the exclusion zones between the first and second solvation shells were narrower; for Na+, the exclusion zone extending in the region 3.1–3.3 Å is not depleted, while for Mg2+, the exclusion zone extending in the region 2.7–3.4 Å is depleted. Conversely, the widths of the second solvation shells were wider, 3.3–5.2 Å for Na+ and 3.4–4.9 Å for Mg2+. Thus, our results may assist in identification of Na+ for Mg2+ ions in crystallographic structures, particularly in the presence of solvating water molecules, and may help in avoiding misidentifications.

Na+ and Mg2+ Interaction with the Guanine Base

The carbonyl oxygen of guanine is a binding site with strong affinity for metal ions, which has been determined to exhibit an exclusion zone (i.e., in this range, the corresponding ion is not observed) at conical angles from 160 to 180° in crystallographic structures.[66] In the crystal phase, the C=O···Na+/Mg2+ angle values were reported to be within the range from 100 to 160°.[66] Remarkably, in solution, the average conical angles were manifestly different for each ion type, with the maximum probability at 108° for Na+ versus 144° for Mg2+, with only 4.4% overlap of the corresponding angle distribution histograms (Figure ). This effectively makes Na+ and Mg2+ distinguishable in solution based on their binding angles to carbonyl oxygen of guanine. The exclusion zone for Na+ is <90° and 125–180° (Figure a), whereas for Mg2+, it is <120° and 170–180° (Figure b).

Figure 5

Histograms of the conical angles of the Na+ ion (a) and Mg2+ ion (b) relative to the C=O group in the guanine base in water solution. The illustrated exclusion zones were determined in an extensive crystallographic PDB survey by Leonarski et al.(66)

Conformational Effects on the Phosphoribose Backbone

Folding the phosphoribose backbone is important for both the structural stability and functional versatility of RNAs.[68,69] To study the effects of metal ions on the RNA backbone conformation, we have computed two-dimensional scattergrams (Figure ) of the angles between the adjacent C–O bonds and the plane formed by the phosphodiester bond and estimated the percentage of time spent by the RNA backbone in gg, gt, tg, or tt conformations (Table S2). Consistent with their dynamic solvation shells, Na+ ions were able to stabilize the RNA backbone in the gg conformation only when directly bound to the phosphate groups (Figure a) but allowed significant transition into the tg conformation when the Na+ ions were solvated or bound to the guanine base, respectively, for 2.3 or 7.1% of the simulated time period (Figures b,c and Table S2). In contrast, consistent with their stable solvation shells, Mg2+ ions were able to stabilize the RNA backbone in the gg conformation both when directly bound to the phosphate groups and when the interaction was mediated through a water molecule (Figures d,e and Table S2). Interestingly, the direct interaction of the Mg2+ ion with the carbonyl group of guanine was accompanied with the indirect interaction with the amino group of cytosine through the water molecule during most of the simulation, which may have caused rare transitions into tg or gt conformations (Table S2). These results support the proposed functional role of Mg2+ ions as potent stabilizers of the RNA secondary structure.[22,70,71]

Figure 6

Scattergrams of the angles between the plane formed by the phosphodiester bond and each of the C1–O or C2–O bonds for (a) CC 2Na+ simulation with two bound sodium ions, (b) CC 2Na+ simulation with two solvated sodium ions, (c) GC 2Na+ simulation with two bound sodium ions, (d) CC Mg2+ simulation with one bound magnesium ion, (e) CC Mg2+ simulation with one solvated magnesium ion, and (f) GC Mg2+ simulation with one bound magnesium ion.

Conclusions

Density functional dynamic simulations of ssRNA interacting with Na+ or Mg2+ countercations in water solution have identified high-affinity binding sites, which coordinated the metal ions for extended periods of time of over 100 ps. While water molecules in the solvation shells of Na+ ions exchange easily, the Mg2+ ions demonstrated lower mobility, much more stable solvation shells, and therefore a stronger stabilization effect upon the phosphoribose backbone of RNA. Our simulations also suggested that magnesium may not only affect the folding of the ssRNA by interacting simultaneously with two distant phosphate groups as suggested earlier[29,60] but also mediate the interaction between nucleobases in these molecules by coordinating simultaneously to both bases via the magnesium ion and one of the water molecules solvating it. This result is consistent with the observation that Mg2+ ions facilitate RNA folding and participate in hairpin formation.[72] It also corroborates experimental evidence that regions with higher GC content are likely to have a more stable secondary RNA structure.[73−75]

The interactions of Na+ or Mg2+ ions with the guanine base in water solution exhibited differential distributions of the conical angle between the ionic bond and the carbonyl group, indicating distinct types of hybridization of the latter as a possible cause of the modified interaction of guanine in RNA secondary or tertiary structures for varying Na+ and Mg2+ electrolyte composition.[76,77] In addition, the typical ranges of the angles C–O–Na/Mg obtained in our simulations in water media suggest somewhat different regions compared to the ranges reported for the crystal structures.[66] These results support the interpretation that the differential effects of Na+ and Mg2+ ions upon the RNA structure may have dynamic origin and arise due to different abilities of the two types of ions to coordinate oxygen atoms, support their inner hydration shells, and move in the water solvent.[31,71,78] The observed 100 ps-long periods of thermal stability of the inner hydration shell of Mg2+ ions, with brief subpicosecond exchange of water molecules, further suggest that the Mg2+-RNA interaction may also involve purely quantum phenomena such as tunneling through potential barriers, in addition to the classical over-the-barrier thermal transition.[79,80] This means that quantum tunneling may speed up certain conformational transitions that would otherwise need higher classical driving potentials and extended periods of thermal agitation.[81,82]

The versatility of the RNA structure and function is intertwined with the chemical fragility of RNA.[83−85] The ability of divalent cations such as Mg2+ to stabilize the RNA structure is pertinent to molecular biology research and synthesis of novel therapeutic agents, including RNA vaccines.[86,87] Thus, the identification of the Mg2+-RNA binding sites and elucidation of the involvement of the solvent in structural stabilization/destabilization could be utilized for improvement of RNA-based medical treatment and prophylaxis of disease.

Classical MD force fields, unlike the ab initio methods, give results dependent on the used parameterization. Properties of the studied systems like coordination number and the mode of interaction with the phosphate groups appear to differ with the used classical force field, especially for the Na+ ions. In this respect, ab initio calculations using explicit water provide more reliable results and can be used for proper parameterization of the classical MD methods.

Data Availability

The data underlying this article will be shared on reasonable request to the corresponding author.