Influence of Sequence and Covalent Modifications on Yeast tRNA Dynamics.

Modified nucleotides are prevalent in tRNA. Experimental studies reveal that these covalent modifications play an important role in tuning tRNA function. In this study, molecular dynamics (MD) simulations were used to investigate how modifications alter tRNA dynamics. The X-ray crystal structures of tRNA(Asp), tRNA(Phe), and tRNA(iMet), both with and without modifications, were used as initial structures for 333 ns explicit solvent MD simulations with AMBER. For each tRNA molecule, three independent trajectory calculations were performed, giving an aggregate of 6 μs of total MD across six molecules. The global root-mean-square deviations (RMSD) of atomic positions show that modifications only introduce significant rigidity to the global structure of tRNA(Phe). Interestingly, RMSDs of the anticodon stem-loop (ASL) suggest that modified tRNA has a more rigid structure compared to the unmodified tRNA in this domain. The anticodon RMSDs of the modified tRNAs, however, are higher than those of corresponding unmodified tRNAs. These findings suggest that the rigidity of the anticodon stem-loop is finely tuned by modifications, where rigidity in the anticodon arm is essential for tRNA translocation in the ribosome, and flexibility of the anticodon is important for codon recognition. Sugar pucker and water residence time of pseudouridines in modified tRNAs and corresponding uridines in unmodified tRNAs were assessed, and the results reinforce that pseudouridine favors the 3'-endo conformation and has a higher tendency to interact with water. Principal component analysis (PCA) was used to examine correlated motions in tRNA. Additionally, covariance overlaps of PCAs were compared for trajectories of the same molecule and between trajectories of modified and unmodified tRNAs. The comparison suggests that modifications alter the correlated motions. For the anticodon bases, the extent of stacking was compared between modified and unmodified molecules, and only unmodified tRNA(Asp) has significantly higher percentage of stacking time. Overall, the simulations reveal that the effect of covalent modification on tRNA dynamics is not simple, with modifications increasing flexibility in some regions of the structure and increasing rigidity in other regions.

Introduction

Naturally occurring covalent modifications to the standard RNA chemistry are prevalent in tRNA. In 561 sequenced tRNAs, from a wide range of organisms, modifications are detected on 11.9% of the residues.[1−3] In the yeast Saccharomyces cerevisiae, 16.2% of the residues of the 28 unique sequenced cytoplasmic tRNA species hold modifications. The range is from 7 to 17 modifications per tRNA.[1] With their frequent occurrence in tRNAs, modifications are recognized as an important device for tuning tRNA structure and function.[4−7]

tRNAs are usually 74–95 nucleotides long, and their sequences (primary structure) typically allow extensive base pairing (secondary structure) and noncanonical interactions (tertiary structure). tRNA secondary structure is nearly universally arranged in a cloverleaf shape, as first realized by Holley et al.[8] A common numbering scheme for tRNA nucleotides is used, which relies on the position of a base relative to the canonical secondary structure. On the basis of the cloverleaf structure, tRNA is composed of four regions: the anticodon stem loop (ASL, nucleotides 27–43), acceptor stem (nucleotides 1–7 and 66–76), D stem loop (nucleotides 10–25), and T stem loop containing the conserved TΨC sequence (nucleotides 49–65). The first tertiary structure was determined by X-ray crystallography studies in 1973.[9,10] The L-shaped tertiary structure is formed by two coaxial stacks of helices in the cloverleaf secondary structure. The acceptor stem and T stem coaxially stack and form one of the arms, and the other arm is formed by coaxial stacking of the D stem and anticodon stem.

The structure and conformation of several unmodified yeast tRNAs have been studied with T7 polymerase transcription. It is shown that fully modified native tRNA and unmodified tRNA have similar overall structure.[11−16] Melting experiments reveal that both modified and unmodified tRNAs are stabilized by addition of Mg2+. Unmodified tRNA, however, has more plasticity than modified tRNA at physiological Mg2+ concentration and is less stable than modified tRNA in the absence of Mg2+ or at low Mg2+ concentration.[15] A nuclease probing study found a greater accessibility of the D- and T-loops of the transcribed tRNA, indicating that the absence of modification may disrupt the tertiary structure of the D–T loop corner in tRNA.[15]

The anticodon loop is the domain that directly interacts with mRNA and the ribosome. Therefore, alteration to tRNA structure at this location by modification directly changes the interaction between tRNA and other partners of translation. This loop is a prominent location for modifications, especially at positions of 34 and 37, where 34, 35, and 36 are the positions of the anticodon nucleotides. Modifications of uridine at position 34 and purine at position 37 are ubiquitous, and only seven tRNAs in Escherichia coli are read by tRNAs without modifications at N34 or R37.[17]

