Structural Stability of Tandemly Occurring Noncanonical Basepairs

14 Oct 2010 - Sukanya Halder and Dhananjay Bhattacharyya*. Biophysics DiVision, Saha Institute of Nuclear Physics, 1/AF, Bidhannagar, Kolkatas700 064,...
0 downloads 0 Views 4MB Size
14028

J. Phys. Chem. B 2010, 114, 14028–14040

Structural Stability of Tandemly Occurring Noncanonical Basepairs within Double Helical Fragments: Molecular Dynamics Studies of Functional RNA Sukanya Halder and Dhananjay Bhattacharyya* Biophysics DiVision, Saha Institute of Nuclear Physics, 1/AF, Bidhannagar, Kolkatas700 064, India ReceiVed: March 30, 2010; ReVised Manuscript ReceiVed: August 19, 2010

Noncanonical basepairs have gained importance over the past few years because of their various functions in RNA biochemistry. These basepairs appear quite frequently in different double helical stems also, but whether they are stabilized by contextual pressure or act as seeds of folding is not yet clear. We have used all-atom molecular dynamics simulations to characterize the stability and functional features of a few noncanonical basepairs within two RNA double helical fragments obtained from ribosome crystal structures. It is anticipated that the noncanonical basepairs would open up spontaneously if they had appeared due to contextual pressure. However, we have found from MD simulations that the noncanonical basepairs occurring in tandem at the central regions of double helical stretches are quite stable. Analysis of basepairing parameters carried out in terms of H-bonding edge-specific axis system indicates that dynamics of the noncanonical basepairs are very similar to those of the canonical ones. The stacking parameters for dinucleotide steps consisting of noncanonical basepairs are rather unusual, but the variability patterns indicate their significant stability. The stacking freeenergy values as presumed from the distributions of structural parameters also appear to be similar for both canonical and noncanonical basepair steps. Introduction Since the description and nomenclature of different noncanonical basepairing types by Leontis and Westhof in 2001,1 studies of many RNA structures solved by NMR and X-ray crystallographic methods have indicated a large number and variety of such basepairs. Various studies have shown that many of these noncanonical basepairs are as stable as the canonical ones,2-5 i.e., Watson-Crick and wobble G:U pairs. A range of thermodynamic studies also indicates reduced but significant stability of these basepairs, although these studies are incapable of detecting the proper basepairing patterns.6-16 Noncanonical basepairing geometries involve edge-to-edge interactions between potential hydrogen-bond donor and acceptor atoms of the nucleotide bases. As indicated earlier,1 all the purine and pyrimidine bases can use three edges for hydrogen bonding: Watson-Crick edge (W), Hoogsteen or C-H edge (H), and sugar edge (S) (Figure 1). These basepairs can be in cis (C) or trans (T) geometries, depending on the relative orientations of their ribose sugars about the hydrogen-bonding pseudoaxes, giving rise to enormous theoretical possibilities.1 Throughout this paper, we have used the nomenclature system for noncanonical basepairs described by Das et al.17 In most of the cases, basepairings involve polar hydrogen bonds mediated by N-H · · · O/N and/or O-H · · · O/N type of directed interactions. Nowadays, nonpolar interactions like C-H · · · O/N are considered to be important in basepair formation. The C8-H8 bonds in purine bases and C2-H2 of adenine can be considered for such basepairing interactions. The C5-H5 bonds of pyrimidine are even more polar according to both AMBER18 and CHARMM19 force fields. In addition to the bases interacting directly, the 2′-OH groups of RNA can also form hydrogen bond with another base. Recent studies have shown that an amino * To whom correspondence should be addressed. Phone: +91-33-2337 0379, Ext: 3506/3408. Fax: +91-33-2337 4637. E-mail: dhananjay. [email protected].

group (-NH2) nitrogen can act as a hydrogen-bond acceptor due to partial pyramidalization and lone pair electrons on the nitrogen atom.20 On the other hand, DNA double helices lack in such variety in basepairing patterns due to the absence of the 2′-OH group in sugar moiety and methylation at the C5position of thymine and, in a large number of cases, of cytosine as well.21 Structural studies of natural and synthetic RNAs have revealed many occurrences of unusual basepairing within RNA duplex.22-28 The small tetraloops, such as GNRA and UNGC, are quite common structural motifs where the last basepair of the helix is generally a noncanonical basepair, such as G:A S:HT in the GNRA tetraloop. In rRNAs, an unexpected high incidence of G:A/A:G and G:A/A:A motifs in a sheared conformation can be observed.29 Examples of other unusual basepairs, like U:U W:WC and A:C W:+C, flanked by regular Watson-Crick or wobble pairs are also very common. Such perturbations in regular RNA helices by noncanonical basepairing motifs are supposed to be functionally important in adopting unusual structures or introducing anchoring sites for metals or proteins.23,27,30,31 They may have a crucial role to play in basesequence specific interactions with proteins or other ligands as they make more functional groups exposed in major or minor grooves of RNA, sometimes by widening the major groove using purine-purine basepairing or by exposing Watson-Crick edges with a maximum number of functional groups.24,28,32 Obviously, the question remains whether these noncanonical basepairs have emerged in the large macromolecules due to constraints of the remaining sections. Had these noncanonical basepairs arose in RNA due to contextual pressure or chance, these would be unstable in absence of any such pressure. Noncanonical basepairs come into play in another important area of modern biology, namely prediction of RNA structures. Unlike DNA, RNA can adopt varieties of structural elements. These are primarily double-stranded regions, which are formed by the RNA chain folding back on itself, or unpaired loops.

10.1021/jp102835t  2010 American Chemical Society Published on Web 10/14/2010

Molecular Dynamics Studies of Functional RNA

Figure 1. H-bonding patterns of noncanonical basepairs: (a) A:A s:hT, (b) A:U H:WT, and (c) A:G H:ST. Potential hydrogen bonds detected by BPFind are indicated by black dotted lines. C1′ atoms are shown as spheres.

Thus, prediction of RNA secondary structures fundamentally leads to prediction of double helical regions in the polymeric sequence. Detection of such double helical structures from base sequence is also essential to determine siRNA or miRNA sequences, which are particularly important for therapeutic applications. The main forces that hold the three-dimensional organization of folded RNA structure are basepairing through hydrogen bonds and stacking interactions of the aromatic rings from nucleotide bases, and the inclusion of noncanonical basepairs increases the number of potential interactions enormously. Theoretically, a total of 156 basepairing types are possible for RNA, whereas DNA molecules show only two. Considering two possible basepairing in DNA, we get 10 unique dinucleotide steps. Upon inclusion of the possibility of G:U wobble pair in RNA, this number increases to 21. It is needless

