Î³ Backbone Conformations in RNA and the Accuracy

Article pubs.acs.org/JPCB

Noncanonical α/γ Backbone Conformations in RNA and the Accuracy of Their Description by the AMBER Force Field Marie Zgarbová,† Petr Jurečka,*,† Pavel Banás,̌ † Marek Havrila,‡ Jiří Šponer,†,‡ and Michal Otyepka† †

Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacky University, 17. listopadu 12, 77146 Olomouc, Czech Republic ‡ Institute of Biophysics, Academy of Sciences of the Czech Republic, Královopolská 135, 612 65 Brno, Czech Republic S Supporting Information *

ABSTRACT: The sugar−phosphate backbone of RNA can exist in diverse rotameric substates, giving RNA molecules enormous conformational variability. The most frequent noncanonical backbone conformation in RNA is α/γ = t/t, which is derived from the canonical backbone by a crankshaft motion and largely preserves the standard geometry of the RNA duplex. A similar conformation also exists in DNA, where it has been extensively studied and shown to be involved in DNA−protein interactions. However, the function of the α/γ = t/t conformation in RNA is poorly understood. Here, we present molecular dynamics simulations of several prototypical RNA structures obtained from X-ray and NMR experiments, including canonical and mismatched RNA duplexes, UUCG and GAGA tetraloops, Loop E, the sarcin−ricin loop, a parallel guanine quadruplex, and a viral pseudoknot. The stability of various noncanonical α/γ backbone conformations was analyzed with two AMBER force fields, ff99bsc0χOL3 and ff99bsc0χOL3 with the recent εζOL1 and βOL1 corrections for DNA. Although some α/γ substates were stable with seemingly well-described equilibria, many were unstable in our simulations. Notably, the most frequent noncanonical conformer α/γ = t/t was unstable in both tested force fields. Possible reasons for this instability are discussed. Our work reveals a potentially important artifact in RNA force fields and highlights a need for further force field refinement.

■

INTRODUCTION RNA is a very versatile molecule that can adopt many different geometrical motifs and form complex structures such as RNA riboswitches and aptamers,1 ribozymes,2 or the ribosome.3 The variability of RNA motifs and interaction patterns stems from the interplay of hydrogen bonding and stacking between nucleobases, as well as hydrogen bonding between nucleobases, backbone phosphate groups, and sugar moieties. An important role is also played by the sugar−phosphate backbone, which may prefer or penalize certain nucleobase arrangements. Thus, the backbone can contribute decisively to the relative stability of certain RNA motifs and co-determine the RNA’s folding and interaction patterns. Although DNA structures exhibit less overall variability, the conformation of the sugar−phosphate backbone in DNA also has important structural consequences. The backbone’s conformation can be described in terms of five dihedral angles, α, β, γ, ε, ζ, along with the sugar pucker P and glycosidic angle χ (Figure 1). In regular RNA or DNA duplexes with canonical Watson−Crick pairing, these angles have “canonical” values.4 However, all of them can take other, noncanonical values corresponding to alternative minima in their torsion profiles. Naturally occurring backbone conformers from X-ray databases have been studied extensively5−7 and comprehensively classified for both RNA8 and DNA.9 For DNA, the most important conformations are those corresponding to the B-DNA, A-DNA, and Z-DNA families. For instance, the second most frequent conformer in the DNA X-ray © 2017 American Chemical Society

Figure 1. Important dihedral angles in the nucleic acid backbone.

database after the canonical BI-DNA is BII-DNA.9 BII-DNA differs from BI-DNA in the values of the ε and ζ backbone angles (BI: ε/ζ = t/g−, BII: ε/ζ = g−/t). The BI/BII equilibrium has been studied in great depth experimentally and theoretically (see, e.g., ref 10) not least because it seems to be important in certain protein−DNA interactions.11,12 Received: January 10, 2017 Revised: February 22, 2017 Published: March 14, 2017 2420

DOI: 10.1021/acs.jpcb.7b00262 J. Phys. Chem. B 2017, 121, 2420−2433

Article

The Journal of Physical Chemistry B

following sets of AMBER biomolecular force fields: ff10, ff12SB, and ff14SB (internal AMBER code abbreviations used for the combined set of protein and nucleic acid force fields). This force field is referred to as χOL3 henceforth. We also performed simulations using χOL3 with the εζOL134 and βOL135 modifications, which were derived using a methodology including conformation-dependent solvation effects.36 This combination is henceforth referred to as χOL3 + εζβOL1. The εζOL1 and βOL1 modifications were originally developed for DNA simulations and significantly improved the description of B-DNA, Z-DNA, and guanine quadruplexes.34,35,37 One of the aims of this work is to test the performance of the χOL3 + εζβOL1 force field in RNA simulations. The α/γ backbone substates of B-DNA have been thoroughly studied using an earlier version of the AMBER force field (ff99),15 revealing that the noncanonical α/γ = g+/t state was overly stabilized by ff99.13 Consequently, undesirable deformations of the B-DNA helices occurred in the MD simulations, accompanied by a loss of helical twist and unwinding of the duplex on a time scale of tens of nanoseconds. This misbehavior was eliminated by adopting the bsc0 modification of the α/γ backbone torsions,16 enabling stable simulations of B-DNA. For the purpose of later discussion, we note that, to a large extent, this modification appears to have succeeded by destabilizing the γ = t conformation with respect to the γ = g+ energy minimum.16 Frequent α/γ transitions were also observed in RNA simulations with ff99.38 However, these transitions were of a slightly different character than those in DNA (α/γ = t/t in RNA vs g+/t in B-DNA). Because the populations of the nonnative substates were clearly too high (almost 30% in canonical A-RNA duplexes), we applied the bsc0 correction in the RNA simulations. Bsc0 almost completely eliminated the non-native α/γ = t/t structures, so we concluded that the bsc0 correction is also beneficial in RNA simulations, although penalization of the α/γ = t/t substate might be a little excessive.17,21 Today, bsc0 is the recommended α/γ potential for both DNA and RNA simulations and, as noted above, along with our glycosidic χOL3 correction,21 is used for RNA in all recent versions of the AMBER biomolecular force field sets. Note that for RNA, the χOL3 correction is necessary to achieve basic stability of the simulations.21 Bsc0 could be dispensable, but it seems to improve RNA simulations significantly when compared with the original α/γ parameterization.21,39 However, the stability of various α/γ backbone substates in RNA has never been systematically assessed using these force fields. The following section discusses the stability and behavior of selected RNA structures in MD simulations with the χOL3 and χOL3 + εζβOL1 force fields. We initially present the results for each structure separately and then summarize the general trends observed in the simulations, the stability of individual α/ γ combinations, and discuss potential reasons for the instability of the α/γ = t/t conformers.

Forty-six discrete RNA backbone conformations, or suites, were identified by Richardson et al. in their X-ray database analysis.8 The dominant conformation is A-RNA (henceforth referred to as conformation 1a, according to Richardson et al.). The second most frequent substate is derived from 1a by a concerted crankshaft-like motion of the α and γ backbone angles from the canonical α/γ = g−/g+ state to the α/γ = t/t state (1c). The analogous noncanonical conformation of DNA (α/γ = g+/t) has been studied in depth13−16 and is believed to be an important substate in protein−DNA complexes.15 The 1c substate in RNA has received much less attention and its structural roles and prevalence are less well understood than those of its DNA counterpart despite its status as the most frequent noncanonical RNA backbone conformation. To fill this gap in our knowledge, we have studied a collection of noncanonical α/γ RNA substates in X-ray and NMR structures of smaller naked RNA molecules. Although we are especially interested in the 1c α/γ = t/t conformation, we also briefly review other possible α/γ combinations. The α and γ angles can each adopt three distinct conformations: gauche+ (g+), trans (t), and gauche− (g−), giving rise to 9 possible α/γ combinations.14 Most of these combinations are represented in the X-ray data for both free RNA and protein−RNA complexes, and many are present in the consensus RNA conformers defined by Richardson et al.8 It has been suggested that noncanonical α/γ conformations play important structural roles in RNA.8 For instance, a change in a single dihedral angle can change the strand direction, potentially causing upstream or downstream bases to be separated or brought closer together, facilitating loop formation. The 1c conformer features simultaneous changes in two dihedral angles (relative to the canonical conformation) that counteract one-another to some degree such that the strand direction is not changed. Consequently, the 1c conformer is readily incorporated into the A-RNA double helix. The most obvious geometric consequence of its incorporation is that the neighboring phosphate groups separated by the α/γ torsions move about 1 Å further apart (the P−P distance increases from about 6 to about 7 Å). This makes it possible to reliably identify the presence of 1c in X-ray datasets because the positions of heavy phosphates are determined relatively well. The situation in solution is more complicated because accurate determination of the α and γ torsions is not a simple task in NMR experiments (for more details, see discussion of the 2DD2 dynamics in the Results and Discussion section). The presence of 1c in the ARNA double helix also reduces its helical twist and widens the major groove.17 Aside from these subtle but non-negligible geometric changes, the presence of 1c does not disrupt the overall geometry of the A-RNA duplex. For molecular modeling of biomolecules, it is important to know how well various RNA conformations are described by the available molecular mechanics RNA force fields. Two major nonpolarizable force field families are used for RNA simulations: the AMBER family based on the Cornell et al. force field,16,18−21 and the CHARMM force field.22−25 In addition, the new polarizable Drude-2013 force field for nucleic acids is undergoing rapid development.26−28 This work focuses on the AMBER family, which is commonly used in RNA simulations. AMBER force fields are under constant development29−33 and several alternative testing versions are currently available.30,32,33 Most of our results are obtained with the currently recommended ff99bsc0χOL3 version. This force field has been consistently used as the RNA force field in the