There is also another way that modification can alter tRNA performance, and that is by tuning the dynamics of the tRNA. Most of the modifications that exist in the T loop and D loop are believed to contribute to dynamics.[18] The absence of the m1A58 modification from initiator tRNA (tRNAiMet) results in its degradation.[19,20] This degradation might be due to a weakening of a tertiary interaction between D and T loops that is unique to tRNAiMet.[20,21] In addition, modification can apply its influence via stabilizing effects that restrict conformational flexibility of tRNA on the angstrom scale.[18] This feature is often observed in NMR studies and molecular dynamics simulations.[22−25] For example, with a methyl group added to the 2′ hydroxyl moiety, 2′-O-methylated nucleotides increase base stacking.[26] In addition, 2′-O-methylation significantly stabilizes the 3′-endo sugar pucker conformation in pyrimidines because of the steric repulsion between the 2-carbonyl group, the 2′-O-methyl group and the 3′-phosphate group in the C2′-endo form.[26]

Pseudouridylation also changes dynamics. It is initiated by the cleavage of the N1–C1′ glycosidic bond, followed by a 180° rotation along the N3–C6 axis of the base and rejoining with the sugar by the formation of a C5–C1′ bond.[27] With this altered base-ribose linkage, an additional imino group is available as a hydrogen bond donor, while the functional groups on the Watson–Crick base pairing face remain intact. Molecular dynamics simulations and experimental data indicate that the pseudouridines (Ψs), in ASL and TΨC loops, are involved in a water mediated base-to-backbone interactions, where the additional (Ψ)N1–H imino group is used to establish an N1–H···O hydrogen bond with water. With this interaction, the conformational sampling of the Ψ nucleotide is reduced.[28−33] In the tRNALys, 3, ASL loop, the conserved Ψ39 increases the melting temperature of the ASL.[34]

With the power to investigate structural and dynamical information on macro-molecular structure in atomic detail, there is a rich history of application of molecular dynamics to study tRNA molecules. Harvey et al. reported the first investigations of the intramolecular dynamics of a tRNAPhe molecule by computer simulation.[35] In their model, solvation was implicit and the effect of the counterions was approximated by reducing the atomic charges on the phosphate groups, and the main features of the crystal structure were preserved in 12 ps of dynamics. Auffinger et al. reported a series of studies on the dynamics of a fully solvated tRNAAsp molecule and anticodon hairpin, in explicit SPC/E water.[22,33,36,37] For the whole tRNAAsp, the simulation reached 500 ps.[22] The MD simulation showed that Ψ32 stabilizes a water molecule, bridging the two adjacent C31 and Ψ32 anionic oxygen atoms through an N1–H···O hydrogen bond. McCrate et al. reported an MD simulation study on the human tRNALys anticodon stem-loop to elucidate roles of tRNA modified bases in mRNA recognition. The simulations were performed on the anticodon stem-loop using eight distinct combinations of modified bases including an unmodified ASL molecule. Each molecule was simulated for 4 ns, except the wild type molecule, which was run to 40 ns to ensure the simulation was stable. This study concluded that the ms2t6 modification at position 37 is required for maintenance of the canonical anticodon stair-stepped conformation and reduces solvent accessibility of U36, while ms2t6A37 generates hydrogen bonds across the loop, which may prevent U36 from rotating into solution and hence stabilizes the loop. In addition, Ψ39 was shown to stabilize tRNA structure through a water-mediated hydrogen-bonding network.[38]

Although experimental studies have been carried out on tRNA modification to get a better understanding of its function, there are few reports about the structural mechanism of modification for tuning tRNA activity. By comparing MD simulations between tRNA molecules with or without modification, this study investigated the structural and dynamical features to understand the effect of tRNA modification on the structure and dynamics of tRNA. There are three yeast tRNAs with structures solved by crystallography, where the tRNA was crystallized in the absence of other components of translation: tRNAAsp,[39] tRNAPhe,[40] and tRNAiMet.[21] The tRNAPhe structure used in this study is the updated and revised structure, which has a higher resolution than the prior structure.[9,10,40]

Here, the roles of modifications of yeast tRNAAsp, tRNAPhe, and tRNAiMet on the dynamical properties of the tRNAs using three separate simulations for each tRNA, both with and without modifications. Across the three simulations, a total of 1 μs of aggregate sampling was performed for each molecule. The set of three simulations for each molecule was used to evaluate the variability of observations. For differences in conformations and dynamics observed between molecules, the statistical significance was tested.

Methods

Initial Structures

