Binding interactions between the core central domain of 16S rRNA and the ribosomal protein S15 determined by molecular dynamics simulations

Abstract

The goal of the current study is to utilize molecular dynamic (MD) simulations to investigate the dynamic behavior of 16S rRNA in the presence and absence of S15 and to identify the binding interactions between these two molecules. The simulations show that: (i) 16S rRNA remains in a highly folded structure when it is bound to S15; (ii) in the absence of S15, 16S rRNA significantly alters its conformation and transiently forms conformations that are similar to the bound structure that make it available for binding with S15; (iii) the unbound rRNA spends the majority of its time in extended conformations. The formation of the extended conformations is a result of the molecule reaching a lower electrostatic energy and the formation of the highly folded, crystal‐like conformation is a result of achieving a lower solvation energy. In addition, our MD simulations show that 16S rRNA and S15 bind across the major groove of helix 22 (H22) via electrostatic interactions. The negatively charged phosphate groups of G658, U740, G741 and G742 bind to the positively charged S15 residues Lys7, Arg34 and Arg37. The current study provides a dynamic view of the binding of 16S rRNA with S15.

Received October 1, 2002; Revised and Accepted November 14, 2002

INTRODUCTION

Investigation of RNA–protein binding mechanisms is of great interest ( 1 – 4 ). Ribosomal subunits contain abundant rRNA–protein interactions ( 5 – 10 ). One example is the assembly of the central domain of the 30S ribosomal subunit, which is initiated by the binding of the core central domain of 16S rRNA with the ribosomal protein S15 ( 9 , 11 – 13 ). The binding sites of 16S rRNA with S15 are located at two three‐helix junctions (3HJ): the lower 3HJ composed of helices 20, 21 (H20 and H21) and the lower portion of helix 22 (H22) and the upper 3HJ composed of the upper portion of H22, H23a and H23b ( 9 , 14 ).

Gel mobility and transient electric birefringence studies have demonstrated that the core central domain of 16S rRNA is highly folded when it is bound to S15 ( 12 , 15 ). In the folded conformation, H22 and H21 stack coaxially and H20 forms an acute angle with H22 ( 15 ). The binding with S15 induces a conformation of the upper 3HJ of 16S rRNA that binds with ribosomal proteins S6 and S18 ( 9 , 11 – 13 ). In the absence of S15, extended conformations of 16S rRNA are observed in which the inter‐helical angle between H20 and H22 increases ( 11 – 13 ). A ‘structure‐capture’ mechanism for 16S rRNA–S15 binding has been proposed whereby the rRNA alternates between the extended and coaxially stacked conformations and the rRNA is recognized by S15 when it transiently adopts the coaxially stacked conformation ( 15 ).

Site‐directed mutagenesis studies have demonstrated that the core central domain of 16S rRNA containing 104 nt is essential for the binding of S15, S6 and S18 ( 11 – 13 , 15 – 17 ). The crystal structure containing 84 nt of 16S rRNA and these three ribosomal proteins shows that the nucleotides that bind to S15 are located in three regions of the rRNA (Fig. 1 ) ( 9 ). In the lower 3HJ, 6 nt (C656, G657, G750, U751, G752 and C754) contact Gln27, Thr21, Gly22, Thr24 and Tyr68 of S15, respectively. These binding residues of S15 are located in or near the loop between the first (α ₁ ) and second (α ₂ ) α‐helices. In the upper portion of H22, 5 nt (G666, G667, G668, C739 and U740) contact residues Asp48, His45 and Ser51. These residues are located in the loop between α ₂ and α ₃ of S15. In addition, 3 nt (A728, A729 and G730) in the tetraloop terminating H23a contact Arg53 of S15 ( 9 ).

Molecular dynamics (MD) simulations are utilized in the current study to investigate the dynamic behavior of the core central domain of 16S rRNA in the presence and absence of S15. Both the Particle‐Mesh‐Ewald (PME) method and the Generalized Born (GB/SA) model were used for computing electrostatic solvation energies ( 18 – 20 ). Using the GB/SA model speeds up searches for molecular conformations, particularly for a large molecular system ( 21 – 28 ). The GB model has been successfully employed in a number of MD simulations for nucleic acids ( 24 – 26 , 29 , 30 ).

Our simulations show that the 16S rRNA–S15 complex assumes a conformation similar to the crystal structure during an equilibration period. In the absence of S15, 16S rRNA transiently forms conformations similar to the bound structure that may make it available for binding with S15. The unbound rRNA spends the majority of its time in extended conformations. The changing conformations of the unbound 16S rRNA supports the ‘structure‐capture’ mechanism for 16S rRNA–S15 binding ( 31 ), and these results provide detailed views of the dynamic features of this rRNA in the extended and folded conformations. Significantly, our MD simulations identify electrostatic binding interactions between G658 and Lys7 of S15 on the lower side of the major groove of H22 and between G742 and Arg34, G741 and Arg34, G741 and Arg37, and U740 and Arg37 on the upper side of the major groove. These interactions play a key role in retaining H21 and H22 of 16S RNA in a coaxially stacked conformation, thereby preventing the large fluctuations of the rRNA.

