Despite the limited number of ribonucleotide residues, RNA molecules are capable to fold into a wide variety of secondary and tertiary structures. This structural diversity is based on the capability of RNA to form networks of stacked nucleobases together with several hydrogen bonds between donor and acceptor functional groups not only in the nucleobase but also in the ribose moiety and along the phosphodiester backbone. The variability in structure is matched by considerable diversity in biological functions. Besides RNA's encoding and translational activity, enzymatic and regulatory functions have been found in biological systems as well as in artificial systems evolved in SELEX experiments. These functions cannot be explained solely based on characteristics of static structures but dynamic properties appear to play an important role in maintaining the functional diversity of RNA (1). RNA dynamics span a broad range of time scales from picoseconds where vibrations and angular fluctuations occur up to seconds where catalytic function and global refolding take place (2–4).

By NMR heteronuclear relaxation measurements, insight can be obtained into the pico- to nanosecond motions (5) by analyzing the relaxation rates R₁, R_1rho and the steady-state heteronuclear NOE of nitrogen and carbon nuclei following the model-free formalism (6–8). Such information has recently been augmented by data stemming from residual dipolar couplings in RNA, where RNA helices have been specifically elongated to shift the time scale of the overall rotational tumbling (9). While analysis of relaxation data for ¹H–¹⁵N sites is widely used for proteins to gain information on the backbone dynamics, the dynamical investigation of nucleic acids from protonated ¹⁵N nuclei is quite limited. Only for base-paired guanine and uridine imino groups, relaxation data can be retrieved, which restricts the analysis mostly to the helical stems and to the nucleobase moieties of RNA. In contrast, ¹³C relaxation has proven to be a strong alternative since data can be obtained for all nucleobase and ribose moieties by different strategies (10–15) to fulfill the model-free requirement of an isolated I-S spin system. Recently, Bax and co-workers (16) investigated protein dynamics and demonstrated that also for nonisolated spin systems found in uniformly ¹³C-labeled proteins or RNA, accurate relaxation data can be obtained and analyzed up to an overall correlation time in the low nanosecond timescale. These reports initiated detailed dynamical analyses of RNA molecules (17,18), which revealed that also parameters such as C–H bond lengths and ¹³C chemical shielding anisotropies need to be known accurately to derive correct information about the amplitude of local motions in RNA.

To compare the findings of NMR relaxation experiments with microscopic models, it is instructive to perform complementary molecular dynamics (MD) studies of RNAs. By using an all-atom force field for solute and solvent, classical MD simulations directly provide information at the atomic level on inter- and intramolecular motions (19,20). In protein research, the combination of NMR and MD investigations has been widely used (21–28) and recently, also comprehensive descriptions of fast conformational dynamics of small RNA and DNA systems have been published (29–32, Duchardt,E., Nilsson,L. and Schleucher,J., submitted for publication). While NMR results may serve as benchmark data to study the accuracy of the MD description, the MD trajectory in turn provides a microscopic model for the underlying conformational dynamics.

To investigate the correlation of structure, dynamics and function, RNA hairpins with YNMG loop motifs (Y = C or U, N = any, M = A or C) have been chosen for the present investigation. The most stable RNA tetraloop with a UUCG sequence closed by a Watson–Crick CG base pair has been intensely studied structurally and dynamically both, by NMR (17,18,33,34,86) and MD (30,35–38). This tetraloop motif is highly abundant in nature, stabilizes structures and serves as folding initiation site (secondary structural elements of the cUUCGg tetraloop are summarized in Table S12 in the Supplementary Data). No ligands are known to bind with high affinity to this tetraloop sequence. In contrast, a structurally highly similar hairpin with a CACG loop sequence closed by an UG wobble base pair is more labile with a melting temperature decreased by more than 15 K (39). In Coxsackievirus B3, this tetraloop interacts directly with the protein proteinase 3C as essential component of the viral replication initiation complex. These differing properties cannot be explained by the surprisingly small structural variations (secondary structural elements of the uCACGg tetraloop are summarized in Table S13 in the Supplementary Data). Our approach of investigating temperature-dependent NMR-derived order parameters S² has been previously applied to the study of (local) protein folding and unfolding (40–45).