■

METHODS Model Molecules. Several smaller RNA molecules were chosen. The selected molecules all had well-resolved X-ray or NMR structures with no ligands or bound proteins (Table 1, Figure 2). We focused primarily on models containing the noncanonical α/γ = t/t substate, preferably in the 1c suite. However, these molecules contained multiple other noncanonical α/γ substates, which we also considered in our analysis. For comparative purposes, we also included two RNA 2421


Article


structures were modified to facilitate simulation. A short GCrich duplex decamer, r(GCACCGUUGG)2, was excised from the full GC-rich duplex (PDB ID 1QC0) to yield the structure denoted 1QC0′ (the duplex was shortened to allow comparison with our previous simulations of the same structure). Similarly, a GAGA tetraloop was excised from the X-ray structure of SRL (PDB ID 1Q9A43) and capped with two GC pairs (r(CGCGAGAGCG)) to suppress fraying.50 The stem of the sarcin−ricin loop was shortened by two noncanonical residues at each end to avoid fraying, and the simulations of the PDB ID 483D system were performed using the A molecule as the starting structure. Magnesium ions were kept at their crystallographic positions in the starting structures for the simulations of Loop E (only the B position of the doubly occupied 4MgA/4MgB structure was retained51) and the BYWV pseudoknot. Two 5′-terminal residues were removed from the BWYV pseudoknot structure, as indicated in Figure 2, to reduce the conformational space to sample, and the simulations of PDB ID 1L2X were performed using the B molecule as the starting structure. The structure of the r(5′GGU GAA GGCU/PCCG AAG CCG5′) mismatched duplex with the AA sheared pair was taken from the NMR structure associated with PDB ID 2DD2, from which dangling residues were removed.

Table 1. List of Studied RNA Structures system duplex 1QC0′ duplex 1RNA UUCG tetraloop GAGA tetraloop SRL pseudoknot loop E duplex 2DD2 quadruplex aptamer

PDB ID

method (resolution in Å)

GC-rich A-RNA duplex AU-rich A-RNA duplex UNCG tetraloop

1QC0′40 1RNA41 2KOC42

X-ray (1.55) X-ray (2.25) NMR

GNRA tetraloop

1Q9A43

X-ray (1.09)

sarcin−ricin loop BYWV viral pseudoknot internal loop mismatched GAA/AAG duplex parallel G quadruplex anti-NFκB RNA aptamer

483D44 1L2X45 354D46 2DD247

X-ray (1.11) X-ray (1.25) X-ray (1.5) NMR

4XK048 2JWV49

X-ray (1.08) NMR

description

duplexes, one AU-rich and the other containing both AU and GC pairs. Three independent 1 μs simulations were run for each system and each force field, except for the first two duplexes in Table 1, for which only one 1 μs simulation was performed. The initial structures for our simulations were taken from the PDB entries with the IDs listed in Table 1. Some of the

Figure 2. Schematic depictions and numbering of test structures. The numbering is the same as in the starting PDB structures. The studied residues and α/γ conformations are shown in color; residues removed from the starting structures are shown in gray. 2422


Article

The Journal of Physical Chemistry B Table 2. Overview of α/γ Conformations in the X-ray Dataset Compiled Using the Richardson et al. Data8a

a Numbers of occurrences are given in italics and suites containing the specified conformations are shown in bold. The PDB IDs of test structures from this work that contain each conformation (and its location within the structure) are shown in regular typeface.