MATERIALS AND METHODS

The starting coordinates of the rRNA and S15 were taken from the crystal structure with a resolution of 2.6 Å (Protein Data Bank entry code 1G1X) ( 9 ). In this rRNA fragment, the native sequences of H20 and H21 were truncated and H21 was capped by adding a GAAA tetraloop. The sequence numbering followed that of Escherichia coli ( 13 , 16 ). Hydrogen atoms were added to the crystal structure using the AMBER6/xleap editing program. The simulations were performed using the Cornell force field ( 32 ) and the Sander module in the molecular simulation software package AMBER6 ( 19 , 20 ). AMBER6/Carnal was used to generate the RMSDs of the trajectories.

The electrostatic interaction energies were obtained using both the PME method and the GB/SA model. The simulations using the PME method and the GB/SA model were performed on the unbound 16S rRNA and a simulation using the GB/SA model was performed on the 16S rRNA–S15 complex.

Simulations using the GB/SA model

All parameters for the GB/SA calculations were taken from those examined by Tsui and Case ( 29 ) using a salt concentration of 1.0 M. In addition, the non‐electrostatic contribution to the solvation free energy was computed using γ· SASA , where γ = 0.005 kcal/mol/Å ² and the solvent‐accessible surface area SASA was calculated using the algorithm of Weiser et al . ( 33 ) that was implemented in AMBER6. The simulations were performed on SGI Origin 2100 and 2000 computers. The production runs followed the same simulation protocol and lasted for 5 ns for both the unbound 84 nt rRNA and the rRNA–S15 complex.

The starting structure was first subjected to 200 steps of steepest descent minimization. This was followed by 50 ps of MD to heat the structure from 0 to 300 K with a 5 kcal/mol/Å ² harmonic constraint on the system relative to the starting structure. At 300 K, the heated structure was then subjected to three 50 ps constrained MD equilibrations with continuously reduced constraints from 2 kcal/mol/Å ² , then to 0.5 and finally to 0.1 kcal/mol/Å ² . Finally, the production run for equilibration was performed with constant temperature.

The time step for all simulations using the GB/SA model was 1 fs. The heat bath coupling constant used was 0.1 ps for heating and 2.0 ps for equilibrations. The cut‐off distance for non‐bonded interactions was 15 Å.

Simulation using the PME method

Sodium ions were added to neutralize the charges of the rRNA. In addition, sodium and chloride ions (40 for each) were added to make an ionic concentration of ∼0.1 M. The rRNA and ions were placed in water boxes with water shells of at least 10 Å using the TIP3P water model.

The systems were minimized for 200 steps using the steepest descent minimization algorithm, followed by 40 ps of dynamics to increase the temperature from 100 to 300 K with a 20 kcal/mol/Å ² constraint on the rRNA. This was then followed by 40 ps of equilibration with the same constraint but only on the solute. This was followed by 200 steps of minimization. After these steps, the water box size was reduced. The new water box contained 12 105 water molecules and the water box sizes varied from initially 84 × 66 × 94 to finally 77 × 60 × 84 Å ³ after the system was initially equilibrated. Using the new water box size, the previous minimization and heating dynamics were repeated. At 300 K, the system was equilibrated with gradually decreased constraints on the rRNA starting from 2 kcal/mol/Å ² , then 0.5 and finally 0.1 kcal/mol Å ² , with 40 ps for each of the equilibrations. The production runs had released constraints on the rRNA and were performed under constant pressure for 5 ns. SHAKE was applied to the bonds involving hydrogen atoms ( 34 ). The time step was 2 fs. The cut‐off distance for non‐bonded interactions was 9 Å.

Calculation of RMSD from B‐factors

The B‐factors were converted to RMSD (RMSD _cry ) by the relationship ( 35 – 37 ):

RMSD _cry = (3 B /8π ² ) ^½ .

The RMSDs of the simulation were derived from 30 snapshots of the equilibration of the rRNA–S15 complex taken over the last 3 ns of the equilibration at 100 ps time intervals. The RMSDs for the phosphate atoms in the backbone of the rRNA and the C4′ atoms that link the phosphate groups and the sugars were chosen to be compared with the RMSD _cry values for each nucleotide. The RMSDs of all phosphate atoms and C4′ atoms of the rRNA from the 30 snapshots relative to the average structure of the equilibration were calculated. Then the RMSDs were averaged over the snapshots, and finally were averaged over the phosphate and C4′ atoms for each nucleotide. These calculations produced the average RMSD of each nucleotide. The RMSD _cry values of the phosphate atom and the C4′ atom were averaged for each nucleotide.