In this work, we investigate the dynamical properties of both hairpin types by NMR ¹³C relaxation and MD simulation. We evaluate in detail parameters such as carbon–hydrogen bond length and ¹³C-chemical shift anisotropies, which significantly affect the interpretation of subnanosecond dynamics by the model-free analysis. The hairpins dynamics were recorded at a temperature close to the melting start of the less stable uCACGg hairpin to investigate the differences of local and global dynamics of the two tetraloops. MD simulations suggest a similar pathway of temperature-induced unfolding of the two hairpins. Interestingly, at high temperatures, nonnative stacking interactions are revealed in one of the four main conformations.

1

2

Further, the CC-dipolar contribution to the relaxation of the carbons C₆ and C_1′ from the adjacent carbon nuclei becomes more relevant with an increase of the overall correlation time τ_c (51):

3

For the smaller 14-mer RNA a minimum error of 2% for the R₁ and R₂ relaxation rates and 3% for the ¹³C–{¹H} steady-state NOE values were used. These minimum errors were increased for the larger 30mer RNA to 3% for the R₁ and R₂ rates and 5% for the ¹³C–{¹H} steady-state NOE. For some of the relaxation rates, the experimentally derived errors were larger and then in these cases used.

For the analysis of the relaxation parameters the program Modelfree 4.15 by Palmer and co-workers (52) was used. The further residue-specific parameters (C–H bond lengths and ¹³C chemical shift anisotropy) are discussed and optimized in this paper. Model selection was according to d’Auvergne and Gooley (53) and Wright and co-workers (54) with a Bayesian Information Criteria (BIC). The procedure is discussed in detail in the Supplementary Data. The global parameters, the overall rotational correlation time τ_c and the diffusion anisotropy were first determined by a hydrodynamic calculation performed with the program hydronmr 5a (55). These values were then optimized by the program Modelfree using the relaxation data.

A twin range cut-off was used for the Lennard–Jones interactions that is interactions between atoms within 1.0 nm were evaluated every step, while interactions between atoms within 1.4 nm were evaluated every five steps. The particle mesh Ewald method (61) was employed to treat Coulomb interactions, using a switching distance of 1.0 nm. Constant pressure p and temperature T were maintained by weakly coupling the system to an external bath at one bar and at the selected temperature, using the Berendsen barostat and thermostat, respectively (62). The RNA, the ions and the solvent were independently coupled to the temperature bath with a coupling time of 0.1 ps. The pressure coupling time was 0.5 ps and the isothermal compressibility 4.5 × 10⁻⁵ bar⁻¹. The bond distances and the bond angle of the solvent water were constrained using the SETTLE algorithm (63). All other bond distances were constrained using the LINCS algorithm (64). A leapfrog integrator with an integration time step of 2 fs was used.

The starting structure of 14-mer cUUCGg hairpin was taken from simulations previously performed at 298 K (32). The starting structure of the 14-mer uCACGg hairpin was taken from the Brookhaven Protein Data Bank (PDB structure 1RFR.pdb) (39). Fifty-nanosecond classical MD simulations have been performed at 298 K, 317 K and 360 K.

Analysis of the trajectories was performed with tools from the GROMACS package and with modified versions of them. To define the presence of a hydrogen bond, an acceptor–donor distance smaller than 0.35 nm and an acceptor–hydrogen-donor angle larger than 150°degrees was requested. To define base–base hydrogen bond interactions, the presence of the N–H···N hydrogen bond between Watson–Crick base pair and N–H···O hydrogen bonds between wobble base pair is used as indicator. Two bases were considered as stacked when their center of mass separation is within 0.15 nm of that seen in the folded structure (simulation at 298 K). Pairs with larger separations are considered as broken (65).

To obtain NMR order parameters, the corresponding internal correlation functions were calculated from Equation (4) (66):

4

5

