Ion Condensation onto Ribozyme is Site-Specific and Fold-Dependent

TL;DR

This study reveals that ion condensation onto the Azoarcus ribozyme is highly specific, site-dependent, and fold-dependent, contrasting with traditional theories of uniform ion binding, and is influenced by the ribozyme's architecture and shape fluctuations.

Contribution

The paper demonstrates through simulations that ion binding to RNA is specific and depends on nucleotide position and folding state, challenging the notion of non-specific ion interactions.

Findings

Ion condensation is highly site-specific and fold-dependent.

Regions with high ion coordination are present even without tertiary interactions.

Ion binding patterns differ markedly from those on rigid charged rods.

Abstract

The highly charged RNA molecules, with each phosphate carrying a single negative charge, cannot fold into well-defined architectures with tertiary interactions, in the absence of ions. For ribozymes, divalent cations are known to be more efficient than monovalent ions in driving them to a compact state although Mg$^{2+}$ ions are needed for catalytic activities. Therefore, how ions interact with RNA is relevant in understanding RNA folding. It is often thought that most of the ions are territorially and non-specifically bound to the RNA, as predicted by the counterion condensation (CIC) theory. Here, we show using simulations of ${\it Azoarcus}$ ribozyme, based on an accurate coarse-grained Three Site Interaction (TIS) model, with explicit divalent and monovalent cations, that ion condensation is highly specific and depends on the nucleotide position. The regions with high coordination between the phosphate groups and the divalent cations are discernible even at very low Mg$^{2+}$ concentrations when the ribozyme does not form tertiary interactions. Surprisingly, these regions also contain the secondary structural elements that nucleate subsequently in the self-assembly of RNA, implying that ion condensation is determined by the architecture of the folded state. These results are in sharp contrast to interactions of ions (monovalent and divalent) with rigid charged rods in which ion condensation is uniform and position independent. The differences are explained in terms of the dramatic non-monotonic shape fluctuations in the ribozyme as it folds with increasing Mg$^{2+}$ or Ca$^{2+}$ concentration.

Figures

b in (Denesyuk15NatChem)). In other words, the midpoint concentration, $C_{\mathrm{m}}$, of the transition measured by $R_{\textrm{g}}$ or other parameters is the (appropriately-weighted) mean of such individual [Mg2+] at which distinct structural elements achieve the native conformation.

Equations

$$U_{\text{Ex}}={\epsilon }_{ij}\left[{\left(\frac{1.6}{r_{ij}+1.6-D_{ij}}\right)}^{12}-2{\left(\frac{1.6}{r_{ij}+1.6-D_{ij}}\right)}^6+1\right],r_{ij}\le D_{ij}$$

$$U_{\mathrm{Ele}}=\frac{q_i{q}_j}{4\pi {\epsilon }_0\epsilon r_{ij}}$$

$$c_i^{\ast }=\frac{1}{N_{\mathrm{A}}{V}_{\mathrm{c}}}{\int }_0^{r_{\mathrm{c}}}{\rho }_i(r)4\pi r^2\mathrm{d}r,$$

$$\Delta c_i^{\ast }=\frac{{⟨{\left(c_i^{\ast }-⟨c_i^{\ast }⟩\right)}^2⟩}^{\text{nicefrac}12}}{⟨c_i^{\ast }⟩},$$

$$\overline{⟨c^{\ast }⟩}=\frac{1}{N}\sum _{i=1}^N⟨c_i^{\ast }⟩,$$

$${\sigma }^2=\frac{1}{N}\sum _{i=1}^N{\left(⟨c_i^{\ast }⟩-\overline{⟨c^{\ast }⟩}\right)}^2.$$

$$Q_m(R)=z_m{\int }_0^R{\rho }_m(R^{\prime })2\pi R^{\prime }\mathrm{d}{R}^{\prime },$$

$$Q_{\mathrm{T}}(R)={\rho }_{{}_{1D}}+Q_{\mathrm{M}{\mathrm{g}}^{2+}}(R)+Q_{{\mathrm{K}}^{+}}(R)+Q_{\mathrm{C}{\mathrm{l}}^{-}}(R).$$

Full text

\corrauthor

[*]dave.thirumalai@gmail.com \papertypeArticle

Introduction

RNA molecules are negatively charged polyelectrolytes (PEs) with each molecule carrying a total charge equal to $-Ne$ where $N$ is the number of nucleotides and $e$ is the elementary charge. The negative charges are localized on the phosphate groups. Because of the large electrostatic repulsion between the phosphate groups, RNA can only fold if the effective charge on the phosphate group is reduced or renormalized (Sun17ARB, Tao18BJ, Pollack11Biopolymers, Kirmizialtin12BJ, Denesyuk13JPCB), which we know occurs by the counterion condensation (CIC) mechanism (Oosawa57JPolymerSci, Manning69JCP1, Oosawa1971Book, Alexander84JCP). In this process, described decades ago in remarkable studies by Oosawa and Manning (Oosawa57JPolymerSci, Manning69JCP1), the condensed ions are localized in an apparent volume occupied by the PE or any highly charged macroion. As a consequence, the effective charge on the negatively charged moieties in the PE (or the phosphate groups in DNA and RNA) is reduced from $-e$ to $-\beta e$ with $\beta<1$ (Oosawa57JPolymerSci, Manning69JCP1, Oosawa1971Book). The renormalized charge parameter, $\beta$, can be calculated for PEs with regular shapes (Oosawa1971Book). For example, the condition for ion condensation for a rod-like macroion is determined by the Manning parameter $\xi=\ell_{\mathrm{B}}/b>\nicefrac{{1}}{{\lvert z\rvert}}$, where $\ell_{\mathrm{B}}$ is the Bjerrum length, $b$ is the mean distance between charges on the polyelectrolyte, and $z$ is the valence of the counterion. The Bjerrum length is the distance at which the electrostatic energy of two univalent charges is equal to the thermal energy, $\ell_{\mathrm{B}}=e^{2}\left(4\pi\varepsilon_{0}\varepsilon k_{\textrm{B}}T\right)^{-1}$, where $\varepsilon_{0}$ is the vacuum permittivity, $\varepsilon$ is the dielectric constant of the solvent, and $k_{\textrm{B}}T$ is the thermal energy. For highly charged rods, ion condensation occurs when $\xi\geq\nicefrac{{1}}{{\lvert z\rvert}}$, with the transition being sharp for long thin rods (high aspect ratio). A similar condition may be obtained for spherical macroions. In this case, ion condensation does occur with a fraction localized in the vicinity of the sphere, but unlike the situation of the cylindrical macroion, there is no critical value of $\xi$ at which counterion condensation is initiated (Oosawa1971Book).