tRNAAsp, tRNAPhe, and tRNAiMet were chosen for this study. Initial coordinates for modified tRNAs were obtained from the X-ray crystal structures available in the Protein Data Bank: tRNAAsp (PDB code 3TRA, resolution = 3.0 Å),[39] tRNAPhe (PDB code 1EHZ, resolution = 1.93 Å),[40] and tRNAiMet (PDB code 1YFG, resolution = 3.0 Å).[21] The initial coordinates for unmodified tRNAs were derived from the corresponding wild type tRNA crystal structures by deleting and replacing the modification group with hydrogen, or by rearranging the position of atoms in the base ring for converting pseudouridine to uridine. Mg2+ ions located in the crystal structures of tRNAAsp and tRNAPhe were not included in the simulations because Mg2+ was reported to distort simulations.[41] This is because of previous limitations in the parametrization of van der Waals parameters for Mg2+ and the inherent difficulties in representing such a highly charged ion with pairwise models.

Force Field Parameters

All MD simulations were performed with the AMBER ff99 force field.[42−44] In this study, there are 17 nonstandard residues in total, including 14 modified residues, m1A, t6A, Ar(p), Cm, m5C, m7G, m22G, Gm, m1G, m2G, yW, m5U, D, and Ψ, and there are three terminal 5′-phosphorylated residues, 5′-phosphate adenosine, 5′-phosphate guanosine, and 5′-phosphate uridine. To be consistent with the ff99 force field, partial charges for nonstandard bases were derived using the RESP approach[45] with charges calculated based on a RESP fit to an HF/6-31G*[46−49] electrostatic potential.

Atom types were assigned for all nonstandard residues. Coordinates for the 17 nonstandard residues were excised from tRNAAsp, tRNAPhe, and tRNAiMet, and saved as PDB files. These coordinates were then used by the ANTECHAMBER and ATOMTYPE tools to assign the atoms for the residues with either AMBER or GAFF format atom types.[50−52] Missing parameters (see Supporting Information Table 1) were manually assigned based upon analogy to the GAFF force field. Additionally, parameters were assigned based on unpublished tRNAPhe parameters archived at the online AMBER parameter database.[53] The parameter library data for nonstandard residues is provided as Supporting Information. A table of parameter assignments and their sources is provided as Table 1 in the Supporting Information.

MD Simulations

tRNA molecules were solvated with a 10 Å isometric box of TIP3P water,[54] such that no solute atom was less than 10 Å from any box edge. Na+ ions were then added to neutralize the system and a 0.1 M solution of NaCl was created by adding an adequate number of Na+ and Cl– ions based on the box volume.

The resulting solvated system was energy-minimized in two steps. First, the tRNA molecule was fixed using a positional restraint on each of the tRNA atoms to its starting structure, and the position of the water and ions was allowed to change using 500 steps of steepest descent minimization, followed by 500 steps of conjugate gradient minimization. The whole system was then energy-minimized for 2500 steps, including 1500 steps of steepest descent minimization and 1000 steps of conjugate gradient minimization.

MD trajectories were calculated using the PMEMD program from AMBER 10 or 11.[43,44] The SHAKE[55] algorithm was applied to all bonds involving hydrogen atoms. The time step was 2 fs. Particle mesh Ewald[56] was used to calculate electrostatics with a direct space cutoff of 10 Å. Constant pressure (1 atm) and temperature (300 K) were maintained during the simulations using the Berendsen barostat[57] and Langevin thermostat,[58] respectively.

For each molecule, three independent simulations were run. Each trajectory of the same molecule started from the same initial configuration, but was run with different random number seeds for the initial velocity and for the Langevin thermostat.[58] Each system was initially heated to 300 K over a 20 ps NVT simulation, followed by 1 ns of NPT equilibration prior to production runs. Each simulation was then run at 300 K for 333 ns.[59] Snapshots were written to disk every 2500 steps (5 ps).

Data Analysis

Trajectories were analyzed using the PTRAJ and CPPTRAJ module of AMBER 11,[44,60] tools from LOOS package, and custom programs built with the LOOS library.[61] Trajectories and structures were visualized using VMD[62] and PYMOL.[63]

Ribose puckering pseudorotation angles (P) and amplitudes (A) were measured following the definitions of Altona and Sundaralingam,[64] where C3′-endo references as roughly P = 0° and C2′-endo references as roughly P = 180°.

Root mean squared deviations (RMSDs) were calculated using PTRAJ. Depending on the analysis, these were calculated either with all-atoms or for backbone heavy atoms (P, O5′, C5′, C4′, C3′, and O3′). These were mass weighted. When a subset of atoms were selected, such as in the regional RMSD calculations, the structure alignment was across only the selected atoms.