The average RMSD of the C _α atom of each residue of S15 from the 30 snapshots was also calculated. The RMSD _cry values of the C _α atoms were converted using the B‐factors of the corresponding residues.

RESULTS

Binding with S15 stabilizes conformation of 16S rRNA

In order to investigate the dynamic features of the core central domain of 16S rRNA in the presence of S15 and the structural stability of the rRNA–S15 complex, we performed an equilibration simulation on the rRNA–S15 complex. The conformation of the complex during equilibration is similar to its crystal structure. An average structure of the rRNA–S15 complex was derived from the last 3 ns of the equilibration (Fig. 2 ). The average RMSD of all the 172 residues of the complex relative to the mean structure of this simulation is 2.4 Å. The rRNA remains in the coaxially stacked conformation throughout the equilibration and the relative orientation of H20 and H22 remains similar to that found in the crystal structure. The equilibration does not significantly change the overall structure of the complex.

In order to evaluate the reliability of the simulation, we compared the atomic thermal motions of the rRNA and S15 as determined from the simulations (RMSDs) and from the crystallographic B‐factors (RMSD _cry ). Overall, the average RMSDs and RMSD _cry values of 16S rRNA show close concordance (Fig. 3 A). The average RMSD and the RMSD _cry of each residue of S15 show an overall good correlation (Fig. 3 B) except that residues 12–17 in the first α‐helix of S15 show higher mobility in the simulation than in the crystal. The high mobility of these residues may be due to the release of the complex from the packed crystal in which the α ₁ helix is in close proximity to the S6 protein in another monomer which includes the 16S rRNA, S15, S6 and S18.

The overall agreement between the RMSDs derived from the simulation and the RMSD _cry values derived from crystallography indicates that the current simulation reasonably reproduces the atomic thermal motions as described by crystallography.

Conformation of 16S rRNA fluctuates in the absence of S15

In order to examine the dynamic features of the core central domain of 16S rRNA when S15 is removed from the crystal complex, equilibration simulations were performed on the rRNA after the release of S15. The GB/SA simulations show that the conformation of the unbound rRNA fluctuates significantly during the equilibration rather than remaining in the crystal structure (Fig. 4 ). When the RMSD of the unbound rRNA is small, it is close to the average RMSD of the rRNA derived from the equilibration performed on the 16S rRNA–S15 complex. When S15 is not present, the 16S rRNA only transiently forms a conformation with a small RMSD relative to its crystal structure.

Similarly, significant conformational changes of the 16S rRNA are also found in the simulation when using explicit solvent with the PME method. The RMSD fluctuates extensively with a maximum of 7 Å. The conformation corresponding to this large RMSD shows a major increase in the angle between H20 and H22, which corresponds to a loss of the coaxial stacking between H21 and H22 (Fig. 5 A). In addition, the highly folded conformation also re‐occurs. Therefore, the conformational changes of the rRNA are similar regardless of the simulation approaches. However, the simulation using the GB/SA model allows a more rapid change of the conformations over a given time period than when using the PME model. Therefore, the GB/SA trajectory may provide better conformational sampling. It should also be noticed that the simulations using salt concentrations of 0.1 and 1.0 M for the GB/SA model produced similar results. We use the trajectory generated by the GB/SA model (0.1 M salt) for the following analyses.

An average structure derived from the fluctuating conformations may not represent any real conformation, but it can be used to illustrate the dominant characteristics of the unbound 16S rRNA. The average structure was derived from the last 1 ns of the trajectory (Fig. 5 B). This structure differs from the fully folded structure of the rRNA in the crystal by an RMSD of 6.0 Å, and it shows that the most significant differences occur in H22. The structure of H22 is the most mobile region compared with other portions of the molecule during the equilibration. When the RMSD is large, H22 is in a significantly different conformation from its crystal structure. Helices 20 and 21 are relatively less mobile during the equilibration. The inter‐helical angle between H20 and H22 increases due to the conformational change of H22. Here, the angle ϕ, defined as the angle between lines AB and AC (depicted in Fig. 7 b for later discussion), is used to illustrate the relative orientation between H20 and H22. The three points designated as A, B and C were chosen at the centers of mass of three nucleotide groups as follows: A, A753, C754 and G587; B, U758, U759, A582 and A583; C, G742 and A663. Angle ϕ increases to 114° in the average MD structure from 80° in the crystal structure. The increase in angle ϕ corresponds to a loss of the coaxial stack between H22 and H21. The inter‐helical angle between H20 and H21 in the simulation remains similar to that observed in the coaxially stacked structure.

