Molecular Dynamics Simulations Reveal an Interplay between SHAPE Reagent Binding and RNA Flexibility.

The function of RNA molecules usually depends on their overall fold and on the presence of specific structural motifs. Chemical probing methods are routinely used in combination with nearest-neighbor models to determine RNA secondary structure. Among the available methods, SHAPE is relevant due to its capability to probe all RNA nucleotides and the possibility to be used in vivo. However, the structural determinants for SHAPE reactivity and its mechanism of reaction are still unclear. Here molecular dynamics simulations and enhanced sampling techniques are used to predict the accessibility of nucleotide analogs and larger RNA structural motifs to SHAPE reagents. We show that local RNA reconformations are crucial in allowing reagents to reach the 2'-OH group of a particular nucleotide and that sugar pucker is a major structural factor influencing SHAPE reactivity.

Ribonucleic acid (RNA) is involved in a large number of functional processes in cells.[1] RNA function depends not only on its sequence and secondary structure but also on its specific 3D fold and structural dynamics.[2,3] Thus there is a longstanding interest in developing tools that can be used to predict or directly determine RNA structure. Chemical probing techniques are popular methods used to determine RNA secondary structure by distinguishing among base-paired and unpaired nucleotides. Hydroxyl radical footprinting, dimethylsulfate probing, selective 2′-hydroxyl acylation analyzed by primer extension (SHAPE), and in-line probing are used most frequently.[4] The number of available data sets obtained by chemical probing techniques has exploded in the last years and has already surpassed by two orders of magnitude the number of RNA 3D structures obtained with traditional X-ray or NMR methods.[5] Statistics from the freely available RMDB database[6] shows that almost half of the chemical probing data sets were obtained by using SHAPE protocols. SHAPE was designed as a method providing a sequence-independent and time-effective analysis of RNAs with single nucleotide resolution.[7,8] Chemically, 2′-hydroxyl (2′-OH) groups are modified by electrophilic molecules (SHAPE reagents) during the acylation reaction, resulting in adducts that are subsequently detected as stops in primer extension by gel electrophoresis or with mutational profiling.[9−11] Recent protocols allow one to analyze massive structures,[12−14] have success in probing RNA in vivo,[15−21] and can also detect motifs participating in RNA–protein interactions.[13,14] Additionally, SHAPE data serve in guiding structure models and improving the accuracy of popular secondary structure prediction tools[22−30] and have been fruitfully combined with a plethora of computational methods and algorithms.[9,19,28,31−33]

Flexible nucleotides are presumed to be more reactive toward SHAPE reagents because they were shown to sample multiple conformations, where few of those could enhance the reactivity of 2′-OH group.[4,8] Thus certain efforts focused on the relationship between SHAPE reactivity and specific geometry of RNA residues.[31,34−36] In particular, Bindewald and coworkers found that SHAPE reactivities are correlated with structural properties of nucleotides, that is, base pairing (cis-Watson–Crick/Watson–Crick) and base stacking.[31] McGinnis and coworkers compared reactivities of small nucleotide analogs with ∼1500 different nucleotides from the small ribosomal subunit and identified three important structural parameters: (i) sugar-pucker of the ribose ring, (ii) conformation of the adjacent phosphate moiety, and (iii) presence of a RNA functional group within hydrogen-bond (H-bond) distance from 2′-OH group.[34] More specifically, the sugar pucker of the ribose ring was identified as an important structural factor[37] and nucleotides with C2′-endo sugar pucker appeared to be more reactive toward SHAPE reagents.[34] Interestingly, transient pucker reconformations of RNA nucleotides can also be observed in NMR studies.[38] Finally, SHAPE data sets were correlated with structural fluctuations obtained from both molecular dynamics (MD) simulations and elastic network models (ENMs), where the best fits surprisingly considered neither 2′-OH groups nor phosphate moieties but rather fluctuations in the distance between consecutive nucleobases.[36] Overall, there is a general consensus that nucleotide flexibility and dynamics are important in yielding SHAPE modifications, but a microscopic interpretation of this statement is totally missing. Indeed, to the best of our knowledge, neither SHAPE reagent interactions with RNA nor atomistic details of the interplay between reagent binding and local nucleotide dynamics have been analyzed in detail.