The analysis of the residence time of water molecules within a coordination shell around a selected set of atoms was conducted in two steps by adapting the method of Impey et al.[65] First, each water was represented by its oxygen atom. The radius of the sphere for each solute heavy atom was set to be 3.5 Å, and a water oxygen center inside this radius was considered to be in contact. For each water, the total number of frames making contact with a heavy atom and the number of times the water came into contact with a heavy atom were determined. Then the residence time was computed by dividing the total contact time (the product of contacting frames and 5 ps time between frames) with the total number of times there were consecutive contacts.

A tool to analyze base stacking was developed with the LOOS library based on the criteria described by previous studies that the distance between the center of mass (COM) of two bases was ≤4 Å and that the angle between the two normal vectors of the bases was ≤20°.[38,66]

Principal component analysis (PCA) was computed with the big-svd program of LOOS by selecting backbone heavy atoms: P, O5′, C5′, C4′, C3′, and O3′.[67] The averaged structure of each simulation was computed by the averager program of LOOS. Then the averaged structure was used as reference to generate the correlated motion trajectory with the enmovie program of LOOS by individually selecting the first three principal components, with displacement scaled by a factor of 20. The resulting movies for each of the first three principal components for each molecule were generated by VMD, and the movies are available as Supporting Information in MPEG-1 format. The all-atom PCA was also computed (excluding the modification atoms), and the covariance of overlap was calculated by the coverlap program of LOOS to assess the similarity between trajectories.[68]

For each sequence, the three independent trajectories allowed a comparison of means for quantities between molecules with a test of significance. Error estimates, when provided, are standard errors using the three trajectories as separate measurements. P values were calculated with a two-tailed t test using the R package.[69] Alpha, the type I error rate, was set at 0.05 for assessing significance.

Results

Overview of Trajectories

To gauge the stability of the simulations and to look for differences in the magnitude of dynamics, the root mean squared deviations (RMSD) as compared to the crystal structure were calculated as a function of time for all heavy atoms in the backbone (Figure 1). The RMSDs support the fact that each molecule retained an ordered structure close to its starting structure. Modified tRNAPhe had a significantly lower mean RMSD compared to unmodified tRNAPhe (P < 0.05). On average, the RMSD of modified tRNAPhe was 3.67 ± 0.041 Å, and unmodified tRNAPhe was 4.32 ± 0.095 Å. No significant differences were detected by comparing the average RMSDs of modified and unmodified tRNAAsp or tRNAiMet. The average RMSDs of modified tRNAAsp and unmodified tRNAAsp were 5.45 ± 0.35 and 5.53 ± 0.17 Å, respectively. The average RMSDs of modified tRNAiMet and unmodified tRNAiMet were 4.61 ± 0.20 and 4.64 ± 0.25 Å, respectively.

Figure 1

Backbone heavy atom RMSD as a function of time for tRNA molecules. The atoms included are P, O5′, C5′, C4′, C3′, and O3′. The left column shows unmodified tRNAs, and the right column shows corresponding modified tRNAs. From top to bottom, tRNAAsp, tRNAPhe, and tRNAiMet plots are shown, and each panel has three independent trajectories plotted in red, yellow, and blue. The average RMSDs of modified tRNAAsp and unmodified tRNAAsp are 5.45 ± 0.35 and 5.53 ± 0.17 Å. The average RMSDs of modified tRNAPhe and unmodified tRNAPhe are 3.67 ± 0.041 and 4.32 ± 0.095 Å. The average RMSDs of modified tRNAiMet and unmodified tRNAiMet are 4.61 ± 0.20 and 4.64 ± 0.25 Å.

RMSD Analysis of Anticodon Stem Loop Domain

Dynamical features of the folding domains may not be apparent in the global RMSD and thus regional RMSDs were calculated. The RMSDs of the anticodon stem loop (ASL) domain were measured by selecting the heavy atoms of the backbone from residue 27 to residue 43, and the differences between modified and unmodified tRNAs were observed (Figure 2). tRNAAsp and tRNAPhe with modifications had significantly lower average RMSDs than corresponding tRNAs without modifications (P < 0.05). The modified tRNAAsp had a mean RMSD of 2.44 ± 0. 030 Å, and unmodified tRNAAsp had a mean of 4.00 ± 0.32 Å. For tRNAPhe, the mean RMSD for the modified molecule was 2.07 ± 0.21 Å, and mean for the unmodified molecule was 3.84 ± 0.12 Å. tRNAiMet, in contrast, was apparently more dynamic for the modified tRNA, having a mean RMSD of 3.70 ± 0.75 Å, while unmodified tRNAiMet had a mean RMSD of 2.90 ± 0.23 Å. This result, however, might be due to crystal packing artifacts in the crystal structure of tRNAiMet, where the base at position 37 is not stacking with the bases in the anticodon as observed in other tRNA structures. This is discussed below in greater detail.