The starting structures were neutralized with K+ ions and the ion concentration was then adjusted to 150 mM using the Joung and Cheatham K+ and Cl− ion parameters.52,53 An octahedral SPC/E54 water box was used to solvate the molecular systems, with a buffer size of 12.0 Å. All simulations were performed using the AMBER 1455 program suite with simulated NPT conditions (1 bar, 298 K, tautp = 1.0 ps, taup = 1.0 ps), hydrogen mass repartitioning, a 4 fs time step,56,57 and a 9 Å direct space nonbonded cutoff. SHAKE constraints were applied to bonds to hydrogen atoms with the default tolerance (1 × 10−5 Å). PME was used with default grid settings (1 Å) and default tolerance (1 × 10−5). PMEMD for CUDA was used in all simulations.58 Initial equilibration was performed as described elsewhere50 and coordinates were stored every 10 ps. Trajectories were analyzed using the cpptraj59 module of AMBER. Two force field variants were tested, χOL3 (ff99bsc0χOL3) and χOL3 + εζβOL1, a combination of χOL3 with the ε/ζOL134 and βOL135 modifications, which were originally developed for DNA simulations (see the Introduction). Parameters for the N3protonated C8+ residue in the BWYV pseudoknot were taken from an earlier publication.60 The structures of Loop E and the BYWV pseudoknot contain magnesium ions, which were modeled with the parameters of Allnér et al.61 The total simulation time for the analyses reported in this work was 52 μs. The gauche+, trans, and gauche− backbone angle conformations are abbreviated as g+, t, and g−, respectively. The C2′- and C3′-endo sugar puckering parameters are abbreviated as 2′ and 3′, respectively, throughout the text. We retained the residue numbering of the PDB files used as starting structures. In addition, we adopted the backbone suite nomenclature of Richardson et al.;8 suite names are quoted in bold in the text (for instance, the canonical RNA backbone is denoted 1a). Backbone angles and suites were analyzed using the Dangle and Suitename programs (both available from http://kinemage. biochem.duke.edu). However, MD simulations produce rather wide distributions for some dihedral angles, which could cause some MD snapshots to not be assigned to the nearest suite despite broadly satisfying its membership criteria and belonging to a continuous distribution associated with that suite. Therefore, for the purposes of this article, we define the g+, t, and g− substates such that unless stated otherwise, angles in the intervals 0−120°, 120−240°, and 240−360° are classified as

g+, trans, and g−, respectively. In some cases a backbone conformation could not be assigned to any of the consensus suites; this usually (but not exclusively) occurred during the MD simulations rather than when analyzing the starting structures. If only one dihedral angle deviated from the suite definition, we indicated this by appending the divergent dihedral angle and an exclamation mark onto the suite symbol. For instance, a structure in the 1z suite with the sugar pucker flipped from C2′-endo to C3′-endo would be denoted 1z!P (see the Discussion of the UUCG tetraloop below). In cases with multiple dihedral angle mismatches, the structure was considered nonclassified and so was assigned to a unique suite denoted nc1, nc2, and so forth. Rare conformations that existed for less than 10% of the simulation time during the MD simulations were ignored in most cases, especially for systems that had at least two more highly populated conformers, as was the case for SRL (see below).

■

RESULTS AND DISCUSSION Noncanonical α/γ backbone conformations are ubiquitous in RNA structures. However, some are much more common than others. Table 2 shows the frequency of occurrence of various α/γ combinations in the X-ray data studied by Richardson et al., 8 along with the consensus suites featuring these combinations. Most RNA structures are based on the canonical α/γ = g−/g+ conformation. The most populated noncanonical conformation is the “crankshaft” 1c suite, which accounts for almost 5% of the backbone substates investigated by Richardson et al. It is present in some stable hairpin tetraloops, such as the GAGA tetraloop; sometimes it accommodates a noncanonical base pair,62,63 but most often it simply appears in a canonical A-RNA double helix, which indicates that it is easily energetically accessible. The α/γ = g+/g+ and t/g+ conformations are also quite common, and even the less populated g+/t and g−/t conformations appear in some important and widely studied structures such as stable tetraloops, Loop E, and the sarcin−ricin loop. However, conformations in which γ = g− are rare. In this work, we focused mainly on the α/γ = t/t conformations. However, because our structures also contain some other noncanonical α/γ backbone substates, we analyzed these as well. Table 2 also specifies the PDB IDs of structures from our test set that feature the listed α/γ substates along with their locations within the structure and suitename to which they belong. The 2423


Article


Figure 3. RMSDs of simulated structures (for residues included in the RMSD calculations, see Table S1 in the Supporting Information). Colors correspond to independent simulations MD1 (black), MD2 (red), and MD3 (green).

refinements, remain far from flawless.64 Thus, the discussion below focuses on the stability of individual backbone conformations during the stable parts of the MD simulations, before any significant structural rearrangement occurred. Canonical Duplexes. The RNA duplexes simulated in this study were the 14-mer r(U(AU)6A) 1RNA containing only AU pairs and the decamer r(GCACCGUUGG)2 1QC0′, which contains mainly GC pairs. Whereas the 1QC0′ structure is purely canonical A-RNA, the 1RNA structure contains two noncanonical α/γ = t/t substates (belonging to the 1c suite) at residues A11 and U24. It should be noted that these residues form intermolecular hydrogen bonding contacts, which may stabilize the unusual geometry.41 Thus, we do not expect the 1c conformer to be fully stable in our solution MD simulations. However, these noncanonical α/γ substates were previously found to be unstable in simulations of the 1RNA crystal with the ff99bsc0χOL3 force field,65 indicating that the stability of the experimentally observed α/γ = t/t substate in this structure may

following section describes the behavior of each test structure in our MD simulations. The individual simulations revealed quite variable behaviors. Some of the tested RNA structures were essentially stable in our MD simulations, maintaining relatively small RMSD values relative to the initial structure (Figure 3) even when some of the native backbone conformations were lost (see below). However, some structures underwent significant rearrangements, such as the S-turn of SRL in simulations with χOL3 + εζβOL1, or the internal loop of the anti-NFκB aptamer, which was unstable in simulations with both force fields. Note that some parts of the structures were excluded when calculating RMSD values, namely, those that were expected to be highly flexible, such as bulged-out residues or flexible loops (for residues included in the RMSD calculations, see Table S1). As such, it is quite likely that at least some of the larger deviations from the experimental structures are due to some deficiency of the existing force fields, which, despite all of the recent 2424


Article

The Journal of Physical Chemistry B be significantly underestimated by the tested force fields. In our simulations, this 1c conformation was predicted to be unstable in solution by both the χOL3 and χOL3 + εζβOL1 force fields; it was consistently converted into the canonical 1a α/γ = g−/g+ conformation within one simulated nanosecond (1c and 1a suites are depicted in Figure S1). As noted above, this is not surprising as the examined substate was not stable, even in its natural crystal environment. The overall shape of the duplex was maintained by both the χOL3 and χOL3 + εζβOL1 force fields, as indicated by the relatively low RMSD values for these structures. As reported previously, the helical parameters of both duplexes were also well within the range expected for A-RNA.35 However, whereas the χOL3 force field maintained the canonical backbone angles almost entirely, χOL3 + εζβOL1 permitted certain backbone reconformations next to the ends of both duplexes (see Figures S2 and S3). These reconformations could be linked to fraying of the ends,50 which was more pronounced with the χOL3 + εζβOL1 force field than with χOL3. We also observed another noncanonical backbone substate near the ends of the 1RNA and 1QC0′ duplexes in the χOL3 + εζβOL1 simulations: the sugar pucker parameter and the χ and ζ angles flipped from the canonical P/χ/ζ = C3′/anti/g− conformation to the noncanonical C2′/high-anti/t conformation. This was accompanied by smaller shifts of the β and ε angles (see Figures S2 and S3). Although the noncanonical states were only observed near the ends of the duplexes, these results might suggest that the ARNA duplexes are less stable in the χOL3 + εζβOL1 simulations than in χOL3 simulations. Although this instability may be unimportant in most applications, it is apparent that the χOL3 force field describes canonical A-RNA duplexes better. UUCG Tetraloop. The UUCG tetraloop is one of the most thermodynamically stable tetraloops, occurs very frequently in RNA structures, and has been widely studied by MD simulations.30,66−68 However, it has been suggested that the UUCG fold is unstable in long MD simulations with the χOL3 force field69,70 (and any other existing force field70), that is, that it is not a global free energy minimum on the folding landscape. This suggestion was recently confirmed by an accurate evaluation of the UUCG folding energy (+30 kJ/mol).71 In unbiased simulations, UUCG unfolding usually begins with tetraloop reconformation, which occurs on the microsecond time scale.69 With the χOL3 + εζβOL1 force field, a similar loop reconformation occurs on an even shorter time scale (see Figure S4 and the RMSD data presented in Figure 3), possibly because of the lower conformational barriers introduced by the εζOL1 correction (see also our previously published ε/ζ profiles34). Nevertheless, the simulation times before unfolding begins are long enough to evaluate the stability of the relevant noncanonical suites in this tetraloop. The discussion below focuses exclusively on the stable period of the simulations, before the tetraloop reconformation begins. Two noncanonical α/γ conformations are found in the UUCG loop region. The first is α/γ = t/g+ on U7 (UUCG), accompanied by a C2′-endo sugar pucker on U7 (suite 1z, Figure S1). Although the U7 residue is exposed, that is, it is outside the loop structure and not stacked or hydrogen bonded with any loop residues, according to the high-resolution NMR experiment,42 the 1z suite is a native and a stable conformation in solution. The 1z substate was stable in MD simulations with χOL3, but with χOL3 + εζβOL1, the sugar pucker of U7 was flipping to C3′-endo for most of the stable part of the simulation, yielding a suite denoted 1z!P (see Tables 3 and 4

Table 3. Backbone Conformations of the UUCG Tetraloop in the Initial NMR Structure (2KOC) and the Stable Parts of Our MD Simulations χOL3

exp.

χOL3 + εζβOL1

location

suite

suite

%

suite

%

U6pU7

1z

1z

100

U7pC8

2[

2[

100

C8pG9

6n

6n

100

1z 1z!P 2[ nc1 6n

34 66 29 64 100

Table 4. Backbone Conformational Suites in the UUCG Tetraloopa

a

Noncanonical values are highlighted in orange.

and Figures S1 and S5). The second noncanonical α/γ conformation is g+/t at G9 (UUCG) in the 6n suite, which contains a noncanonical C2′-endo sugar pucker at C8 together with ε = g+ and ζ = g− substates. This conformation was stable in all our MD simulations before loop reconformation. GAGA Tetraloop. GAGA belongs to the GNRA family of tetraloops, which occur frequently in diverse RNAs and participate in protein−RNA and RNA−RNA interactions.43 We studied an isolated tetraloop with a short three-base-pair canonical stem, r(CGCGAGAGCG). The GAGA tetraloop in the SRL structure (see below) is identical to this one, and in our simulations it behaves in virtually the same way as the isolated GAGA loop. Two noncanonical α/γ substates were observed in GAGA (Tables 5 and 6 and Figure S1). The first is Table 5. Backbone Conformations of the GAGA Tetraloop in the Initial X-ray Structure (1Q9A) and Our MD Simulations χOL3

X-ray

χOL3 + εζβOL1

location

suite

suite

%

suite

%

G2659pA2660 A2662pG2663

1g 1c

1g 1a

∼100 ∼100

1g 1a

∼100 ∼100

in the base step G2659pA2660, which is classified as belonging to the 1g GAGA U-turn suite with an α/γ = t/g+ conformation. The other is the crankshaft α/γ = t/t (1c) conformation in the A2662pG2663 step between the loop and stem regions. These Table 6. Backbone Conformational Suites Observed in the GAGA Tetraloopa

a

2425

Noncanonical values are highlighted in orange. DOI: 10.1021/acs.jpcb.7b00262 J. Phys. Chem. B 2017, 121, 2420−2433

Article


Many features of the S-turn were poorly described in our MD simulations. Among other things, the S-turn was unstable in two out of three χOL3 + εζβOL1 simulations because the residue G2655 in the GpU platform lost its interactions with the phosphate group of A2665 on the opposite strand and bulged out into solution after 210 ns (MD1) and 535 ns (MD3); see the RMSD values for SRL in Figures 3 and S1. The analysis below considers only the parts of the simulations before such reconformations occurred. Figure 4 shows the time development of selected backbone angles and backbone suites (listed at the top of each graph) for one of three MD simulations for each force field (the remaining data are presented in Figure S7). The suites are colored in shades of green if consistent with the initial X-ray structure and red if they diverge appreciably from the X-ray suite. For clarity, states existing for less than 10% of the simulation time and unclassified conformations not listed in Table 7 are shown in gray. The topmost row of Figure 4 shows results for the U2653pA2654 step, which adopts a conformation belonging to the 5z suite but is also observed in a nonclassified (nc2) conformation in the crystal structure. The 5z geometry was only approximately reproduced in our MD simulations; whereas α ∼ 164° in 5z,8 we observed values of around 130° or even 120°. Because these values are between those assigned to the 5z suite and the “wannabe” suite 5p with α ∼ 75° (see Richardson et al.8), we denoted this conformation 5z/5p. However, because the departure from the canonical 5z conformation was small and inaccuracies in empirical potentials often shift the maxima of angle distributions, we colored the dominant 5z/5p conformation green to indicate relative closeness to the native state. The minor nc2 conformation is also close to that observed in the X-ray structure (Table 7), so both the χOL3 and χOL3 + εζβOL1 force fields provide a reasonably good overall description of this backbone region containing the α/γ = t/g+ conformation. Similarly, the 4s and nc3 suites at the A2654pG2655 step are also described quite well, although the values of the β and ζ angles are slightly shifted with both force fields (see Table 7 and Figures 4 and S7). We do, however, note that the 4s suite containing the γ = t conformation is less populated in our MD simulations (21 and 40% with χOL3 and χOL3 + εζβOL1, respectively) than in the crystal structures (51, 47, and 100% in 483D, 1Q9A, and 3DVZ, respectively). The last part of the S-turn is the G2655pU2656 step. In the X-ray structure, this step is assigned to the #a suite and has a canonical α/γ = g−/g+ conformation; we included it in our analysis to determine its effect on the Sturn’s stability. The #a suite is unstable in the χOL3 force field, and the residues in question spend most of their time in the 6g conformation (Table 7; for clarity only the two most populated conformations are listed). When the G2655pU2656 step is in the 6g conformation, a native base−phosphate interaction between the U2656 phosphate and G2664 residue of the opposite strand is broken. However, #a is fully stable in simulations with χOL3 + εζβOL1 and the native base−phosphate interaction is retained. Possibly, this is due to stabilization of the ζ = t region by the εζOL1 potential.34 The last motif with the noncanonical α/γ = g−/t conformation is the 1e suite conformation at the G2664pA2665 step opposite the S-turn. This motif is only partially stable in both force fields, accounting for 77 and 56% of the simulation time in χOL3 and χOL3 + εζβOL1, respectively. In the simulations, these residues occasionally adopt a nonnative &a conformation that lacks the native γ = t conformation associated with the 1e suite.

conformations are consistently observed in GAGA X-ray structures. In our MD simulations, GAGA exhibited very small RMSD values relative to its X-ray structure. However, only one of the two noncanonical substates, 1g with α/γ = t/g+ at G2659pA2660, was stable (see Figure S6). The 1c suite conformation at A2662pG2663 was unstable, being converted into the canonical 1a conformation during the first few nanoseconds in all our simulations. Sarcin−Ricin Loop. The sarcin−ricin loop is an intricate internal loop with a characteristic S-turn backbone topology, and it is capped with a GAGA tetraloop in our model molecule. The GAGA tetraloop behaves like the isolated GAGA tetraloop discussed above, so we focus mainly on the S-turn, which has been thoroughly studied before.72−74 The S-turn is considered to be one of the most extreme RNA backbone topologies. The backbone in this region features two sharp bends, giving it an “S” shape and necessitating the adoption of noncanonical backbone substates at residues U2653, A2654, G2655, and U2656 (see Figure 2). The conformations of these four residues in the X-ray structure with PDB ID 3DVZ can be denoted U26535zA26544sG2655#aU2656, but some crystals (483D, 1Q9A) exhibit double occupancies at the U2653pA2654 and A2654pG2655 steps (see Tables 7 and 8 and Figure S1). Because the two Table 7. Backbone Conformations of SRL in the Initial X-ray Structure (483D) and Our MD Simulations χOL3

X-ray

χOL3 + εζβOL1

location

suite

%

suite

%

suite

%

U2653pA2654

G2655pU2656

5z nc2 4s nc3 #a

66 34 51 49 100

85 14 53 40 92

1e

100

61 38 56 21 43 40 77 21

5z/5p nc2 nc3 4s #a

G2664pA2665

5z/5p nc2 nc3 4s 6g #a 1e &a

1e &a

56 35

A2654pG2655

Table 8. Backbone Conformational Suites in SRLa

a


possible structures have similar occupancy numbers in crystals (Table 7 only shows occupations for 483D), we expected that the backbone might also visit two different conformations in the MD simulations. Another noncanonical backbone conformation (belonging to the 1e suite) is found at residue A2665 in the strand opposite the S-turn, which is bound to U2656 of the S-turn’s GpU platform. This conformer exhibits the 1e suite’s characteristic α/γ = g−/t substate but differs from the standard 1e suite in that it has an unusually low β value of around 80°. 2426


Article


Figure 4. Time development of selected dihedral angles in the sarcin−ricin loop and the associated conformational suites (top of each graph). Suites consistent with the initial X-ray structure are shown in shades of green; those inconsistent with the X-ray structure are shown in red. States observed for less than 10% of the simulation time and unclassified conformations not listed in Table 7 are shown in gray.

which are strongly bound to the stem structures: C8+ is hydrogen bonded to G12 and A9 stacks with the next residue, C10. In our MD simulations, Loop 1 and the canonical stem regions, Stems 1 and 2, were stable with both force fields. However, the bulged-out residues U13 and G19 moved freely in solution, and Loop 2 also underwent substantial reconformation, as would be expected after the loss of crystal contacts. Consequently, most of the noncanonical backbone substates observed in these regions vanished during the first nanoseconds of our MD simulations (Figure S8). In addition, a noncanonical unclassified substate with an α/γ = g−/t conformation at

BYWV Pseudoknot. The crystal structure of the BYWV pseudoknot (PDB ID 1L2X) features several noncanonical α/γ substates. However, many of these are probably induced by extensive intermolecular interactions in the crystal lattice or are present in the regions corresponding to flexible loops. In particular, residues GTP1 and G2 (not included in our model) and U13 and G19 (see Figure 2) are bulged-out and strongly hydrogen bonded and/or stacked with neighboring molecules in the crystal lattice. In addition, residues A20 to A25, which comprise Loop 2, are only weakly bound to Stem 1. Conversely, Loop 1 contains only two residues, both of 2427


Article


Table 12. Backbone Conformational Suites in Loop Ea

residue A25 (A24pA25 step), which is only weakly bound to Stem 1, disappeared very quickly in all of the simulations. Finally, the 1c substate at the terminal residue G28 was also unstable in our MD simulations. However, as a terminal residue, it is expected to be quite flexible. The situation with residues C8+, C17, and A23 is different. C8+ and C17 are strongly bound to the stable Stem 1. All three nucleobases, and the bases immediately upstream and downstream, are part of the stable pseudoknot structure and form no contacts with neighboring molecules in the X-ray structure. We therefore expected these conformations to be stable. However, all three 1c substates were converted to the canonical 1a conformation during the first few nanoseconds of MD simulations with both tested force fields (Tables 9 and 10

a

354D structure contains magnesium cations that occupy the central part of the major groove near the GpA steps (see Methods). However, none of these ions interact directly with the GpA steps. Although the GpA steps can return back to their initial 1e geometry, they do so only very rarely for a negligible portion of the full simulation time. As such, the 1e substates in Loop E are unstable in our MD simulations. 2DD2 NMR Duplex. The 2DD2 structures of the NMR ensemble contain multiple α/γ = t/t conformations in the GAA triplet of noncanonical base pairs in the middle of the 5′GGU GAA GGCU/PCCG AAG CCG5′ duplex. However, a detailed analysis of the 27 lowest energy structures deposited in the PDB database revealed that canonical α/γ = g−/g+ combinations are also strongly populated in the alternative structures, as well as other conformers (Table 13). This may

Table 9. Backbone Conformations of the BYWV Pseudoknot in the Initial X-ray Structure (1L2X) and Our MD Simulations χOL3

X-ray

χOL3 + εζβOL1

location

suite

suite

%

suite

%

G7pC8+ G16pC17 C22pA23 C8+pA9

1c 1c 1c 5d

1a 1a 1a 5d

∼100 >90 ∼100 ∼100

1a 1a 1a 5d

∼100 >90 >90 ∼100

Table 13. α/γ Conformers of the 27 Lowest Energy NMR Models of the Mismatched 2DD2 Duplex α/γ in 27 lowest NMR models

Table 10. Backbone Conformational Suites in the BYWV Pseudoknota

a


Table 11. Backbone Conformations of Loop E in the Initial X-ray Structure (354D) and Our MD Simulations χOL3

χOL3 + εζβOL1

location

suite

suite

%

suite

%

G72pA73, G98pA99

1e

1a

∼100

1a

∼100

location

α/γ in model 1

t/t

g−/g+

g+/t

U3pG4 G4pA5 A5pA6 A15pA16 G14pA15 C13pG14

t/t (1c) t/t (1t) t/t (nc4) t/t (1c) g+/t (nc5) t/t (1t)

11 4 2 11 8 14

16 18 5 4 10 12

20 12 8 1

t/g+

g−/t

5

1

indicate either the dynamic nature of the backbone or uncertainty in the NMR data interpretation. We should note that determination of the α/γ conformation from NMR measurements is difficult. The γ angle can be determined from the 3J couplings, but these are not always measured and were not available for 2DD2. The α angle is even more problematic. Although explicit information on the α angle can in principle be obtained through NMR spectroscopy,76 it is difficult to measure and thus is usually unavailable, as is the case for this structure. Moreover, the available NOE data are not sufficient to determine the γ angle for the 2DD2 structure. Consequently, no direct information on the α/γ combinations of this duplex could be gathered from the primary NMR data, and their assignment is uncertain. Unfortunately, this is rather common with RNA NMR structures. Thus, although it would be exceptionally useful to have solution reference data for the RNA backbone structure, the limitations of the NMR techniques do not allow for regular determination of all backbone dihedral angles, except for highly sophisticated measurements, such as those applied to elucidation of the UUCG structure by Nozinovic et al.42 Partial information on the α/γ backbone substates in noncanonical RNAs can also be found in several works of Turner and co-workers.63,77 Our initial MD simulations of the 2DD2 duplex used the lowest energy model as the starting geometry. All γ = t substates on the GAA triplet transitioned to the canonical α/γ = g−/g+ 1a conformation during the first few nanoseconds in three independent simulations, and remained stable for the

and Figure S8). Interestingly, another noncanonical α/γ = g+/g+ substate assigned to suite 5d on residue A9 (C8pA9 step, Figure S1) remained stable in all simulations. Although the A9 residue is not involved in hydrogen bonding within the pseudoknot structure, and despite its participation in intermolecular hydrogen bonding interactions, it appears to be well stabilized by stacking with the preceding C8+ residue and thus represents a stable feature of this pseudoknot structure. Loop E. Bacterial 5S rRNA Loop E is a widely studied51,73,75 duplexlike structure with multiple non-WC base pairs (Figure 2). The backbone adopts mostly canonical A-RNA conformations, except for the GA steps G72pA73 and G98pA99, in which the guanines are involved in trans Hoogsteen/sugar edge (sheared) GA pairs, and the adenines form trans WC/ Hoogsteen pairs and adopt a 1e conformation with γ = t and β values in the g+ region (Figure S1). This conformation can be easily accommodated in the A-RNA duplex structure. In MD simulations with both the χOL3 and χOL3 + εζβOL1 force fields, these residues transitioned from the 1e suite to 1a within the first few nanoseconds (Figure S9 and Tables 11 and 12). The

X-ray


2428


Article


above a G18U13 Wobble pair. In our simulations, the 1g suite conformation in the GUAA tetraloop is stable in all trajectories and force fields, like the similar conformation in the GAGA tetraloop. However, the 1c crankshaft conformation disappeared within the first few nanoseconds of all of our simulations (see Tables 16 and 17 and Figure S12). In contrast to the case

remainder of the simulated period (see Figure S10). Thus, the simulated structure is pure A-DNA and does not exhibit dynamic interconversions between different backbone substates. If the backbone of 2DD2 was dynamic, as suggested by the number of structures consistent with NMR data, observation of the pure canonical form in the MD simulations would support our hypothesis of an excessively destabilized γ = t substate (see below). However, we prefer not to overemphasize this result because of the above-explained ambiguity in the NMR data. Another interesting feature of the 2DD2 duplex is that the C2′-pucker on residue A5, which was identified using NMR Jcoupling data, was fully stable in the χOL3 + εζβOL1 force field but completely unstable in the χOL3 force field. This indicates that χOL3 + εζβOL1 may, in some cases, provide a more realistic description of certain noncanonical features of RNA molecules, as was also reported for the catalytic center of the hairpin ribozyme.78 Parallel Guanine Quadruplex. The stems of parallel RNA quadruplexes usually have canonical backbone conformations. However, noncanonical conformations are present in the 4XK0 crystal near the ends of the quadruplex stems, in the G1pG2 base pair steps. There are four identical substates, one in each of the identical strands (Figure 2), which could almost be assigned to the 4n suite except that their ε angles are around −90°, whereas the standard value for the 4n suite is −133° (Figure S1). Consequently, these substates are denoted 4n!ε. Because both G1 and G2 are firmly bound in the stem of the quadruplex and not involved in any intermolecular contacts, we expected this noncanonical substate to be stable in our MD simulations. However, in practice, it disappeared within the first few nanoseconds, being converted to a nonclassified substate (nc6) with χOL3 and a different nonclassified substate (nc7) with χOL3 + εζβOL1 (see Tables 14 and 15 and Figure S11). Thus, although the noncanonical α/γ combination remains unchanged, the precise backbone conformation is unstable in both force fields.

Table 16. Backbone Conformations of the RNA Aptamer in the Initial X-ray Structure (2JWV) and Our MD Simulations

χOL3

a

suite

suite

%

suite

%

G1pG2

4n!ε

nc6

∼100

nc7

∼100

%

suite

%

G14pU15 G18pG19

1g 1c

1g 1a

∼100 ∼100

1g 1a

∼100 ∼100


■

DISCUSSION We have analyzed the capability of contemporary MD simulations to reproduce the noncanonical α/γ backbone substates present in experimental RNA structures. Only some of the noncanonical α/γ substates in our test structures were stable in simulations with the tested force fields. The stabilities of individual substates, sorted according to their α/γ combinations, are listed in Table 18. The 2DD2 results are not shown because backbone conformation was not reliably determined in the initial NMR structure. Because similar results were obtained with both force fields, we only show those obtained with χOL3, which yielded more stable simulations overall with our test set of RNA molecules. First, we should note that not all of the native backbone conformations in our test set were expected to be stable in MD simulations. Some of the X-ray structures exhibit intermolecular contacts that may stabilize certain conformations, which would be unstable in a molecule moving freely in solution. An example is the pseudoknot structure 1L2X, in which residues U13 and G19 are involved in extensive hydrogen bonding and stacking with surrounding molecules in the crystal lattice but have no stabilizing contacts with their parent pseudoknot. After the loss

Table 15. Backbone Conformational Suites in the Guanine Quadruplexa

a

suite

of the 2DD2 mismatched duplex, the 1c conformation is present in all 10 of the structural models that have been derived from NMR data, but was unstable in our MD simulations. Still, we should reemphasize that this α/γ backbone conformation was not derived directly from the primary NMR data, and even in this case, it remains a result of the refinement protocol. Even though the remaining part of the anti-NFκB RNA aptamer contains no other noncanonical α/γ substates and we did not investigate it in detail, we should note that it was not fully stable in any of our MD simulations. The region around the cross-strand base stacking of residues G8 and G23 in the internal loop underwent significant structural rearrangements in simulations with both tested force fields. However, because no α/γ reconformations are involved in this process and the region is relatively far from the investigated backbone conformations, we did not analyze this structural change further.

χOL3 + εζβOL1

location

χOL3 + εζβOL1

suite

Table 17. Backbone Conformational Suites in the RNA Aptamera

Table 14. Backbone Conformations of the Guanine Quadruplex in the Initial X-ray Structure (4XK0) and Our MD Simulations X-ray

χOL3

NMR location


NMR Structure of the RNA Aptamer. The NMR structure of the anti-NFκB RNA aptamer features two noncanonical α/γ substates. One is located in the GUAA tetraloop on the G14pU15 step with an α/γ = t/g+ conformation belonging to the 1g suite (Figure S1). The other belongs to the 1c suite and is located in the stem at residue G19, which is part of a canonical G19C12 WC pair that is located immediately 2429


Article

The Journal of Physical Chemistry B Table 18. Stability of α/γ Conformational Substates in MD Simulations with the χOL3 Force Fielda

crankshaft conformer occurs very frequently and is found in various types of RNA structures.8 As noted above, it is also easily distinguishable in X-ray structures because it pushes the neighboring phosphates about 1 Å further apart than they are in the canonical 1a conformation, which is clearly visible in X-ray densities. We would therefore expect this conformer to be preserved (or at least significantly populated) in the MD simulations of at least some structures. Strikingly, as shown in Table 18, none of the simulated 1c conformers remained stable in our simulations. Instead, 1c always transitioned to the canonical 1a within the first few nanoseconds. Although the 1c conformer reappeared a few times in multiple simulations, it always disappeared almost immediately. This indicates that the transition between 1a and 1c occurs easily on a 1 μs time scale. The total population of this substate in MD simulations was much lower than the population of almost ∼5% seen in the Xray database (see above). Thus, the 1c conformer is kinetically accessible but is thermodynamically unstable in the MD simulations. It seems very likely that this substate is described poorly by current AMBER force fields. To better understand the instability of the 1c conformers, we also analyzed other noncanonical α/γ substates in our test set. We initially focused on substates in which one of the α and γ angles is in the trans conformation, and the other is either g+ or g−. Whereas only two of the five substates featuring the γ = t conformation were fully stable (see Table 18), all five conformations featuring the α = t and γ = g+ conformation were stable, as were the two additional α/γ = g+/g+ conformations. The instability of the 1c (α/γ = t/t) conformations is thus probably due to the γ dihedral potential. As mentioned in the Introduction, the χOL3 force field uses the bsc0 correction, which was introduced to remove excessive α/γ = g+/t transitions in DNA and primarily operates by destabilizing the γ = t part of the potential. It is tempting to suggest that this destabilization was excessive and is responsible for the unrealistically low stability of noncanonical conformations with γ = t in simulations conducted with χOL3 and the χOL3 + εζβOL1 force field, which use the same bsc0 α/γ potential. This indicates a need for further refinement of the α/γ potential. We should note, however, that other force field terms, such as van der Waals or electrostatic, may contribute to instability of the noncanonical configurations as well. Also, we would like to emphasize that in spite of the above mentioned shortcomings, the bsc0 correction provides a much improved description of the A-RNA duplexes compared to that of the original ff99 potential.21,39 This is mainly because population of the α/γ = t/t substates is unnaturally high in simulations without the bsc0 α/γ correction even with modified glycosidic torsions and also for the ff99 + χOL3 combination (unpublished results). It remains to be seen if it is possible to reduce penalization of the γ−trans states without negatively affecting some other structural features in the simulations. We also tested the effect of augmenting the current RNA χOL3 force field by the additional εζβOL1 refinements, which were originally derived for DNA. This variant has been abbreviated as the χOL3 + εζβOL1 force field. χOL3 + εζβOL1 provided stable simulations of canonical A-RNA stems, although the extent of end-fraying in the AU-rich structure 1RNA and the GC-rich duplex 1QC0′ was more extensive than with χOL3. The structures of the GAGA tetraloop, Loop E, parallel quadruplex and mismatched duplex were generally stable with this force field, exhibiting similar behavior to that seen with χOL3. Improvements were observed for the C2′-endo

a

Less reliable results in brackets. bInvolved in crystal contacts. Unstable even in crystal environment.65 dLower population than expected.

c

of their intermolecular contacts, these residues would be expected to sample a range of conformations while they remain bulged out of the main RNA structure and exposed to the solvent. This can be considered natural behavior, and we will not draw any conclusion from the instability of these conformations. In Table 18, these substates are put in parentheses to indicate that they are less reliable. Also, the NMR structures may provide a less reliable reference (with an exception of the excellent UUCG structure42 in which backbone angles were determined explicitly) and thus we did not include the 2DD2 duplex results in Table 18. The noncanonical α/γ backbone combination was not derived from the primary NMR data for the 2JWV aptamer structure; however, the deposited structural ensemble is more consistent in this case. Further, there are noncanonical α/γ substates that we expected to be stable in our simulations, such as those found repeatedly in well-resolved X-ray structures that are not influenced by intermolecular contacts. In particular, the 1c 2430


Article

The Journal of Physical Chemistry B sugar pucker in the mismatched 2DD2 duplex and the #a suite in the sarcin−ricin loop, possibly due to stabilization of the ζ = t region by the εζOL1 potential. However, in several cases the χOL3 + εζβOL1 potential produced significant deviations from the native geometry. In SRL, the G2655 residue of the GpU platform in the S-turn shifted from its native location in two out of three MD simulations and bulged out into the solvent. Similarly, the loop structure of the UUCG tetraloop was lost in all three 1 μs χOL3 + εζβOL1 simulations because the guanine of the GU pair bulged out of the tetraloop structure. Although the UUCG tetraloop also lost its accurate backbone structure with the χOL3 force field (the guanine of the GU pair did not bulge out of the tetraloop structure on our time scale), its reconformation was faster with χOL3 + εζβOL1. The faster reconformation could possibly be caused by reduced rotational barriers due to the εζOL1 potential of the later force field. The region around the cross-strand base stacking of residues G8 and G23 in the internal loop of the anti-NFκB RNA aptamer also undergoes significant structural rearrangements with χOL3 + εζβOL1, but the same thing happens in simulations with the χOL3 force field. Some of the instabilities observed with χOL3 + εζβOL1 could be of kinetic origin due to faster reconformations resulting from the use of the εζOL1 potential, which could lead to quicker shifts to more stable (as predicted by the force field) non-native conformations. However, we cannot exclude the possibility that the χOL3 + εζβ OL1 modifications also thermodynamically stabilize some non-native states compared to χOL3. Thus, after all of the tests performed, we do not recommend the use of the χOL3 + εζβOL1 potential for RNA simulations, although the εζβOL1 refinements are very useful for DNA simulations. Instead, we recommend staying with the current χOL3 RNA force field.

would like to emphasize that we do not suggest excluding the bsc0 α/γ refinement from the current RNA force fields because it provides significant improvements as compared with the original ff99 α/γ potential.21,39 Finding a better solution for the γ-potential may be very challenging and may require concerted modifications of some other force field terms. The εζOL1 and βOL1 modifications implemented in the extended χOL3 + εζβOL1 force field version yield stable simulations for canonical and noncanonical duplexes, the GAGA tetraloop, the parallel RNA quadruplex, and Loop E, and even improvements in the description of some noncanonical backbone geometries (see also the previously reported improvements achieved for the reactive center of the hairpin ribozyme with εζOL1).78 However, we observed notable deviations from the native X-ray and NMR structures of several systems with the χOL3 + εζβOL1 force field, including the S-turn of the sarcin−ricin loop, the anti-NFκB RNA aptamer, and the UUCG tetraloop. Some of these systems are also unstable with the χOL3 force field, albeit to a lesser extent and on a longer time scale. The more rapid undesired structural changes observed with χOL3 + εζβOL1 can potentially be explained as a consequence of lowered barriers resulting from the use of the εζOL1 potential as compared to the εζ potential of ff99bsc0.34,35 Although the εζOL1 and βOL1 modifications are definitely very useful for DNA simulations, based on our results, we do not currently recommend using them for RNA simulations. The exception could be cases where noncanonical ε, ζ, and β values are known to be present or are expected to play an important role.

■

ASSOCIATED CONTENT

* Supporting Information

■

S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b00262. Time series of backbone dihedral angles for noncanonical backbone substates of simulated systems (PDF)

CONCLUSIONS Rotameric conformational substates involving noncanonical α/γ backbone dihedral angle combinations were analyzed in several prototypical RNA structures, including highly stable tetraloops, an S-turn, Loop E, a parallel quadruplex, and a pseudoknot. Our focus was mainly on the α/γ = t/t (1c) conformation, which is the most populated noncanonical α/γ substate in crystal structures of RNA molecules. Molecular dynamics simulations of the selected NMR structures and wellresolved X-ray structures were performed using two force fields, the currently recommended ff99bsc0χOL3 (χOL3) and a combination of χOL3 with our recent modifications for DNA, χOL3 + εζβOL1. Although the results obtained with these two AMBER force fields were not identical, several common features were evident. In general, many of the native noncanonical α/γ substates were not stable in either force field. Strikingly, none of the simulated α/γ = t/t backbone substates were stable in any of our simulations. All native 1c suite conformations were converted into 1a conformations very rapidly, typically within a few nanoseconds. This may highlight a deficiency in the tested force fields. By analyzing the stability of various α/γ conformations, we identified a possible problem in the bsc0 correction’s description of the γ potential, which seems to excessively destabilize the γ = t region. In keeping with this hypothesis, it has previously been suggested that the bsc0 correction excessively destabilizes this region in Z-DNA, whose backbone conformation is also described poorly by current force fields.34,35 These results reveal a potentially important artifact in current AMBER force fields that should be considered in future force field refinements. Nevertheless, we

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +420 585 63 4760. ORCID

Petr Jurečka: 0000-0002-3741-3672 Michal Otyepka: 0000-0002-1066-5677 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by grant 14-29874P (M.Z.) from the Grant Agency of the Czech Republic. Further funding was provided by project LO1305 of the Ministry of Education, Youth and Sports of the Czech Republic (P.J., P.B., J.S., and M.O.).

■

REFERENCES

(1) Tucker, B. J.; Breaker, R. R. Riboswitches as Versatile Gene Control Elements. Curr. Opin. Struct. Biol. 2005, 15, 342−348. (2) Serganov, A.; Patel, D. J. Ribozymes, Riboswitches and Beyond: Regulation of Gene Expression without Proteins. Nat. Rev. Genet. 2007, 8, 776−790. (3) Steitz, T. A. A Structural Understanding of the Dynamic Ribosome Machine. Nat. Rev. Mol. Cell Biol. 2008, 9, 242−253.

2431


Article

The Journal of Physical Chemistry B (4) Saenger, W. Principles of Nucleic Acid Structure; Springer-Verlag: New York, 1984; pp 9−28. (5) Schneider, B.; Moravek, Z.; Berman, H. M. RNA Conformational Classes. Nucleic Acids Res. 2004, 32, 1666−1677. (6) Murray, L. J. W.; Arendall, W. B.; Richardson, D. C.; Richardson, J. S. RNA Backbone Is Rotameric. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 13904−13909. (7) Hershkovitz, E.; Tannenbaum, E.; Howerton, S. B.; Sheth, A.; Tannenbaum, A.; Williams, L. D. Automated Identification of RNA Conformational Motifs: Theory and Application to the HM LSU 23S rRNA. Nucleic Acids Res. 2003, 31, 6249−6257. (8) Richardson, J. S.; Schneider, B.; Murray, L. W.; Kapral, G. J.; Immormino, R. M.; Headd, J. J.; Richardson, D. C.; Ham, D.; Hershkovits, E.; Williams, L. D.; et al. RNA Backbone: Consensus AllAngle Conformers and Modular String Nomenclature (an RNA Ontology Consortium Contribution). RNA 2008, 14, 465−481. (9) Č ech, P.; Kukal, J.; Cerny, J.; Schneider, B.; Svozil, D. Automatic Workflow for the Classification of Local DNA Conformations. BMC Bioinf. 2013, 14, 205. (10) Ben Imeddourene, A.; Elbahnsi, A.; Gueroult, M.; Oguey, C.; Foloppe, N.; Hartmann, B. Simulations Meet Experiment to Reveal New Insights into DNA Intrinsic Mechanics. PLoS Comput. Biol. 2015, 11, No. e1004631. (11) Djuranovic, D.; Hartmann, B. DNA Fine Structure and Dynamics in Crystals and in Solution: The Impact of BI/BII Backbone Conformations. Biopolymers 2004, 73, 356−368. (12) Svozil, D.; Kalina, J.; Omelka, M.; Schneider, B. DNA Conformations and Their Sequence Preferences. Nucleic Acids Res. 2008, 36, 3690−3706. (13) Várnai, P.; Zakrzewska, K. DNA and Its Counterions: A Molecular Dynamics Study. Nucleic Acids Res. 2004, 32, 4269−4280. (14) Svozil, D.; Sponer, J. E.; Marchan, I.; Perez, A.; Cheatham, T. E.; Forti, F.; Luque, F. J.; Orozco, M.; Sponer, J. Geometrical and Electronic Structure Variability of the Sugar-Phosphate Backbone in Nucleic Acids. J. Phys. Chem. B 2008, 112, 8188−8197. (15) Várnai, P.; Djuranovic, D.; Lavery, R.; Hartmann, B. Alpha/ Gamma Transitions in the B-DNA Backbone. Nucleic Acids Res. 2002, 30, 5398−5406. (16) Pérez, A.; Marchan, I.; Svozil, D.; Sponer, J.; Cheatham, T. E.; Laughton, C. A.; Orozco, M. Refinenement of the Amber Force Field for Nucleic Acids: Improving the Description of Alpha/Gamma Conformers. Biophys. J. 2007, 92, 3817−3829. (17) Besseová, I.; Otyepka, M.; Reblova, K.; Sponer, J. Dependence of A-RNA Simulations on the Choice of the Force Field and Salt Strength. Phys. Chem. Chem. Phys. 2009, 11, 10701−10711. (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. A 2nd Generation Force-Field for the Simulation of Proteins, Nucleic-Acids, and Organic-Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (19) Cheatham, T. E., III; Cieplak, P.; Kollman, P. A. A Modified Version of the Cornell Et Al. Force Field with Improved Sugar Pucker Phases and Helical Repeat. J. Biomol. Struct. Dyn. 1999, 16, 845−862. (20) Wang, J. M.; Cieplak, P.; Kollman, P. A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21, 1049−1074. (21) Zgarbová, M.; Otyepka, M.; Sponer, J.; Mladek, A.; Banas, P.; Cheatham, T. E., 3rd; Jurecka, P. Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. J. Chem. Theory Comput. 2011, 7, 2886−2902. (22) Foloppe, N.; MacKerell, A. D. All-Atom Empirical Force Field for Nucleic Acids: I. Parameter Optimization Based on Small Molecule and Condensed Phase Macromolecular Target Data. J. Comput. Chem. 2000, 21, 86−104. (23) MacKerell, A. D.; Banavali, N. K. All-Atom Empirical Force Field for Nucleic Acids: II. Application to Molecular Dynamics

Simulations of DNA and Rna in Solution. J. Comput. Chem. 2000, 21, 105−120. (24) Hart, K.; Foloppe, N.; Baker, C. M.; Denning, E. J.; Nilsson, L.; MacKerell, A. D. Optimization of the Charmm Additive Force Field for DNA: Improved Treatment of the BI/BII Conformational Equilibrium. J. Chem. Theory Comput. 2012, 8, 348−362. (25) Denning, E. J.; Priyakumar, U. D.; Nilsson, L.; Mackerell, A. D. Impact of 2′-Hydroxyl Sampling on the Conformational Properties of Rna: Update of the Charmm All-Atom Additive Force Field for RNA. J. Comput. Chem. 2011, 32, 1929−1943. (26) Savelyev, A.; MacKerell, A. D. Balancing the Interactions of Ions, Water, and DNA in the Drude Polarizable Force Field. J. Phys. Chem. B 2014, 118, 6742−6757. (27) Savelyev, A.; MacKerell, A. D. All- Atom Polarizable Force Field for DNA Based on the Classical Drude Oscillator Model. J. Comput. Chem. 2014, 35, 1219−1239. (28) Lemkul, J. A.; Huang, J.; Roux, B.; MacKerell, A. D. An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications. Chem. Rev. 2016, 116, 4983−5013. (29) Gil-Ley, A.; Bottaro, S.; Bussi, G. Empirical Corrections to the Amber RNA Force Field with Target Metadynamics. J. Chem. Theory Comput. 2016, 12, 2790−2798. (30) Chen, A. A.; Garcia, A. E. High-Resolution Reversible Folding of Hyperstable RNA Tetraloops Using Molecular Dynamics Simulations. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 16820−16825. (31) Cesari, A.; Gil-Ley, A.; Bussi, G. Combining Simulations and Solution Experiments as a Paradigm for RNA Force Field Refinement. J. Chem. Theory Comput. 2016, 12, 6192−6200. (32) Yildirim, I.; Stern, H. A.; Kennedy, S. D.; Tubbs, J. D.; Turner, D. H. Reparameterization of RNA Chi Torsion Parameters for the Amber Force Field and Comparison to NMR Spectra for Cytidine and Uridine. J. Chem. Theory Comput. 2010, 6, 1520−1531. (33) Yildirim, I.; Kennedy, S. D.; Stern, H. A.; Hart, J. M.; Kierzek, R.; Turner, D. H. Revision of Amber Torsional Parameters for RNA Improves Free Energy Predictions for Tetramer Duplexes with GC and iGiC Base Pairs. J. Chem. Theory Comput. 2012, 8, 172−181. (34) Zgarbová, M.; Luque, F. J.; Sponer, J.; Cheatham, T. E., 3rd; Otyepka, M.; Jurecka, P. Toward Improved Description of DNA Backbone: Revisiting Epsilon and Zeta Torsion Force Field Parameters. J. Chem. Theory Comput. 2013, 9, 2339−2354. (35) Zgarbová, M.; Sponer, J.; Otyepka, M.; Cheatham, T. E., 3rd; Galindo-Murillo, R.; Jurecka, P. Refinement of the Sugar-Phosphate Backbone Torsion Beta for Amber Force Fields Improves the Description of Z- and B-DNA. J. Chem. Theory Comput. 2015, 11, 5723−5736. (36) Zgarbová, M.; Luque, F. J.; Sponer, J.; Otyepka, M.; Jurecka, P. A Novel Approach for Deriving Force Field Torsion Angle Parameters Accounting for Conformation-Dependent Solvation Effects. J. Chem. Theory Comput. 2012, 8, 3232−3242. (37) Galindo-Murillo, R.; Robertson, J. C.; Zgarbova, M.; Sponer, J.; Otyepka, M.; Jurecka, P.; Cheatham, T. E., 3rd Assessing the Current State of Amber Force Field Modifications for DNA. J. Chem. Theory Comput. 2016, 12, 4114. (38) Réblová, K.; Lankas, F.; Razga, F.; Krasovska, M. V.; Koca, J.; Sponer, J. Structure, Dynamics, and Elasticity of Free 16s Rrna Helix 44 Studied by Molecular Dynamics Simulations. Biopolymers 2006, 82, 504−520. (39) Bešsě ová, I.; Banas, P.; Kuhrova, P.; Kosinova, P.; Otyepka, M.; Sponer, J. Simulations of A-RNA Duplexes. The Effect of Sequence, Solute Force Field, Water Model, and Salt Concentration. J. Phys. Chem. B 2012, 116, 9899−9916. (40) Klosterman, P. S.; Shah, S. A.; Steitz, T. A. Crystal Structures of Two Plasmid Copy Control Related RNA Duplexes: An 18 Base Pair Duplex at 1.20 A Resolution and a 19 Base Pair Duplex at 1.55 A Resolution. Biochemistry 1999, 38, 14784−14792. (41) Dockbregeon, A. C.; Chevrier, B.; Podjarny, A.; Johnson, J.; Debear, J. S.; Gough, G. R.; Gilham, P. T.; Moras, D. Crystallographic 2432


Article

The Journal of Physical Chemistry B Structure of an RNA Helix - [U(UA)6A]2. J. Mol. Biol. 1989, 209, 459−474. (42) Nozinovic, S.; Furtig, B.; Jonker, H. R. A.; Richter, C.; Schwalbe, H. High-Resolution NMR Structure of an RNA Model System: The 14-mer cUUCGg Tetraloop Hairpin RNA. Nucleic Acids Res. 2010, 38, 683−694. (43) Correll, C. C.; Swinger, K. Common and Distinctive Features of GNRA Tetraloops Based on a GUAA Tetraloop Structure at 1.4 Angstrom Resolution. RNA 2003, 9, 355−363. (44) Correll, C. C.; Wool, I. G.; Munishkin, A. The Two Faces of the Escherichia Coli 23 S rRNA Sarcin/Ricin Domain: The Structure at 1.11 Angstrom Resolution. J. Mol. Biol. 1999, 292, 275−287. (45) Egli, M.; Minasov, G.; Su, L.; Rich, A. Metal Ions and Flexibility in a Viral RNA Pseudoknot at Atomic Resolution. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 4302−4307. (46) Correll, C. C.; Freeborn, B.; Moore, P. B.; Steitz, T. A. Metals, Motifs, and Recognition in the Crystal Structure of a 5s Rrna Domain. Cell 1997, 91, 705−712. (47) Chen, G.; Kennedy, S. D.; Qiao, J.; Krugh, T. R.; Turner, D. H. An Alternating Sheared AA Pair and Elements of Stability for a Single Sheared Purine-Purine Pair Flanked by Sheared GA Pairs in RNA. Biochemistry 2006, 45, 6889−6903. (48) Chen, M. C.; Murat, P.; Abecassis, K.; Ferre-D’Amare, A. R.; Balasubramanian, S. Insights into the Mechanism of a G-QuadruplexUnwinding Deah-Box Helicase. Nucleic Acids Res. 2015, 43, 2223− 2231. (49) Reiter, N. J.; Maher, L. J.; Butcher, S. E. DNA Mimicry by a High-Affinity Anti-Nf-KappaB RNA Aptamer. Nucleic Acids Res. 2008, 36, 1227−1236. (50) Zgarbová, M.; Otyepka, M.; Sponer, J.; Lankas, F.; Jurecka, P. Base Pair Fraying in Molecular Dynamics Simulations of DNA and RNA. J. Chem. Theory Comput. 2014, 10, 3177−3189. (51) Auffinger, P.; Bielecki, L.; Westhof, E. The Mg2+ Binding Sites of the 5S rRNA Loop E Motif as Investigated by Molecular Dynamics Simulations. Chem. Biol. 2003, 10, 551−561. (52) Joung, I. S.; Cheatham, T. E. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112, 9020−9041. (53) Joung, I. S.; Cheatham, T. E. Molecular Dynamics Simulations of the Dynamic and Energetic Properties of Alkali and Halide Ions Using Water-Model-Specific Ion Parameters. J. Phys. Chem. B 2009, 113, 13279−13290. (54) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (55) Case, D. A.; Babin, V.; Berryman, J. T.; Betz, R. M.; Cai, Q.; Cerutti, D. S.; Cheatham, T. E.; Darden, T. A.; Duke, R. E.; Gohlke, H.; et al. Amber 14; University of California: San Francisco, 2014. (56) Feenstra, K. A.; Hess, B.; Berendsen, H. J. C. Improving Efficiency of Large Time-Scale Molecular Dynamics Simulations of Hydrogen-Rich Systems. J. Comput. Chem. 1999, 20, 786−798. (57) Hopkins, C. W.; Le Grand, S.; Walker, R. C.; Roitberg, A. E. Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning. J. Chem. Theory Comput. 2015, 11, 1864−1874. (58) Götz, A. W.; Williamson, M. J.; Xu, D.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with Amber on GPUs. 1. Generalized Born. J. Chem. Theory Comput. 2012, 8, 1542−1555. (59) Roe, D. R.; Cheatham, T. E., 3rd Ptraj and Cpptraj: Software for Processing and Analysis of Molecular Dynamics Trajectory Data. J. Chem. Theory Comput. 2013, 9, 3084−3095. (60) Banás,̌ P.; Sklenovsky, P.; Wedekind, J. E.; Sponer, J.; Otyepka, M. Molecular Mechanism of Preq(1) Riboswitch Action: A Molecular Dynamics Study. J. Phys. Chem. B 2012, 116, 12721−12734. (61) Allnér, O.; Nilsson, L.; Villa, A. Magnesium Ion-Water Coordination and Exchange in Biomolecular Simulations. J. Chem. Theory Comput. 2012, 8, 1493−1502. (62) Pan, B. C.; Mitra, S. N.; Sundaralingam, M. Crystal Structure of an RNA 16-mer Duplex r(GCAGAGUUAAAUCUGC)(2) with

Nonadjacent G(syn).A(+)(anti) Mispairs. Biochemistry 1999, 38, 2826−2831. (63) Wu, M.; Turner, D. H. Solution Structure of (rGCGGACGC) (2) by Two-Dimensional NMR and the Iterative Relaxation Matrix Approach. Biochemistry 1996, 35, 9677−9689. (64) Šponer, J.; Krepl, M.; Banas, P.; Kuhrova, P.; Zgarbova, M.; Jurecka, P.; Havrila, M.; Otyepka, M. How to Understand Atomistic Molecular Dynamics Simulations of RNA and Protein-RNA Complexes? Wiley Interdiscip. Rev.: RNA 2016, DOI: 10.1002/ wrna.1405. (65) Liu, C. M.; Janowski, P. A.; Case, D. A. All-Atom Crystal Simulations of DNA and RNA Duplexes. Biochim. Biophys. Acta, Gen. Subj. 2015, 1850, 1059−1071. (66) Garcia, A. E.; Paschek, D. Simulation of the Pressure and Temperature Folding/Unfolding Equilibrium of a Small RNA Hairpin. J. Am. Chem. Soc. 2008, 130, 815−817. (67) Miner, J. C.; Chen, A. A.; Garcia, A. E. Free-Energy Landscape of a Hyperstable RNA Tetraloop. Proc. Natl. Acad. Sci. U.S.A. 2016, 113, 6665−6670. (68) Banas, P.; Hollas, D.; Zgarbova, M.; Jurecka, P.; Orozco, M.; Cheatham, T. E.; Sponer, J.; Otyepka, M. Performance of Molecular Mechanics Force Fields for RNA Simulations: Stability of UUCG and GNRA Hairpins. J. Chem. Theory Comput. 2010, 6, 3836−3849. (69) Haldar, S.; Kuhrova, P.; Banas, P.; Spiwok, V.; Sponer, J.; Hobza, P.; Otyepka, M. Insights into Stability and Folding of GNRA and UNCG Tetra Loops Revealed by Microsecond Molecular Dynamics and Well-Tempered Metadynamics. J. Chem. Theory Comput. 2015, 11, 3866−3877. (70) Bergonzo, C.; Henriksen, N. M.; Roe, D. R.; Cheatham, T. E. Highly Sampled Tetranucleotide and Tetraloop Motifs Enable Evaluation of Common Rna Force Fields. RNA 2015, 21, 1578−1590. (71) Bottaro, S.; Banas, P.; Sponer, J.; Bussi, G. Free Energy Landscape of GAGA and UUCG RNA Tetraloops. J. Phys. Chem. Lett. 2016, 7, 4032−4038. (72) Spacková, N.; Sponer, J. Molecular Dynamics Simulations of Sarcin-Ricin rRNA Motif. Nucleic Acids Res. 2006, 34, 697−708. (73) Kruse, H.; Havrila, M.; Sponer, J. Qm Computations on Complete Nucleic Acids Building Blocks: Analysis of the Sarcin-Ricin RNA Motif Using DFT-D3, HF-3c, PM6-D3H, and MM Approaches. J. Chem. Theory Comput. 2014, 10, 2615−2629. (74) Havrila, M.; Reblova, K.; Zirbel, C. L.; Leontis, N. B.; Sponer, J. Isosteric and Nonisosteric Base Pairs in RNA Motifs: Molecular Dynamics and Bioinformatics Study of the Sarcin-Ricin Internal Loop. J. Phys. Chem. B 2013, 117, 14302−14319. (75) Réblová, K.; Spackova, N.; Stefl, R.; Csaszar, K.; Koca, J.; Leontis, N. B.; Sponer, J. Non-Watson-Crick Basepairing and Hydration in Rna Motifs: Molecular Dynamics of 5s rRNA Loop E. Biophys. J. 2003, 84, 3564−3582. (76) Nozinovic, S.; Richter, C.; Rinnenthal, J.; Furtig, B.; DuchardtFerner, E.; Weigand, J. E.; Schwalbe, H. Quantitative 2d and 3d Gamma-HCP Experiments for the Determination of the Angles Alpha and Zeta in the Phosphodiester Backbone of Oligonucleotides. J. Am. Chem. Soc. 2010, 132, 10318−10329. (77) Hammond, N. B.; Tolbert, B. S.; Kierzek, R.; Turner, D. H.; Kennedy, S. D. RNA Internal Loops with Tandem AG Pairs: The Structure of the 5′GAGU/3′UGAG Loop Can Be Dramatically Different from Others, Including 5′AAGU/3′UGAA. Biochemistry 2010, 49, 5817−5827. (78) Mlýnský, V.; Kuhrova, P.; Zgarbova, M.; Jurecka, P.; Walter, N. G.; Otyepka, M.; Sponer, J.; Banas, P. Reactive Conformation of the Active Site in the Hairpin Ribozyme Achieved by Molecular Dynamics Simulations with Epsilon/Zeta Force Field Reparametrizations. J. Phys. Chem. B 2015, 119, 4220−4229.

2433


Î³ Backbone Conformations in RNA and the Accuracy

Recommend Documents