There are two implications of these well-known aspects of the CIC theory, summarized above, that are relevant to the current study. (1) Monovalent ions ($\lvert z\rvert=1$) are not as efficient in charge renormalization as higher-valent ions ($\lvert z\rvert>1$). Based on this consideration, we calculated that about 90% of the phosphate charge in Tetrahymena ribozyme is neutralized in the presence of divalent and trivalent cations (HeilmanMiller01JMBeq). (2) The condensation mechanism depends critically on the shape of the macroion. As noted above, the theoretical predictions are different for spheres and rods. Therefore, we expect that the CIC mechanism ought to be different for RNA, which undergoes substantial conformational fluctuations during the ion-driven folding process. The shape of Azoarcus ribozyme that is used as an example here, changes as the concentration of ions (Mg2+) is increased from a low to a high value. This ribozyme, and many other RNA molecules, are not spherical even when folded (Hyeon06JCP). Moreover, even at low Mg2+ concentrations, Azoarcus ribozyme does not adopt globally cylindrical structures although individual helices may be rigid enough to be pictured as stiff small cylinders.

Because of the irregular shapes of RNA molecules at all ion concentrations, ion-RNA interactions, which are often thought to involve diffuse ions that are territorially (within a volume of RNA where the electrostatic potential is substantial) but non-specifically bound, require scrutiny. Motivated by the considerations described above, we answer the following question here: How does CIC occur in such molecules whose shapes are not only irregular but also change as they fold? We are able to quantitatively answer this question because of the development of an accurate coarse-grained force field, which has been used to quantitatively predict the thermodynamic properties of Azoarcus ribozyme and other RNA molecules (Denesyuk15NatChem) in the presence of K+, Mg2+ and Ca2+.

There is an additional reason for undertaking this investigation. The physically appealing CIC mechanism qualitatively explains RNA folding both from the stability and kinetic perspective (HeilmanMiller01JMBeq). The mechanism is based on the basic assumption that condensed ions are diffuse and freely move within the volume surrounding the RNA where the electrostatic potential is the largest. Indeed, the condition for ion condensation is obtained by assuming that the mobile but bound ions are in equilibrium with free ions (the chemical potential of these two species are identical) that explore the region outside the sphere of influence of the macroion (Oosawa1971Book). Such an assumption does qualitatively explain the success of even the $z$-dependent folding kinetics of Tetrahymena ribozyme (HeilmanMiller01JMBkin, Thirumalai01AnnRevPhysChem). However, our recent studies (Denesyuk15NatChem) suggest that the assumption of uniform condensation, often made in rationalizing electrostatic effects in RNA, may be qualitatively incorrect. In other words, nucleotide-specific interactions between ions (especially Ca2+ and Mg2+ ions) provide the key in quantitatively describing charge renormalization in ribozyme folding.

In order to fully understand the nature and consequence of ion condensation onto RNA, we undertook two types of simulations. We first probed the interaction of divalent cations with small finite-sized rigid rods of different linear charge density, expressed as $b=n/L$ where $n$ is the total number of charges and $L$ is the cylinder length. It follows that $b$ is the mean distance between charges. However, in our simulations the distance between the charges ($b$) is a constant. We also investigated the condensation of divalent cations onto Azoarcus ribozyme using an accurate model introduced recently (Denesyuk15NatChem). By comparing the results of these two systems, we elaborate on the importance of the architecture of the macroion (rod versus ribozyme) on the ion condensation mechanism. Our major finding is that the shape of the RNA matters for counterion condensation. In other words, ion-condensation to RNA depends exquisitely on the folded RNA architecture, and hence cannot be treated as a non-specific interaction that is used to understand CIC for macroion with definite simple shapes. A surprise in our study is that the specificity of ion binding plays a vital role in driving RNA folding even at very low divalent ion concentrations.

Model and Methods

Rigid rod simulations

A 5 nm-long rigid cylindrical rod is placed at the center of 15 nm-wide cubic box (Fig. 1a and b). The diameter of the rod is 0.42 nm, which is the size of the phosphate bead in the coarse-grained RNA model described below. Technically, the rod is constructed by overlapping spheres that are separated by 0.1 nm. The rod is completely rigid and fixed at the center of the simulation box, whereas ions are allowed to move freely around the rod. We use periodic boundary conditions. Three different charge densities on the rod, $\rho_{{}_{\mathrm{1D}}}=$ $-1$, $-2$, and $-5$ $e\,\textrm{nm}^{-1}$ were simulated, which translates to the constant separation distance ($b$) between neighboring charges as 1.0, 0.5, and 0.2 nm, respectively. For comparison, the Bjerrum length ($\ell_{\mathrm{B}}$) is 0.73 nm at the simulation temperature, $T=37^{\circ}\textrm{C}$. The interactions between the charges on the rod and the ions are the same as in the coarse-grained RNA-ion model (Denesyuk15NatChem). Ions (K+, Cl*-*, and Mg2+) are modeled as spheres and interact with each other and with the charged rod via the excluded volume potential, which is identical to the one used in the RNA model,

[TABLE]

where $r$ is the distance of the two beads, $D_{ij}=R_{i}+R_{j}$, and $\varepsilon_{ij}=\sqrt{\varepsilon_{i}\varepsilon_{j}}$. The values of $R_{i}$ and $\varepsilon_{i}$ for the ions and the RNA sites are tabulated in Table S1 in the Supplementary Information (SI). The $R_{i}$ and $\varepsilon$ for the charged rod are the same as the phosphate site in the RNA model. If $r_{ij}>D_{ij}$ then $U_{\textrm{EX}}=0$. For any pair of charged beads, the standard Coulomb potential is used,

