NMR and MD studies of the temperature-dependent dynamics of RNA YNMG-tetraloops.

In a combined NMR/MD study, the temperature-dependent changes in the conformation of two members of the RNA YNMG-tetraloop motif (cUUCGg and uCACGg) have been investigated at temperatures of 298, 317 and 325 K. The two members have considerable different thermal stability and biological functions. In order to address these differences, the combined NMR/MD study was performed. The large temperature range represents a challenge for both, NMR relaxation analysis (consistent choice of effective bond length and CSA parameter) and all-atom MD simulation with explicit solvent (necessity to rescale the temperature). A convincing agreement of experiment and theory is found. Employing a principle component analysis of the MD trajectories, the conformational distribution of both hairpins at various temperatures is investigated. The ground state conformation and dynamics of the two tetraloops are indeed found to be very similar. Furthermore, both systems are initially destabilized by a loss of the stacking interactions between the first and the third nucleobase in the loop region. While the global fold is still preserved, this initiation of unfolding is already observed at 317 K for the uCACGg hairpin but at a significantly higher temperature for the cUUCGg hairpin.

INTRODUCTION

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 R1, R 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 1H–15N sites is widely used for proteins to gain information on the backbone dynamics, the dynamical investigation of nucleic acids from protonated 15N 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, 13C 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 13C-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 13C 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 S2 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 13C relaxation and MD simulation. We evaluate in detail parameters such as carbon–hydrogen bond length and 13C-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.

METHODS

Sample preparation

The 14-nt UUCG-hairpin (5′-PO3−-PO2−-PO2−-GGCACUUCGGUGCC-3′; bold residues indicate the YNMG loop-residues) was purchased from Silantes GmbH (Munich, Germany) as a uniformly 13C/15N-labeled RNA. The concentration of the NMR sample was 0.7 mM in 20 mM KH2PO4/K2HPO4, pH 6.4, 0.4 mM EDTA and 10% v/v D2O. For the investigation of the 30-nt CACG-hairpin (5′-PO4−-GGCACUCUGGUAUCACGGUACCUUUGUGUC-3′; bold residues indicate the loop-residues) two samples were synthesized by in vitro transcription with T7-RNA polymerase with a linearized plasmid DNA as template and purified as described previously (46). One sample contained uniformly 13C/15N-labeled guanosine and cytidine residues with a RNA concentration of 0.6 mM, while the second was synthesized with uniformly 13C/15N-labeled adenosine and uridine residues at a final concentration of 1.2 mM. The buffer composed 10 mM KH2PO4/K2HPO4 adjusted to pH 6.2, 40 mM KCl, 0.2 mM EDTA and 99.99% v/v D2O.

NMR spectroscopy

All NMR experiments were carried out on a 600 MHz Bruker (Rheinstetten, Germany) spectrometer equipped with a 5 mm 1H {13C/15N} Z-Grad TXI probe. 13C R, R and 13C–{1H} steady-state NOE data were obtained using 13C modifications of the Bruker standard pulse sequences for 15N relaxation measurements (hsqct1etf3gpsi3d, hsqctretf3gpsi3d.2 and invnoef3gpsi) (47–50). For the 14 nt cUUCGg-hairpin, data were collected as published before (17). The 30 nt uCACGg-hairpin data were measured with slight changes compared to the 14 nt RNA. The carrier frequency was set to 137 and 89 p.p.m. for the aromatic and the ribose carbons with spectral widths of 10 and 12 p.p.m., respectively. During the carbon evolution, off-resonant carbon Q3 pulses (512 µs) were applied with an offset of −7000 Hz and −5000 Hz in order to suppress the 1J(C5, C6) and 1J(C1′, C2′) couplings, respectively. Further, long-range homonuclear couplings were assumed to be negligible. For the R1 and R1rho measurements, a relaxation delay of 2 s was applied between each of the 32 scans for all 104 complex points in the indirect dimension. R1- and R1rho-subspectra with varying relaxation delays were acquired in one interleaved pseudo-3D experiment. For the acquisition of the R1 relaxation rates, the relaxation delay τM was set to 10, 50, 100, 200, 400, 700 ms, 1 and 1.5 s, while spectra with τM = 50 and 400 ms were recorded twice for error determination. R rates were determined as described previously (17) with a spin-lock field of 3.6 kHz at an offset of 2000 Hz and relaxation delays τM of 12, 24, 36, 48, 64, 80, 104, 128 ms. Duplicate measurements were carried out for τM = 24 and 80 ms. The 13C–{1H} steady-state NOE experiments were recorded with 8 scans for 128 t-increments in an interleaved manner, with alternating proton-presaturated and nonpresaturated spectra. During the premeasurement delay of 6 s, on-resonance presaturation was applied for 3 s subsequent to a 3 s relaxation delay. To avoid heating differences between the proton-presaturated and the nonpresaturated spectra, an off-resonance presaturation was applied during the blank experiment. The interleaved spectra were separated by a Bruker standard macro. All measurements were repeated four times for error determination.