Figure 2

Backbone heavy atom RMSD as a function of time for tRNA anticodon stem loops (nucleotides 27–43). The left column shows unmodified tRNAs, and the right column shows corresponding modified tRNAs. From top to bottom, tRNAAsp, tRNAPhe, and tRNAiMet plots are provided, and each panel has three independent trajectories plotted in red, yellow, and blue, where the trajectory colors correspond to those in Figure 1. The average RMSDs of modified tRNAAsp and unmodified tRNAAsp are 2.44 ± 0. 030 and 4.00 ± 0.32 Å. The average RMSDs of modified tRNAPhe and unmodified tRNAPhe are 2.07 ± 0.21 and 3.84 ± 0.12 Å. The average RMSDs of modified tRNAiMet and unmodified tRNAiMet are 3.70 ± 0.75 and 2.90 ± 0.23 Å.

To further examine how modification affects tRNA dynamics regarding function in translation, the three anticodon-residues were studied with all-atom RMSD calculations (Figure 3). To have an equivalent and unbiased atom selection for the comparison between modified and unmodified tRNAs, the methyl group at position 34 of modified tRNAPhe was not included in the calculation. Interestingly, in contrast to the ASL domain RMSDs, comparison between anticodon all-atom RMSDs of modified and unmodified tRNAs showed that modifications tended to introduce flexibility to the anticodon, although the differences were not statistically significant for simulations of this length. On average, the RMSD of modified tRNAAsp was 3.81 ± 0.68 Å, whereas unmodified tRNAAsp was 3.09 ± 0.52 Å. For tRNAPhe, the mean RMSD for modified is 2.77 ± 0.24 Å, and it was slightly greater than the mean of unmodified molecule, which was 2.57 ± 0.82 Å. This is due to the fact that one trajectory (trajectory 3, blue in Figures 1–3) of the modified tRNAPhe stayed close to the starting configuration during the simulation, and the mean RMSD for that trajectory is 1.00 Å. The RMSDs of the modified tRNAiMet and the unmodified tRNAiMet were 3.40 ± 0.59 and 2.67 ± 0.20 Å, respectively.

Figure 3

All-atom RMSD as a function of time for tRNA anticodon loops. The left column shows unmodified tRNAs, and the right column shows corresponding modified tRNAs. From top to bottom, tRNAAsp, tRNAPhe, and tRNAiMet plots are listed, and each panel has three independent trajectories plotted in red, yellow and blue, where the trajectory colors correspond to those in Figure 1. The average RMSDs of modified tRNAAsp and unmodified tRNAAsp are 3.81 ± 0.68 and 3.09 ± 0.52 Å. The average RMSDs of modified tRNAPhe and unmodified tRNAPhe are 2.77 ± 0.24 and 2.57 ± 0.82 Å. The average RMSDs of modified tRNAiMet and unmodified tRNAiMet are 3.40 ± 0.59 and 2.67 ± 0.20 Å.

Sugar Pucker of Pseudouridine

To quantitatively assess how pseudouridine (Ψ) selectively populates the 3′-endo sugar pucker, which was observed in previous experimental studies,[32] the sugar conformation was measured for all the Ψ residues of modified tRNAs by calculating the pseudorotation angle for the ribose.[64] Ψ13, Ψ32, and Ψ55 of tRNAAsp and Ψ39 and Ψ55 of tRNAPhe were compared with corresponding uridines of the unmodified tRNAs. Sugar pucker pseudorotation angles as a function of time for tRNAAsp and tRNAPhe Ψs and corresponding uridines are presented in Supporting Information Figure S1. The plot shows that all the Ψ residues maintained 3′-endo conformations with only rare frames of transitions to 2′-endo for Ψ13 of tRNAAsp. The sugar in U13 and U32 of unmodified tRNAAsp and U39 of unmodified tRNAPhe, however, sampled both the C2′-endo and C3′-endo pucker conformations during the simulations, as shown in Supporting InformationFigure S1. One trajectory of unmodified tRNAAsp stayed in the 2′-endo conformation for the full length of the trajectory. For uridines at position 55 of tRNAAsp and tRNAPhe, although there were no notable transitions from 3′-endo to 2′-endo pucker observed, the plot shows that Ψ residues at this position had a more rigid sugar pucker conformation compared to uridine. These findings support that Ψ favors the 3′-endo conformation in RNA.

Water Residence Time around Ψ