To compare the conformational fluctuations of the hairpin at different temperatures, we performed a principal component analysis (67–69) of the trajectories at 298, 317 and 360 K, where only atoms of residues C5-G10/U13-G18 were considered. Preceding the analysis, each conformation was translated and rotated (through a least root mean square fit) to give the best fit to the averaged structure. For the cUUCGg hairpin, the first two (of in total 579) eigenvectors describe 57 and 47% of the total fluctuations at 298 and 360 K, respectively. For the uCACGg hairpin, 46 and 44% of the total fluctuations are covered.

Recently, ab initio calculations were performed by Case to incorporate zero-point motion averaging of the bonds in proteins (73) and in nucleic acids (Case,D., personal communication). The corresponding bond lengths for nucleotides are summarized in Table 1. The incorporation of the zero-point motion averaged bond lengths into the model-free analysis results in increased order parameters S² close to or even exceeding the maximum value of unity. These unphysical results (S² > 1) can be remedied for the base vectors C6H6 and C8H8 by correction of the carbon chemical shielding anisotropy (¹³C CSA).

Table 1. C–H bond lengths used in earlier studies for relaxation analysis (R) and zero-point motion averaged bond lengths which give an effective bond length reff

In 2006, ¹³C-CSA values have been reported both, for a helical A-form RNA segment and a helical B-form DNA dodecamer (74). These CSA values differ significantly from the previous values (Table 2) with changes up to 16%. Trantirek and co-workers (75) calculated the conformational dependence of cross-correlated relaxation rates across the glycosidic torsion angle χ of deoxynucleosides. From their analysis, a dependence of the ¹³C-CSA values on the sugar pucker mode was proposed (Brumovska,E., Sychrovsky,V., Vokacova,Z., Sponer,J., Schneider,B. and Trantirek,L., manuscript in preparation). While for the aromatic carbons the influence of the sugar pucker mode on the CSA is negligible, the dependence on the torsion angle χ is significant. For a purine nucleotide in syn-conformation (χ-angle of around 60°), the absolute value of the carbon C8 CSA is reduced by 10 to 13 p.p.m. For pyrimidines, the difference of 20 to 30 p.p.m. for the C6 CSA is even more pronounced. These results are of particular interest for the analysis of the hairpins studied here, since the guanosine residue in the forth position of an YNMG-loop adopts a syn-conformation. In our previous analysis, this guanosine residue appeared as the second most flexible moiety while the residue to which guanosine is base paired to at the first position of the loop was as rigid as the stem residues (17). In addition, the order parameters derived from ¹³C relaxation rates were different from those derived from a ¹⁵N relaxation analysis (34). Taking into account the new predictions of a χ-value dependence of the CSA (Table 2) alleviates the differences between the previous analyses and yields fully consistent data for the guanine residue (Figure S1 in the Supplementary Data).

Table 2. Chemical shift anisotropies (CSA) of aromatic carbon nuclei in nucleobases

For the ribose carbon C1′, no CSA values have been published so far that take a zero-point motion averaged bond length into account. The published CSA range from 29 to 60 p.p.m. and were determined by DFT calculations on ribose and deoxyribose model compounds (76,77), solid-state measurements on isolated nucleotides (78) or liquid-crystal NMR on an A-form helical RNA (79). These measurements showed no differences for the individual nucleotides. For our analysis, we used the previous parameters with a bond length of 1.09 Å and an absolute value of the chemical shielding anisotropy of 45 p.p.m.

Hydrodynamic calculations using the program hydroNMR (80) of the 14 mer RNA resulted in starting values of the global correlation time τ_c ≈ 2.30 and a diffusion anisotropy D_/ ≈ 1.50. The parameters were further optimized using the program Modelfree (52) during the fitting procedure. This optimization results in a correlation time of 2.27 ns, which lies between our previously published values of 2.18 and 2.44 ns, as a result of separately fitting the carbons in the nucleobases and the ribose moieties. The diffusion anisotropy of 1.48 is close to the value calculated by hydroNMR.