J. Phys. Chem. B, Vol. 114, No. 44, 2010 14029 to say that the cumulative figure becomes huge when all the other basepairing possibilities are taken into account. Unfortunately, none of the RNA structure prediction algorithms utilize these beyond the 21 stacking sequences. Systematic characterization of these basepairs and basepair stacks is essential for proper understanding of their contribution toward RNA stability. Though there have been both classical and quantum chemical studies on the structure and energetics of different noncanonical basepairs, little can be found out regarding the basepair steps containing them. Nanosecond-scale MD simulations with explicit solvent model have been successfully employed earlier in different issues to supplement the limitations of experimental methods. This technique has been proved to be very useful in structural studies of proteins, as well as nucleic acids. It has been shown earlier that basepair dynamics are generally picoseconds events4 and, thus, are detectable in nanosecond time scales. The terminal opening and bubble formation are possibly slower processes. There have been many cases where B-form double helices of DNA have been computationally investigated to determine their sequence-dependent structure, deformability, elasticity, stacking interactions, or basepairing characteristics.33-37 Molecular dynamics simulation studies on structure and dynamics of noncanonical basepairs situated at termini of RNA double helices or at the beginning of hairpin loops are also reported,38-42 but the number of such studies, which focus on the structural variations in RNA duplex introduced by noncanonical basepairs, is very limited.43,44 For this present study, we have picked out all detected doublehelical stretches from a nonredundant set of 110 RNA structures available on April 28, 2010.45 We find that there are several noncanonical basepairs present within these helices, though few are centrally located (see Table S1, Supporting Information). Among them, motifs like (C:A W:+C)::(A:C +:WC) or (C:C W:+C)::(A:A s:hT)::(A:G H:ST) have been avoided as the force field parameters for protonated basepairs are not standardized (here “::” represents stacking between the basepairs). The analysis shows that (G:A S:HT)::(A:G H:ST), (A:A s:hT)::(A:U H:WT)::(A:G H:ST), and (G:A S:HT)::(A:G H:ST)::(A:G H:ST) are the most common among various naturally occurring motifs. We have chosen the first two from 23S rRNA and 16S rRNA structures as representatives. One of the selected motifs, (G:A S:HT)::(A:G H:ST), is present within helix 44 in 16S as well as helix 25 of 23S rRNA structures of a variety of bacterial species, e.g., Thermus thermophilus, Escherichia coli, Deinococcus radiodurans, and Haloarcula marismortui. The other motif occurring within helix 24 of a small ribosomal subunit (16S rRNA), namely, (A:A s:hT)::(A:U H:WT)::(A:G H:ST), is a common feature in E. coli and T. thermophilus. We performed molecular dynamics simulation for 20 ns for each of the double-helical fragments with water and charge-neutralizing counterions at 300 K to mimic the physiological environment as to the extent possible. We have also carried out MD for an ideal A-form double-helical stretch of comparable length containing Watson-Crick basepairing only, to compare the effect of noncanonical motifs on helix stability. We have further analyzed the trajectories of basepair and local doublet parameters to understand their structure and flexibility, which can, in principle, predict their stability and free energy. Materials and Methods Coordinates of the double helices have been obtained from Protein Data Bank (PDB46). We have chosen the noncanonical motifs because of their high frequency of occurrence: (A:A

14030

J. Phys. Chem. B, Vol. 114, No. 44, 2010

Halder and Bhattacharyya

TABLE 1: Description of the Helices and RMSD Values source

sequence

1N32

PDB ID 1N32

5′-GCAAACCGG-3′ 3′-UGAUGGGCC-5′

2AW4

PDB ID 2AW4

RNA11

fiber model

5′-GUGGGAGCACG-3′ 3′-CGUCAGUGUGC-5′ 5′-GAACGCUGGAG-3′ 3′-CUUGCGACCUC-5′

noncanonical motif

residues nos. in crystal

no. of Na+ ions

ion conc. (M)

A:A s:hT A:U H:WT A:G H:ST G:A S:HT A:G H:ST

778-786 796-804

16

1.382 × 10-4

533-543 550-560

20

1.598 × 10-4

20

1.898 × 10-4

s:hT)::(A:U H:WT)::(A:G H:ST) from 1N32 and (G:A S:HT)::(A:G H:ST) from 2AW4 (Table 1). The initial structure of the 11-bp long RNA double helix (RNA11) containing only Watson-Crick basepairs (5′-GAACGCUGGAG-3′) is generated using the fiberdiffraction model based on helical symmetry.47 We have also analyzed structures of the other PDB files belonging to the families45 in which 1N32 and 2AW4 exist. The 1N32 family members are 1FJG, 1FKA, 1HNW, 1HNX, 1HNZ, 1HR0, 1I94, 1IBK, 1IBL, 1IBM, 1J5E, 1N32, 1N33, 1XMO, 1XMQ, 1XNQ, 1XNR, 2E5L, 2HHH, 2J00, 2J02, 2UU9, 2UUA, 2UUB, 2UUC, 2UXB, 2UXC, 2UXD, 2VQE, 2VQF, 2WDG, 2WDH, 2WDK, 2WDM, 2WH1, 2WH3, 2ZM6, 3D5A, 3D5C, 3F1E, 3F1G, 3HUW, 3HUY, 3KNH, 3KNJ, 3KNL, and 3KNN, and the 2AW4 family members are 1HC8, 1VS6, 1VS8, 2AW4, 2AWB, 2BH2, 2I2T, 2I2V, 2QAM, 2QAO, 2QBE, 2QBG, 2QOZ, 2QP1, 3DF2, 3DF4, 3I1N, 3I1P. In all the simulations, we have used the CHARMM27 force field19 to describe the biomolecular interactions, as this force field does not have any bias toward any particular form of nucleic acids.48 Each of the double-helical stretches has been solvated in an orthorhombic water box containing TIP3P water molecules in such a way that there remains enough solvent layer (∼15 Å) surrounding the double helices in all directions to reduce interaction between the solute molecules in adjacent cells. The systems contain Na+ ions in exact number required to maintain the electroneutrality of the whole system. We have avoided incorporation of a large number of ions in order to alleviate the polarization effects. The positioning of these counterions is generated by Monte Carlo simulation in the absence of water, only considering the interaction between the ions and the double helix concerned.49 These solvated systems have been subjected to minimization and, finally, MD simulation under no force constraint using periodic boundary conditions. We have used particle mesh Ewald (PME50) method for treating the long-range electrostatic interactions with a κ-value of 0.35. A force-switch method is applied for nonbonded interactions with a 12 Å cutoff. The systems are then slowly heated from 0 to 300 K in 30 ps. All these calculations are done using CHARMM.51 A constant pressure-temperature MD simulation (CPT Dynamics) is carried out by NAMD software52,53 using an extended system algorithm without coupling with a 1 fs integration time step. Translational and rotational movements of the center of mass have been removed at an interval of 5 ps. A production run of 20 ns has been carried out at 300 K and 1 atm pressure (piston mass of 400 amu). Trajectory analyses are performed using CHARMM and the intra- and interbasepair parameters have been calculated by NUPARM software54,55 at an interval of 1 ps. The basepairing patterns during the MD run are determined by BPFind software,17 with a cutoff distance of 3.8 Å and cutoff pseudoangle of 120°. For most of the basepairs, the intrinsic standard deviation values were reported earlier from our laboratory4 but that for the A:A s:hT basepair was not calculated. We have calculated

canonical

normal modes of vibrations, following earlier calculations using the B3LYP/6-31G** basis set56,57 by GAMESS.58 The frequencies (νj), force-constants (k), and intrinsic standard deviations (σcalc) corresponding to vibrations about the six degrees of freedom are calculated as described earlier4 for this basepair using eqs 1-4.

k ) 4π2c2ν¯ µ

for translational motions

(1)

for rotational motions

(2)

or,

k ) 4π2c2ν¯ I

where k is the force constant, c is velocity of light, νj is normalmode frequency in cm-1, µ is reduced mass, and I is the reduced moment of inertia about a specific axis along which the basepair vibration takes place. These force constants are used to predict the probability distribution F(R) (eq 3)

F(R)/F(R0) ) e-1/2k(R

- R0)2/kBT

(3)

where R0 is the equilibrium value of parameter R, kB is Boltzmann’s constant, and T is the absolute temperature. The standard deviation value corresponds to the half-width at halfmaxima of this distribution (eq 4).

σcalc ) √2ln 2kBT/k

(4)