Pseudouridylation of uridine replaces the N1 atom from the N1–C1′ glycosidic bond with the C5 atom, which provides an additional imino nitrogen for hydrogen bonding. Experimental and MD simulation studies showed that the (Ψ)N1–H increases the chance to form hydrogen bonds to water.[28−33,38] To examine the hydrophilicity shift, water residence times around Ψ and the corresponding uridine were calculated. The investigated residues were Ψ13, Ψ32 and Ψ55 of tRNAAsp and Ψ39 and Ψ55 of tRNAPhe. The averaged distributions of water residence time for Ψ or U base were computed from the three trajectories of each molecule. Comparing the distribution of water residence time between Ψ and U bases indicated that a higher population of water molecules tends to be present at the short residence time range (shorter than 10 ps) for U (Supporting Information Figure S2). For Ψ, however, the frequencies at the long residence time range (longer than 30 ps) are higher (Figure 4). Moreover, in the long time range, the difference of the residence time between Ψ and U bases was more remarkable for tRNAAsp than it was for tRNAPhe.

Figure 4

Averaged distribution of U or Ψ base water residence time for the water molecules in the simulation system. The distributions are averaged from three trajectories of each molecule, and ± a single standard error is displayed as error bars. To focus on longer residence time range, the time axes in the plots are truncated from 20 to 600 ps. Plots for the full time range (5 to 600 ps) are shown in Supporting Information Figure S2. Panel A: From top to bottom, plots for bases at position 13, 32, and 55 of tRNAAsp are listed, with unmodified tRNAAsp plotted in blue and modified tRNAAsp plotted in yellow. Panel B: From top to bottom, plots for bases at position 39 and 55 tRNAPhe are listed; unmodified tRNAPhe is plotted in blue, and modified tRNAPhe is plotted in yellow.

Base Stacking Analysis in the Anticodon Loop

The conformation of anticodon bases in solution is crucial for tRNA interaction with mRNA. The bases in the anticodon loop, positions 34–37, are stacked in the crystal structures of all three tRNAs studied here, with the notable exception that tRNAiMet has t6A37 not stacked. In addition, the bases at position 37 of all three tRNAs are modified, and the modifications at this position facilitate codon-anticodon interaction and maintain translational fidelity.[17] Therefore, the fraction of time that 34 stacks on 35, 35 stacks on 36, and 36 stacks on 37 were determined for all the trajectories (Figure 5) to elucidate the dynamics of stacking. During the simulation, for the stacking of 34 and 35, modified tRNAAsp (0.033 ± 0.014) had a significantly (P < 0.05) lower fraction of time stacking than unmodified tRNAAsp (0.105 ± 0.018). Although no significant difference was detected except the stacking of 34 and 35 in tRNAAsp, in general, modified tRNAs showed lower fraction of time stacking when compared to corresponding unmodified tRNAs. Only one exception was found for tRNAAsp at 36 and 37 stacking, modified tRNAAsp (0.408 ± 0.018) showed a higher fraction of time stacking than unmodified tRNAAsp (0.273 ± 0.211). These findings correspond to what was observed with anticodon RMSDs, i.e. anticodon residues had more flexibility in conformation in modified tRNAs than unmodified tRNAs

Figure 5

Fraction of stacking conformation during simulations for anticodon loop bases. Three stacks are shown for each molecule. From left to right are unmodified tRNAAsp, modified tRNAAsp, unmodified tRNAPhe, modified tRNAPhe, unmodified tRNAiMet, and modified tRNAiMet. In each panel, base 34 to 35 (blue), base 35 to 36 (orange), and base 36 to 37 (yellow) stacking conformation fractions are displayed.

Principal Component Analysis

To study correlated motions of the tRNA structure, principal component analysis (PCA) was performed.[70,71] In this study, the heavy atoms of the backbone were used for the analysis and all three trajectories for each molecule were analyzed together to facilitate the comparison of modified and unmodified tRNA. The first principal component represents the largest correlated atomic fluctuation in the sampled conformations. The second principal component represents the next largest correlated atomic fluctuation orthogonal to the axis of the previous mode, and so on. For all the molecules in this study, the amplitudes of the first three principal components were comparably close. Therefore, in order to summarize the major motions, the first three principal components were examined. Furthermore, to avoid losing the correlations within each principal component, the first three principal components were extracted and visualized individually. The motions observed in animations of the first three principal components essentially involve the movements of the two L-shaped arms relative to each other. The L-shaped tertiary structure provides the flexibility that allows the two arms to move relative to each other in three directions, which are bending the arms as hinge-bending motion in the plane of the arms, bending the arms out of the plane of the arms and compressing along the axis of the helix as a spring-like motion.