To reveal the effects of SHAPE reagent on RNA structural dynamics, we report here a computational analysis of reagent interactions with diverse RNA systems, that is, the RNA core of the signal recognition particle (SRP–RNA), GNRA tetraloop, and two RNA nucleotide analogs. Classical as well as enhanced sampling MD simulations revealed that reagent is mostly stabilized in conformations required for the acylation reaction by stacking interaction with ribose sugar rings. Reaction rates are dependent on the dynamics of the whole RNA motif and on the ability of each nucleotide to accommodate the reagent in close proximity of its 2′-OH group.

The SRP–RNA (Figure ) containing several unpaired as well as Watson–Crick (WC) and noncanonical (NC) base-paired nucleotides was used to investigate the accessibility of nucleotides toward fast acting benzoyl cyanide (BzCN) SHAPE reagent. We mostly used BzCN in our simulations because interactions between complex RNA systems and larger reagents, for example, N-methylisatoic anhydride (NMIA), appeared to be more complicated (see the Supporting Information (SI) part 3 for details). Starting structure was taken from the protein–RNA complex determined by X-ray crystallography.[39] Structural dynamics of individual nucleotides was first quantified by running a 500 ns long classical MD simulation and some enhanced sampling simulations.[40] Namely, bias-exchange MD simulations[41] were applied to enforce C3′-endo/C2′-endo sugar-pucker flips, and a similar approach[41,42] was then used to investigate the interplay between reagent binding and RNA reconformation to overcome binding free-energy barriers and bring the reactive carbon of the BzCN reagent toward the 2′-OH group of each SRP–RNA nucleotide. Results were correlated against a recent experiment using BzCN reagent for investigating cotranscriptional folding of SRP–RNA.[43] SHAPE data sets from this experiment are freely available and were retrieved from the RMDB database.[6] In addition, we performed other tests involving smaller RNA systems: (i) a gccGUAAggc GNRA tetraloop with BzCN SHAPE reagent and (ii) two nucleotide analogs, that is, 3′-5′-cyclic-adenosine monophosphate (cAMP) and 3′-5′-cyclic-cytosine monophosphate (cCMP), for which we used larger (NMIA) SHAPE reagent to directly compare results with available experimental data.[7] All simulations were performed using GROMACS[44] and PLUMED[45] (see the SI part 2 for details). We note that the simulation setup of SRP–RNA motif contained soft restrains that were required to (i) compensate for previously reported force-field deficiencies[46−50] and to (ii) avoid unfeasibly long simulations in an attempt to capture the whole structural ensemble upon multiple RNA unfolding/refolding events. Restraints were applied using an RNA-dedicated metric (ϵRMSD)[49,51] (see the SI parts 1−7 for details).

Figure 1

Tertiary and secondary structures of SRP–RNA motif. The left panel shows the tertiary structure used as starting point for MD simulations (PDB ID 1DUL).[39] Nucleotides from GNRA tetraloop, first and last noncanonical base pairs from the symmetric loop, all nucleotides from the asymmetric internal loop, and U6-A36 closing base pair from the asymmetric loop are labeled (unpaired nucleotides in capital letters). Structural model of small BzCN reagent is displayed (not to scale) on the left side with the reactive carbon (C*) highlighted. The right panel shows the corresponding secondary structure model, where the coloring scheme for each nucleotide is based on experimentally measured SHAPE reactivity.[43]