[TABLE]

where $q_{i}$ is the charge of the $i$th bead. The dielectric constant, $\varepsilon$, is temperature dependent and is 74 at the simulation temperature $T=37^{\circ}\mathrm{C}$. A series of Langevin dynamics simulations were performed for $2\times 10^{7}$ iteration time steps, in which ion positions are saved every 100 steps. Thus, in total $2\times 10^{6}$ snapshots are used for analyses for each setup.

We ought to emphasize that the results of rigid rod simulations for different purposes have been reported numerous times in the literature (see for example (GronbechJensen97PRL, Borukhov97PRL, Deserno00Macromolecules)) although new phenomena continue to be found in (Cha18PRL). The simulations presented here merely serve as a control in order to illustrate the dramatically different findings for the Azoarcus ribozyme.

Simulations of group I intron RNA

Because the model and simulation methods of the ribozyme were described in detail previously (Denesyuk15NatChem), we give only a brief description here. We used a coarse-grained Three Interaction Site (TIS) RNA model that has three interaction sites, for each nucleotide, corresponding to phosphate, sugar and base moieties (Hyeon05PNAS). All the ions are explicitly treated whereas water is modeled implicitly using the temperature-dependent dielectric constant. We considered concentration of Mg2+ ions in the range from 0 to 30 mM in the presence of either 50 mM or 12 mM of K+ that is typically used in the Tris buffer (Behrouzi12Cell). We added an appropriate number of anions ($\textrm{Cl}^{-}$) to ensure that the entire system is charge neutral. All the simulations were conducted using a periodic-boundary cubic box with each side being 35 nm. We performed a series of Langevin dynamics simulations at 37${}^{\circ}\textrm{C}$ in order to efficiently sample conformations of the system containing the ribozyme and ions (Honeycutt92Biopolymers). For each salt concentration, several independent simulations were started from the crystal structure and trajectories were generated for sufficiently long time to ensure that equilibrium is reached before we obtained snapshots for analyses. For low Mg2+/Ca2+ concentrations, at which the ribozyme is unfolded state, we generated seven trajectories with 7.5 $\mu$s of equilibration followed by 15–22.5 $\mu$s of production runs. The length of these simulations are long enough to observe multiple unfolding/folding events in the low-friction Langevin dynamics. For higher Mg2+/Ca2+ concentrations, at which the ribozyme is predominantly in the folded conformations, we ran three trajectories with length 7.5 $\mu$s and discarded the first $10^{5}$ steps of each trajectory. Elsewhere (Denesyuk15NatChem) we provide more detailed descriptions of the model and the simulation protocol.

The crystal structure of Azoarcus group I intron (PDB 1U6B (AzoXray04N), Fig. 1c and d) was used as the reference structure of the native conformation. The numbering of nucleotides from G12 through G206 is consistent with the crystal structure  (AzoXray04N) and the first nucleotide of the exon sequence is labeled 207.

Data analyses

To quantify ion condensation onto the RNA and the rigid rod, we calculated the local ion concentration, a fingerprint of ion interaction with the macroion, around the $i$th phosphate (or point charges on the rod) defined as,

[TABLE]

where $\rho_{i}(r)$ is the number of ions per unit volume (number density) of the ion at the distance $r$ from the $i$th phosphate, $V_{c}$ is the spherical volume of radius $r_{c}$, $N_{\mathrm{A}}$ is the Avogadro’s number to represent $c_{i}^{\ast}$ in molar units. We use the Bjerrum length as the cutoff distance, $r_{\textrm{c}}=\ell_{\mathrm{B}}=0.73\,\textrm{nm}$ because at that distance the favorable Coulomb attraction between a cation and an anion ($z=1$) is $k_{\mathrm{B}}T$.

We also calculated the fluctuations in the local ion concentration as the normalized root-mean-square deviation from the mean,

[TABLE]

where $\left\langle c_{i}^{\ast}\right\rangle$ is the ensemble average of the local ion concentration around the $i$th nucleotide. Thus, $c_{i}^{\ast}-\left\langle c_{i}^{\ast}\right\rangle$ is an instantaneous deviation from the ensemble average. The average ion concentration is,

[TABLE]

where $N$ is the number of nucleotides. Note that the upper bar indicates an average over the sequence, whereas the bracket indicates the ensemble average. The extent of deviation of ion concentration fluctuation from $\left\langle c^{\ast}\right\rangle$ along the sequence was calculated using,

[TABLE]

In the results section, for simplicity, we omit the brackets for the ensemble average. Thus, all quantities are averaged over an ensemble of configurations in the system (ions and the charged RNA or the rod) unless otherwise stated.

For charged rods, we also computed the total charge per unit length due to contributions from each ionic species. The effective charge densities around the rod were calculated as the amount of accumulated charge per unit length for the ion species $m$,

[TABLE]

where $\rho_{m}(R)$ is the number density of species $m$ at the radial distance $R$ from the axis of the cylinder, and $z_{m}$ is its valence. Note that $Q$ can be computed for each ion Mg2+, K+, and Cl*-* independently. Finally, the total effective charge per unit length is the sum of the rod charge and contributions from those ions,

[TABLE]

Note that we use $R$ for the radial distance from the center of the rods in the cylindrical coordinate, whereas $r$ is for the distance from the point charges on the rod in the spherical coordinate (see Fig. 1a).

Results

Charged cylindrical rods