The conformational change of H22 can be localized to individual nucleotides by comparing changes in their RMSDs. The fluctuation of the RMSD of the single bulge nucleotide C748 shows the most significant change relative to the crystal structure among the nucleotides of the rRNA. The average RMSD of the nucleotides is ∼0.1 Å when measured individually, but the average RMSD of C748 is ∼0.4 Å, and this correlates with the fluctuation of the rRNA during the equilibration. The high mobility of nucleotide C748 accompanies the changes in the backbone geometry of the lower portion of H22 (Fig. 5 C). This base does not stack between base pairs G657·C749 and G658·C747 in the crystal structure. The high mobility of the unstable single bulge, nucleotide C748, appears to be primarily responsible for changes in the backbone conformation of H22.

Influenced by the relaxation of the irregular base arrangement of the lower portion of H22, the width of the major groove of H22 fluctuates significantly. In order to show these changes, two snapshots (named A and B) were chosen from the equilibration performed on the unbound rRNA (Fig. 4 ). Snapshot A has a relatively small RMSD relative to the crystal structure (Fig. 5 ). It shows a highly folded conformation with a compact H22 and a narrow major groove. Snapshot B shows H22 with an extended backbone geometry and a wider groove with a large RMSD (Fig. 5 ). The distance between G657:P and U743:P was used here to characterize the width of the major groove. This distance fluctuates during the equilibration (Fig. 6 ) and the fluctuations are mostly correlated with the fluctuations in the RMSD of the unbound rRNA (Fig. 4 ). The groove is narrow in snapshot A, with the distance between G657:P and U743:P being ∼12 Å. The groove is much wider in snapshot B, with the distance being ∼18 Å (Fig. 6 ). This conformational change causes an increase in the radius of gyration derived from the apical four nucleotides of H20 and the internal loop (A663–A665, G741 and G742). The radius of gyration is 18 Å in snapshot A and 25 Å in snapshot B.

In order to characterize the backbone geometry of the major groove, an RNA fragment with A‐form geometry was constructed using the program INSIGHTII (Fig. 5 G). In the coaxially stacked conformation, the major groove is narrower than that of the standard A‐form RNA and G657:P and U743:P are located in close proximity across the groove. In the extended conformation, the backbone geometry of the major groove is close to that of the A‐form RNA (Fig. 5 H). When the phosphates of G657 and U743 are in close proximity, a high electrostatic repulsion occurs.

Unbound 16S rRNA spends the majority of its time (>90%) in extended conformations and these conformations can be significantly different from the crystal‐like conformation. This probably explains why an extended conformation was observed experimentally in the absence of S15 ( 11 , 12 , 15 – 17 ).

The changing conformations of the unbound 16S rRNA include a structure that is similar to the bound rRNA as shown in snapshot A. The structure of H22 in snapshot A shows high similarity to the average structure of H22 in the bound rRNA (Fig. 5 F). Thus, the transiently and spontaneously formed crystal‐like conformation of the unbound rRNA may be available for binding with S15.

Causes of the conformational changes of the unbound 16S rRNA

In order to examine the potential causes for the dynamic behavior of the unbound 16S rRNA, the potential energies of the highly folded, crystal‐like conformation of snapshot A and the extended conformation of snapshot B were compared. The total potential energies of snapshots A and B are –18639 and –18668 kcal/mol, respectively, therefore the potential energy of the extended conformation is only marginally lower (–29 kcal/mol) than that of snapshot A (Table 1 ). However, the electrostatic energy of the extended conformation (–3039 kcal/mol) is significantly lower than that of the highly folded conformation (–1086 kcal/mol). Therefore, the formation of an extended conformation of the 16S rRNA results in a lower electrostatic energy. In contrast, the solvation energy of the extended conformation of snapshot B is significantly lower (–18923 kcal/mol) compared to that of the highly folded conformation of snapshot A (–20892 kcal/mol). Therefore, the formation of the highly folded conformation of snapshot A results in a lower solvation energy. The decrease in the electrostatic energy from snapshot A to B, –1953 kcal/mol, is approximately counterbalanced by the increase in the solvation energy, +1969 kcal/mol, during the conformational change. The highly folded conformation and the extended conformation are alternately formed and neither conformation is dominantly stable.

Structural reorientation of the upper 3HJ

Experiments have demonstrated that the formation of a coaxial stack at the lower 3HJ induces a conformation in the upper 3HJ that is suitable for binding to S18 and S6 ( 13 , 15 ). Only a partial upper 3HJ is simulated in the current study since the coordinates of H23b were not included in the 84 nt crystal structure (Fig. 1 ). The simulation performed on the unbound rRNA shows a change in the conformation of this portion of the upper 3HJ when the fully folded crystal structure moves to the extended conformation.