As a preliminary step, we analyzed the X-ray structure and the plain MD trajectory using previously proposed approaches, namely, computing solvent accessibility of 2′-OH and fluctuations of base–base distances. Results are reported in the SI (part 7, Figures S1–S3) and are in agreement with previously reported analyses of different RNA molecules.[7,34,36,52] We then analyzed RNA structural dynamics, in particular, monitoring the transient C2′-endo sugar puckers. Nucleotides from SRP–RNA were enforced to undergo C3′-endo/C2′-endo sugar-pucker flips during enhanced sampling simulations. The C2′-endo population analysis revealed that the data are correlated with experimental reactivity (Pearson’s linear R value of ∼0.5, Figure ) and the sugar pucker could be considered as a possible structural factor for distinguishing among reactivity patterns (see the SI part 7 for details). To explicitly take into account the RNA–reagent dynamical interplay, we performed enhanced sampling simulations and modeled the propensity of the reagent to reach a reactive conformation. We estimated relative binding constants between BzCN reagent and 2′-OH groups of nucleotides from SRP–RNA. The results from these bias-exchange-like simulations revealed rather poor correlation (∼0.3) against experimental SHAPE reactivities (see Figure S5 in the SI). However, the correlation could be enhanced up to ∼0.6 by combining the two above-mentioned procedures, that is, by applying the sugar-pucker population analysis to the results from enhanced sampling simulations (Figure ). An equivalent analysis made on a simulation where RNA was kept completely frozen lead to zero correlation, unambiguously revealing that RNA flexibility is a fundamental determinant for SHAPE reactivity (Figure ). In a further attempt to increase the correlation, we analyzed the presence of nearby RNA groups with the ability to function as general bases in deprotonation of 2′-OH,[34] but the overall correlation was not improved. Interestingly, the outputs from two different MD-based approaches involving a completely different analysis, that is, (i) monitoring fluctuations of distance between consecutive nucleobases[36] in classical MD simulations and (ii) combining bias-exchange-like simulations with sugar-pucker analysis, appear to capture the same phenomena as their results are highly correlated against each other (R ∼ 0.8, see Figure S6 in the SI). In summary, all of the successful approaches discussed here explicitly take into account RNA dynamics, supporting the fact that flexibility is crucial for SHAPE reactivity.

Figure 2

Correlations between experimental and calculated SHAPE reactivities of nucleotides from SRP–RNA. Top (A,B) and bottom (C,D) panels show results from simulations with frozen (X-ray like) RNA and flexible RNA, respectively. On the left side (A,C), results based on RNA-only structural analysis are shown, that is, the distribution of C2′-endo puckers in the X-ray structure and their respective populations from enhanced sampling simulations (see Figure S4 in the SI for the comparison between classical MD and enhanced sampling simulations). On the right side (B,D), both structural analysis of sugar puckers and relative binding rate constants from bias-exchange-like simulations with BzCN reagent were used. See Figure S5 in the SI for the results from enhanced sampling simulations not modified by additional sugar-pucker population analysis. Computed and experimental reactivities were derived as a free energy (as logarithm of SHAPE reactivity) and subsequently normalized against one particular WC base-paired nucleotide (see the SI part 4 for details).

We notice that experimental reactivities from four independent experiments are available,[43] but attempts to correlate our predictions against different experimental data sets revealed minimal differences. Hence, only R values against the first data set are discussed in the main text (see Table S1 in the SI for all raw data sets and selected cross-correlations). Experimental data sets are significantly correlated against each other (R values from 0.73 to 0.98, see the SI part 5 for details), confirming that the obtained reactivities are reproducible. However, the R values indicate that any prediction on the order of 0.75 or better would be practically undistinguishable from experimental data sets.

Our results from bias-exchange-like simulations with the SRP–RNA system might be influenced by the choice of the X-ray structure used as a starting point for MD simulations, which was taken from a protein–RNA complex.[39] Other X-ray structures of SRP–RNA are available, but they are either in complex with proteins as well or display sequences that are different from the one used in SHAPE experiment (especially in less structured parts). The protein is interacting with the symmetrical loop of SRP–RNA stabilizing five NC base pairs (from C14-A31 to U18-A27, see Figure ).[39] NMR data in solution are, however, available for only a small part of SRP–RNA (terminal RNA hairpin) and suggest that NC duplex region may prefer different structural order.[53] In particular, considering the X-ray structure, the classical 500 ns long MD simulation (see the SI part 2 and 7 for details), as well as all enhanced sampling simulations, one of the reported nuclear Overhausser effect (NOE) contacts[53] involving adenine H2 protons, that is, A15(H2)...G29(H1′), is never formed, and the other (C19(H6)...A27(H2)) is often violated during both classical and enhanced sampling simulations. In addition, NMR data suggest that U13, A23, A24, and C30 should transiently sample C2′-endo pucker,[53] but we did not detect any repuckering of these nucleotides either in classical MD or during bias-exchange-like simulations with BzCN reagent. The bias-exchange simulation enforcing pucker flips revealed only minimal populations of C2′-endo pucker for C13, A23, A24, and C30 nucleotides (probing them as nonreactive, Table S1 in the SI). Thus our choice to initialize the simulations from the X-ray structure could explain the disagreements in reactivity between prediction and experiment that were observed in NC duplex and tetraloop regions (Figure S7 in the SI). On the contrary, we propose that discrepancies between our calculations and experimental SHAPE data could be used to identify dynamical regions, where a single X-ray structure might not be representative. This is probably the case for the asymmetric loop (NMR data not available), which contains a closing U6-A36 base pair (Figure ) and is probably dynamic, as suggested by both experimental[43] and computed SHAPE reactivities (see Figures S7 and S8 in the SI).