We first describe condensations of ions onto charged rods with varying spacing between fixed charges, $b=-e/\rho_{{}_{\mathrm{1D}}}$, where $\rho_{{}_{\mathrm{1D}}}$ is the charge density of the rod. Based on the CIC theory, we expect that ions should condense uniformly along the rod. The condition for condensation in the CIC theory is that $\xi>\nicefrac{{1}}{{\lvert z\rvert}}$, where $z$ is the valence of the ion species. In the absence of multivalent salt, $\xi>1$ and thus we expect ions should condense on rods when $\rho_{{}_{\mathrm{1D}}}=$ $-5\,e/\mathrm{nm}$ ($\xi=3.6$) and $\rho_{{}_{\mathrm{1D}}}=$ $-2\,e/\mathrm{nm}$ ($\xi=1.5$) among the three values of charge densities we examine here. This expectation was borne out in the simulations even though the aspect ratio ($\approx$ 11.9) of the rod was relatively small. In Fig. 2(a-c), the cumulative charges, $Q$ (Eq. 7), around the rod are plotted for each ion species and the net charge. The condensation of monovalent ions appears as a sharp increase in $Q_{\mathrm{K}^{+}}$ around the vicinity of the rod ($R\sim 5\,\mathrm{nm}$), with a subsequent reduction of the net negative charge at longer distances (thick lines in gray). The amount of condensation is larger for the rod with $\rho_{{}_{\mathrm{1D}}}=$ $-5\,e/\mathrm{nm}$ than for one with $\rho_{{}_{\mathrm{1D}}}=$ $-2\,e/\mathrm{nm}$. In the case of the lowest charge density ($-1$ $e/\mathrm{nm}$ corresponding $\xi=0.7$ ), condensation is not observed as predicted by the CIC theory (Fig. 2c).

The condition for condensation, $\xi>\nicefrac{{1}}{{\lvert z\rvert}}$, also implies that the charged rod must be surrounded to a greater degree by $\textrm{Mg}^{2+}$ than $\textrm{K}^{+}$. This prediction was also borne out in the simulations. In Fig. 2 (d and e), Mg2+ ions condense in the immediate vicinity of the rods taking over the role of K+ in the charge renormalization. The replacement of K+ by Mg2+ in the vicinity of the rod is expected based on the counterion release mechanism (leads to an entropy of the whole system increase), which is another well-known consequence of the CIC theory. The condensation of cations, and the subsequent net reduction of negative charge on the rod, are more efficient in the presence of Mg2+ than when $z=1$ (compare panels (a) and (d) in Fig. 2). Since the condensation condition for Mg2+ is $\xi>0.5$, a distinct increase of $Q_{\mathrm{Mg}^{2+}}$ is visible even at the lowest density case, $\rho_{{}_{\mathrm{1D}}}=$ $-1$ $e/\mathrm{nm}$ ($\xi=0.7$ ) albeit the extent of condensation is not as great as found in the two rod systems with higher density (Fig. 2f).

The results quoted above using the control simulations are entirely in accord with the CIC theory predictions. The extremely close correspondence between the CIC theory and simulations for charged rods motivated us to make several assessments of ion condensation phenomena in ribozyme folding, in which we use the same explicit-ion model. In order to make quantitative comparisons between the rod and RNA folding simulations, we need to devise a quantity to precisely measure the ion condensation. The radial distributions in cylindrical coordinates are not suitable for flexible RNA chains, for which the persistence length (1 $\sim$ 2 nm for a single strand) and the Debye screening length (1.3 nm in the simulation conditions) are similar. With this in mind, we defined local ion concentrations ($c^{\ast}$) and the associated fluctuations (Eqs. 3-4) in the spherical coordinate for each point charge. The results of the rod simulations in terms of $c^{\ast}$ are shown in Fig. 3 (a and b) for the three rod charge densities in the presence of 5 mM MgCl2 and 50 mM KCl. It is clear that the density of condensed ions around the charges in the rod is uniform; there is no position dependence along the rod. From panels (c) and (d) in Fig. 3, the fluctuation in $c^{\ast}$ is the lowest in the case of Mg2+ condensed for the highest charge density rod. Both lowering the linear charge density and reducing the valence of cations (Mg2+ $\rightarrow$ K+) make the condensed state less favorable, thus increasing the magnitude of charge fluctuations. These quantitative descriptions for the charged rods set the stage for analyzing the completely contrasting behavior in ion-ribozyme interactions.

[Mg2+]-dependent folding of Azoarcus ribozyme

Before proceeding to examine the roles of ions in RNA folding, we summarize briefly the [Mg2+]-dependent folding of the Azoarcus ribozyme, which has been investigated in a previous study using the elaborate version of the TIS model (Denesyuk15NatChem). The folding process of the ribozyme (Fig. 1) strongly depends on the concentration of $\textrm{Mg}^{2+}$. In Fig. 4a, the dependence of the radius of gyration ($R_{\textrm{g}}$) on [Mg2+] is shown in gray lines for two different conditions of monovalent salt concentration, [K+] = 12 and 50 mM. The global compaction of the RNA appears as a single transition in $R_{\textrm{g}}$ occurring at $C_{\textrm{m}}\approx 1.5$ mM, where we designate $C_{\textrm{m}}$ as the midpoint concentration of [Mg2+] separating the folded and unfolded states. Comparison of 12 mM K+ simulation data and SAXS experiments at the equivalent condition (Behrouzi12Cell) (squares in Fig. 4) shows excellent agreement, thus establishing that our model captures monovalent- and divalent-ion dependence of folding of Azoarcus ribozyme. In the absence of Mg2+, the key tertiary interactions (marked in red in Fig. 1c) are disrupted and only some of the helical domains, P2, P4, P5, and P8 are stably formed. As [Mg2+] increases, other helices and tertiary interactions form in multiple steps. At [Mg2+] $\sim 1$ mM, tertiary interactions, TH, G site, and SE form, followed by folding of other peripheral interactions, TL2-TR8 and TL9-TR5, at [Mg2+] $>\,2$ mM. The formations of these tertiary interactions lead to a compact intermediate state with $\left.\langle R_{\textrm{g}}\rangle\right.\sim 3.5$ nm. Some interesting cooperativity and anti-cooperativity, such as cooperative formation of P3 helix and SE, were found and discussed in the previous study (Denesyuk15NatChem). The complete native structure, including the active site, is spontaneously formed at [Mg2+] $>\,4$ mM. In subsequent parts of this paper, we report detailed analyses of ion distributions around the RNA. Comparisons with the interaction of divalent ions with highly charged rods illustrate the importance of ribozyme architecture on ion–RNA interactions.