Results and Discussion 1. Rmsd Value and Energy Fluctuation. In order to study the overall behavior of the simulations, we have first considered the issue of convergence. We have calculated root-mean-square deviations of the trajectories with respect to energy-minimized structures of the double helices. In all the cases, the rms distances are seen to oscillate around some average values without any tendency to become larger. These stable trajectories from initial stages of simulations reflect that the structures do not undergo significant conformational alterations. This is perhaps expected, as the initial models are taken from crystal structures, which have all the sequence-directed effects and special structural features to accommodate the noncanonical basepairs. Nevertheless, this signifies that the structures remain stable without the rest of the RNA or protein components of ribosomes from which these double helices are fragmented out. Similar events of early equilibrium have been observed previously for many DNA simulations.49,59-61 For the two helices, 1N32 and 2AW4, we have calculated the rmsd values both for the complete fragments as well as considering only the central regions consisting of noncanonical basepairing motifs flanked by one canonical basepair, i.e., standard A:U, G:C pairs (W:

Molecular Dynamics Studies of Functional RNA

J. Phys. Chem. B, Vol. 114, No. 44, 2010 14031

TABLE 2: Backbone Dihedrals and the Phase Angle Values (deg) for Noncanonical Motifsa base 1N32-A780(s)

1N32-A781(H)

1N32-A782(H)

1N32-G800((S)

1N32-U801(W)

1N32-A802(h)

2AW4-G537(S)

2AW4-A538(H)

2AW4-G555(S)

2AW4-A556(H)

RNA11

R

β

γ

δ

ε

ζ

χ

phase

-77.3(14.8) -81.6 -73.9(15.4) -64.6(9.7) 30.4 34.7(18.6); 166.3(8.2)b -90.2(18.1) -102.8 -72.0(21.6) -75.4(17.4) -64.8 -64.9(8.4) -76.9(11.7) -73.2 -71.9(10.8) -81.4(13.0) -58.0 -63.8(33.6) -81.1(20.5) -80.3 -76.8(11.3) -66.9(16.9) -49.4 -48.4(10.9) -72.8(17.5) 137.0 132.6(54.4) -65.3(10.8) -63.4 -62.7(10.9) -78.9

-181.1(11.1) -190.9 -188.3(5.7) -171.7(10.5) -149.4 -187.9(38.6)c -185.4(9.2) -175.5 -186.1(8.5) -184.5(11.1) -185.2 -184.7(7.5) -182.7(10.1) -190.4 -193.4(8.1) -181.7(9.7) -192.9 -185.6(9.4) -180.9(11.9) -180.9 -176.9(8.1) -183.6(12.1) -222.5 -199.5(10.8) -193.7(18.2) -98.9 -122.5(18.9) -185.4(9.6) -213.0 -208.4(8.2) -188.2

58.8(9.4) 69.2 62.4(13.2) 57.0(11.5) -69.6 -75.6(5.8); 168.9(5.2)b 73.3(8.6) 58.5 52.8(7.6) 60.4(9.3) 54.1 53.6(4.9) 60.4(10.3) 51.2 52.2(8.0) 60.8(9.8) 55.2 61.9(17.4) 61.2(10.6) 62.7 56.0(6.0) 62.1(10.2) 55.0 52.9(9.5) 63.2(10.1) 171.4 163.5(28.0) 62.6(9.3) 52.6 54.4(8.5) 63.6

76.9(5.2) 92.0 87.0(3.7) 69.7(5.2) 85.5 88.7(3.5) 73.9(5.3) 82.5 80.8(2.2) 78.6(5.3) 86.9 87.5(1.8) 77.6(4.9) 80.0 81.3(4.7) 78.1(5.0) 86.3 85.8(2.1) 85.5(22.5) 84.3 83.6(1.4) 79.5(5.3) 79.3 81.4(1.2) 77.5(5.1) 85.6 88.9(4.8) 80.5(5.1) 84.5 86.2(2.1) 77.8

-179.8(6.9) -154.8 -166.6(16.8) -147.4(17.6) -170.3 -132.7(25.6) -138.5(24.7) -139.1 -136.5(6.3) -173.8(8.1) -150.9 -151.0(7.8) -153.2(11.7) -153.0 -157.3(6.0) -131.4(16.7) -128.8 -141.7(8.8) -176.9(16.7) -138.5 -169.5(11.8) -140.7(23.2) -151.7 -144.9(8.5) -181.7(7.6) -116.0 -136.4(15.4) -127.4(15.2) -137.2 -130.6(11.2) -151.6

-96.1(7.1) -121.5 -104.9(19.9) -26.4(16.9) 0.6 -32.4(19.7) -37.0(30.6) -53.5 -53.8(5.1) -89.9(7.9) -102.4 -106.4(7.8) -64.5(9.5) -78.7 -73.7(5.8) -46.1(16.0) -57.2 -47.5(5.7) -109.9(25.5) -131.3 -111.5(9.5) -35.5(27.6) -68.2 -68.3(5.1) -98.7(8.5) -155.2 -153.4(6.0) -47.5(15.8) -51.3 -51.0(8.0) -63.3

-150.5(7.8) -151.8 -144.4(9.6) -163.5(8.2) -160.1 -167.4(9.3) -169.8(7.2) -160.4 -160.9(8.5) -151.4(6.7) -156.5 -153.0(4.6) -163.1(6.4) -157.4 -159.6(6.5) -163.2(6.7) -161.9 -158.2(7.4) -155.3(13.8) -157.5 -155.0(3.5) -167.2(7.7) -159.7 -159.9(4.5) -148.6(8.7) -181.9 -171.0(13.8) -164.6(7.4) -173.9 -174.8(3.3) -157.9

11.5(7.7) 27.6 17.6(6.0) 22.9(6.8) 5.8 6.1(6.6) 16.7(8.2) 17.7 18.3(3.2) 10.3(7.9) 8.0 9.7(3.5) 12.2(7.3) 17.1 20.4(5.4) 8.7(6.9) 10.2 13.0(5.3) 27.8(56.5) 15.7 12.4(3.2) 6.6(7.3) 32.8 18.4(6.5) 8.3(7.4) 29.5 36.6(12.2) 6.9(7.3) 17.2 15.0(5.2) 14.5

a Standard deviations are given within parentheses; average values have been calculated leaving the first 2 ns runtime. Values present in the second and third lines give backbone dihedrals as seen in initial crystal structures and an average from the crystal structure family of 1N32 and 2AW4. b For this nucleotide, we observe a bimodal distribution in R-γ space, corresponding to (g+,g-) and (t,t)-conformations. Average values and standard deviations for both are reported here. c A larger range of β torsion in the trans region only leads to the higher standard deviation.

WC) or wobble G:U pairs (W:WC), on either side (bases 779-783 and 536-539 along with their complementary residues in 1N32 and 2AW4, respectively). The average values of structural fluctuation in Cartesian coordinate for the noncanonical motifs in these two simulations are 0.85 Å (0.09) and 1.45 Å (0.22), respectively (standard deviation values are given within parentheses; see Table S2 and Figure S1 in the Supporting Information). Both the potential and kinetic energy components of the two helical systems display a stable nature; the overall double helical structures are also maintained throughout the time propagation with moderate terminal melting. Such terminal melting phenomena are the inherent nature of nucleic acids in solutions also and have been frequently observed previously.62-64 Our analysis of hydration or ion-binding indicates the presence of a rigid sodium ion in the major groove of 1N32, near the negatively charged phosphate group of C784, but no other rigid water or ion could be located near the RNA. 2. Backbone Conformation within Noncanonical Motifs. A detailed study of different backbone dihedrals (R, β, γ, δ, ε, ζ) and glycosidic torsion angle (χ) along with the sugar pucker pseudorotation phase angle (P) has been performed in order to get a general idea about the conformational changes introduced in the sugar-phosphate backbone during the MD run (Table 2). The average values for most of the torsion angles obtained from our simulations lie within allowed regions in conformational space. The data also corroborates with the crystal structure analysis for A-form DNA helices.65 The angles involved in phosphodiester linkages, R and ζ, show the usual g--conformation, while β and ε assume the sterically favorable transorientation. The R and γ torsions are seen to be anticorrelated, with the (g-,g+)-combination being most popular. The pseudorotation phase angle, P, describes the sugar puckering most

effectively. Its values by and large remain in the C3′-endo regions for most of the sugars. A closer look at the sugar-phosphate backbone conformation suggests that whenever a base uses its Hoogsteen or sugar edge for H-bonding, backbone dihedrals show a deviation from the usual values. This distortion is overcome either by the variability being distributed over all the backbone torsions or by maintaining correlation between a pair of them. A little deviation in ζ-angle is observed for the noncanonical basepairs using their Hoogsteen or sugar edge for H-bonding. For example, the average ζ-values are -27° and -37° for the A781 and A782 of 1N32, where the bases use their Hoogsteen edges. Similar is the case for A538 of 2AW4 with an average ζ of -38°. On the other hand, use of the sugar edge results in more negative value of ζ compared to the normal g--orientation, as observed for A780 (-96°) and G800 in 1N32 (-90°) and G537 (-118°) and G555 (-99°) in 2AW4 (Table 2). While we find that the MD-averages of backbone dihedrals maintain a regular pattern, the torsion angle values obtained from initial crystal structures sometimes fall in different domains. The ζ angle for A781 in the crystal structure (1N32) indicates an unfavorable eclipsed conformation around the O3′-P bond with a value of 0.6° (Table 2). Similar instances of eclipsed conformation around the O3′-P bond have been found in many other crystal structures of T. thermophilus 16S rRNA (1N32 family), e.g., 1FJG. Though some crystal structures of this group show gauche-conformation for this dihedral, the general trend of eclipsed or near eclipsed conformation lowers the magnitude of the average ζ for this ensemble to -32.4°. We observe a magnesium ion with broken coordination with this phosphate group, which also coordinates with another part of the RNA chain. Presumably, a strain arising due to orientation

14032

J. Phys. Chem. B, Vol. 114, No. 44, 2010

of the nearby chain causes an anomalous ζ-value in some of the structures of this family. During simulation, it alters to a more favorable staggered form within the first ∼400 ps in absence of the folded environment and retains regular features for the remaining production phase. Likewise, anomalous values are observed in 2AW4 for G555 with R ) 137°, β ) -99°, γ ) 171°, and ζ ) -155°, although G555 does not form any special contact in the crystal structure. It should be noted that these variations are not restricted to noncanonical basepairs, as a few Watson-Crick basepairs also exhibit torsion angle values outside the normal range. For example, ζ for G543 takes up values like -3.7° and 5.4° for crystal structures 3DF4 and 3I1P, respectively, in the 2AW4 family, though the other members display the usual g--conformation. Usually, such backbone variations reorient the phosphate groups, which is presumably required in the compact folding of RNA double helix with the rest of the large macromolecule and proteins in ribosome crystal structures. During simulations, the conformational strains are released, as the simulated RNA duplex remains surrounded by a uniform distribution of water and counterions only. Almost all the 3′-terminal bases show sugar repuckering from the C3′-endo to the C2′-endo conformation for significant periods of time during the simulations. The most significant variation has been observed for G543 in 2AW4 (26.8%) (Figure 2b), whereas U804 in 1N32 and C560 in 2AW4 show transitions to C2′-endo puckering to a lesser extent, 6.5% and 4.5%, respectively. In 2AW4, the longest span through which C2′endo puckering is retained for G543 is ∼2.5 ns near the end of simulation time. This appears to be an inherent characteristics of most A-form helices, as the canonical helix RNA11 also shows C2′-endo puckering at either or both the 3′-termini. Crystal structures of several A-DNA helices, such as 28DN, 2A7E, 2D95, 2PKV, 317D, 339D, also display similar tendency to attain C2′-endo conformations at either or both of their 3′ends. Moreover, this seems to be independent of force-field parameters, as simulations carried out in AMBER force-field also show similar conformational transitions of the sugar moieties (unpublished results). It has been suggested before that the absence of a phosphate group at the 3′-termini allows the sugar moiety to conform to C2′-endo-like puckering, as no steric clash appears between the 2′-OH group and the next phosphate group or base. The G537 of 2AW4, involved in G:A S:HTpairing, also undergoes repuckering to the C2′-endo conformation (∼20.5%) for almost 4.3 ns time-span, although it is situated near the center of the double helix (Figure 2c). Concomitant to this, the backbone takes up a (g-,t)-conformation in the ε-ζ space with the glycosidic angle varying in high anti region, which is typical for BII-form helices. Purine-purine mismatches in B-DNA structures are also known to have similar preference for BII-like backbone phosphate conformation.66,67 These features may appear to avoid steric clash with the 3′-phosphate group and to maintain a stable basepair. However, other A:G pairs lying within noncanonical motifs do not show the alternate form. Limited sampling allows us to visualize this transition for one of the G:A S:HT pair for ∼4 ns, indicating this to be a less likely event. Presumably longer simulation may show such features for the other G:A S:HT basepairs also. 3. Analysis of the Canonical Basepairs. Tables 3 and 4 give the average values of intrabasepair (buckle, open angle, propeller twist, stagger, shear, stretch) and interbasepair stacking (tilt, roll, twist, shift, slide, rise) parameters for the two RNA double helices with noncanonical motifs. The definitions and nomenclatures of these parameters can be obtained from the reports of IUPAC/IUB conventions.68,69 It is evident from the data that