Figure 1. Order Parameters S2 of the 14 nt cUUCGg hairpin (left of the diagrams) extracted from NMR 13C relaxation data (open diamonds) and from a MD trajectory (squares) for the nucleobase vectors C6H6/C8H8 (left) and the ribose vectors C1′H1′ (more ...)

Table 3. Results of the model-free analysis of the 14 nt cUUCGg RNA at 298 K extracted from the NMR relaxation data and the MD trajectory

An increase of the temperature to 317 K affects only the global overall tumbling of the RNA, whereas the internal motions do not change significantly. The diffusion anisotropy of D_/ = 1.47 is maintained and the global correlation time τ_c is decreased to 1.53 ns as expected (the hydroNMR derived value is 1.70 ns). In addition to the faster global motion, also the order parameters decrease slightly at a higher temperature. For the nucleobase moieties, the average decrease is larger and more uniform (0.98 to 0.87 for C8H8 and 0.87 to 0.80 for C6H6) than for the ribose moieties (0.90 to 0.87) where the loop riboses of C8 and G9 even have a slight increase of S² (Figure 1 (mid right) and Table S8 in the Supplementary Data). Of all nucleobase moieties, only U7 is flexible with a fast internal motion and a very small R_ex-term.

Figure 2 shows the order parameters of the uCACGg 30-mer RNA hairpin. In the uCACGg hairpin, the nucleobase of residue A15 at the second loop position is the only flexible moiety (Figure 2 top left). The dynamics of the nucleobase of A15 can be described with a fast motion (<10 ps) and a slow motion with an internal correlation time of 441 ± 121 ps (Table 4). The analyses of the C1′ relaxation rates are more difficult. A number of residues around the loop can hardly be fit, especially at the temperature of 298 K. The order parameters of these residues fluctuate between 0.8 and 1.0. In the loop, the sugar moiety of A15 is the most flexible with a slow motion of 827 ± 445 ps. The other C1′H1′ ribose sites for which fits can be obtained have R_ex-terms between 1.9 and 11.8 s⁻¹ (Table 4).

Figure 2. Order Parameters S2 of the 30-nt (NMR) and 14-nt (MD) uCACGg hairpin (left of the diagrams) extracted from NMR 13C relaxation data (open diamonds) and from a MD trajectory (squares) for the nucleobase vectors C6H6/C8H8 (left) and the ribose vectors C1′H1′ (more ...)

Table 4. Results of the model-free analysis of the 30 nt uCACGg RNA hairpin part G9 to C22 at 298 K extracted from the NMR relaxation data and of the 14 nt uCACGg RNA extracted from the MD trajectory

For the uCACGg hairpin, the changes at higher temperatures become more pronounced. At a temperature of 317 K at which the imino signals of the non-Watson–Crick base pairs can no longer be observed in a ¹⁵N–HSQC (39), no significant changes for the order parameters of the aromatic and anomeric carbons can be observed in the NMR analysis. The global correlation time decreases from 6.55 to 4.15 ns, which is comparable to the predicted decrease in global correlation time by hydroNMR of over 2 ns. For the internal dynamics of the nucleobase moieties, a decrease of the order parameter by <0.03 is observed while the overall profile remains unchanged (Figure 2 mid left). This result suggests that the disappearance of the imino resonances as detected in the ¹⁵N–HSQC is due to an increase of exchange rates with the solvent water at elevated temperatures but not due to an increase in dynamics. The doubling of the R_ex-term to 15.3 s⁻¹ of the C8 atom of the G17 nucleobase in the syn-conformation is significant and in agreement with the time range of the water exchange. The other loop nucleobases C14 and C16 display R_ex-terms of 7.4 s⁻¹ and 3.1 s⁻¹, indicating the appearance of a slow motion. The data for nucleobase of A15 is described by the same motional model with similar order parameters , and an analogous internal correlation time of 404 ± 136 ps (Table S9 in the Supplementary Data).