Radial profiles of ion condensation

We calculated the density profiles of each ion species around phosphate, $\rho_{i}(r)$, where $r$ is the distance from the phosphate of the $i$th nucleotide. We also computed the profiles of phosphate–phosphate distance, which can be used to trace the compaction of the ribozyme. The ensemble-averaged profiles, for given Mg2+ concentrations, are shown in Fig. 5. From the profiles of the phosphate (gray dotted lines in Fig. 5), we can track the folding of the ribozyme as follows. The peak at $r\approx 0.55$ nm is roughly at the distance between two consecutive phosphates, and the smaller second peak corresponds to the location of the second nearest-neighbor phosphates. The amplitudes of these two peaks do not show significant change as [Mg2+] increases. However, the folding process can be discerned by the increase in $\rho_{{}_{\mathrm{P}}}(r)$ in the range of 0.7 $<r<$ 0.9 nm (red arrows in Fig. 5). At low [Mg2+], there is practically no density in $\rho$ in this range ($\rho(r)\approx 0$). As [Mg2+] increases, we find that the amplitude increases, which shows that Mg2+ ions facilitate the approach of non-neighboring phosphates. The increase in $\rho(r)$ in this range coincides with the compaction of the ribozyme and subsequent folding as shown in Fig. 4(a).

In the absence of Mg2+, K+ cations are condensed around the ribozyme, regardless of the monovalent salt concentration, [KCl] = 12 mM or 50 mM (Figs. 5a and S1a). Under these conditions with [Mg2+] $=0$, all tertiary interactions are disrupted but helices, except P3 and P7, form secondary structures (Fig. 1c). Thus, K+ cations must be condensed around the vicinity of phosphates to balance the electrostatic repulsion between negative charges on the phosphates. The distance to condensed K+, $r\approx 0.4$ nm, is closer than the neighboring phosphate. In the presence of Mg2+, $\rho(r)$ for K+ decreases dramatically as the Mg2+ concentration increases (Fig. 5 b-d). This shows, rather vividly, that K+ ions condensed around RNA are replaced with $\textrm{Mg}^{2+}$ ions because divalent cations screen negative charges on the RNA more efficiently than monovalent cations. By replacing two K+ ions by one $\mbox{Mg}^{2+}$, there is a net gain in translational entropy of the ions. This is the counterion release mechanism predicted by the CIC theory.

Fingerprint of local ion concentration $c^{\ast}_{i}$

In order to examine the relationship between ion condensation and Mg2+-driven RNA folding, we calculated the local ion concentration of the $i$th nucleotide $c_{i}^{\ast}$ (Eq. 3). The cutoff distance to define the local volume is chosen to be the Bjerrum length, $r_{c}=0.73$ nm. The range of this distance covers the peaks in the radial distributions of both K+ and Mg2+ (Fig. 5). Accordingly, $c_{i}^{\ast}$ is expected to capture the amount of ion atmosphere around each nucleotide, and such a fingerprint is dependent on the architecture of the folded state of the RNA. The overall ion condensation of the entire molecule, $\overline{c^{\ast}}$, is calculated as the average of the local ion concentration along the RNA sequence (Eq. 5, see Models and Methods).

The changes in $\overline{c^{\ast}}$ as a function of [Mg2+] is shown in Fig. 4(b). As noted previously, more $\textrm{K}^{+}$ cations condense at lower Mg2+ than at higher Mg2+ concentrations. Indeed, $\overline{c_{\textrm{K}^{+}}^{\ast}}$ monotonically decreases as [Mg2+] increases (green lines). The folding midpoint of Mg2+ concentration is around 1.5 mM. Above the midpoint, Mg2+ ions are more dominantly found in the vicinity of the RNA than K+. $\overline{c_{\textrm{Mg}^{2+}}^{\ast}}$ does not vary significantly as a function of K+ concentrations for the entire range of [Mg2+] (Fig. 4b, solid line for [K+] = 50 mM and dotted line for 12 mM). The sharp increase in $\overline{c_{\textrm{Mg}^{2+}}^{\ast}}$ is well correlated with the global compaction of the ribozyme, that is evident as the decrease in $R_{\mathrm{g}}$ in Fig. 4(a).

Ion condensation is site-specific

In Fig. 4c, we show variances, $\sigma^{2}$, of the local ion concentration over the RNA sequence (Eq. 6). This quantity indicates how much of each ion species are localized at particular positions along the sequence at each solution condition. As Mg2+ concentration increases, $\sigma^{2}$ for K+monotonically decreases, although it is much smaller compared to divalent cations for the entire range of [Mg2+]. This means that Mg2+ ions tend to be localized at specific positions more than K+ ions. In particular, $\sigma^{2}$ for Mg2+ sharply increases reaching a maximum around [Mg2+] $\approx 3$ mM. This is because Mg2+ ions are preferentially localized to specific places, which are, surprisingly, determined by the tertiary structure of the ribozyme. Beyond the midpoint $C_{\mathrm{m}}$, the value of $\sigma^{2}$ decreases modestly because the specific sites of higher phosphate densities are occupied by existing localized Mg2+. The newly added Mg2+ ions start contributing to lessening the gap between specific and other non-specific sites.

Fig. 6 shows the distributions of the local ion concentration along the RNA sequence. This figure shows that all non-zero values of Mg2+ condensation is highly specific. The extent of condensation depends on the nucleotide. Even at the lowest [Mg2+] $=0.4$ mM (much lower than $C_{\mathrm{m}}$) there is a significant density of Mg2+ around $i\approx 125$, which covers the triple helix (TH) of Azoarcus ribozyme. Thus, we conclude that Mg2+ binding is highly specific. More importantly, charge neutralization resulting in $\beta<1$, occurs at the lowest [Mg2+] not uniformly, as predicted by the CIC theory, but does so in a manner that initiates folding. This would mean that different regions of this ribozyme order at somewhat distinct Mg2+ concentrations, as illustrated previously (see Fig. 2b in (Denesyuk15NatChem)). In other words, the midpoint concentration, $C_{\mathrm{m}}$, of the transition measured by $R_{\textrm{g}}$ or other parameters is the (appropriately-weighted) mean of such individual [Mg2+] at which distinct structural elements achieve the native conformation.