Halder and Bhattacharyya

Figure 2. Plot of pseudorotation phase angles (P) for the 3′-terminal bases of (a) 1N32 and (b) 2AW4 showing the duration for C2′-endo puckered states. Yellow curves are for the first strand (G543 in 2AW4, G786 in 1N32), black for the second strand (C560 in 2AW4, U804 in 1N32). Plot c presents phase-angle variation for G537 (2AW4). It shows the C2′-endo form to be present for the first ∼4 ns.

the basepairs situated at either terminus of double helices show unusual average values for the structural parameters with larger standard deviations because of their unstable nature. The other canonical basepairs located in the middle show regular behavior in terms of basepair parameters. The regular Watson-Crick basepairs in central regions of the helices have small values of buckle, opening, and stagger; propeller twists are negative and stretch values remain around 3 Å, all with small standard deviations. Shear and open-angles of G:U wobble basepairs are systematically larger due to their hydrogen-bonding criterion. Although the A:T W:WC basepairs in B-DNA show larger negative values of propeller twist as compared to G:C W:WC, we do not observe any such correlation in the simulated RNA

Molecular Dynamics Studies of Functional RNA

J. Phys. Chem. B, Vol. 114, No. 44, 2010 14033

TABLE 3: Intrabasepair Parameters for the Noncanonical Helices 1N32 and 2AW4a basepairs 778G:804U W:WC 779C:803G W:WC 780A:802A s:hT 781A:801U H:WT 782A:800G H:ST 783C:799G W:WC 784C:798G W:WC 785G:797C W:WC 786G:796C W:WC

buckle

opening

6.6(13.0) 4.2 -7.5(8.2) 13.8(9.8) 7.2 -0.5(8.8) -17.9(7.8) -15.3 -23.2(6.4) -10.9(8.8) 0.9 -4.1(7.9) -13.7(8.1) -16.6 -16.2(7.0) 7.4(8.0) -2.4 16.3(7.4) 8.1(9.7) 1.9 6.3(5.8) -12.0(11.6) -6.0 -8.1(9.4) -6.0(14.8) 18.1 -22.4(8.4)

-8.2(17.0) 0.5 0.9(3.4) 1.7(4.6) 0.6 2.9(3.1) 2.2(4.2) 11.1 0.7(5.9) -0.4(4.5) -5.2 0.6(5.3) 6.7(4.4) 14.7 12.3(5.0) 8.0(11.6) 4.5 0.4(3.5) 2.1(6.2) -4.1 4.2(2.6) 4.9(11.8) -2.2 3.4(4.5) 5.3(17.4) -4.4 0.2(3.0)

0.9(21.2) -2.8 1. 5(1.7) 12.0(15.4) 0.1 4.6(2.3) -9.5(10.6) 14.3 -3.4(2.4) -2.8(9.1) 3.4 -3.7(4.0) -2.4(8.3) 33.9 -0.9(1.9) -6.6(8.9) -15.6 1.4(1.7) 5.2(10.1) -8.3 -2.7(3.3) 14.7(14.7) 13.5 3.2(2.6) -12.4(14.5) 5.6 -0.3(1.8) -1.2(12.9) 11.5 1.2(3.9) -37.7(55.2) -12.0 -4.1(2.9)

28.8(38.9) 0.3 -2.9(1.2) 23.2(33.0) -11.3 -3.2(2.6) 26.6(24.2) -1.9 3.5(1.9) 3.4(5.5) -0.9 -0.6(0.6) -7.2(5.3) 20.0 -13.2(2.3) 5.8(6.6) -6.3 14.9(3.1) 31.9(23.8) -10.7 0.2(1.6) 10.2(18.2) 4.7 -0.7(1.0) 8.0(22.2) -0.7 1.2(3.2) 46.9(26.3) -9.3 -0.7(1.2) 30.4(34.8) 15.3 -1.2(2.6)

propeller 1N32

stagger

shear

stretch

-8.1(13.3) -6.6 -9.1(6.7) -14.1(7.3) -18.5 -17.4(7.1) 5.2(5.8) -6.8 4.0(5.0) -7.2(7.0) 0.4 -10.4(7.6) 0.3(6.0) 3.4 -11.2(5.4) -7.9(6.9) -5.1 -15.8(6.6) -10.1(7.5) -23.2 -2.8(6.1) -10.8(10.0) -12.7 -11.3(7.9) -7.2(11.5) 4.6 -8.5(5.6)

-0.2(0.9) -0.1 0.0(0.6) -0.1(0.4) -0.2 -0.2(0.4) 0.5(0.4) 0.2 0.2(0.6) 0.2(0.4) 0.1 0.0(1.0) 0.4(0.4) 0.4 0.0(0.4) 0.0(0.4) 0.0 -0.4(0.5) -0.4(0.4) -0.1 -0.1(0.3) -0.4(0.4) 0.4 -0.4(0.5) -0.2(0.6) 0.4 -0.3(0.3)

-3.1(1.7) -2.8 -2.1(0.3) 0.4(0.4) 0.6 0.2(0.3) 2.3(0.2) 2.7 2.6(0.4) 0.1(0.3) 0.1 0.0(0.7) 2.7(0.3) 2.6 2.8(0.3) 0.6(0.7) 0.8 0.0(0.3) 0.3(0.4) 0.4 0.0(0.4) -0.2(0.6) -0.7 0.0(0.4) -0.5(0.9) 0.3 0.0(0.3)

3.0(0.7) 2.8 2.8(0.1) 2.9(0.1) 2.8 2.8(0.1) 2.9(0.2) 2.9 2.8(0.2) 2.8(0.1) 2.8 2.9(0.3) 3.4(0.2) 3.4 3.3(0.3) 3.0(0.2) 2.9 2.8(0.2) 2.9(0.2) 2.9 2.9(0.1) 3.0(0.2) 2.9 2.8(0.2) 3.1(0.6) 3.0 2.8(0.1)

1.5(16.1) 1.4 -3.2(1.8) 0.4(13.3) 0.9 -3.9(4.4) -3.2(10.1) -16.6 -3.4(2.0) -4.8(8.0) 9.5 -2.9(1.1) 4.3(8.1) 33.9 -0.6(1.1) 6.6(6.5) -1.1 4.8(4.7) -9.2(9.2) -6.1 0.6(2.4) -14.5(11.2) -9.5 -1.0(3.8) -8.4(11.0) -24.5 -3.0(4.5) -3.8(12.5) -1.9 -3.0(4.5) -14.0(29.4) 3.3 -0.8(3.9)

-1.0(1.4) -0.2 0.0(0.2) -0.8(1.1) -0.1 0.0(0.2) -0.5(0.6) -0.2 -0.2(0.2) -0.3(0.4) 0.0 -0.2(0.1) -0.2(0.4) -1.3 0.0(0.2) -0.2(0.4) -0.6 0.1(0.1) 0.0(0.7) 0.0 0.0(0.1) -0.4(0.5) -0.2 0.1(0.2) -0.7(0.9) -0.8 0.0(0.1) -0.1(1.0) -0.2 0.0(0.1) -1.5(2.3) -1.2 -0.2(0.4)

-0.8(2.4) 0.4 -0.1(0.1) 3.0(0.7) 2.7 2.4(0.1) -2.7(0.6) -2.7 -2.3(0.1) -0.4(0.4) -0.3 -0.1(0.1) 2.2(0.3) 1.9 2.5(0.2) 2.0(0.4) 2.0 2.7(0.1) -2.9(0.6) -2.7 -2.4(0.1) 0.4(0.7) 0.4 0.2(0.1) 0.1(0.4) 0.0 0.0(0.1) 3.1(2.5) 0.9 0.3(0.1) -1.8(3.9) -1.5 -0.1(0.1)

3.8(2.0) 3.0 2.9(0.0) 3.2(0.8) 2.9 2.8(0.0) 3.3(0.5) 2.8 2.8(0.1) 2.9(0.1) 2.9 2.9(0.0) 3.4(0.2) 4.2 3.2(0.1) 3.4(0.2) 3.8 3.3(0.0) 3.4(0.5) 2.9 2.9(0.0) 3.0(0.4) 2.9 2.9(0.0) 2.9(0.7) 2.8 2.8(0.1) 3.3(1.4) 3.0 2.9(0.0) 4.2(3.5) 2.8 2.9(0.1)

2AW4 533G:560C W:WC 534U:559G W:WC 535G:558U W:WC 536G:557C W:WC 537G:556A S:HT 538A:555G H:ST 539G:554U W:WC 540C:553G W:WC 541A:552U W:WC 542C:551G W:WC 543G:550C W:WC

a Standard deviation values are given within parentheses; average values have been calculated leaving the first 2 ns runtime to exclude the probable equilibration period. It was seen that equilibration was attained at the early stages of simulation. Values in the second and third lines are those for initial crystal structures and averages from crystal structure families of 1N32 and 2AW4.