A further temperature increase to 325 K does not change the ¹³C–HSQC of the aromatic and of the anomeric region significantly, revealing that the averaged conformation is similar to the determined structure at 298 K. Therefore, a model-free analysis of the relaxation data can be performed (Table S10 in the Supplementary Data). The model-free analysis of the relaxation data at 325 K results in a global correlation time of 3.09 ns, which is exactly the same value as derived from hydroNMR. The average order parameters of all nucleobases drop to ~0.77 and the values obtained for the stem residues are only slightly higher (between 0.78 for the pyrimidine C6H6 and 0.81 for the purine C8H8). For the nucleobase of A15, the order parameters of the slow and the fast motion are very close to the values which are gained at lower temperatures. Only the internal correlation time of the slow motion is reduced by a factor of 2 to 224 ± 162 ps. The surrounding nucleobases of the loop and the loop closing base pair show the largest and most significant changes (Figure 2 bottom). For nucleobases U13, C14 and C16 involved in stacking interaction in the ground state conformation, the order parameters decrease and a similar internal correlation time of ~25 ps is obtained. The nucleobase of G17 exhibits a more dramatic increase of dynamics and becomes the second most flexible moiety. In addition, also a significant R_ex-term (20.9 ± 2.8 s⁻¹) and a fast internal correlation time (13 ± 11 ps) are obtained. The adjacent G18 shows a slower motion of 64 ± 29 ps and has a large decrease in . These observations indicate that the 3′-end of the loop is the primary melting site, while the 5′-end is more stabilized by the stacking of bases U13, C14 and C16.

The model-free analysis of the adenine C2 relaxation rates results in very similar parameters compared to the corresponding adenine C8 data in both RNAs and at all measured temperatures (Tables 3, 4 and S5, S6 and S8–S11 in the Supplementary Data). Although the NMR relaxation data of the C1′ can be fit at 325 K, a wave-like modulation of the order parameters is observed which we do not further investigate here.

The good overall agreement between experimental and calculated order parameters indicates that the underlying theoretical model is suitable to reproduce the temperature dependence of the RNA structure and dynamics and suggests that valid microscopic dynamic information can be extracted from the MD simulations. The MD trajectories were analyzed by a principal component analysis on the coordinates of residues around the loop region to identify the principal motions of the RNA hairpins (see Methods section). Figure 3 presents the resulting population probability along the first two principal components at the three temperatures under consideration, which reveal the coexistence of several conformational states. Figure 4 shows representative conformations of the hairpin for these states.

Figure 3. The cUUCGg (left) and uCACGg (right) hairpin configurations at 298 K (top), 317 K (middle) and at 360 K (bottom) projected on the first two PCA eigenvectors. The PCA analysis was performed on the atoms coordinates of residues C5–G10/U13–G18. (more ...)

Figure 4. Representative conformations of the PCA clusters (Figure 3) from the MD simulations at 298 K, 317 K and 360 K. Dashed lines around residue labels indicate that the residues point into the solvent. Parallel lines between residues indicate base pair interactions (more ...)

For the cUUCGg hairpin, the trajectories at 298 and 317 K are seen to consist of two conformational states (cluster 1 and 2). They reflect a rearrangement of the backbone involving residues U6 and U7. A closer analysis reveals that the observed dihedral transition correlates with the increase of the syn-oriented population of the nucleobase U7 (Figure 5). In contrast, the secondary structural interactions in the loop region are not affected. Upon increase of the temperature to 360 K, the principal component analysis reveals an additional conformational state (cluster 3), which is characterized by an unstacked nucleobase C8 (see Figure 4 bottom left). As expected, the unusually stable cUUCGg RNA hairpin generally does not exhibit major changes of structure and dynamics in the temperature range under consideration.

Figure 5. Dynamic behavior of the loop-out residues U7 (top) and A15 (bottom) at 298 K extracted from the molecular dynamics simulation trajectory revealed by the torsional angle χ.

The situation is different for the uCACGg hairpin. Due to its significantly lower melting temperature, this RNA loop exhibits substantial changes of the local dynamics on the pico- to nanosecond timescale upon an increase of the temperature from 298 to 360 K. While only one conformation is observed at 298 K in the principal component analysis in Figure 3 (cluster 1′), two and four widely spread clusters are populated at 317 K and 360 K, respectively. A first conformational change is observed at 317 K (cluster 2′), where a backbone conformational rearrangement causes the lost of the stacking interactions between the residues C16 and C14, see Figure 4 (mid right).