On the basis of our results, we propose that the modification of a particular RNA nucleotide by SHAPE reagent occurs in a sequence of stages (Figure ). First, reagent reaches the nucleotide and is stabilized in the close proximity of its 2′-OH group mainly by stacking with sugar rings and nucleobases and by some (rather weak) H-bond contacts. We clustered all possible reactive conformations of SRP–RNA nucleotides and identified that stacking to sugar rings is the most frequent in unpaired as well as WC and NC base-paired nucleotides (Figures S9 and S10 in the SI). WC paired nucleotides mostly accommodated BzCN on their sugar ring, whereas unpaired and NC paired rather favored stacking on other sugars or nucleobases. This is not surprising because stacking in canonical A-RNA duplex is mostly saturated.[35] Our limited sample of diverse nucleotides shows that reactive geometries of WC and NC paired nucleotides typically provided a similar number of different clusters, whereas the number of clusters from unpaired nucleotides was slightly (∼1.3 times) higher on average (Figures S9 and S10 in the SI). Notice that identified stacking preference to sugars here could be driven by the specific SHAPE reagent (BzCN) or due to the number of restrains required for the simulations (see the SI part 3 for details). In general, we expect more frequent reagent stacking to nucleobases, larger reconformations of RNA motifs, and other structural changes to be induced, especially by larger SHAPE reagents (like NMIA and its derivative).

Figure 3

Hypothetical mechanism of RNA modification by SHAPE reagent. Schematic profile, where the free energy depends on the distance between reactive carbon (C*) of BzCN reagent and O2′ oxygen from particular nucleotide. Insets represent snapshots of SRP–RNA motif with BzCN in unbound state (U), bound state in close proximity from 2′-OH group (B, highlighted in double-inset), and anticipated reaction adduct after the acylation reaction (P). Bound states are here defined as structures for which the distance between C* and O2′ is <4.0 Å. Water molecules and counterions are not shown for clarity. We find that the barrier for binding (ΔGbind‡) and unbinding is marginal (∼1 kcal/mol, Figure S11) in comparison with the activation barrier of the acylation reaction (estimated ∼20 kcal/mol; see the SI part 4 for details). The activation free energy could be affected by some structural factors, for example, by the sugar pucker of the particular nucleotide in a reaction complex with reagent (in the range of a few kcal/mol favoring C2′-endo sugar pucker in a free energy; Table S2 in the SI).

After the initial binding, a reaction undergoes with a rate proportional to exp(−ΔGreact‡/RT). Considering the acylation reaction producing SHAPE adduct as irreversible, the measured amount of each adduct is directly proportional to the reactivity of the corresponding nucleotide. We here find that a simple model where the reaction barrier (ΔGreact‡) only depends on the sugar conformation such that the C2′-endo pucker is more reactive is sufficient to obtain a high correlation with experiments. The typical barriers that we observe for unbinding are very small (∼1 kcal/mol, see Figure S11 in the SI for calculated free-energy binding profiles of all analyzed SRP–RNA nucleotides) and certainly much lower than the expected ΔGreact‡ of the following chemical reaction (Figure ). The reactivity of a particular nucleotide is thus proportional to the binding rate constant of the SHAPE reagent toward the corresponding 2′-OH group with a correction taking into account the sugar pucker. Consequently, we are able to distinguish reactivity patterns without the need to explicitly simulate the acylation reaction. This makes the presented protocol robust even without usage of computationally demanding ab initio methods or problematic semiempirical potentials that were required for a recent investigation of inline probing experiments.[54] Assuming the typical reagent concentration in SHAPE experiments (∼100 mM)[9,10,43] and our calculated absolute binding constants (from 0.2 to 6.4 M–1) between BzCN and nucleotides in a common structural motif, the roughly estimated ΔGreact‡ should be ∼20 kcal/mol to get one modification per molecule under experimental time scale of minutes, which is a reasonable value for this type of reaction (see the SI part 4 for details). In addition, the estimated absolute binding constants indicate that the reagent binding likely occurs before the activation (deprotonation) of the particular 2′-OH group.