The analysis of the RMSD on the unbound rRNA indicates that the major difference between the upper 3HJ and the crystal structure is in H23a. The distance of the tetraloop from H20 increases as angle ϕ increases. The distance d from the mass centers of the nucleotides U758, U759, A582 and A583 to the nucleotides C726 and G731 increases significantly to 35.9 Å in the extended conformation from 24.1 Å in the crystal structure as angle ϕ increases from 80° in the crystal structure to 114° in the average structure of the rRNA (Fig. 7 ). These changes are primarily caused by the elimination of the interaction of this tetraloop with S15, by the relaxation of the internal loop G718–G721 and by the change in the backbone geometry of the lower portion of H22.

The conformation of the upper portion of H22 (above the internal loop) remains similar to the highly folded structure during the equilibration, with an average RMSD of 1.6 Å. However, a significant change occurs in the large internal loop of H23a (G718–A733) with an average RMSD of 4.4 Å. This internal loop may be further stabilized by the binding of S6 and S18 to the upper 3HJ. Thus the binding of the upper 3HJ to S6 and S18 may not take place until a suitable conformation of the upper 3HJ is induced by the binding of the rRNA with S15. This may explain why the binding of 16S rRNA with S15 acts as an initiator for the integral assembly of the central domain of the 30S subunit ( 9 , 11 – 13 ).

Binding sites of 16S rRNA with S15

Based on the crystal structure ( 9 ), 16S rRNA interacts with S15 in three regions of the rRNA: (i) the lower portion of H22 involving the lower 3HJ (referred to as region 1); (ii) the upper portion of H22 involving the upper 3HJ (referred to as region 2); (iii) the tetraloop terminating H23a (referred to as region 3) (see Fig. 1 ). The dynamic behavior of the intermolecular hydrogen bonds obtained from the equilibration are summarized in Table 2 . The atomic labeling of S15 is assigned using the program INSIGHTII.

16S rRNA–S15 binding in region 1. The binding nucleotides in region 1 are located in the two strands of the lower portion of H22 (Fig. 8 ) ( 9 ). In the crystal, residue Asp20 in the loop between α ₁ and α ₂ of S15 contacts G750, but this interaction becomes weak during the equilibration (hydrogen bond length Asp20:OD1–G750:HO′2 is >3 Å). The contacts of Gly22 to G750, Thr24 to U751, and Gln27 to both C656 and G750 also remain throughout the equilibration. However, during the equilibration, G658 is observed to bind with Lys7 of S15.

G752 contacts Tyr68 via hydrogen bonds G752:H5(1)– Tyr68:OH and G752:H5(2)–Tyr68:OH in the crystal structure, but these bond distances significantly increase to ∼10 Å after the MD equilibration runs for 1.5 ns or longer.

U751 contacts Thr24 and Gln27 via two intermolecular hydrogen bonds Gln27:HE21–U751:O2 and U751:HO′2– Thr24:OG1. The first bond fluctuates during the equilibration and thus results in a weak average bond length of 4.44 Å over the MD equilibration (Table 2 ). The second bond is weakly formed at the beginning of the equilibration, but it becomes stronger as the equilibration time evolves. This bond length fluctuates roughly between 5 and 2 Å, and the bond strengthens (2 Å) near the end of the equilibration.

G750 is involved in three hydrogen bonds with S15. G750:H2(2)–Gly22:O is a strong bond with an average bond length of 2.15 Å (maximum bond length 5.23 Å, mini mum bond length 1.61 Å). This bond is retained during the equilibration and is as strong as it is in the crystal structure. This bond is the strongest one in binding region 1. G750 also contacts Gln27 via two hydrogen bonds, Gln27:HE21– G750:N2 and Gln27:HE22–750:N2. The dynamic features of these two bonds are correlated during the equilibration with average bond lengths of 3.90 and 4.10 Å, respectively. Overall, these two bonds are weak with significant fluctuations in length.

C656 strongly interacts with Gln27 via hydrogen bond C656:HO′2–Gln27:OE1 in the crystal structure with a bond length of ∼2 Å. The equilibration results in large fluctuations of the bond length. However, the bond is strengthened during the final nanosecond of the equilibration.

G657 strongly interacts with Thr21 via hydrogen bond G657:HO′2–Thr21:O in the crystal structure with a bond length of ∼2 Å (Table 2 ). This strong bond is retained only for the first 2 ns of the equilibration, and then the bond length increases. This bond occasionally has a length of ∼2 Å during the second half of the equilibration, but the majority of the conformations show a bond length of ∼3 Å (Fig. 9 ). In addition, the interactions evident in the crystal between G750 and Gly22, and G750 and Asp20, as well as C754 and Tyr68, do not remain strong during the equilibration.