At 360 K, the secondary structure elements in the loop region are gradually lost and the closing base pair of U13 and G18 is also disturbed. The time evolution of the torsion angle χ and the pseudorotational angle P reveal that the motions of ribose and nucleobase are mostly uncorrelated (Supplementary Figure S2 in the Supplementary Data). Only the nucleobase and the ribose moiety of G18 have a high correlation value (0.84). Figure 4 (bottom right) comprises the main conformational states found in the 360 K trajectory of the uCACGg hairpin. In the cluster 3', the hairpin mostly satisfies the structural characteristics of the YNMG-motif with only the hydrogen bond between C14 and G17 missing. The cluster 4' displays that the nucleobase of C16 is not stacked anymore on top of C14 but loops out and potentially has stacking interactions with the neighboring A15. A mostly unfolded structure is observed in the cluster 5'. This cluster dominates the simulation with a population of ~40%. In the cluster 6’, the nucleobases seem rigid since stacking interactions between the neighboring residues C14 and A15 as well as C16 and G17 stabilize the loop. Snapshots are in Figure S3 in the Supplementary Data.

The MD calculated order parameters contain contributions arising from the characteristic fast dynamics for each of these conformations, while the conformational change between the clusters occur at 10 ns timescale, thus are not detected by ¹³C NMR relaxation measurement. The presence of several nonnative stacking interactions involving residues of the loop regions explain the similar flexibility observed for these residues ( ≈ 0.7–0.6 and ≈ 0.75–0.6). In particular, nonnative stacking interactions involving the loop-out residue A15 at the premelting temperature may explain the reduced flexibility of this nucleobase with respect to lower temperatures ( ≈ 0.5–0.4 at 298–317 K) and ( ≈ 0.7 at 360 K).

Compared to protein structures, RNA structural architecture appears to be considerably more modular. Tetraloops are an important structural modular element; they play an important role in RNA structure since they connect single-stranded segments involved in regular A-form helices. Loops are therefore often exposed and involved in biomolecular interactions and have also been postulated to form sites where folding and melting are initiated. The two tetraloops investigated here have similar ground state conformations (Figure 4 top) but differ in their thermal stability. At a temperature of 298 K, the two hairpins also reveal similar conformational dynamics on the pico- to nanosecond timescale (Figures 1 and 2). In both hairpins, increased motions are only found for the looped-out nucleobase at the second loop position (U7 in the cUUCGg hairpin and A15 in the uCACGg hairpin). This observation is reflected in the relative differences of the order parameters between stem and the looped-out nucleobase, which clearly differ for the two hairpins: - ≈ 0.26 (NMR) and 0.34 (MD) versus − ≈ 0.40 (NMR) and 0.44 (MD). The finding is a consequence of the nucleobase motion around the glycosidic bond angle χ (Figure 5). While U7 exhibits frequent anti/syn transitions, its fluctuations in these conformations are rather small (±30°). The anti/syn transitions are on a nanosecond timescale, longer than the overall correlation time of the molecule and therefore cannot be detected by NMR. The glycosidic angle of A15 varies continuously around 0° with an amplitude of 60° and only a slight preference for −60°. These fluctuations are fast enough to be detected by NMR ¹³C relaxation and therefore yield a lower order parameter.