structures. In B-DNA, the larger negative propeller values are attributed to the formation of strong cross-strand hydrogen bonds between successive basepairs, which is impossible in RNA due to the sliding motions. This presumably leads to equivalent and slightly smaller propeller for all the canonical basepairs. It has been previously speculated that A-form helices lack diagonal cross-chain purine clashes as the basepairs are displaced more from the center of the helix along the pseudodyad axis, with larger magnitude of X-displacement or the associated slide. As the steric clash between purines in purine-pyrimidine steps are overcome by large negative slide of around -1.5 Å in

A-form structures, they can easily adopt any roll value. In B-form structures without such usual negative slide, however, a larger negative roll is required to reduce the cross-strand purine clash. This negative slide also increases the overlap between bulky purine bases from opposite strands for pyrimidine-purine steps.70 Thus, pyrimidine-purine steps usually attain larger roll values to avoid the cross-strand purine clash, which is in agreement with Calladines’ rule.70 The average values of roll for the purine-pyrimidine steps are 12.9°, 10.4° for (AC).(GU) in 2AW4 RNA11, respectively, and 5.3° for (GC).(GC) in RNA11, respectively, whereas those for the pyrimidine-purine steps are 14.4°, 11.4° for (CG).(CG) in 1N32 RNA11, respec-

14034

J. Phys. Chem. B, Vol. 114, No. 44, 2010

Halder and Bhattacharyya

TABLE 4: Interbasepair Step Parameters for the Noncanonical Helices 1N32 and 2AW4a base-doublet

tilt

roll

778G:804U::779C:803G

-5.7(3.2) -4.6 -0.1(1.8) 0.2(2.2) 4.8 2.0(2.8) -1.6(4.1) 1.7 -10.6(6.5) 3.5(3.4) 0.2 -5.0(3.4) -1.6(2.0) 1.2 -1.2(1.5) 0.9(2.6) -3.5 -3.4(2.3) -0.0(3.2) 7.2 1.0(1.4) 1.4(3.8) 13.4 1.2(3.7)

8.7(7.1) 14.7 8.3(2.9) 6.8(4.9) 12.2 5.9(4.2) 10.3(5.3) 10.6 -1.7(13.0) 12.2(5.6) 13.9 -10.1(4.6) 11.8(5.3) 8.2 7.6(4.3) 7.6(5.2) 8.4 12.2(5.4) 14.4(7.5) 0.2 8.4(3.0) 7.8(6.2) -8.2 9.9(4.7)

twist

shift

slide

rise

26.0(5.0) 27.6 25.3(1.8) 18.0(3.2) 21.1 12.1(1.5) 67.5(2.4) 61.1 91.1(4.7) 49.9(2.2) 58.8 38.0(2.4) 15.5(3.3) 18.0 15.1(1.4) 30.0(3.6) 27.4 28.1(1.7) 33.1(4.0) 23.6 28.1(1.9) 30.5(4.8) 36.4 34.2(1.4)

-0.4(0.6) -0.8 -0.2(0.2) -0.4(0.4) -0.3 -0.6(0.3) -1.3(0.3) -1.3 -0.2(1.5) 1.2(0.3) 1.4 -1.9(0.5) 0.2(0.5) 0.5 0.3(0.2) 0.0(0.5) -1.6 -0.8(0.2) 0.1(0.9) 2.4 0.7(0.3) -0.0(0.7) 1.2 0.1(0.3)

-2.0(0.6) -2.0 -1.1(0.2) -3.1(0.3) -2.8 -0.4(0.2) 0.7(0.3) 0.6 -1.3(1.0) -1.2(0.3) -1.3 1.9(0.4) -2.8(0.3) -2.6 -0.7(0.2) -2.1(0.4) -1.6 -2.2(0.2) -1.8(0.5) -1.7 -2.2(0.2) -1.9(0.5) -1.6 -2.3(0.2)

3.2(0.2) 3.1 3.1(0.1) 3.9(0.3) 3.6 4.2(0.3) 2.9(0.3) 2.9 2.6(0.3) 3.1(0.3) 3.2 3.4(0.2) 3.0(0.3) 3.1 2.8(0.3) 3.2(0.2) 3.1 3.4(0.1) 3.5(0.2) 3.4 3.4(0.2) 3.2(0.2) 3.5 3.2(0.1)

12.7(11.3) 33.3 27.7(1.3) 41.1(4.9) 43.0 46.5(1.0) 18.2(5.2) 23.9 25.8(0.8) 11.9(3.6) 2.5 13.4(1.2) 70.8(4.2) 77.2 80.5(2.0) 19.7(4.1) 25.1 11.0(1.6) 17.1(6.3) 24.3 24.4(1.2) 32.7(4.3) 35.3 30.8(1.3) 14.0(10.2) 22.4 28.6(1.9) 24.9(31.4) 39.5 38.3(2.0)

0.1(1.1) -1.0 -0.4(0.2) 0.3(0.7) 1.4 0.8(0.1) -1.2(0.7) -0.5 0.3(0.2) -0.6(0.4) 0.6 0.3(0.2) -0.3(0.6) -0.5 1.1(0.2) 0.4(0.6) 0.3 -0.8(0.2) -1.1(0.7) -0.3 -0.8(0.2) -0.2(0.9) -0.8 -0.2(0.2) 1.4(1.0) 1.0 0.4(0.3) -0.0(0.9) 1.5 0.9(0.2)

-2.8(0.9) -1.8 -1.0(0.1) -2.0(1.0) -1.1 -2.1(0.2) -2.6(0.4) -2.0 -1.8(0.2) -3.0(0.4) -1.9 -0.6(0.2) 0.6(0.3) 1.2 -2.0(0.2) -2.7(0.5) -2.8 -0.4(0.1) -2.3(0.8) -1.4 -1.2(0.1) -1.9(0.6) -1.4 -1.9(0.2) -2.1(0.8) -1.3 -1.7(0.2) -0.6(4.0) -1.1 -2.0(0.3)

3.4(0.3) 3.4 3.3(0.1) 3.8(0.5) 3.2 2.9(0.1) 3.4(0.2) 3.2 3.2(0.1) 3.5(0.3) 3.0 3.3(0.1) 3.3(0.4) 4.1 3.1(0.2) 3.4(0.3) 3.2 3.6(0.1) 3.3(0.2) 3.4 3.4(0.1) 3.6(0.3) 3.7 3.1(0.1) 3.3(0.3) 3.2 3.3(0.1) 3.7(2.2) 3.6 3.2(0.1)

1N32

779C:803G::780A:802A 780A:802A::781A:801U 781A:801U::782A:800G 782A:800G::783C:799G 783C:799G::784C:798G 784C:798G::785G:797C 785G:797C::786G:796C

533G:560C::534U:559G 534U:559G::535G:558U 535G:558U::536G:557C 536G:557C::537G:556A 537G:556A::538A:555G 538A:555G::539G:554U 539G:554U::540C:553G 540C:553G::541A:552U 541A:552U::542C:551G 542C:551G::543G:550C

-1.7(4.1) 0.2 0.1(0.6) 1.4(4.6) 6.2 0.9(1.3) -1.5(2.6) 3.6 0.9(1.4) 0.3(2.1) 4.6 1.8(0.9) -1.9(6.2) -14.5 9.4(1.8) -2.0(2.3) -3.2 -1.9(0.9) -1.7(2.5) -2.9 -2.5(1.7) -0.2(3.9) 1.5 -1.8(0.8) 0.5(3.1) 7.1 0.5(1.4) 5.9(13.5) 2.8 4.4(3.1)

2AW4 10.4(7.5) 5.5 2.8(2.7) 3.8(9.2) 11.8 19.3(2.4) 4.5(5.7) 14.8 7.1(2.7) 8.7(5.4) 27.3 12.1(2.4) 0.3(4.6) -2.4 0.1(2.3) 7.2(5.8) 9.1 10.2(2.1) 9.0(5.6) 1.4 5.6(2.6) 13.7(9.4) 9.7 13.6(2.3) 10.4(7.3) 20.7 6.1(2.5) 4.8(19.6) 5.0 9.9(5.0)

a Standard deviation values are given within parentheses; average values have been calculated leaving the first 2 ns runtime to exclude the probable equilibration period. It was seen that equilibration was attained at the early stages of simulation. Values in the second and third lines are those for initial crystal structures and averages from crystal structure families of 1N32 and 2AW4.

tively, and 13.7°, 12.8° for (CA).(UG) in RNA11 and 2AW4, respectively (Tables 4 and S3, Supporting Information). 4. Noncanonical Basepairing MotifssStructure and Stability. Among the six intrabasepair parameters, buckle, propeller twist, and stagger reflect the overall nonplanarity of a basepair, whereas opening, shear, and stretch define the H-bonding pattern and proximity of the bases, thus describing H-bonding strength. In our case, the noncanonical basepairs are situated near the

central regions of double helices, and hence the out-of-plane movements, in the direction of buckle, propeller and stagger, are likely to be restricted because of the stacking interaction with adjacent basepairs on either side. Earlier quantum chemical studies on different noncanonical basepairs occurring in functional RNA had shown that the movements for most of them is favorable along the direction of propeller and buckle, but rather restricted for open, shear, and, in particular, stretch, as they are

Molecular Dynamics Studies of Functional RNA