The equilibration time course results reveal that a strong charge interaction can occur between the negatively charged phosphate group of G658 and the positively charged residue Lys7 of S15. This forms the hydrogen bond Lys7:HZ3– G658:O1P with an average bond length of 2.79 Å. The maximum bond length of 7.68 Å occurs at the beginning of the equilibration, but it significantly decreases with a minimum bond length of 1.58 Å during the equilibration.

16S rRNA–S15 binding in region 2. In binding region 2, the rRNA–S15 interactions in the crystal structure occur between C739 and His41, between G666 and Ser51, between U740 and Ser51 and between G667 and Asp48, as well as between G668 and His45. The simulation shows that the first two interacting pairs do not retain strong association during the equilibration, but U740, G667 and G668 do retain a stable contact with S15. In addition, the equilibration results indicate that the negatively charged phosphate groups of G741 and G742 interact with the two positively charged residues, Arg34 and Arg37. These binding details are now described more completely.

Located in the upper portion of H22, G667 and G668 interact with His45 and Asp48 that are in the loop between helices α ₂ and α ₃ of S15. G667 contacts the side chain of Asp48 via two hydrogen bonds, G667:HO′2–Asp48:OD2 and G667:H2(2)–Asp48:OD2, and the side chain of Ser51 via hydrogen bond G667:H2(2)–Ser51:OG. These bonds have relatively large fluctuations (deviations of 1.15 and 0.82 Å) (Table 2 ). G667 occasionally forms a bond with Ser51, so that G667 weakly interacts with S15. G668:H2(2) strongly binds to His45:NE2 with an average bond distance of 2.18 Å during the equilibration, although a maximum bond distance of ∼9 Å occurs at the beginning of the MD equilibration.

On one side of the major groove of H22, the negatively charged oxygen atoms in the phosphate groups of U740, G741 and G742 interact with the side chains of the positively charged residues Args34 and Arg37 (Fig. 9 ). These hydrogen bonds strengthen during the equilibration. In particular, hydrogen bond G742:O2P–Arg34:HH22 is very strong with an average bond length of 2.28 Å (Table 2 ). These contacts, via the electrostatic interactions, are induced during the equilibration. In the crystal structure, the interactions between G741, G742 and Arg34, Arg37 are quite weak. These identified interactions may play a crucial role in stabilizing the upper part of H22 in a crystal‐like conformation.

16S rRNA–S15 binding in region 3. The GAAG tetraloop terminating H23a interacts with S15 in the crystal structure; A728 contacts Arg53 via the hydrogen bond A728:O1P– Arg53:HH12. A729 and G730 also weakly interact with His50 via the hydrogen bonds A729:H6(2)–His50:NE2m, G730:H1–Hig50:NE2 and G730:O6–Hig50:HD2. Our simulation shows that overall this tetraloop does not interact strongly with S15. Only the hydrogen bond A728:O1P– Arg53:HH12 is weakly retained with an average bond length of ∼3 Å during the last 3 ns of the equilibration. The hydrogen bond A729:H6(2)–His50:NE2 forms occasionally. The contact between G730 and His50 is basically lost during the equilibration (the average bond length is ∼6 Å). The simulation indicates that the interaction between the GAAG tetraloop with S15 may locally restrict the orientation of the GAAG tetraloop but does not play a crucial role in maintaining the bound complex of 16S rRNA and S15.

DISCUSSION

The simulation results are consistent with experimental findings

The dynamic characteristics of the core central domain of 16S rRNA generated using the MD simulations are similar to experimental descriptions. Experiments have demonstrated that a fully folded structure of the rRNA showed an acute inter‐helical angle between H20 and H22, and this angle increases to ∼120° in the absence of S15 ( 11 , 12 , 38 ). Similar conformations are produced by the current simulations for the bound and unbound rRNA. In the presence of S15, the bound rRNA remains in a highly folded, crystal‐like conformation. In the absence of S15, angle ϕ (Fig. 7 ) increases to 114° in the average structure of the unbound rRNA, compared to 80° in the crystal structure. These changes in the conformation of 16S rRNA were observed in the simulations using both the PME method and the GB/SA model. In addition, in order to examine the dependency of the fluctuating behavior of this rRNA on the starting structures of the simulations, we performed a simulation using an intermediate snapshot as the starting structure for the GB/SA calculation (data not shown). This intermediate structure also generates fluctuating conformations with the RMSD ranging between 2.5 and 6.6 Å relative to the crystal structure.