By comparison of MD and NMR data, we can delineate the overall thermal unfolding pathway of the two hairpins and identify hot spots where unfolding is initiated at a temperature where the global fold is still maintained. When the temperature is increased, the loss of the stacking interactions between the first and the third nucleobase in the loop region is the first event that destabilizes the hairpin. In the MD simulations of the cUUCGg hairpin, it takes a temperature of 360 K for this to happen. At this temperature, the C8 nucleobase partially unstacks from the nucleobase U6 and extrudes into the solution (see cluster 3 in Figure 4), which is reflected in a S² decrease of U7 and C8. For the uCACGg hairpin, these changes already occur at a temperature of 317 K (see cluster 2′ in Figure 4). In both sequences, it is a cytidine stacking on the 5′-end of the loop, which has an intraresidual hydrogen bond to its ribose and a cross-strand hydrogen bond to the phosphate backbone between the first and second loop residue. Less favorable stacking interactions of the nucleobase C16 on the cytidine C14 and/or a slight deviation of the backbone along the top loop residues A15 and C16 which is observable in the structural calculations of the uCACGg hairpin (39) and the cUUCGg hairpin (Nozinovic,S., Fürtig,B., Jonker,H., Richter,C. and Schwalbe,H., manuscript in preparation) may contribute to the decreased stability of the uCACGg hairpin.

A further increase of the temperature to 325 K in the experiment and to 360 K in the simulation destabilizes all residues of the uCACGg hairpin (Figure 4 bottom right). None of the hydrogen bonds are detectable anymore in the NMR spectra (data not shown). The MD trajectory, too, displays enlarged distances and unfavorable arrangements between the hydrogen bond donors and acceptors in the loop. Moreover, the adjacent wobble base pair (U13–G18) is strongly destabilized in comparison to the canonical base pair (C5–G10) in the cUUCGg hairpin at this elevated temperature.

In summary, we have presented a detailed investigation of the thermal unfolding of two RNA tetraloops, which exhibit different stabilities although they are closely related in their ground state conformation. The large temperature range considered represents a challenge for both NMR and MD. On one hand, MD simulations are known to yield too high melting temperatures, as standard biomolecular force fields have been parameterized to reproduce room-temperature properties. To obtain the physically correct simulation temperature, the temperature was therefore rescaled according to a previously performed replica exchange MD simulation. For the NMR, on the other hand, the large temperature range represents a stringent test whether the parameters (effective bond length, CSA) have been chosen correctly to facilitate an appropriate and physically meaningful model-free analysis over the entire temperature range. It has been found that more recent data suggest a general increase for ¹³C CSAs of the aromatic sites in the nucleobase and in addition the conformation around the glycosidic torsion angle χ (syn versus anti) needs to be considered. Such correction has a pronounced effect on the CSA that propagates into the order parameter analysis. Since the experimental determination of ¹³C CSA is dependent on the choice of the C–H bond length, the consistent treatment of zero-motion averaged bond lengths is necessary. Following the recommendations in the protein field (73), we and others (74) utilize longer bond lengths. These adjustments now allow determination and interpretation of experimental data that translates into a precision for of a few percentages when using the model-free formalism.

Employing the thus established methods, we have presented a detailed temperature dependent NMR/MD study of the melting of the most abundant loop-motif in RNA. It has been found for both systems that the loss of the stacking interactions between the first and the third nucleobase in the loop region is the first event that destabilizes the hairpin. While the global fold is still preserved, this hot spot initiating unfolding was observed at 317 K for the uCACGg hairpin but only at an MD temperature of 360 K for the more stable cUUCGg hairpin. Increasing the temperature further, both NMR and MD elucidate the premelting of the uCACGg hairpin, which is accompanied by a loss of specific hydrogen bonds and stacking interactions. The study demonstrates the potential of joint MD/NMR studies to assess the dynamic properties of RNA secondary structure, which is ultimately required for a microscopic understanding of the various functions of RNA systems.

Supplementary Data are available at NAR Online.

We would like to thank Oliver Ohlenschläger, Matthias Görlach for providing us with a sample of the 30mer RNA and David Case, Ad Bax, Lucas Trantirek for sharing their data prior to publication. The work was supported by the Deutsche Forschungsgemeinschaft (DFG) (SFB 579: ‘RNA-Ligand-Interactions’), the state of Hesse [Center for Biomolecular Magnetic Resonance (BMRZ)] and the Fonds der Chemischen Industrie (H.S. and G.S.). Funding to pay the Open Access publication charges for this article was provided by the DFG-funded SFB 579: ‘RNA-Ligand-Interactions’.

Conflict of interest statement. None declared.