In the absence of or at low Mg2+ concentrations, the amount of condensed K+ depends significantly on positions along the RNA sequence. At low Mg2+ concentrations, the ribozyme is unfolded with some helices intact. To see if there is any relationship between K+ condensation and such residual structure, we compare the distribution of $c^{\ast}$ with the fraction of secondary structure formation in Fig. 7(a). The profile of condensed K+ (green in Fig. 7a) has several characteristic minima that correspond to loops in the stem-loop structures of P2, P5, P6, P8, and P9. This indicates that the distribution of K+ reflects the shape of RNA that is mostly due to residual secondary structures at low [Mg2+].

At high Mg2+ concentration ([Mg2+] $>1.5$ mM), the localization of condensed Mg2+ ions depends on the tertiary structures of the RNA (Fig. 6c and d). The local concentration of Mg2+ correlates with the local density of phosphates (compare blue and gray lines in Fig. 7b). Under such conditions, all the native secondary structures are formed with almost unit probability, but we do not find any secondary-structure dependence of condensed Mg2+.

The relationship between ion condensation and local density of phosphates is vividly illustrated in Fig. 8(a, b). At low [Mg2+], the ribozyme is mostly unfolded, predominantly containing only secondary structures. Because there are no stable tertiary interactions, $c_{\mathrm{P}}^{\ast}$ values are small (Fig. 8a, horizontal axis. $c_{\mathrm{P}}^{\ast}\approx 3\,\textrm{M}$). Accordingly, $c^{\ast}$ for Mg2+ are relatively low and K+ are predominantly condensed. However, there are six nucleotides that have higher $c^{\ast}$ for both P and Mg2+ (see the rectangle of $c_{\mathrm{P}}^{\ast}>3.5\,\mathrm{M}$ and $c_{\mathrm{Mg^{2+}}}^{\ast}>0.5\,\mathrm{M}$ in Fig. 8a). We identified all those nucleotides to be ones involving the central triple helix (TH) formation; nucleotides 50, 51 and from 124 to 127. In accordance with the previous study (Denesyuk15NatChem), TH starts to form around [Mg2+] = 0.4 mM ($<C_{m}$). The data in Fig. 8 indicates that there are Mg2+ ions specifically bound to the sites around TH. In contrast, at high [Mg2+], the entire ribozyme is folded and $c_{\mathrm{P}}^{\ast}$ can be as high as \sim$$7 M depending on the local architecture (Fig. 8b). The amount of condensed Mg2+ is well correlated with the local concentration of phosphates. There is a weak anti-correlation between K+ condensation and local phosphate concentration.

Fluctuations of condensed ions

RNA folding is a stochastic process driven by thermal fluctuations and the solution conditions. Consequently, it would be interesting to see how these condensed ions fluctuate at different stages of folding. In Fig. S3, normalized fluctuations of the local ion concentrations, $\Delta c^{\ast}$ (Eq. 4), for each nucleotide are shown for the several solution conditions. Below or around the midpoint Mg2+ concentration, where the ribozyme folding is not completed, fluctuations of condensed K+ ions and Mg2+ have similar position dependence. This correlation is more clearly seen in the logarithm plot of $\Delta c_{\mathrm{Mg}}^{\ast}$ vs $\Delta c_{\mathrm{K}}^{\ast}$ in Fig. 8(c). Comparing Figs. 6 and S3, we also find that nucleotide positions that have high $c^{\ast}$ values also have small $\Delta c^{\ast}$. This indicates that the specific positions, where more ions are condensed, bind those ions also more tightly. Once the ions reach the “equilibrium” positions, they tend to be localized.

The correlation between the fluctuations of $\Delta c_{\mathrm{Mg}}^{\ast}$ and $\Delta c_{\mathrm{K}}^{\ast}$ changes dramatically after the midpoint, [Mg2+] $>1.5$ mM (Fig. S3d and e). There is anti-correlation between $\log\Delta c_{\mathrm{Mg}}^{\ast}$ and $\log\Delta c_{\mathrm{K}}^{\ast}$ (Fig. 8d). As discussed in previous sections, there are specific positions where more Mg2+ ions are condensed (Fig. 6) beyond the midpoint of [Mg2+]. At these positions, there are smaller fluctuations in $c_{\mathrm{Mg}}^{\ast}$ (smaller $\Delta c_{\mathrm{Mg}}^{\ast}$) but higher fluctuations in $c_{\mathrm{K}}^{\ast}$ (larger $\Delta c_{\mathrm{K}}^{\ast}$), reflecting strong condensation of Mg2+.

Ion condensation in Calcium-driven folding

We have hitherto focused on the condensation of K+ in the absence and presence of Mg2+, and the condensation of Mg2+ itself. It is known that Ca2+ could also facilitate $Azoarcus$ ribozyme folding although the folded state is not catalytically active and is modestly less compact than it is in the presence of Mg2+ (Behrouzi12Cell, Denesyuk15NatChem). We repeated the same analyses in the case of Ca2+-dependent folding instead of Mg2+. Both average local concentration ($\overline{c^{\ast}}$) and its variance ($\sigma^{2}$) show similar trends as the case of Mg2+ (Fig. 4b and c). At higher solution concentration ([Mg/Ca2+] $\gtrsim 1$ mM), the amount of condensation given by $\overline{c^{\ast}}$ is slightly less in the case of Ca2+.

Fig. S4 shows the sequence dependence of the local concentration ($c^{\ast}$) and fluctuation ($\Delta c^{\ast}$) in the case of Ca2+. Comparing them with Fig. 6 and S2, the condensation profiles of Ca2+ and K+ bear remarkable resemblances to the Mg2+ case. At the highest concentration ([Mg/Ca2+] $=30$ mM), the local concentration of Ca2+ has the same peak positions as Mg2+, but has $\sim$20 % lower values. Comparison of the fingerprint of ion concentration profiles between Mg2+ and Ca2+ clearly illustrates the efficacy of Mg2+ of creating a compact and functionally competent ribozyme.