No striking difference was found from the correlated motions between unmodified and modified tRNAs. In the first principal component, the type of hinge-bending motion was predominant among all 18 trajectories; with three trajectories as exceptions: one trajectory of unmodified tRNAAsp (trajectory 2, yellow in Figures 1–3) displayed a combination of motions, with bending of the anticodon arm out of the plane and an acceptor arm spring-like motion; one trajectory of modified tRNAAsp (trajectory 1, red in Figures 1–3) displayed a motion bending the arms out of the plane; and one modified tRNAPhe (trajectory 1, red in Figures 1–3) had a linear spring-like motion for the acceptor arm and a bending planar motion for the anticodon arm. The predominant in-plane motion between the two arms of the tRNA was previously noted in normal-mode analysis.[72] In the second components, the hinge-bending motion decreases and loses its predominance as in the first principal component, and the molecules showed a combination of all three types of motion for the two arms in the L. Interestingly, in the third principal component, the occurrence of the spring-like motion increased noticeably as compared to the second principal component.

To assess the similarity in correlated motions between trajectories, the covariance overlap of the principal components was determined.[73] When the two trajectories are for identical molecules, the covariance overlap tests the convergence of sampling. The covariance overlap is one only if the sampled subspaces are identical. It is zero when the sampled subspaces are completely orthogonal. When the two trajectories are of different molecules, the covariance overlap measures the similarity in correlated motions. Another set of PCA analysis was conducted to measure the covariance overlap between trajectories, where all atoms except those on modifications were selected for the calculation. All-atom selection gives the best coverage of the molecule’s dynamics, and the exclusion of modification atoms allows the covariance overlap to be calculated between modified and unmodified molecules. The three trajectories of each molecule were used to calculate the averaged covariance overlap by taking the mean of covariance overlaps over three pairs of trajectories. Additionally, the averaged covariance overlap between modified and unmodified tRNAs was measured from the nine pairs of modified-unmodified tRNA trajectories. The results showed that the averaged covariance overlaps of trajectories for modified and unmodified tRNA with itself were significantly higher than between unmodified and modified tRNA. This suggests that the correlated motions of unmodified tRNAs are different than those of the corresponding modified tRNAs (Figure 6).

Figure 6

Averaged covariance overlap between PCA of backbone heavy atoms. From left to right, three sets of bars represent tRNAAsp, tRNAPhe, and tRNAiMet, respectively. The blue bar is the averaged covariance between trajectories of unmodified molecule, the orange bar is the averaged covariance between trajectories of modified molecules, and the yellow bar is the averaged covariance between trajectories of unmodified molecule and modified molecule. Mean values are plotted, and a single standard deviation is shown as error bars.

Discussion and Conclusions

Solvated yeast tRNAAsp, tRNAPhe, and tRNAiMet in the presence and absence of modifications were studied by MD simulations with an aggregate sampling of one microsecond across three simulations. A key question addressed in this study was whether modifications of tRNA resulted in a significant alteration in the structure and dynamics of the molecule. Both functional and structural studies have demonstrated that modified residues are involved in changing tRNA structure and result in tuning the tRNA activity.[4−7,11−16,26]

The MD simulations of all three tRNAs, both modified and unmodified, were stable. Previous experimental studies showed that modified residues can increase the rigidity and thermal stability of the structure of tRNA or have exactly the opposite effect.[18,26] On the basis of the global RMSD calculations for these simulations, however, no significant rigidity was introduced by modifications. There are two possibilities that could contribute to this apparent inconsistency. First, a total of 1 μs simulation time may not have been long enough to fully sample tRNA conformations. For example, Förster resonance energy transfer studies on tRNA fluctuation during translocation suggest that the lifetime of tRNA fluctuation in solution is on the order of a hundred milliseconds.[74−76] Second, different modifications may have different effects on tRNA structural stability. Ashraf et al. observed from a UV-monitored thermal denaturation studies that the stem modification, m5C40, improved RNA thermal stability, but had a deleterious effect on ribosomal binding. In contrast, the loop modification, m1G37, enhanced ribosome binding, but dramatically decreased thermal stability.[77] In this study, modifications were not studied in subsets of combinations, and hence opposite stability effects of modifications could cancel each other.

In the tRNA crystal structures, tRNAAsp and tRNAPhe both have Mg2+ ions located at Mg2+ binding sites. In this study, however, Mg2+ ions were not added to the simulations because of the missing location information on tRNAiMet and the fact that parameters of Mg2+ were reported to be inaccurate.[41] The trajectories of tRNAAsp and tRNAPhe showed that Na+ ions consistently reside in the groove of tRNA molecules, including the Mg2+ binding sites, and this suggests that a Na+ ion may act as a Mg2+ surrogate. As in the cases of tRNAAsp and tRNAPhe, stable dynamics of the whole tRNA in the absence of a Mg2+ ion were observed.[36,41] Evidently, the tertiary structure of tRNA is stabilized by Mg2+ accumulating in the locations of highly negative electrostatic potentials, and previous experimental studies have revealed that the presence of Mg2+ alters the local structure and mobilities.[78,79] There is evidence that tRNA core modifications improve the binding of Mg2+ ions.[80,81] It would be interesting to investigate this question in subsequent MD simulation studies using, for example, recently reported Mg2+ parameters.[82]