J. Phys. Chem. B, Vol. 114, No. 44, 2010 14035

Figure 3. Open angle variation for (a) 780A:802A s:hT (1N32), (b) 781A:801U H:WT (1N32), (c) 782A:800G H:ST (1N32), (d) 537G:556A S:HT (2AW4), and (e) 538A:555G H:ST (2AW4). The open angle for the 781A:801U-basepair of 1N32 lies outside the normal range for the first ∼400 ps.

associated with distortion of hydrogen bonding between the bases.4 The vibration along stretch direction was found to be happening at high frequency, indicating larger energy cost. This is true for our data also (Table 3). Thus, the noncanonical motifss(780A:802A s:hT)::(781A:801U H:WT)::(782A:800G H:ST) in 1N32 and (537G:556A S:HT)::(538A:555G H:ST) in 2AW4sexhibit high stability in terms of basepair parameters within the double helices. The plots in Figure 3 do not show much fluctuation of open-angle, and average values are close to normal. Even the C-H · · · O-mediated basepair, 780A:802A s:hT of 1N32, does not display any tendency to open up or become distorted, indicating high significance of C-H · · · O/Nlike nonpolar interactions in stabilizing the biomolecular structures. There is, however, a possibility of three hydrogen bonds in this basepairsC2-H2 · · · N7, N3 · · · H6-N6, and O2′ · · · H6-N6srendering significant stability to the system.

Although this is a less stable energy-minimized basepair in terms of quantum chemical interaction energy alone,4 its stacking with respect to its neighbors probably keeps it in a stable equilibrium geometry. It has been previously speculated that basepairs involving sugar edge are usually nonplanar in order to avoid a steric clash involving the bulky ribose sugar moiety.4,55 Here also, we find that the average values of buckle for A:A s:hT and A:G H:ST basepairs in 1N32 and 2AW4 are in general larger than their canonical counterparts. The values of shear as well are higher than usual for these two types, as they use their sugar edges for H-bonding to the paired bases. The propeller twist assumes large positive values for the trans-basepairs, as reported earlier.55 The near-universal tendency of negative propeller twist does not remain true for the trans basepairs, and these are positive more often than negative. We notice very small standard deviation values for stretch and shear for these

14036

J. Phys. Chem. B, Vol. 114, No. 44, 2010

Halder and Bhattacharyya

TABLE 5: Standard Deviation Values of Intrabasepair Parameters As Calculated from Normal-Mode Analysis for the Noncanonical Basepairs Present in 1N32 and 2AW4a basepair

buckle

open angle

propeller

stagger

shear

stretch

A:A s:hT A:U H:WT A:G H:ST

9.2; 7.8 10.6; 8.8 9.7; 8.1 -; 8.3 -; 8.9

-; 4.2 4.8; 4.5 4.8; 4.4 -; 5.3 -; 6.6

16.0; 5.8 13.9; 7.0 10.3; 6.0 -; 8.1 -; 6.5

1.0; 0.4 0.5; 0.4 0.6; 0.4 -; 0.4 -; 0.4

-; 0.2 0.4; 0.3 0.4; 0.3 -; 0.3 -; 0.4

-; 0.2 -; 0.1 -; 0.2 -; 0.2 -; 0.2

a

First set of values are those calculated from normal-mode frequency and the latter ones are from simulation data. For A:G H:ST, data from both 1N32 and 2AW4 are reported in different lines, respectively.

noncanonical basepairs, indicating their structural rigidity and high stability during simulations, even if the geometry may appear distorted with respect to canonical basepairs. Roy et al.4 had earlier reported intrinsic standard deviations of different basepairs from normal-mode analysis using density functional theory calculations. In their attempt to compare these with respect to the standard deviations calculated from crystal structure data, they noted that the intrinsic standard deviations (σ) are always smaller than the corresponding observed ones. The crystallographic data, however, were taken from diverse environments and it appeared likely. In our samples, the variability of basepair orientation parameters are from identical environments and hence more similar to the isolated basepairs in vacuum. The frequencies and standard deviations of A:A s:hT basepair were not reported earlier and hence we have calculated that following the reported protocol. A comparison of the intrinsic standard deviations (from ref 4 and calculated for A:A s:hT) and the standard deviations calculated from MD simulations for different basepairs are shown in Table 5. We note better agreement between the intrinsic flexibilities of the basepair parameters with respect to those calculated from simulations. The exception is for propeller twist, which shows smaller variability in double helical environment. In both cases, the inplane movements are more restricted compared to those along out-of-plane directions, e.g., buckle, propeller, and stagger, indicating high stability of these basepairs against deformation or rupture of H-bonds. We always notice that the normal mode vibrations along stretch direction are coupled with bond stretching and bond bending motions. Thus, the frequencies for stretch vibrations could not be confirmed to calculate the force constants and expected standard deviations. Likewise, A:A s:hT shows high rigidity along open angle and shear. The intrinsic variability for these movements appears to be high, as they are coupled with energetically less favorable bond stretching and bond bending motions. This is further supported by small standard deviation values obtained from simulation data (Table 5). However, the intrinsic fluctuations are in general larger than that obtained from simulations; this is predictable, as the basepairs remain more constrained within the double-helical stretches. In the case of the standard Watson-Crick basepairs, the values of basepair parameters calculated by NUPARM are in very good agreement with those obtained from 3DNA.71 However, there are some gross differences for the noncanonical basepairs. Small magnitudes of parameters in NUPARM definition indicate stable basepairing even for the noncanonical ones, while the 3DNA values for them are unusually large. This is because 3DNA uses one and the same basepair axis system for both canonical and noncanonical basepairs; but NUPARM follows an edge-specific axis system that is related to the basepairing pattern.55 Due to the choice of appropriate axis system, the true nature of deformation is evident from NUPARM

Figure 4. View of buckle (top panel) and propeller (bottom panel) angle for 780A:802A s:hT in 1N32. The views are taken along the basepair long axis and short axis for buckle and propeller, respectively. The A780 base is shown in the X-Z plane to highlight the movements of the A802 base in the X-Y plane.

values, which was a recommendation of the IUPAC/IUB convention.69 For example, a true large buckle value has been calculated by NUPARM for the basepair 780A:802A s:hT in 1N32 crystal structure (-19.1°) (Figure 4a). This basepair has a small propeller motion of 3.3°, which can be estimated from the perpendicular view (Figure 4b). However, 3DNA indicates this basepair to be highly propeller twisted (-19.1°) with a small buckle (1.3°). Deviations are also prominent along the direction of shear and open angle (Table S4, Supporting Information). For example, 3DNA suggests 781A:801U of 1N32, which is a H:WT basepair, to be present in opened-up state in crystal with an opening of -92.9°, whereas its open-angle value as calculated by NUPARM is 3.5°, indicating a stable basepaired structure. In accordance with this, the basepair remains stable throughout the simulation. The local doublet parameters for the dinucleotide steps containing noncanonical basepairs have also been analyzed. As there is no previous report on stacking parameters of noncanonical steps, we have compared the values with the initial crystal structures only. Both sets of values are in good agreement, except for 780A:802A::781A:801U step in 1N32. For this basepair step, the values of tilt, shift, and slide in the crystal lie outside the range of MD structures. It may be mentioned that BPFind17 results indicate that the 781A:801U H:WT basepair in 1N32 is part of a triplet where C1514 of the 16S chain is the third base, being situated within another double helix. Thus, 780A:802A::781A:801U in its initial structure makes important contact for molecular recognition, inducing an unusual conformation in the crystal, represented by a large open-angle and shear value (Figure 3b). Absence of such pressure of another helix through base triplet formation and

Molecular Dynamics Studies of Functional RNA

J. Phys. Chem. B, Vol. 114, No. 44, 2010 14037

Figure 5. Twist in the noncanonical regions. Structures a and c show the increased stacking at the junction of canonical and noncanonical regions, while structures b and d show near-perpendicular arrangement of long axes within a doublet with twist values close to 90°: (a) 779C:803G::780A:802A (1N32), (b) 780A:802A::781A:801U (1N32), (c) 536G:557C::537G:556A (2AW4), and (d) 537G:556A::538A:555G (2AW4). The long axes are drawn by joining C6-C8 (Y:R) or C8-C8 (R:R) atoms. C1′-atoms, shown as spheres, are labeled with residue number; the top bases are shown in CPK color.