Discussion

The central result of this study is that the nature of ion condensation onto Azoarcus ribozyme differs drastically from theoretical predictions for regularly shaped polyanions. For example, the CIC theory shows that ions condense uniformly onto highly charged rods. The positively charged ions are localized at the positions of the charges on the rod. In sharp contrast, in Azoarcus ribozyme, and presumably in other RNA molecules as well, the condensation is highly non-uniform and nucleotide specific, and critically depends on the structure of the RNA. We note in passing that the prevailing view in much of the RNA literature is that most of the ions are uniformly condensed onto RNA although crystal structures have shown that there are site-specific localizations of Mg2+. A natural explanation of our finding is that the electrostatic potential around the phosphate groups is non-uniform, with the folded structure with fully formed tertiary interactions containing pockets of highly negative potentials. Consequently, ions are specifically drawn, with higher probability, to regions of large negative electrostatic potential with smaller ions approaching these regions easily. Indeed, this was first shown in our previous study using coarse-grained and atomic-detailed simulations (see Fig. S11 in (Denesyuk15NatChem)). However, what is surprising in this study is that the ions sense these regions even at extremely low divalent ion concentrations when there is no tertiary interaction present.

Link to the shape fluctuations in the ribozyme

The non-uniform and site-specific ion-RNA interactions can be succinctly illustrated by considering the ion concentration-dependent changes in the shape parameters, asphericity $\Delta$ ($0\leq\Delta\leq 1$) and prolateness $S$ ($-0.25\leq S\leq 2$) (Aronovitz86JP, Dima04JPCB), which are readily computed from the inertia tensor that is related to the square of the radius of gyration. For reference, it is worth pointing out both $\Delta$ and $S$ are unity for rods, whereas $\Delta=S=0$ for spheres. For long self-avoiding polymers, numerical results show that $\Delta=0.55$ and $S=0.92$ (a prolate ellipsoid) (Honeycutt89JCP). In contrast to the values of the shape parameters for rods, the changes in $\left<\Delta\right>$ and $\left<S\right>$ as a function of [Mg2+] for the ribozyme, at the two values of K+ concentration, exhibit unusual non-monotonic dependence as [Mg2+] is varied (inset in Fig. 9), which is not reflected in $R_{\mathrm{g}}$. Fig. 9 also shows the distributions $P(\Delta)$ and $P(S)$ at three Mg2+ concentrations. Surprisingly, the widths of $P(\Delta)$ and $P(S)$, which depend non-monotonically on [Mg2+], are relatively large. Therefore, shape fluctuations in the ribozyme are substantial. Two comments are worth making: (1) The values of $\Delta$ and $S$ depend greatly on the K+ concentration, which is also reflected in the $R_{\mathrm{g}}$ as a function of Mg2+ concentration. (2) It is is interesting that the minima in $\left<\Delta\right>$ and $\left<S\right>$ are at [Mg2+] $\approx 1.5$ mM, the midpoint of the folding transition, where the ribozyme almost resembles a sphere. At high and low Mg2+ concentrations, the values of $\left<\Delta\right>$ show that globally the ribozyme may be pictured as a prolate ellipsoid. It also implies that this ribozyme and others are not as densely packed as proteins (Hyeon06JCP). The irregular shape leads to the high specificity of ion–RNA interactions. It is likely that this finding would apply to any RNA molecules whose shape has grooves and whose surface are irregular. It is most interesting that there is a calculable fingerprint for ion-RNA interactions that is a reflection of the folded structure. A clear implication is that, even at very low ion concentration, binding is highly specific. This interpretation differs from the usual assertion that the stability of RNA is determined largely by diffuse divalent cations (Soto07Biochem).

Correlation between ions at different RNA sites

We computed the correlations of Mg2+/K+ binding to two different nucleotides (Fig. S5). The correlations for a pair of nucleotides ($i$, $j$) was calculated as the Pearson’s coefficient of the number of Mg2+/K+ ions bound to site $i$ and site $j$ for the ensemble at a given Mg2+ and K+ concentrations. Overall, we did not find strong correlations involving Mg2+ or K+, except that some signatures appear as consequences of formations of secondary and tertiary structures. For example, there are strong correlations near the diagonal region because of helix domains P5, P6, and P8, as indicated by the red circles (Fig. S5a). In addition, there are high correlations at the interfaces of P3/P7 (see Figure 4 of Ref. (Denesyuk15NatChem)), and P2/P8. The observed high correlations arise because these nucleotides are close in space and share the same Mg2+ ion when the secondary/tertiary interactions form. Interestingly, at the midpoint concentration [Mg2+] $\sim 1.5$ mM, there are weak anti-correlations between P2 and P6a, and P2 and J8/7 (blue circles). Helix P2 is located at a peripheral part of the ribozyme, whereas P6a and J7/8 are close to the core region. These anti-correlations are a consequence of anti-cooperative folding of tertiary interactions, arising from topological frustration (Denesyuk15NatChem). Around the midpoint concentration of [Mg2+], there are still not enough amount of Mg2+ ions to simultaneously satisfy all the folding elements. Thus, Mg2+ binding to P2 reduces the probability of Mg2+ binding to the core region. The correlation maps between K+-K+ ions are less instructive.

Flexible polyelectrolytes