Although the mean global RMSDs do not show differences that suggest modifications alter tRNA structural rigidity, the covariance overlap calculation shows that sampled subspaces overlap within either modified tRNA or unmodified tRNA molecules to a significantly greater extent than between modified and unmodified tRNA. This suggests that there are differences in correlated motions in the presence or absence of modifications.

With the PCA analysis, the lowest principal components revealed how the two arms in the L-shaped tertiary structure moved. The arms can either bend in directions in the plane of the arms, or bend out of the plane of the arms, or, in addition, the arms can compress as springs. The hinge-bending motion was predominantly found among all tRNAs in the first mode, which may indicate that hinge-bending is the efficient way for tRNA to “walk” from one site to another. On the other hand, with a similar amplitude, the compressing spring motion of the anticodon arm allows the head of tRNA to rotate, which was revealed in previous studies of tRNA translocation.[83]

In this study, statistical significance was estimated for the observed differences when comparing molecules. The set of three simulations for each molecule provided a level of statistical power to assess significance. Other, subtler, differences in the dynamics between molecules cannot be shown to be significant with only three simulations of this duration. Additional simulations could be run in the future to provide additional power to differentiate the dynamics between molecules on the microsecond time scale.

The properties of Ψ were assessed for the tRNAAsp and tRNAPhe based on experimental studies that showed Ψ enhances the 3′-endo sugar pucker in RNA oligonucleotides and the released N1 atom might increase the capability of hydrogen bond formation. According to the sampled sugar pucker conformations, the simulations supported that pseudouridylation does cause Ψ to favor the 3′-endo conformation. Additionally, the water residence time (Figure 4) also demonstrated that Ψ tends to interact with water more tightly than uridine.

The comparison of ASL regional RMSDs between modified and unmodified tRNAs indicates that modifications introduce rigidity into the loop region, with the exception of tRNAiMet, which is more dynamic with modifications. The inconsistency of tRNAiMet ASL behavior in the simulation might be due to artifacts in its initial structure. The crystal structure of tRNAiMet was not defined clearly at the anticodon loop because of the low resolution.[21] Additionally, the anticodon loop was involved in lattice contacts, where the base of residue 37 stacked with the same base from the neighboring molecules, which resulted in a conformation different from that found in other tRNA structures.[21] Conformational flexibility of the anticodon was observed throughout the trajectories, suggesting that the residue 37 might be exploring the conformational space to find a stacking conformation, which was not found in 333 ns.

The ASL domain is a crucial part of tRNA because it recognizes the translating codon within the ribosome and the modified residues are believed to be involved in regulating the fidelity of translation. Conformational flexibility and rigidity are important for the anticodon arm to interact with mRNA during translation. The high frequency of modification in the ASL could be important for tuning this balance between flexibility and rigidity. The ASL domain RMSD showed that modifications increased the rigidity of the backbone of the ASL domain (Figure 2), which was mostly an influence on the loop region. It was observed (by RMSD) that the stem region for both modified and unmodified tRNAs are stiff and difficult to differentiate (data not shown). The all-atom RMSD calculation of the anticodon, however, showed that modifications tended to increase the dynamics of anticodon residues (Figure 3). Taken together, these findings suggest an arm and hand model for the ASL. The ASL domain of tRNA works as an arm. The stem is stiff for translocation during translation, and the loop domain is a palm that is structured, but can move in relationship to the helix, as noted in previous coarse-grained studies,[72,84] to be able to position the anticodon. The anticodon residues then work as fingers to recognize the right codon. Flexibility is the key feature for the anticodon, and modification strengthens this property by introducing dynamics.

Although it was observed that some modified residues showed distinct properties as compared to unmodified residues, in this study only local tRNA features were found significantly altered. One explanation could be that the tRNA molecule works as an adaptor, interacting with other molecular devices, such as aminoacyl-tRNA synthetases, ribosomes, and mRNAs. The modified nucleosides may alter the interface between tRNA and its interactive components, such as the interactions studied for tRNA interacting with the ribosome.[85] Changes in dynamics caused by interactions of modifications with other components would be difficult to capture in the scope of this study. Additional studies that simulated modified and unmodified tRNA interacting with other components of translation machinery would be interesting and might elucidate further roles for covalent modifications.