Mg2+-mediated phosphate-phosphate interaction (as discussed in section 2) relaxes the structure during the first ∼400 ps to a regular one, which persists for rest of the simulation, and this step behaves normally, as reflected by the small variation of step parameters. The general feature of large positive roll values in A-form double helices is followed by the noncanonical basepairs also. On the other hand, usual negative slide around -1.5 Å for all dinucleotide steps in A-form structures are not maintained in the basepair steps containing noncanonical basepairs. Slide is seen to vary widely from small positive values to even as large as -3 Å for the steps containing one or more noncanonical basepairs. We have also observed that the average twist angles are rather unusual for some of the steps involving noncanonical basepairs; e.g., we found 11.9° twist for G:C W:WC::G:A S:HT of 2AW4 (Table 4). The smaller twist angles at the interface of noncanonical and canonical pairings arise from increased stacking interactions between the basepairs, which tend to stabilize the double helix. For a dinucleotide step, where two adjacent bases of one strand use Hoogsteen edge and sugar edge successively, as is seen in the structure of 1N32 [780A(s):802A(H)::781A(H): 801U(W)] and 2AW4 [537G(S):556A(H)::538A(S):555G(H)] (Figure 5), the basepair long axis, which lies along the vector joining C1′-C1′ atoms,54 is aligned nearly perpendicular to each other, giving rise to high twist values. This also indicates the large extent of cross-strand purine-purine stacking contributing to the helix stability. Tandem A:G H:ST pairs within duplex regions are known to be stabilized through such cross-strand stacking forces. It is clear from Figure 5 that such high twist values are not artifacts. Crystal structure analysis suggests that intrabasepair parameters for G:U wobble pairs are similar to canonical Watson-Crick basepairs, except for higher values of shear than that of the standard A:U or G:C-pairs.55,72 The larger shear can also affect the twist values of the corresponding dinucleotide steps, as one of the bases moves leading to a rotation of the basepair long axis (Figure 6a). For the 534U:559G::535G:558U step in 2AW4, we have large shearing motions of two consecutive basepairs, but in opposite directions, e.g., 3.1 and -3.0 Å. This high value

Figure 6. The deviation of a G:U basepair from normal G:C is shown in (a) by superposing the two structures [G:C in CPK, G:U in brown]. (b) shows the 534U:559G::535G:558U step in 2AW4, where the twist is ∼40°. This arises from high ∆shear, one with positive and other with negative values. The long axes are drawn by joining C6-C8 (Y:R) or C8-C8 (R:R) atoms. C1′ atoms are shown as spheres.

of ∆shear (-6.1 Å) between the successive basepairs makes the corresponding twist angle higher (40.0°) than the ideal values for A-DNA or A-RNA (∼32°) (Figure 6b). Similarly, positive ∆shear values of 2.3 and 3.3 Å reduce the twist angles of 535G:558U::536G:557C (19.5°) and 539G:554U::540C:553G (19.3°) in 2AW4, respectively (Table 4). A very similar nature is found in all structures belonging to the 2AW4 family. It can be speculated that a ∆shear value of 6 Å can reduce twist of two ideally stacked G:U basepairs in Gi:Uj::Ui+1:Gj-1 step around 10° without any deformation. However, this correlated change of twist angle with shearing motion of the consecutive basepairs is not reflected by 3DNA results. Though the values for other

14038

J. Phys. Chem. B, Vol. 114, No. 44, 2010

Halder and Bhattacharyya

step parametersstilt, roll, shift, slide, and risesare in good agreement for NUPARM and 3DNA, the values of twist show just the reverse trend. Whereas the twist for (UG).(UG) step is higher than a Watson-Crick doublet and less for (GG).(CU) and (GC).(GU), as shown pictorially (Figure 6b), 3DNA calculates a twist angle of ∼20° in the former case and ∼40° for the latter ones (Table 6). Apart from the value of shear, we find that, for some of the G:U W:WC pairs, the magnitudes of open angle deviate from standard Watson-Crick basepairs, even when the basepairs are not located at termini (Table 3) compared to the normal open angle value of ∼0° for G:U W:WC basepairs, as reported earlier. A detailed analysis indicates that such large values of average open angle arise from structural alterations occurring within the basepairs. We further note that even if the average values of local doublet parameters involving noncanonical basepairs are sometimes unusual, the basepair steps with noncanonical pairings have very small standard deviation values (Table 4). As the stacking free energies are related to the probability distribution of local doublet parameters (eq 2 of Samanta et al.73), Figure 7 indicates that these noncanonical steps have stacking free energy values comparable to those between canonical ones. Similar distributions are observed for the other noncanonical basepair steps also. Calculation of these free-energy components for all the basepair steps is a gigantic task and beyond scope of the present work. Furthermore, such an attempt would require absolute values of stacking interaction energies calculated by advanced quantum chemical methods, which is also tedious for such large systems containing hundreds of electrons. Nevertheless, our simulation results highlight free energy values of six unique dinucleotide steps containing noncanonical basepairs, namely, C:G::A:A, A:A::A:U, A:U::A:G, A:G::C:G, G:A::A:G, and A:G::G:U, which are equivalent to the Watson-Crick basepaired regions. Conclusion In our present study, we have employed all-atom MD simulation technique with explicit solvent model to comprehend the structural stability and role of noncanonical basepairs present within RNA double helices, our long-term goal being detailed characterization of different non-Watson-Crick basepairings. As an initial step, we have chosen two double helical fragments from rRNA with noncanonical basepairs at their central regions. We have observed that noncanonical basepairs that use the Hoogsteen or sugar edges for hydrogen bonding, at tandem, show significant stability in double helices. Thus it can be speculated that tandem noncanonical motifs occurring within double helices are stable enough to act as seeds of helix formation, though several other noncanonical motifs are needed to be analyzed for verification. Properties of the other tandem noncanonical motifs would be presented elsewhere. Further studies on the stability of single noncanonical basepairs with Hoogsteen or sugar edges, situated within duplex regions, can reveal interesting facts about their structural stability. As

Figure 7. Normalized frequency distribution for (a) twist and (b) roll for canonical and noncanonical steps. Black represents 782A:800G:: 783C:799G in 1N32 and gray is for a canonical step of RNA11. The equivalence of width for the distributions indicate their similarity in free energy despite the location of the mean values.

structures of such double helices are scarce yet, one needs to do a special model building study prior to their structural analyses. We have pointed out that the helix 1N32 makes a tertiary contactsthrough noncanonical pairingswith another helical fragment from a different part of the long RNA chain. Such interactions are vital for the three-dimensional organization of the macromolecules. In addition to this, we observe that both the helices 1N32 and 2AW4 are part of a stem-loop motif and reside at the surface of crystal structures, indicating their importance in binding of the ribosomal section with protein moieties. Again, occurrence of the noncanonical motifs on crystal surface gives evidence of their spontaneous origin without any contextual pressure. We further note that the dinucleotide steps containing noncanonical pairs are quite rigid in terms of relative orientation between successive basepairs, indicating they have high free-energy of stacking. This signifies their possible large contribution toward RNA folding by acting as seeds of helix formation and high stability that should not

TABLE 6: Comparison of NUPARM and 3DNA Results for Dinucleotide Steps Containing G:U-Basepairs As Obtained from Crystal Structures of 23S rRNA of E. coli (2AW4 family) (UG).(UG) (GG).(CU) (GC).(GU)

NUPARM 3DNA NUPARM 3DNA NUPARM 3DNA

tilt

roll

twist

shift

slide

rise

0.9(1.3) 1.9(2.2) 0.9(1.4) -0.1(2.8) -2.5(1.7) -3.5(2.8)

19.3(2.4) 18.3(2.4) 7.1(2.7) 7.6(2.5) 5.6(2.6) 5.1(2.6)

46.5(1.0) 20.4(1.4) 25.8(0.8) 37.8(0.8) 24.4(1.2) 39.0(1.1)

0.8(0.1) 1.1(0.3) 0.3(0.2) 0.3(0.2) -0.8(0.2) -0.7(0.2)

-2.1(0.2) -2.2(0.2) -1.8(0.2) -1.7(0.2) -1.2(0.1) -1.4(0.1)

2.9(0.1) 2.9(0.2) 3.2(0.1) 3.3(0.1) 3.4(0.1) 3.2(0.1)