The original study of Weeks and coworkers identified that nucleotide analogs with 3′-phosphodiester moiety constrained away from the 2′-OH group by 3′-5′-linkage have one of the highest reaction rates with NMIA reagent among various nucleotide analogs[7] and are thus denoted as “hyper-reactive”.[34] We used two of these analogs (cAMP and cCMP) for additional simulations to test if our designed simulation setup was able to detect possible differences in reactivity among different nucleotides, that is purines and pyrimidines. We analyzed NMIA reagent binding toward cAMP and cCMP and estimated absolute binding constants of 6.6 and 2.9 M–1 for cAMP and cCMP, respectively (see the SI part 4 for details). In other words, NMIA reagent established a measurably stronger interaction with cAMP (here used as an analog of purine nucleotides) than cCMP (as corresponding analog of pyrimidine nucleotides). A similar observation was reported from an experimental analysis of large data sets, where cytosines were shown to be the least reactive among nucleotides with relative reactivities almost two times lower than adenines.[55] Those differences could partly originate from the different pKa of 2′-OH groups between purine and pyrimidine residues.[55,56] On the basis of our results, we speculate that the difference could be rather explained by the better ability of purine nucleotides to stabilize SHAPE reagents in the close proximity of 2′-OH groups, namely, by stronger stacking interaction. However, this systematic difference in reactivity between purine and pyrimidine nucleotides could be hidden or even reversed when analyzing individual structural motifs.[8] For instance, considering results for SRP–RNA, the average reactivity of pyrimidine nucleotides is ∼1.5 times higher than that of purines, both in our prediction and in the experimental data set.[43] Indeed, other effects such as stacking interactions with other residues and ribose repuckering might be more important than nucleobase identity.

In summary, we here provide a computational attempt to characterize the interactions between RNA nucleotides and SHAPE reagents. We developed and applied a novel protocol to predict reactivity patterns of small RNA structural elements as well as larger RNA motifs, estimated the reagent’s binding constants, and compared results with experimental SHAPE data sets. Our analysis shows that RNA flexibility is crucial to allow the reagent to reach the 2′-OH group, which requires some local RNA reconformations. Indeed, among all of the tested methods, only those where RNA dynamics was explicitly or implicitly (as in the ENM approach) included were capable of reaching a significant correlation with experiments. The explicit simulation of RNA–reagent binding allowed us to reach important structural insight on the SHAPE modification of nucleotides. Sugar rings and nucleobases are directing the reagent to the close proximity of 2′-OH group by the stacking interaction, and moieties from neighboring and other residues located further away along the RNA chain are frequently involved in stabilizing the reagent in the reactive conformation. The sugar pucker of the ribose ring is also an important structural factor that can determine the reactivity of particular nucleotides. Ultimately, our data show that the probability to observe the reagent in the binding position at the usually employed concentrations is much higher than the probability to observe a deprotonated 2′-OH group, which is required for the following acylation reaction to happen. We speculate that the protocol developed here could be employed to compute the reactivity patterns of predicted structures associated with sequences for which SHAPE data are available and, subsequently, to rank those putative structures by their agreement with experimental data. If the capability to discriminate decoy structures will be confirmed, then this procedure may allow the direct utilization of SHAPE data in 3D structure prediction. In addition, because many RNA chemical probing methods are based on the idea of using small molecules to probe equivalent sites located on different nucleotides, we anticipate that the introduced protocol may be suitably adapted to predict reactivities obtained using other experimental techniques.