Data analysis

Spectra were processed and analyzed with Topspin 1.3 (Bruker, Germany). R1 and R relaxation decays were fit from peak heights to monoexponential two-parameter functions with the program Sparky 3.112 (University of California, San Francisco, USA). The measured Rmeas rates contain spin-lock offset and spin-lock power dependent contributions of R1 and R2 following Equation (1): where θ is the angle of the effective spin-lock field with the B0 field for each nucleus defined by Equation (2): with υ being the spin-lock field strength and Ω the offset of the nucleus resonance from the spin-lock carrier, both in Hertz. The influence of homonuclear Hartmann–Hahn magnetization transfer on the R1rho rates can be neglected for the nuclei investigated (17).

Further, the CC-dipolar contribution to the relaxation of the carbons C6 and C1′ from the adjacent carbon nuclei becomes more relevant with an increase of the overall correlation time τc (51): in which ρCi is the dipolar contribution to the relaxation of the 13C nucleus arising from the attached 1H or 13C nucleus; γi is the gyromagnetic ratio of the attached nucleus; rCi is the distance to the attached nucleus; ωC is the carbon Larmor frequency and τc the overall correlation time. But also for the 30 nt RNA with a τc at 298 K of about 6 ns the contributions are still below 1% and therefore are neglected in the model-free analysis.

For the smaller 14-mer RNA a minimum error of 2% for the R1 and R2 relaxation rates and 3% for the 13C–{1H} steady-state NOE values were used. These minimum errors were increased for the larger 30mer RNA to 3% for the R1 and R2 rates and 5% for the 13C–{1H} 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 13C 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.

MD simulations

The MD simulations were performed using the GROMACS suite of programs (version 3.2) (56,57). The AMBER force field (parm98) (58,59) was employed to describe the 14-mer cUUCGg and uCACGg RNA hairpins. The hairpin was placed in a rhombic dodecahedron box (edge length ∼5 nm), which was subsequently filled with TIP3P water molecules (60). To neutralize the system, 13 sodium ions were placed randomly in the simulation box.

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−5 bar−1. 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): where is a unit vector pointing along the C–H bond, is the second Legendre polynomial, and denotes an equilibrium average over the MD trajectory from which molecular rotation and translation was removed. Subsequently, the first 100 ps of the correlation functions were fit to the Lipari–Szabo expression (6): where S2 is the order parameter and τ is the effective (or internal) correlation time for the C–H dipole. Using the first nanoseconds of C(t) hardly changes the results (32).

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.

RESULTS

Analysis of relaxation rates and parameter optimization for the model-free analysis