Molecular Dynamics Studies of Functional RNA be ignored while predicting RNA secondary structure. RNA double helices are more prone to melt from termini, yet they are always rather short. It has been found very often that these short helical stretches have capping noncanonical basepairs, especially before the hairpin loops and at the junctions of continuous stacks. Detailed characterization of the noncanonical basepairs requires clear understanding about the dynamic stabilities of these helices. Acknowledgment. We acknowledge financial support from Dept. of Atomic Energy and Dept. of Biotechnology, Govt. of India. We are grateful to Prof. Manju Bansal, IISc, Bangalore, India, for important suggestions. We sincerely acknowledge computational facilities at Bioinformatics Resources and Application Facility at Center for Development of Advanced Computing, Pune and Centre for Applied Mathematics and Computational Sciences at Saha Institute of Nuclear Physics, Kolkata, India. We thank Parijat Majumder of Biophysics Division, Saha Institute of Nuclear Physics, Kolkata, India, and Somdutta Saha of West Bengal University of Technology, Kolkata, India, for their technical support in the work. Supporting Information Available: A detailed list of noncanonical motifs along with structural parameters for the simulation RNA11 and a comparison between structural parameters of noncanonical basepairs calculated by NUPARM and 3DNA. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Leontis, N. B.; Westhof, E. RNA 2001, 7, 499. (2) Hobza, P.; Sponer, J. Chem. ReV. 1999, 99, 3247. (3) Mladek, A.; Sharma, P.; Mitra, A.; Bhattacharyya, D.; Sponer, J.; Sponer, J. E. J. Phys. Chem. B 2009, 113, 1743. (4) Roy, A.; Panigrahi, S.; Bhattacharyya, M.; Bhattacharyya, D. J. Phys. Chem. B 2008, 112, 3786. (5) Sponer, J.; Jurecka, P.; Hobza, P. J. Am. Chem. Soc. 2004, 126, 10142. (6) Allawi, H. T.; SantaLucia, J., Jr. Biochemistry 1998, 37, 9435. (7) Allawi, H. T.; SantaLucia, J., Jr. Nucleic Acids Res. 1998, 26, 2694. (8) Allawi, H. T.; SantaLucia, J., Jr. Biochemistry 1998, 37, 2170. (9) Davis, A. R.; Znosko, B. M. Biochemistry 2007, 46, 13425. (10) Li, Y.; Zon, G.; Wilson, W. D. Biochemistry 1991, 30, 7566. (11) Morse, S. E.; Draper, D. E. Nucleic Acids Res. 1995, 23, 302. (12) Peyret, N.; Seneviratne, P. A.; Allawi, H. T.; SantaLucia, J., Jr. Biochemistry 1999, 38, 3468. (13) SantaLucia, J., Jr.; Kierzek, R.; Turner, D. H. Biochemistry 1991, 30, 8242. (14) Sugimoto, N.; Nakano, M.; Nakano, S. Biochemistry 2000, 39, 11270. (15) Turner, D. H. Curr. Opin. Struct. Biol. 1996, 6, 299. (16) Hammond, N. B.; Tolbert, B. S.; Kierzek, R.; Turner, D. H.; Kennedy, S. D. Biochemistry 2010, 49, 5817. (17) Das, J.; Mukherjee, S.; Mitra, A.; Bhattacharyya, D. J. Biomol. Struct. Dynam. 2006, 24, 149. (18) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (19) Foloppe, N.; MacKerell, A. D. J. Comput. Chem. 2000, 21, 86. (20) Sponer, J.; Mokdad, A.; Sponer, J. E.; Spackova, N.; Leszczynski, J.; Leontis, N. B. J. Mol. Biol. 2003, 330, 967. (21) Jeltsch, A. Chembiochem 2002, 3, 275. (22) Baeyens, K. J.; Debondt, H. L.; Holbrook, S. R. Nat. Struct. Biol. 1995, 2, 56. (23) Baeyens, K. J.; DeBondt, H. L.; Pardi, A.; Holbrook, S. R. Proc. Natl. Acad. Sci. U.S.A. 1996, 93, 12851. (24) Battiste, J. L.; Mao, H. Y.; Rao, N. S.; Tan, R. Y.; Muhandiram, D. R.; Kay, L. E.; Frankel, A. D.; Williamson, J. R. Science 1996, 273, 1547. (25) Cate, J. H.; Gooding, A. R.; Podell, E.; Zhou, K. H.; Golden, B. L.; Kundrot, C. E.; Cech, T. R.; Doudna, J. A. Science 1996, 273, 1678. (26) Lietzke, S. E.; Barnes, C. L.; Berglund, J. A.; Kundrot, C. E. Structure 1996, 4, 917.

J. Phys. Chem. B, Vol. 114, No. 44, 2010 14039 (27) Pley, H. W.; Flaherty, K. M.; McKay, D. B. Nature 1994, 372, 68. (28) Shen, L. X.; Cai, Z. P.; Tinoco, I. FASEB J. 1995, 9, 1023. (29) Gautheret, D.; Konings, D.; Gutell, R. R. J. Mol. Biol. 1994, 242, 1. (30) Cate, J. H.; Doudna, J. A. Structure 1996, 4, 1221. (31) Scott, W. G.; Finch, J. T.; Klug, A. Cell 1995, 81, 991. (32) Walczak, R.; Carbon, P.; Krol, A. RNA 1998, 4, 74. (33) Beveridge, D. L.; Barreiro, G.; Byun, K. S.; Case, D. A.; Cheatham, T. E.; Dixit, S. B.; Giudice, E.; Lankas, F.; Lavery, R.; Maddocks, J. H.; Osman, R.; Seibert, E.; Sklenar, H.; Stoll, G.; Thayer, K. M.; Varnai, P.; Young, M. A. Biophys. J. 2004, 87, 3799. (34) Gardiner, E. J.; Hunter, C. A.; Packer, M. J.; Palmer, D. S.; Willett, P. J. Mol. Biol. 2003, 332, 1025. (35) Madhumalar, A.; Bansal, M. J. Biomol. Struct. Dynam. 2005, 23, 13. (36) Young, M. A.; Ravishanker, G.; Beveridge, D. L. Biophys. J. 1997, 73, 2313. (37) Dixit, S. B.; Beveridge, D. L.; Case, D. A.; Cheatham, T. E.; Giudice, E.; Lankas, F.; Lavery, R.; Maddocks, J. H.; Osman, R.; Sklenar, H.; Thayer, K. M.; Varnai, P. Biophys. J. 2005, 89, 3721. (38) Du, Z.; Yu, J.; Andino, R.; James, T. L. Biochemistry 2003, 42, 4373. (39) Menger, M.; Eckstein, F.; Porschke, D. Biochemistry 2000, 39, 4500. (40) Proctor, D. J.; Ma, H.; Kierzek, E.; Kierzek, R.; Gruebele, M.; Bevilacqua, C. Biochemistry 2004, 43, 14004. (41) Sorin, E. J.; Engelhardt, M. A.; Herschlag, D.; Pande, V. S. J. Mol. Biol. 2002, 317, 493. (42) Villa, A.; Widjajakusuma, E.; Stock, G. J. Phys. Chem. B 2008, 112, 134. (43) Spackova, N.; Sponer, J. Nucleic Acids Res. 2006, 34, 697. (44) Ditzler, M. A.; Otyepka, M.; Sponer, J.; Walter, N. G. Acc. Chem. Res. 2010, 43, 40. (45) Ray, S. S.; Halder, S.; Bhattacharyya, D. Conference Proceedings of International Conference on Frontiers of Interface Between Statistics and Sciences; 2009; p 724. (46) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. Nucleic Acids Res. 2000, 28, 235. (47) Arnott, S.; Hukins, D. W. L.; Dover, S. D. Biochem. Biophys. Res. Commun. 1972, 48, 1392. (48) MacKerell, A. D.; Banavali, N. K. J. Comput. Chem. 2000, 21, 105. (49) Samanta, S.; Mukherjee, S.; Chakrabarti, J.; Bhattacharyya, D. J. Chem. Phys. 2009, 130, 115103. (50) Darden, T.; Perera, L.; Li, L. P.; Pedersen, L. Struct. Fold. Des. 1999, 7, R55. (51) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (52) Kale, L.; Skeel, R.; Bhandarkar, M.; Brunner, R.; Gursoy, A.; Krawetz, N.; Phillips, J.; Shinozaki, A.; Varadarajan, K.; Schulten, K. J. Comput. Phys. 1999, 151, 283. (53) Nelson, M.; Humphrey, W.; Kufrin, R.; Gursoy, A.; Dalke, A.; Kale, L.; Skeel, R.; Schulten, K. Comput. Phys. Commun. 1995, 91, 111. (54) Bansal, M.; Bhattacharyya, D.; Ravi, B. Comput. Appl. Biosci. 1995, 11, 281. (55) Mukherjee, S.; Bansal, M.; Bhattacharyya, D. J. Comput. Mol. Des. 2006, 20, 629. (56) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (57) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (58) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (59) Bevan, D. R.; Li, L. P.; Pedersen, L. G.; Darden, T. A. Biophys. J. 2000, 78, 668. (60) Cheng, Y. H.; Korolev, N.; Nordenskiold, L. Nucleic Acids Res. 2006, 34, 686. (61) Korolev, N.; Lyubartsev, A. P.; Laaksonen, A.; Nordenskiold, L. Nucleic Acids Res. 2003, 31, 5971. (62) Pan, Y. P.; MacKerell, A. D. Nucleic Acids Res. 2003, 31, 7131. (63) Reblova, K.; Lankas, F.; Razga, F.; Krasovska, M. V.; Koca, J.; Sponer, J. Biopolymers 2006, 82, 504. (64) Romanowska, J.; Setny, P.; Trylska, J. J. Phys. Chem. B 2008, 112, 15227. (65) Schneider, B.; Neidle, S.; Berman, H. M. Biopolymers 1997, 42, 113. (66) Chou, S. H.; Zhu, L. M.; Reid, B. R. J. Mol. Biol. 1997, 267, 1055. (67) Skelly, J. V.; Edwards, K. J.; Jenkins, T. C.; Neidle, S. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 804. (68) Dickerson, R. E.; Bansal, M.; Calladine, C. R.; Diekmann, S.; Hunter, W. N.; Kennard, O.; et al. EMBO J. 1989, 8, 1.

14040

J. Phys. Chem. B, Vol. 114, No. 44, 2010

(69) Olson, W. K.; Bansal, M.; Burley, S. K.; Dickerson, R. E.; Gerstein, M.; Harvey, S. C.; Heinemann, U.; Lu, X. J.; Neidle, S.; Shakked, Z.; Sklenar, H.; Suzuki, M.; Tung, C. S.; Westhof, E.; Wolberger, C.; Berman, H. M. J. Mol. Biol. 2001, 313, 229. (70) Calladine, C. R. J. Mol. Biol. 1982, 161, 343. (71) Lu, X. J.; Olson, W. K. Nucleic Acids Res. 2003, 31, 5108.

Halder and Bhattacharyya (72) Olson, W. K.; Esguerra, M.; Xin, Y. R.; Lu, X. J. Methods 2009, 47, 177. (73) Samanta, S.; Bhattacharyya, D.; Chakrabarti, J. J. Biomol. Struct. Dyn. 2010, 27, 429.

JP102835T