Binding interactions between 16S rRNA and S15

One of the interesting issues associated with 16S rRNA is the understanding of the mechanism by which the rRNA and S15 recognize and bind to each other. The ‘structure‐capture’ mechanism proposes that the unbound 16S rRNA alternates between extended conformations and a highly folded, crystal‐like conformation and that the latter conformation is recognized by and binds to S15 ( 31 ). Both the extended conformations and the highly folded conformations of the rRNA have been observed experimentally ( 38 ). The recognition of the rRNA by S15 may be critically determined by a structural matching between these two molecules. The current simulations provide a dynamic view of 16S rRNA conformations in the unbound state and provide support for the ‘structure‐capture’ binding mechanism.

Our results suggest that a decrease in the electrostatic energy concomitant with an increase in the solvation energy of the unbound 16S rRNA, and vice versa, results in the formation of both the highly folded, crystal‐like conformation and the extended conformations. The changes in the conformations occur primarily in the backbone geometry of H22. This geometry is critical for the binding of the rRNA to S15. These interactions include the binding of G658 to Lys7, U740 to Arg37 and G741 to Arg34 and Arg37, as well as G742 to Arg34. These interactions are the result of electrostatic interactions between the negative charges of the nucleotide phosphate groups and the positive charges of the amino acid residues of S15. As shown in Figure 0 A, as the unbound rRNA transiently forms a highly folded conformation, the intermolecular interactions between G658 and Lys7 take place on one side of the major groove and the interactions between U740, G741 and G742 and Arg34 and Arg37 occur on the other side of the groove. These interactions would be expected to maintain the stability of the major groove in the binding geometry.

In an extended conformation of the rRNA, the major groove of H22 is widened. The respective distances between nucleotides G741 and G742 and residues Arg34 and Arg37 increase. If Lys7 were able to bind to G658 on the lower side of the groove, then the widened groove would not allow Arg34 and Atg37 to contact G741 and G742 on the upper side (Fig. 0 B). Therefore, the increased distance does not confer the required geometry for the binding. Thus, the backbone geometry of H22 crucially determines whether or not the binding of the rRNA with S15 can take place.

In summary, the current MD simulations provide a dynamic view of the binding interactions between 16S rRNA and S15. The binding may take place when 16S rRNA is folded in an available conformation, and our results identify new binding sites between rRNA and S15. 16S rRNA can form this required conformation spontaneously but transiently in the absence of S15.

ACKNOWLEDGEMENTS

We thank Dr Jacob Maizel for his helpful discussions and the Laboratory of Experimental and Computational Biology for their support. We appreciate the computer support from the National Cancer Institute’s Advanced Biomedical Computing Center.

Figure 1. Sequence and secondary structure of the core domain of 16S rRNA in the crystal structure. Nucleotide numbering follows the E.coli sequence. Those nucleotides in lower case are not in the native sequence of 16S rRNA. The nucleotides involved in hydrogen bonds with S15 in the complex in three binding regions (BR1–BR3) are in red ( 9 ).

Figure 2. ( A ) A stereo view of the crystal structure of S15 binding with the core central domain of 16S rRNA from T.thermophilus ( 9 ). Helices 20, 21 22 and 23a are blue, yellow, green and pink, respectively, and the ribbon of S15 is red. ( B ) A stereo view of the average structure of the rRNA–S15 complex. The RNA (blue) is displayed with a ribbon representing its backbone with the heavy atoms of the nucleotides. S15 (light blue) is displayed only with a ribbon representing the backbone. The backbone segments of the binding nucleotides are red, while the binding residues of S15 are yellow.

Figure 3. Comparisons of the converted crystallographic B‐factors, RMSD _cry values (red) with the simulation RMSDs (blue). ( A ) The average RMSDs of phosphate atoms and C4′ atoms of the 16S rRNA. The numbering of the nucleotides is continuous in the plot, but the rRNA has two separate strands. The arrow points to the end of the first strand (A582–A675). ( B ) The average RMSDs of the C _α atoms of S15.

Figure 4. The RMSDs of the bound rRNA (blue) and the unbound rRNA (pink) relative to the starting structures of their equilibrations using the GB/SA model. The green circles are the approximate locations of snapshots A and B that have relatively small and large RMSDs, respectively, relative to the crystal. The RMSD curve has been smoothed by averaging the data over 20 ps.