The modulations of the dipolar and the chemical shift anisotropy fields are the two main components influencing the relaxation of NMR-active nuclei. The model-free analysis of 13C relaxation data, in turn, depends critically on the value of the C–H bond length and the carbon chemical shielding anisotropy (5,17). These two parameters are fundamental, since they weigh the two dominant contributions to the 13C relaxation and their determination is often mutually dependent. In addition, it is often assumed that the chemical shift anisotropies do not depend on conformation. Studies on proteins, however, have shown a clear conformational dependence of the 15N CSAs (70–72).

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 S2 close to or even exceeding the maximum value of unity. These unphysical results (S2 > 1) can be remedied for the base vectors C6H6 and C8H8 by correction of the carbon chemical shielding anisotropy (13C 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

	bond	R [Å]	reff [Å]
Ribose	C1′–H1′	1.095	1.118
Pyrimidine	C6–H6	1.082	1.102
Purine	C8–H8	1.079	1.102
	C2–H2	1.079	1.102

In 2006, 13C-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 13C-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 13C relaxation rates were different from those derived from a 15N 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

Carbon nucleus	Nucleobase	CSA [p.p.m.]
		Formerly used	New
C8	Adenine	−134	−144
	Guanine (anti-conformation)	−134	−133
	Guanine (syn-conformation)		−122
C6	Cytidine	−179	−208
	Uridine	−179	−208
C2	Adenine		−168

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.

Temperature dependence of S2

In order to study the temperature dependence of the conformational dynamics of both RNA hairpins, NMR relaxation measurements have been performed at room temperature (298 K), at a temperature (325 K) close to the melting point of the uCACGg hairpin, and at an intermediate temperature (317 K). To compare these results with corresponding MD simulations, one needs to take into account that MD melting studies typically overestimate melting temperatures (81–84). Besides sampling problems, this temperature shift is most likely related to deficiencies of the force field, since standard biomolecular force fields are parameterized to reproduce properties at room temperature. In fact, a recent replica exchange MD simulation study on 14-mer RNA hairpins (85) resulted in melting temperatures of 425 K and 397 K for cUUCGg and uCACGg hairpins, respectively. To achieve premelting conditions for uCACGg hairpin in the MD simulation, we therefore employed a temperature of 360 K instead of the experimental temperature of 325 K. As to our knowledge, all current empirical force fields are plagued by this problem; it remains to be shown by future theoretical models whether the simple temperature rescaling yields the correct structures and dynamics.

Experiment

Figure 1 shows the NMR-derived order parameters of the cUUCGg 14-mer hairpin for the nucleobase vectors C6H6/C8H8 (left) and the ribose vectors C1′H1′ (right). The results of the model-free analysis at 298 K are given in Table 3 and in Figure 1 (top). The nucleobase of loop residue U7 is by far the most flexible moiety in the entire RNA molecule. This observation is in agreement with the high-resolution structures (33; Nozinovic,S., Fürtig,B., Jonker,H., Richter,C. and Schwalbe,H., manuscript in preparation) that reveal no stabilizing interactions like hydrogen bonds or stacking for this nucleotide. Also, the loop riboses of U7, C8 and G9 display slightly increased dynamics compared to the rigid structured helical stem. In the high-resolution NMR solution structure (Nozinovic,S., Fürtig,B., Jonker,H., Richter,C. and Schwalbe,H., manuscript in preparation), the riboses U7 and C8 are in C2′-endo conformations [pseudorotation phase P (U7) = 163° and P (C8) = 176°] and the ribose of G9 has a syn-conformation of its nucleobase.

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′ (right) at temperatures of 298 K (top), 317 K (middle) and 360 K (bottom). The loop region is shaded in gray. For the highest temperature of 360 K, NMR relaxation data are not available due to stability reasons.

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

	C6H6/C8H8/C2H2a					C1′H1′b
	NMR				MD	NMR			MD
	S2	Sf2	τe [ps]	Rex [s−1]	S2	S2	τe [ps]	Rex [s−1]	S2
G1					0.89	0.88 ± 0.06		13.97 ± 1.24	0.86
G2	0.97 ± 0.02			0.59 ± 0.38	0.92	0.93 ± 0.02		2.62 ± 0.38	0.91
C3	0.91 ± 0.01				0.91	0.92 ± 0.02		0.87 ± 0.33	0.93
A4	0.97 ± 0.01				0.92	0.92 ± 0.02		1.12 ± 0.33	0.93
	0.93 ± 0.02			2.8 ± 1.10	0.92
C5	0.89 ± 0.02		82 ± 44		0.91	0.90 ± 0.02		1.28 ± 0.31	0.93
U6	0.92 ± 0.01				0.89	0.92 ± 0.02		0.51 ± 0.31	0.91
U7	0.68 ± 0.03	0.76 ± 0.02	311 ± 204		0.57	0.87 ± 0.02		3.67 ± 0.36	0.82
C8	0.83 ± 0.02			0.81 ± 0.47	0.82	0.86 ± 0.02		4.81 ± 0.38	0.86
G9	0.95 ± 0.02			1.50 ± 0.37	0.89	0.89 ± 0.02		1.36 ± 0.33	0.87
G10	1.00 ± 0.01				0.92	0.94 ± 0.02		0.56 ± 0.35	0.93
U11	0.94 ± 0.01				0.91	0.91 ± 0.02		1.26 ± 0.32	0.93
G12	0.99 ± 0.01				0.93	0.90 ± 0.02		2.95 ± 0.34	0.93
C13	0.92 ± 0.01				0.91	0.92 ± 0.02		0.80 ± 0.31	0.92
C14	0.87 ± 0.02		36 ± 29		0.90	0.84 ± 0.03	36 ± 15	0.72 ± 0.33	0.85

aThe vector C6H6 of the pyrimidines, the vector C8H8 of the purines and the adenine C2H2 vectors (italics) were analyzed with a bond length of 1.102 Å and a residue-specific carbon CSA (Table 2).

bFor the C1′H1′-vector a bond length of 1.09 Å and a carbon CSA of 45 p.p.m. were used.

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 S2 (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 Rex-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 Rex-terms between 1.9 and 11.8 s−1 (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′ (right) at temperatures of 298 K (top), of 317 K (middle) and at the onset of the hairpin melting at 325 K in NMR and 360 K in MD, respectively. The loop region is shaded in gray.

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

	C6H6/C8H8/C2H2					C1′H1′
	NMR				MD	NMR				MD
	S2	Sf2	τe [ps]	Rex [s−1]	S2	S2	Sf2	τe [ps]	Rex [s−1]	S2
G9	0.99 ± 0.03			5.5 ± 1.8	0.89	1.0 ± 0.0			11.8 ± 0.0	0.87
G10	0.97 ± 0.03				0.92	0.97 ± 0.02				0.92
U11	0.92 ± 0.04			2.9 ± 2.8	0.92	0.94 ± 0.02				0.93
A12	0.98 ± 0.02				0.93	0.89 ± 0.02				0.93
	0.93 ± 0.02				0.92
U13	0.92 ± 0.02				0.91	0.96 ± 0.02				0.93
C14	0.94 ± 0.03				0.89	0.92 ± 0.03			4.4 ± 1.6	0.92
A15	0.58 ± 0.02	0.70 ± 0.02	441 ± 121		0.49	0.76 ± 0.06	0.86 ± 0.06	827 ± 445		0.88
	0.57 ± 0.02	0.68 ± 0.02	635 ± 233		0.43
C16	0.86 ± 0.03			6.7 ± 2.2	0.87	0.94 ± 0.04			2.8 ± 1.6	0.89
G17	1.0 ± 0.0			8.9 ± 0.0	0.89	1.0 ± 0.0			2.3 ± 0.0	0.87
G18	0.99 ± 0.02				0.92	0.85 ± 0.03			7.4 ± 1.2	0.92
U19	0.90 ± 0.02				0.92	1.0 ± 0.0			1.9 ± 0.0	0.93
A20	0.96 ± 0.03				0.92	0.86 ± 0.03			3.9 ± 1.2	0.92
	0.93 ± 0.02				0.91
C21	0.89 ± 0.03			5.8 ± 2.4	0.90	0.97 ± 0.02				0.90
C22	0.88 ± 0.03			7.4 ± 2.9	0.89	0.95 ± 0.02				0.85

Parameters used are the same as in Table 3.

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 15N–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 15N–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 Rex-term to 15.3 s−1 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 Rex-terms of 7.4 s−1 and 3.1 s−1, 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 13C–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 Rex-term (20.9 ± 2.8 s−1) 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.

MD simulation

The comparison of measured and calculated order parameters of the cUUCGg 14 mer RNA hairpin shown in Figure 1 reveals an excellent overall agreement of experiment and theory. At both the temperatures, 298 and 317 K, the deviations between and are within the estimated systematic errors of ±0.05. When going from 298 to 317 K, the experimental order parameters show an average decrease of ≈0.1, while the calculated values of S2 change only marginally. Even at 360 K, only a modest decrease of the MD-derived order parameters is observed. However, particularly the central loop residues U7 and C8 become more flexible. Figure 2 compares the order parameters obtained for the 30 nt (NMR) and 14 nt (MD) uCACGg hairpin. Again, a convincing overall agreement of experiment and MD calculation is observed. In particular, the calculations nicely match the temperature-dependent changes of the nucleobase measurements.

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. Shown are the projections every 20 ps for five different time intervals: 0–10 ns (black); 10–20 (red); 20–30 ns (green); 30–40 ns (blue) and 40–50 (beige).

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 and thick lines stacking interactions between nucleobases.

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 13C 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).

DISCUSSION

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 13C 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.

CONCLUSIONS

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 13C 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 13C 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

Supplementary Data are available at NAR Online.