Just as is the case for RNA, which does adopt a well-defined folded structure at high ion concentrations, the shapes of flexible PEs are also highly variable. Therefore, it follows that the CIC theory should not be quantitative in predicting the nature of ion condensation. Indeed, early simulations showed (Lee01Macromolecules) that ion condensation is non-uniform in the presence of divalent and trivalent cations for PEs in poor solvents. Indeed, in certain cases, the divalent cations and the monomers are packed as in a crystalline arrangement, which would be inconsistent with the CIC theory. Using theory, Muthukumar (Muthukumar04JCP) has shown that, even in monovalent salts ($z$=1), the condensation process in flexible PEs is non-uniform and could be understood by the variations in the local dielectric mismatches. A more recent simulation study (Chremos16SoftMatter) has also arrived at the conclusion that the condensation of monovalent ions depends crucially on the shape of the PE. In particular, they showed that non-uniform condensation is determined by the molecular topology. Thus, the conclusion that shapes determine ion condensation holds for both synthetic PEs and RNA. The latter has additional features (for example stacking interactions and hydrogen bonding), which make RNA different from a synthetic. In particular, the ion-induced shape changes in PEs are not as spectacular as they are in RNA.

Conclusion

The Oosawa-Manning theory of ion condensation, which assumes that there is an equilibrium between bound and free ions, is the basis of our understanding of many thermodynamics properties of highly charged polyelectrolytes. As shown multiple times previously and in the simulations here, the predictions of the CIC theory are quantitatively correct if the charged polyelectrolytes adopt regular shapes, such as rods and spheres. This is not the case in RNA, which although is a polyanion, does not adopt a simple shape but undergoes large shape transitions as it folds. Our simulations of Azoarcus ribozyme in the presence of Mg2+ or Ca2+ in a buffer containing K+ ions show that ion binding is highly specific. Surprisingly, the specificity of binding depends on the folded structure of the RNA. In other words, in this ribozyme with a complex fold, the charge neutralization of phosphates does not occur uniformly as predicted by the CIC theory for rod-like macroions. There are specific regions where condensation occurs even at very low concentrations of the divalent cations. Although we have arrived at this unexpected finding using simulations of one ribozyme, we expect that the main conclusion should be applicable for other RNA molecules as well.

We conclude with the following additional comments. (1) We find, perhaps not surprisingly, that monovalent ions unbind (termed as counterion release process) when RNA is titrated with divalent ions (Fig. 6). This validates one of the predictions of the CIC theory. The release of K+, however, occurs from specific nucleotide locations on the RNA, which is not anticipated in the CIC theory. (2) It is unclear to us if our findings, especially the exquisite specificity of the predicted nucleotide-dependent binding of ions to RNA, could be tested in experiments because of the highly correlated and many-body nature of ion condensation. Multicolor FRET experiments may be used to label regions of high affinity for divalent cations. The conformational changes as a consequence of ion binding and unbinding at these locations could be used, most likely in conjunction with simulations, to determine the architecture-dependent interaction of ions with RNA. (3) Our results show that charge renormalization of phosphates occurs in a highly heterogeneous manner and not by the diffuse binding picture of ion-RNA interactions. The fingerprints of ion localizations show that the extent of charge neutralization depends on the nucleotide position – finding that has to be incorporated in calculating free energy changes of RNA as it folds. (4) The clear implication of our results is that, even when using coarse-grained models, ions (especially those with valence greater than unity) should be modeled explicitly to reflect their sizes and charge densities. It is possible that electrostatic interactions due to monovalent ions can be accounted for by using effective interactions (Denesyuk18JPCB). However, in order to obtain accurate thermodynamics of RNA folding, divalent ions have to be explicitly modeled (Denesyuk15NatChem, Hayes15PRL).

Author Contributions

NH, NAD, and DT designed the research; NH and NAD performed all the simulations; NH, NAD, and DT analyzed the data; NH and DT wrote the article.

Acknowledgments

We thank Dr. Hung T. Nguyen for useful discussions. This work was supported by a grant from the National Science Foundation (CHE 16-36424) and the Welch Foundation (F-0019) through the Collie–Welch Chair.

Supplementary Material

Figure S1: Radial profiles of local ion concentration, $\rho(r)$, for $\textrm{Mg}^{2+}$, $\textrm{K}^{+}$, $\textrm{Cl}^{-}$, and phosphate (P) at [K+] = 12 mM and various [Mg2+]: (a) 0, (b) 0.4 mM, (c) 1.5 mM, (d) 4.0 mM, and (e) 30 mM. (f) The profile at [K+] = 50 mM and [Mg2+] $=$ 30 mM. Appropriate amount of $\textrm{Cl}^{-}$ were added to the solution for the electroneutrality. Note that the scale for $\rho_{{}_{\text{Mg}^{2+}}}$ and $\rho_{{}_{\textrm{P}}}$ (left axis) is 10 times greater than $\rho_{{}_{\textrm{K}^{+}}}$ and $\rho_{{}_{\textrm{Cl}^{-}}}$ (right axis).

Figure S2: Local ion concentrations, $c_{i}^{\ast}$ ($i$ is the nucleotide position), for $\textrm{Mg}^{2+}$ (blue) and $\textrm{K}^{+}$ (green) at [K+] = 12 mM and various [Mg2+]: (a) 0, (b) 0.4 mM, (c) 1.5 mM, (d) 4.0 mM, and (e) 30 mM. (f) The profile at [K+] = 50 mM and [Mg2+] $=$ 30 mM.

Figure S3: Normalized fluctuations of local ion concentrations, $\Delta c^{\ast}$ (Eq. 4), for Mg2+ (dotted blue) and K+ (solid green), in 12 mM KCl (left column; a–e) and 50 mM KCl (right column; f–j).

Figure S4: Sequence dependence of the local ion concentration (left column, $c^{\ast}$) and its fluctuations (right column, $\Delta c^{\ast}$) in the presence of Ca2+ ions instead of Mg2+. The concentration of Ca2+ is labeled in each panel. [KCl] is fixed at 50 mM.

Figure S5: (a) Correlations between Mg2+ condensations to two different nucleotides. The horizontal and vertical axes are the nucleotide positions $i$ and $j$. Color in each cell ($i$, $j$) represents the Pearson correlation coefficient (scale on the right) of the number of bound ions at site $i$ and $j$. The Mg2+ concentration is 1.5 mM and [K+] = 50 mM. (b) Correlations between Mg2+-K+ condensations is most prominent along the diagonal. This implies that K+ ion is released upon Mg2+ binding.