Figure 5. Illustration of the backbone geometry of H22 in various conformations. The backbones are displayed using ribbons, and the hydrogen atoms are not displayed. ( A ) Snapshots with a small RMSD (orange) and a large RMSD (blue) relative to the crystal structure generated using the PME method. ( B ) The crystal structure of the unbound 16S rRNA (orange) is superposed onto the average structure of the unbound 16S rRNA derived from the equilibration performed on the unbound rRNA during the last 1 ns (blue) using the GB/SA model. ( C ) The fragment from the dashed box of (B) showing the backbone around C748. ( D ) The extended conformation (blue) and the crystal‐like conformation (orange) of H22 taken from the two snapshots A and B of the GB/SA calculation. ( E ) Superposition of H22 in snapshot B (purple) onto the average structure generated from the simulation performed on the 16S rRNA–S15 complex (orange). ( F ) Superposition of the crystal‐like snapshot A (blue) onto the average structure generated from the simulation performed on the 16S rRNA–S15 complex (orange). ( G ) Superposition of the crystal structure (orange) onto a constructed A‐form RNA (green). ( H ) Superposition of the average structure of the 16S rRNA (purple) in the absence of S15 onto the A‐form RNA fragment (green).

Figure 6. Distance between G657:P and U743:P over the equilibration performed on the unbound 16S rRNA. The curve has been smoothed by averaging the data over 30 ps.

Figure 7. Schematic representation of the relative orientations of H20, H21 and H22. The inter‐helical angle between H20 and H22 is measured using ϕ. A, B and C are located at the mass centers of three nucleotide groups (see text).

Figure 8. Binding sites of S15 in helix 22 and helix 23a of the 16S rRNA. The nucleotides contacting S15 in the crystal structure are red. Those binding with S15 only in the MD equilibration are green.

Figure 9. Binding sites of S15 in H22 occurring in the average structure of the complex. ( A ) The backbone of H22 is shown as a ribbon in conjunction with the heavy atoms of the nucleotides binding with S15. ( B ) Ribbon of S15 with atomic representations of the binding residues. Hydrogen atoms are not shown. ( C ) A stereo view of the binding of H22 to S15.

Figure 10. Structure of H22 and its influence on its binding with S15. ( A ) The structure of H22 is taken from the snapshot showing a crystal‐like conformation (left). The orange arrow on the left represent the width of the major groove. G658 is located close to residue Lys7 of S15 (red ribbon at right) and U740–G742 are located close to Arg34 and Arg37. ( B ) The structure of H22 taken from the snapshot showing an extended conformation (left). In this structure, the binding nucleotides U740–G742 are too distant to interact with Arg34 and Arg37 (dashed arrow).

Table 1.

Comparisons of the potential energies of snapshots A and B

	Total	Bond	Angle	DHR	VDW	Elec	EGB
84nt_A	–18 640	1038	1434	1772	–973	–1086	–20825
84nt_B	–18 668	1002	1431	1798	–1005	–3039	–18855

84nt_A and 84nt_B are for the 84 nt rRNA in the conformations of snapshots A and B. The energies (kcal/mol) refer to the total potential (Total), bond energy (Bond), angle energy (Angle), dihedral energy (DHR), van der Waals energy (VDW), electrostatic energy (Elec) and solvation energy (EGB).

Table 2.

The hydrogen bond distances of the binding residues (Å) derived from the MD equilibration simulation

	Hydrogen bonds	Average	Deviation	Max	Min
Binding region 1	Gln27:HE21–U751:O2	4.44	1.38	9.86	1.87
	U751:HO′2–Thr24:OG1	3.87	0.97	8.42	1.65
	G750:H2(2)–Gly22:O	2.15	0.46	5.23	1.61
	Gln27:HE21–G750:N2	3.90	0.92	8.24	2.18
	Gln27:HE22–G750:N2	4.10	0.90	8.43	2.47
	C656:HO′2–Gln27:OE1	3.40	1.16	7.06	1.57
	G657:HO′2–Thr21:O	4.16	1.61	8.25	1.68
	Lys7:HZ3–G658:O1P	2.79	0.79	7.68	1.58
Binding region 2	Lys7:HZ3–G658:O1P	2.69	1.15	6.54	1.42
	G667:H2(2)–Asp48:OD2	3.64	0.64	5.54	1.61
	G667:H2(2)–Ser51:OG	3.64	0.82	8.61	1.81
	G668:H2(2)–His45:NE2	2.18	0.29	9.00	1.75
	Arg37:HH21–U740:O1P	3.41	0.43	5.43	1.71
	Arg37:HH11–G741:O1P	3.34	0.33	5.43	1.82
	Arg34:HH12–G742:O2P	2.90	0.73	5.31	1.70
	Arg34:HH22–G742:O2P	2.28	0.61	4.15	1.56

These hydrogen bonds remain relatively stable when compared with other possible bonds that have large fluctuations and have bond distances >5 Å. The average bond distance (Average), deviations (Deviation), maximum bond distance (Max) and minimum bond distance (Min) are derived from the MD trajectory.