Article pubs.acs.org/JPCB
Specific Recognition Mechanism between RNA and the KH3 Domain of Nova‑2 Protein Qingfen Yu,†,∥ Wei Ye,†,∥ Cheng Jiang,† Ray Luo,*,‡ and Hai-Feng Chen*,†,§ †
State Key Laboratory of Microbial Metabolism, Department of Bioinformatics and Biostatistics, College of Life Sciences and Biotechnology, Shanghai Jiaotong University, 800 Dongchuan Road, Shanghai, 200240, China ‡ Department of Molecular Biology and Biochemistry, University of California, Irvine, California 92697-3900, United States § Shanghai Center for Bioinformation Technology, 1278 Keyuan Road, Shanghai, 200235, China S Supporting Information *
ABSTRACT: The KH3 domain of Nova-2 protein can precisely recognize the sequence-specific target RNA of human glycine receptor α2. However, the recognition mechanism between the protein and its target RNA is still hotly debated. In this study, molecular dynamic simulations in explicit solvent were utilized to understand the recognition mechanism. The structural analysis and the Kolmogorov−Smirnov P-test statistics reveal that the KH3 domain might obey a conformational selection mechanism upon RNA binding. However, the induced fit mechanism could not be completely ruled out. Unfolding kinetics indicates that the folding of RNA and KH3 happens first and then the binding between RNA and KH3 follows. Principle component analysis shows that the invariant Gly-Lys-Gly-Gly loop moves toward to the RNA molecule but the C-terminal domain moves away from the RNA molecule upon binding. These specific dominant motions were hypothesized to stabilize the complex structure. The hydrophobic and hydrogen bonding interactions were found to be the driving forces for the specific recognition, in contrast to the dominant electrostatic interactions for nonspecific recognition.
■
hydrophobic cleft consists of helices α1 and α2 and the edge of strand β2. The invariable Gly-Lys-Gly-Gly loop and the variable segment flank the hydrophobic cleft. The core recognition sequence 5′-UCAC-3′ (nucleotides 12−15) lies upon the cleft and is gripped by the invariant loop and the variable segment. The RNA molecule consists of a standard Aform hairpin with four Watson−Crick base pairs (G1:C20, A2:U19, G3:C18, G4:C17). The complex structure is shown in Figure 1. The previous works have demonstrated that the KH3 domain is sufficient and necessary to recognize the core sequence 5′-UCAY-3′ of the target single-strand RNA, with high affinity and sequence specificity.1−3,5,10,11 Structural comparisons of various KH-type proteins and their target single-strand RNAs reveal some useful binding rules.12,13 A systematic analysis of the RNA motif displays recognition preferences and functions of sequence-specific RNA-binding proteins.14 However, due to the complicated and subtle interactions between proteins and RNAs,15 the recognition mechanism is still poorly understood. To gain better insight into the recognition mechanism of KH-type proteins to their target sequence-specific RNAs, we
INTRODUCTION Nova proteins are pre-mRNA splicing factors that are restricted to express in the central nervous system.1−3 The K homology (KH) domains of these proteins can interact with single-strand RNAs.4 The KH domains include a conserved motif for sequence-specific RNA binding, whose function is RNA recognition, splicing, or metabolism in neurons.1−3,5 The KH motif is constituted of approximately 70 amino acids and includes a characteristic hydrophobic core, an invariant Gly-XX-Gly (X usually represents Arg, Lys, or Gly) loop, and a variable linker segment.6−8 Both Nova-1 and Nova-2 proteins contain three KH motifs (KH1−3).1−3,5,7 The KH3 domain exclusively recognizes the core sequence 5′-UCAY-3′ (Y usually represents Cyt or Ura) of target GABA-A3,9,10 and glycine receptor pre-mRNA.3,10 The defection in the pre-mRNA specific splicing can result in motor incapacitation or death in mice.3,10 The complex crystal structure of the Nova-2 KH3 domain and the stem-loop RNA of human glycine receptor α2 was released in 2000 (pdb code: 1EC6).10 The resolution of the complex is 2.4 Å. There are 87 residues for the KH3 domain and 20 bases for RNA. The KH3 domain of Nova-2 consists of three α helices (α1−3) and three β strands (β1−3). The invariant Gly-Lys-Gly-Gly (amino acids 19−22) loop connects helices α1 and α2. The amino acids from Lys40 to Arg49 form the variable segment between strands β2 and β3. The © 2014 American Chemical Society
Received: August 5, 2014 Revised: October 8, 2014 Published: October 9, 2014 12426
dx.doi.org/10.1021/jp5079289 | J. Phys. Chem. B 2014, 118, 12426−12434
The Journal of Physical Chemistry B
Article
All solvated systems were first minimized by 1000-step steepest descent to remove any structural clash, followed by 20 ps heating up and brief equilibration in the NPT ensemble at 298 K. A Berendsen barostat was used in the NPT ensemble. Pressure was set to 1 bar in the NPT ensemble. The time step was 2 fs with a friction constant of 1 ps−1 in Langevin dynamics. Ten independent trajectories of 20.0 ns each in the NPT ensemble at 298 K were simulated with PMEMD of AMBER 11.18 Then, 10 independent high-temperature trajectories of 20.0 ns each were performed to investigate comformational reorganization pathways for each solvated system in the NVT ensemble at 498 K. Four mutant systems (U12C, C13U, A14G, and C15U) were simulated for five trajectories of 40.0 ns each at 298 K. The simulation conditions for all systems are listed in Table 1. A total of 2 μs trajectories were collected for these solvated systems. It took about 137 500 CPU hours in the inhouse Xeon (3.0 Ghz) cluster. Native contacts of the apo and bound states for KH3 and RNA were monitored to detect the binding-induced reorganization processes. It was found the 18 ns were sufficient to reach the equilibrium state for both apo and bound states at 498 K. Therefore, the first 18 ns (a total 180 ns for each system) were used to monitor the reorganization kinetics. The remaining 2 ns (a total of 20.0 ns for each system) were used to study the reorganized equilibrium state. Principle Component Analysis. A few essential degrees of freedom obtained from principle component analysis (PCA)24 can be used to describe the dominant motions of the high dimensional molecular systems.25,26 PCA assumes that the essential degrees of freedom can be used to identify the slow motions relevant to biological functions.25−28 To eliminate the overall translation and rotation motions, all the structures were superimposed onto a particular reference conformation. The fitted internal motion is represented by a trajectory X(t), where X is a column vector with 3N atomic coordinates and t is the time index representing each snapshot of the trajectory. The correlations of atomic fluctuations are represented by the covariance matrix C (eq 1):
Figure 1. Ribbon representation of the X-ray structure for the third KH domain (KH3) of Nova-2 protein and RNA of glycine receptor α2. Helices α1, α2, and α3 (light blue, cyan, and red), strands β1, β2, and β3 (blue, green, and yellow), the invariant Gly-X-X-Gly loop, and the variable segment of KH3 are labeled. Nucleotides A11 to C15 (blue, cyan, green, yellow, and orange) of RNA are labeled. The N and C terminal domain of KH3 and 5′ and 3′ ends of RNA are labeled.
utilized all-atom molecular dynamics (MD) simulations in explicit solvent to analyze the coupling between binding and conformational reorganization in the complex of the Nova-2 KH3 domain and its target sequence-specific RNA.16,17 At first, we revealed the binding mode between RNA and KH3, and then, the unfolding simulation was used to reveal the dissociation and reorganization kinetics of RNA. The results of unfolding kinetics suggest that the folding of RNA happens first, which is then followed by the binding of RNA. Finally, KS P-test statistics was used to evaluate the significance of induced fit and conformational selection for RNA and KH3.
■
MATERIALS AND METHODS Molecular Dynamic Simulations. The atomic coordinates of the Nova-2 KH3 domain/RNA complex were obtained from X-ray crystal structure (pdb code: 1EC6).10 Point mutations of U12C, C13U, A14G, and C15U were modeled with SYBYL-X software. All hydrogen atoms were added using the LEAP module of AMBER 11.18 Counterions were added to maintain system neutrality. All systems were solvated in a truncated octahedron box of TIP3P waters with a buffer of 10 Å.19 The particle mesh Ewald (PME)20 method was applied to handle long-range electrostatic interactions with the default setting in AMBER 11.18 The parm99SBildn force field21,22 was used to model the KH3 domain and RNA. The SHAKE algorithm23 was employed to constrain bonds including hydrogen atoms.
C = cov(X ) = ⟨(X − ⟨X ⟩)(X − ⟨X ⟩)T ⟩
(1)
where ⟨ ⟩ represents the average over the structures sampled in the MD trajectories. The symmetric covariance matrix C can be diagonalized by the orthogonal coordinate transformation T (eq 2)
Λ = TTCT
(2)
where the elements of Λ represent the eigenvalues and the ith column vector of T is the eigenvector for the ith eigenvalue. The eigenvector represents the coordinate axis in the
Table 1. Simulation Conditions for All Systems system WT
mutant
ions +
waters
water box size (Å3)
temperature (K)
simulation time (ns)
number of trajectories
6730
286 569
298 498 298 498 298 498 298
20 20 20 20 20 20 40
10 10 10 10 10 10 5
bound KH3-RNA
15 Na
Apo-KH3
5 Cl−
6010
252 280
Apo-RNA
20 Na+
4899
199 453
bound U12C bound C13U bound A14G
14 Na+
6688 6898 6673
335 919 346 659 335 557 12427
dx.doi.org/10.1021/jp5079289 | J. Phys. Chem. B 2014, 118, 12426−12434
The Journal of Physical Chemistry B
Article
Figure 2. Variations and landscapes of distance differences for apo-RNA, apo-KH3, and complex: (A) C5′ and Cα atomic fluctuation for RNA and KH3, respectively; (B) α, β, γ, δ, ε, and ζ variation for RNA and Φ/ψ variation for KH3; (C) landscape of distance differences for RNA; (D) landscape of distance differences for KH3.
where Ebond is the molecular mechanics bond energy including bond, angle, and dihedral energies; Evdw is the van der Waals energy contribution; Eelec is the electrostatic energy; GPB and GSA are polar and nonpolar contributions to the solvation energy; T is the absolute temperature; and SS is the solute entropy. All the binding free energies were calculated using the GB and PB models in MMPBSA.py.29 Considering that the equilibration period may affect the result, only the last 10 ns of MD trajectories was used to calculate the free binding energy. Other Data Analysis Details. The tertiary contact assignments were performed with in-house software.30−32 Any two nonadjacent residues are in native contact when their center of mass is less than 6.5 Å apart. Any positively and negatively charged residues are in electrostatic interaction when their centers of mass are less than 11 Å apart.33 A stable interaction/contact is defined as a contact occurring in more than 60% of snapshots in all simulation trajectories. The fraction of native tertiary contact (Qf) and the fraction of native binding contact (Qb) were used to study the reorganization and dissociation kinetics. All fitted kinetics curves were analyzed with Origin 8.5. Analysis of Binding Mechanism. The mechanism of specific binding is one of the most interesting issues in our study. There are two mainstream hypotheses to explain ligand binding: (1) the “induced fit” model,34−36 emphasizing the fact that the conformation of the receptors is induced to change upon the ligand binding (binding first), and (2) the “conformational selection” model,35−43 where the ligand will directly select the most optimal conformation of receptor among several different and competing conformations. If the binding obeys the induced-fit mechanism, significant conformational change must occur in the receptor upon ligand binding,
conformational space, and the eigenvalue corresponds to the mean square deviance along the corresponding eigenvector. The eigenvectors are sorted according to the descend order of eigenvalues. Projection of the structures onto the first few important eigenvectors can be used to display the collective distribution in the conformational space as P = TT(X − ⟨X ⟩)
(3)
Finally, the projections can be transformed back to the Cartesian coordinates for visualizing by X * = TP + ⟨X ⟩
(4)
The detail analysis step for PCA was as follows: The input data were the conformers of bound and apo states at room temperature. PTRAJ was first used to compute eigenvalues and corresponding eigenvectors of the covariance matrices. Next, VMD was used to project three-dimensional structural snapshots along each principle component. Finally, PYMOL was used to visualize the final results of PCA. MM/PBSA Energy Calculation. Free-energy calculation together with MD simulations can provide quantitative predictions of protein−RNA binding energies. The free energy of binding, ∇Gbind, is given by ∇G bind = G P + D − (G P + G D)
where GP+D, GP, and GD are the free energies of the complex, the isolated protein, and RNA, respectively. In the MM/PBSA approach, each free energy term in the previous equation is calculated as G = E bond + Evdw + Eelec + G PB + GSA − TSS 12428
dx.doi.org/10.1021/jp5079289 | J. Phys. Chem. B 2014, 118, 12426−12434
The Journal of Physical Chemistry B
Article
This suggests that the stability of the secondary structures of KH3 does not change much upon RNA binding. To clearly illustrate the conformational difference, the landscapes of distance differences between bound and apo states for RNA and KH3 are shown in Figure 2C and D. The landscapes may reflect the relative conformational changes of RNA and KH3 backbones. The dark red areas show that the differences between nucleotides 7−11 and 13−20 are positive. This indicates that the loop and stem regions are stretched upon KH3 binding. The dark blue areas show that the distance differences are negative between nucleotides 1−5 and 11−15, suggesting that the stem region moves closer to the core recognition sequence. In contrast, most regions of KH3 are not changed dramatically upon binding from the distance differences. The positive differences were found between amino acids 40−50 and 83−87, indicating that the variable segments connecting strands β2 and β3 are stretched away from the Cterminal domain. However, the distance differences between amino acids 12−30 and 80−87 are negative, suggesting that both helices α1 and α2 are close to the C-terminal domain upon RNA binding. Therefore, the C-terminal domain might play a key role in the formation of the KH3/RNA complex. To explore the conformational differences between bound and apo states for KH3 and RNA, the energy landscapes with the reaction coordinate of RMSD and the radius of gyration (Rg) were analyzed, as shown in Figure 3. It was found that bound-RNA is more elongated than apo-RNA with the maximum Rg changed from 16 to 18 Å upon binding (Figure 3A and B). The conformation distribution of bound-KH3 is also more diffusive than that of apo-KH3 (Figure 3C and D). There is one free energy basin for apo-RNA, with Rg between 11.5 and 14 Å and C5′ atomic RMSD between 2 and 6 Å. For the bound-RNA, the free energy basin is at the conformational space with Rg from 15.5 to 16.5 Å and RMSD from 2 to 4 Å. There is one free energy basin for apo-KH3 with Rg between 13 and 14.5 Å. For bound-KH3, the basin is at the conformational space with Rg from 15.5 to 16.5 Å and Cα atomic RMSD from 1 to 3 Å. The average structures corresponding to the free energy minimum for RNA and KH3 are shown in Figure 3E. Structural comparison shows that the loop of bound-RNA is more stretched than that of apoRNA. In summary, bound-KH3 is mostly similar to that of apoKH3 except for the C-terminal domain, which is more extended than apo-KH3, consistent with the experimental observation.10 To study the driving force for these binding-induced conformational adjustments, the electrostatic, hydrophobic, and hydrogen-binding interactions between KH3 and RNA were analyzed. Figure 4A illustrates the stable electrostatic contacts in all 10 independent trajectories. Fifteen electrostatic interactions were found between the positively charged amino acids and negatively charged phosphates of RNA, with populations higher than 60%. The positively charged amino acids, such as Lys20, Lys23, Arg35, Lys40, Lys41, Arg49, Arg51, and Arg80, form salt-bridge interactions with phosphates of C6, C7, U8, A9, U12, C13, A14, and C16 in RNA. It is interesting to note that Lys40, Lys41, and Arg49 are located at the variable loops of the protein. This suggests that the variable loops play a key role in the formation of the complex, which is consistent with the experimental observations.10 The hydrogen bonds are illustrated in Figure 4B. Seven stable hydrogen bonds were found between KH3 (Gly19, Lys20, Gly21, Ile38, and Arg51) and RNA (U12, C13, and A14). Residues Gly19, Lys20, and Gly21 belong to the conserved Gly-Lys-Gly-Gly loop,
especially near the binding site. Otherwise, the binding tends to obey the conformational-selection mechanism. Following this logic, the existence of induced fit could be inferred from the RMSD between a bound structure and its most similar apo structure (which has the lowest RMSD from the bound structure). In the induced-fit model, the relative magnitude of conformational selection could be characterized from the average RMSD between this most similar apo structure and the other apo structures.44 The RMSD for induced fit of 10 trajectories was given as a function of distance from the binding site. Because the RMSD values are not normally distributed, which is not suitable for a t test, a standard two-sample Kolmogorov−Smirnov (KS) test was used to evaluate the statistical significance of the RMSDs’ variations. Recent studies have demonstrated that many binding processes may employ a mixed mechanism: selecting an optimal conformation globally and inducing the local structural adjustment of the ligandbinding site.45−47 Atomic RMSDs of all the bases between aligned bound and apo structures were first calculated for all distance ranges from the binding site. A KS P test was employed to evaluate the significance of the difference between RMSDs of atoms included in a certain distance from the binding site and the atomic RMSDs of the whole structure. A more detailed description on the identification of the binding mechanism can be found in ref 44. The detail analysis step for binding mechanism was as follows: First, the last (fully relaxed) 500 conformers from each of the bound and apo sets of trajectories were extracted. Next, the global RMSD difference was calculated for these conformers between bound and apo states. Then, we ordered these conformer pairs according to the global RMSD and obtained 500 pair conformers for bound and apo states. For these 500 pair conformers, RMSDs were recalculated as a function of distance from the binding site. Finally, the KS test was used to evaluate the significance of the RMSD-distance profiles.
■
RESULTS Binding Mode between the KH3 Domain and RNA. A limited number of trajectories for MD simulations (5−10) are sufficient to capture the average properties of a biomolecule.48 To study the recognition between the Nova-2 KH3 domain (KH3) and its target RNA, 10 independent trajectories of 20.0 ns each for apo-RNA, apo-KH3, and their complex were simulated at room temperature (298 K), respectively. C5′ and Cα atomic variances for apo and bound states are illustrated in Figure 2A. The C5′ fluctuations of bound-RNA are significantly smaller than those of apo-RNA, especially in the region of core recognition sequence 5′-UCAC-3′ (nucleotides 12−15). This indicates that bound-RNA becomes more stable upon KH3 binding. The Cα fluctuations of bound- and apo-KH3 are more similar. This suggests that the structures of bound-KH3 are not significantly different from those of apo-KH3. These results are consistent with the experimental observations.10,13 In general, the angle parameters (α, β, γ, δ, ε, and ζ) are used to describe the secondary structure of RNA. It was found that α, β, γ, δ, ε, and ζ fluctuations for bound-RNA were lower than those for apo-RNA, especially the core recognition sequence (shown in Figure 2B). This indicates that the secondary structure of apo-RNA is more disordered. However, the Φ/ψ variations of bound-KH3 were similar to those of apo-KH3. 12429
dx.doi.org/10.1021/jp5079289 | J. Phys. Chem. B 2014, 118, 12426−12434
The Journal of Physical Chemistry B
Article
Figure 4. Interactions between RNA and KH3. (A) Electrostatic contacts. 1 for U12/Lys41; 2 for A14/Lys41; 3 for C13/Lys41; 4 for U8/Arg49; 5 for U12/Arg51; 6 for C13/Lys20; 7 for A9/Arg49; 8 for C7/Lys41; 9 for U12/Lys20; 10 for C13/Lys23; 11 for U8/Lys41; 12 for U12/Arg80; 13 for C16/Arg35; 14 for C7/Lys40; 15 for C6/ Lys41. (B) Hydrogen bonds. 1 for N3(C13)/NH1(Arg51); 2 for N6(A14)/O(Ile38); 3 for O2(C13)/NH2(Arg51); 4 for O1P(C13)/ N(Gly21); 5 for O2′(U12)/O(Gly19); 6 for N1(A14)/N(Ile38); 7 for O2(U12)/N(Lys20). (C) Hydrophobic contacts. 1 for C13/ Leu18; 2 for A14/Ile38; 3 for C13/Val14; 4 for A14/Leu18; 5 for U12/Ala16; 6 for A14/Leu25.
region. The dominant motions of the complex are shown in Figure 5E and F. The major motion is found in the invariant loop, moving toward the target RNA and away from the Cterminal domain. This shows that the contacts between the invariant loop and the target RNA play key roles in the formation of the complex. The target RNA in the complex has a bending motion at the hairpin structure that is different from the stretching motion in apo-RNA. It implies that the RNA molecule becomes more stable and ordered upon KH3 binding (Figure 5F). Dissociation and Reorganization Kinetics. In order to reveal the recognition process, the unfolding kinetics was employed in this study. The fraction of native tertiary contacts (Qf) and the fraction of native binding contacts (Qb) were used to analyze reorganization and dissociation kinetics. Time evolutions of Qf and Qb for bound and apo states are shown in Figure 6, suggesting that the tertiary reorganization and dissociation can be fitted well by single exponential functions, indicating first-order kinetics in the NVT ensemble at the high temperature condition. The fitted kinetics data are listed in Table 2. The kinetics analysis shows that the dissociation half time is 2.965 ± 0.014 ns. The tertiary reorganization half time is 5.779 ± 0.020 and 7.374 ± 0.027 ns for the complex and bound-KH3, respectively. This indicates that the dissociation of the complex is faster than that of the complex and bound-KH3 reorganization. However, the dissociation is slower than the tertiary reorganization of bound-RNA with 1.003 ± 0.003 ns half time. The half times of tertiary reorganization are 0.816 ± 0.003 and 5.376 ± 0.018 ns for apo-RNA and apo-KH3, respectively. This suggests that the tertiary reorganization time scale of the RNA molecule mostly stays the same upon binding, while the tertiary reorganization time scale of the KH3 domain is changed noticeably upon binding.
Figure 3. Free energy landscapes with respect to Rg and RMSD and average structures corresponding to the energy minimum for apo and bound states of RNA and KH3: (A) apo-RNA; (B) bound-RNA; (C) apo-KH3; (D) bound-KH3; (E) average structures for apo-RNA, bound-RNA, apo-KH3, and bound-KH3.
indicating that the invariant loop is critical in the conformational adjustments. Besides electrostatic contacts and hydrogen bonds, there are also six stable hydrophobic interactions, shown in Figure 4C, such as U12/Ala16, C13/Val14, C13/Leu18, A14/Leu18, A14/Leu25, and A14/Ile38. Hydrophobic residues Val14, Ala16, and Leu18 are located at helix α1 of the broad binding cleft, next to Leu25 and Ile38 in the binding cleft. Remarkably, most electrostatic interactions, hydrogen bonds, and hydrophobic contacts are formed with core nucleotides of U12, C13, and A14. This suggests that the three core nucleotides might be crucial to the recognition between KH3 and RNA, in agreement with the experimental observations.10,11 Principle Component Analysis. To analyze the overall conformational changes upon binding, the principle component analysis was used to examine the dominant motions of RNA, KH3, and their complex. For apo-RNA, the dominant motion as represented by the first two principal components (PC1 and PC2) can be described as a stretch motion of the RNA chain, as shown in Figure 5A and B. Here the first two principal components account for approximate 63% of all the fluctuations. For apo-KH3, there is a dominant swing of the C-terminal region, as shown in Figure 5C and D. In addition, the invariant loop (Gly19-Gly22), the variant segment (Lys40Arg49), and the helix α3 move away from the C-terminal 12430
dx.doi.org/10.1021/jp5079289 | J. Phys. Chem. B 2014, 118, 12426−12434
The Journal of Physical Chemistry B
Article
Figure 5. Dominant motions as represented by the first two principle components (PC1 and PC2) for apo-RNA, apo-KH3, and their complex. The arrows indicate the direction of motion for each residue. Nucleotides 11A to 15C of RNA are colored with blue, cyan, green, yellow, and orange. (A) PC1 for apo-RNA; (B) PC2 for apo-RNA; (C) PC1 for apo-KH3; (D) PC2 for apo-KH3; (E) PC1 for complex; (F) PC2 for complex.
are in good agreement with the structural analysis that the invariant loop forms important interactions with RNA. Furthermore, experimental data indicate that the variable segment (Lys40-Arg49) interacts with C6, C7, U8, A9, C13, and C16,10 which is consistent with our simulations where positive charged amino acids of Lys40, Lys41, and Arg49 form stable electrostatic contacts with C6, C7, U8, A9, U12, C13, A14, and C16. The X-ray structure suggests that U12 forms van der Waals contacts with Ala16,10 which was also found as a stable hydrophobic contact. Structural data also reveal the existence of a complicated interaction network among C13, Glu11, Val14, Leu18, and Arg51.10 Most of these interactions also appear as stable in our simulations. Finally, the experimental data show that A14 interacts directly with amino acids of Leu18, Leu25, and Ile38 in the hydrophobic α/β platform,10 in agreement with our simulations. In addition, experiment demonstrates that the C-terminal amino acids of KH3 are crucial for RNA binding,10,11 which is consistent with our observation that there is a stable electrostatic contact between Arg80 and U12. In summary, these results are in accordance with the experimental observation. Comparison with Nonspecific Sequence RNA Recognition. MD simulations of the RNA/KH3 complex reveal that electrostatic interactions between RNA bases and amino acids
Figure 6. Kinetics fitting of tertiary reorganization (Qf) and tertiary dissociation (Qb) for apo and bound states for RNA and KH3.
■
DISCUSSION Comparison with Experiment. The X-ray structure shows that the invariant loop Gly-Lys-Gly-Gly (residue 19−22) forms important interactions with U12, C13, and A14.10 Two stable hydrogen bonds for U12/Gly19 and C13/Gly21 and two stable electrostatic contacts for U12/Lys20 and C13/Lys20 were found in the room-temperature MD simulations. These results 12431
dx.doi.org/10.1021/jp5079289 | J. Phys. Chem. B 2014, 118, 12426−12434
The Journal of Physical Chemistry B
Article
Table 2. RNA and KH3 Unfolding Kinetics Constantsa τ (ns) apo-RNA bound-RNA apo-KH3 bound-KH3 complex complex a
Qf Qf Qf Qf Qf Qb
0.816 1.003 5.376 7.374 5.779 2.965
± ± ± ± ± ±
A
0.003 0.003 0.018 0.027 0.020 0.014
0.995 0.934 (0.563 (0.655 (0.636 0.858
± ± ± ± ± ±
0.003 0.002 6.465) × 10−4 7.494) × 10−4 7.444) × 10−4 0.002
R2
B (0.021 (0.005 (0.200 (0.153 (0.160 (0.057
± ± ± ± ± ±
3.256) 2.096) 5.282) 8.691) 6.757) 7.118)
× × × × × ×
−4
10 10−4 10−4 10−4 10−4 10−4
0.963 0.986 0.988 0.991 0.988 0.961
All curves are fitted by y = A exp(−t/τ) + B.
Information). However, it should be pointed out that the entropy effect that might be important in this relative affinity analysis was not considered in the MMPBSA calculation due to the difficulty to obtain a reliable estimation. Recognition Mechanism. The unfolding kinetics reveals the recognition process that the folding of RNA and KH3 happens first and then the binding follows. Conformational selection and induced fit are further used to interpret the recognition between RNA and protein.36,53 According to the conformational selection paradigm,54−56 various conformational ensembles explore the free energy landscapes that correspond to diverse stable unbound states in equilibrium. During the binding process, the favorable conformation compatible with binding selectively stabilizes and the population of conformational ensembles shifts toward the stabilizing state. However, the induced fit scenario proposes that the favorable conformation results from significant changes of unbound ensembles upon allosteric binding.57 It is worth pointing out that conformational selection and induced fit models cannot be distinguished in a clear-cut manner.36 Indeed, both mechanisms may be involved in some systems.58,59 To explore the recognition mechanism, the average RMSD deviations of the bound conformation and the globally most similar apo conformation were analyzed as a function of distance from the centroid of the binding interface60 and shown in Figure 7A−C. In the case including the C-tails, the average RMSD gradually increases for KH3. In the case excluding the C-tails, the average RMSD of KH3 significantly decreases. This
are crucial to the complex stability. Such interactions are also common in nonspecific binding by single-stranded RNAs.12,49,50 However, nonspecific binding proteins predominantly form hydrophilic interactions with RNA backbone groups.51 This is different from the RNA/KH3 complex that is also stabilized by hydrophobic contacts. In addition, most hydrogen bonds are formed with the RNA backbone groups in nonspecific binding,52 mostly responsible for the binding stability of nonspecific sequence single-stranded RNAs.12,51,52 In contrast, it is found that the KH3 domain forms hydrogen bonds with both phosphate groups and base moieties of RNA, where the binding to the base moieties of RNA plays key roles in the stability of the complex. Overall, the complex hydrogenbonding network is necessary and important for binding the unique nucleotide in the RNA/KH3 complex.10 Mutation Analysis. Both the X-ray crystal structure10 and MD simulation analysis show that the core sequence from nucleotides U12 to C15 plays key roles in recognition of the KH3 domain. In order to reveal the mutation effects of each nucleotide in the core sequence for binding, MMPBSA was used to calculate the binding free energy for mutants U12C, C13U, A14G, and C15U and the wild type of the complex, as summarized in Table 3. The binding free energy for the wild Table 3. Binding Free Energies of Wild Type and Mutants with Generalized Born and Poisson Boltzmann Models model wild type U12C C13U A14G C15U
GB (kcal/mol)
PB (kcal/mol)
−82.9 −96.4 −81.3 −105.4 −110.6
−80.2 −90.0 −75.9 −101.4 −102.6
± ± ± ± ±
9.7 8.8 17.9 16.3 21.5
± ± ± ± ±
11.8 4.9 12.4 14.2 17.1
type is −80.151 ± 11.836 kcal/mol with the PB model, in comparison with −89.993 ± 4.907 kcal/mol for U12C. This is consistent with the experimental observation that U12 could be replaced by C12.10,11 The binding free energy lost is about 4.236 kcal/mol for C13U, which is in agreement with the experimental observation that C13 is the key base for the stability of the RNA/KH3 complex.10,11 For position 15, the nucleotide is always Cyt or Ura.10,11 Indeed the binding free energy for C15U is −102.575 ± 17.126 kcal/mol. This suggests that C15U mutant should also be favorable for binding of the KH3 domain. However, the binding free energy for A14G is −101.428 ± 14.236 kcal/mol, which is higher than that of the wild type. This is in conflict with the experimental observation that position 14 is entirely specific.10,11 The simulation data can be explained by the following observations: A14G introduces more stable interactions in the RNA/KH3 complexthere are more electrostatic contacts and hydrogen bonds in A14G than in the wild type (shown in Tables S1 and S2, Supporting
Figure 7. Local conformational RMSD differences for 500 pairs bound and the globally most similar apo conformations as a function of distance from the centroid of the binding partner for KH3 including C-tails (A), KH3 excluding C-tails (B), and RNA (C). Statistical significance of the conformational difference in KH3 including C-tails (D), KH3 excluding C-tails (E), and RNA (F). 12432
dx.doi.org/10.1021/jp5079289 | J. Phys. Chem. B 2014, 118, 12426−12434
The Journal of Physical Chemistry B
■
suggests that the area of conformational change for KH3 is focused on the C-tails. However, the RMSD gradually decreases for RNA. This suggests that conformational adjustment in RNA is focused nearby the binding interface and RNA obeys an induced-fit mechanism upon binding to the KH3 domain. In order to investigate the statistical significance of the local conformational deviations, the two-sample Kolmogorov− Smirnov (KS) P test61 was conducted. Note that the KS test, as a nonparametric test, is a good choice for this study because the distributions of magnitudes of atomic deviations do not fit well to any distribution for parametric test. Figure 7C and D illustrated the median of P values and the fraction with P ≤ 0.01 for all 500 pairs of KH3 conformations in each distance bin. The fraction of typical P values was less than 0.4 (even 0.3). This suggests that the local conformational changes are not significant for the KH3 domain, especially when excluding Ctails. According to the hypothesis of the nonparametric test, the results indicate that KH3 might obey a conformational selection mechanism upon binding to RNA. For RNA, the conformational difference was statistically significant up to ∼20 Å away from the binding interface with the median P values typically 0.8. This suggests that the local conformational changes are significant for RNA. Thus, the sequence-specific recognition of RNA follows an induced-fit mechanism on the RNA side. It is consistent with the recognition mechanisms for kink-turn sRNA/L7Ae, RNA/U1A, and other RNA and protein complexes.53,62,63
AUTHOR INFORMATION
Corresponding Authors
*E-mail:
[email protected]. Phone: 1-949-8249528. Fax: 1-9498248551. *E-mail:
[email protected]. Phone: 86-21-34204348. Fax: 86-21-34204348. Author Contributions ∥
(Q.Y. & W.Y.) These authors contributed equally to this work.
Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported by Center for HPC at Shanghai Jiao Tong University, grants from the Ministry of Science and Technology of China (2012CB721003), the National Hightech R&D Program of China (863 Program) (2014AA021502), the National Natural Science Foundation of China (J1210047 and 31271403), the Innovation Program of the Shanghai Education Committee (Grant No. 12ZZ023), and Medical Engineering Cross Fund of Shanghai Jiaotong University (YG2013MS68).
■
REFERENCES
(1) Buckanovich, R. J.; Darnell, R. B. The neuronal RNA binding protein Nova-1 recognizes specific RNA targets in vitro and in vivo. Mol. Cell. Biol. 1997, 17, 3194−3201. (2) Yang, Y. Y. L.; Yin, G. L.; Darnell, R. B. The neuronal RNAbinding protein Nova-2 is implicated as the autoantigen targeted in POMA patients with dementia. Proc. Natl. Acad. Sci. U. S. A. 1998, 95, 13254−13259. (3) Jensen, K. B.; Dredge, B. K.; Stefani, G.; et al. Nova-1 regulates neuron-specific alternative splicing and is essential for neuronal viability. Neuron 2000, 25, 359−371. (4) Siomi, H.; J.Matunis, M.; Michael, W. M.; et al. The pre-mRNA binding K protein contains a novel evolutionarily conserved motif. Nucleic Acids Res. 1993, 21, 1193−1198. (5) Buckanovich, R. J.; Yang, Y. Y. L.; Darnell, R. B. The onconeural antigen Nova-1 is a neuron-specific RNA-binding protein, the activity of which is inhibited by paraneoplastic antibodies. J. Neurosci. 1996, 16, 1114−1122. (6) Musco, G.; Stier, G.; Joseph, C.; et al. Three-Dimensional Structure and Stability of the KH Domain: Molecular Insights into the Fragile X Syndrome. Cell 1996, 85, 237−245. (7) Lewis, H. A.; Chen, H.; Edo, C.; et al. Crystal structures of Nova1 and Nova-2 K-homology RNA-binding domains. Structure 1999, 7, 191−203. (8) Grishin, N. V. KH domain: one motif, two folds. Nucleic Acids Res. 2001, 29, 638−643. (9) Dredge, B. K.; Darnell, R. B. Nova Regulates GABAA Receptor 2 Alternative Splicing via a Distal Downstream UCAU-Rich Intronic Splicing Enhancer. Mol. Cell. Biol. 2003, 23, 4687−4700. (10) Lewis, H. A.; Musunuru, K.; Jensen, K. B.; et al. Specific RNA Binding by a Nova KH Domain: Implications for Paraneoplastic Disease and the Fragile X Syndrome. Cell 2000, 100, 323−332. (11) Jensen, K. B.; Musunuru, K.; Lewis, H. A.; et al. The tetranucleotide UCAY directs the specific recognition of RNA by the Nova K-homology 3 domain. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 5740−5745. (12) Auweter, S. D.; Oberstrass, F. C.; Allain, F. H. T. Sequencespecific binding of single-stranded RNA: is there a code for recognition? Nucleic Acids Res. 2006, 34, 4943−4959. (13) Valverde, R.; Edwards, L.; Regan, L. Structure and function of KH domains. FEBS J. 2008, 275, 2712−2726. (14) Ray, D.; Kazan, H.; Cook, K. B.; et al. A compendium of RNAbinding motifs for decoding gene regulation. Nature 2013, 499, 172− 177.
■
CONCLUSION Molecular dynamics (MD) simulations were used to study the recognition mechanism between sequence-specific binding between the Nova-2 KH3 domain and its target RNA. The unfolding kinetics indicates that the folding of RNA and KH3 happens first and then the binding follows. The average RMSD analysis and Kolomogorov−Smirnov P test statistics further indicated that the KH3 domain may follow a conformationalselection mechanism in its binding to DNA. This is consistent with the binding mode analysis at 298 K that the secondary and tertiary structures of bound-KH3 are similar to those of apoKH3 except for the C-terminal domain. However, RNA may obey an induced-fit mechanism for its recognition to the KH3 domain. Both the loop and stem of bound-RNA become more stable and structured upon binding to the KH3 domain. The principle component analysis shows that the dominant motion between the KH3 C-terminal domain and RNA was essential for the formation of a complex. Furthermore, the invariant GlyLys-Gly-Gly loop plays a key role in stability. This is further confirmed by the stability of hydrogen bonds between the invariant loop and the nucleotides U12, C13, and A14. Upon comparison with nonspecific sequence RNA binding proteins, the sequence-specific binding KH3 domain can form various important contacts, including hydrophobic contacts and hydrogen bonds, with the base moieties of RNA. All of these findings provide insight into the recognition patterns between the KH-type proteins and their target sequence-specific singlestranded RNAs.
■
Article
ASSOCIATED CONTENT
S Supporting Information *
Tables showing the interactions between RNA and KH3 for wild type and A14G mutant. This material is available free of charge via the Internet at http://pubs.acs.org. 12433
dx.doi.org/10.1021/jp5079289 | J. Phys. Chem. B 2014, 118, 12426−12434
The Journal of Physical Chemistry B
Article
(41) Tsai, C. J.; Ma, B.; Nussinov, R. Folding and binding cascades: shifts in energy landscapes. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 9970−9972. (42) Tsai, C. J.; Ma, B.; Sham, Y. Y.; et al. Structured disorder and conformational selection. Proteins 2001, 44, 418−427. (43) Weikl, T. R.; von Deuster, C. Selected-fit versus induced-fit protein binding: kinetic differences and mutational analysis. Proteins 2009, 75, 104−110. (44) Wlodarski, T.; Zagrovic, B. Conformational selection and induced fit mechanism underlie specificity in noncovalent interactions with ubiquitin. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 19346−19351. (45) Bucher, D.; Grant, B. J.; McCammon, J. A. Induced fit or conformational selection? The role of the semi-closed state in the maltose binding protein. Biochemistry 2011, 50, 10530−10539. (46) Anthis, N. J.; Doucleff, M.; Clore, G. M. Transient, sparsely populated compact states of apo and calcium-loaded calmodulin probed by paramagnetic relaxation enhancement: interplay of conformational selection and induced fit. J. Am. Chem. Soc. 2011, 133, 18966−18974. (47) Silva, D. A.; Bowman, G. R.; Sosa-Peinado, A.; et al. A role for both conformational selection and induced fit in ligand binding by the LAO protein. PLoS Comput. Biol. 2011, 7, e1002054. (48) Day, R. Ensemble versus single-molecule protein unfolding. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 13445−13450. (49) Albertini, A. A. V. Crystal structure of the rabies virus nucleoprotein-RNA complex. Science 2006, 313, 360−363. (50) Green, T. J. Structure of the vesicular stomatitis virus nucleoprotein-RNA complex. Science 2006, 313, 357−360. (51) Sengoku, T.; Nureki, O.; Nakamura, A.; et al. Structural basis for RNA unwinding by the DEAD-Box protein drosophila vasa. Cell 2006, 125, 287−300. (52) Tawar, R. G.; Duquerroy, S.; Vonrhein, C.; et al. Crystal structure of a nucleocapsid-like nucleoprotein-rna complex of respiratory syncytial virus. Science 2009, 326, 1279−1283. (53) Leulliot, N.; Varani, G. current topics in RNA-Protein recognition control of specificity and biological function through induced fit and conformational capture. Biochemistry 2001, 40, 7947− 7956. (54) Frauenfelder, H.; Sligar, S. G.; Wolynes, P. G. The energy landscapes and motions of proteins. Science 1991, 254, 1598−1603. (55) Ma, B.; Kumar, S.; Tsai, C.-J.; et al. Folding funnels and binding mechanisms. Protein Eng. 1999, 12, 713−720. (56) Tsai, C.-J.; Kumar, S.; Ma, B.; et al. Folding funnels, binding funnels, and protein function. Protein Sci. 1999, 8, 1181−1190. (57) Koshland, D. E. Application of a theory of enzyme specificity to protein synthesis. Proc. Natl. Acad. Sci. U. S. A. 1958, 44, 98−104. (58) James, L. C.; Tawfik, D. S. Conformational diversity and protein evolution − a 60-year-old hypothesis revisited. Trends Biochem. Sci. 2003, 28, 361−368. (59) Okazaki, K. i.; Takada, S. Dynamic energy landscape view of coupled binding and protein conformational change: Induced-fit versus population-shift mechanisms. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 11182−11187. (60) Wlodarski, T.; Zagrovic, B. Conformational selection and induced fit mechanism underlie specificity in noncovalent interactions with ubiquitin. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 19346−19351. (61) Frank, J.; Massey, J. The kolmogorov-smirnov test for goodness of fit. J. Am. Statist. Assoc. 1951, 46, 68−78. (62) Qin, F.; Chen, Y.; Wu, M.; et al. Induced fit or conformational selection for RNA/U1A folding. RNA 2010, 16, 1053−1061. (63) Ye, W.; Yang, J.; Yu, Q.; et al. Kink turn sRNA folding upon L7Ae binding using molecular dynamics simulations. Phys. Chem. Chem. Phys. 2013, 15, 18510−18522.
(15) Westhof, E.; Fritsch, V. The Endless Subtleties of RNA-Protein Complexes. Structure 2011, 19, 902−903. (16) Onoa, B.; Tinoco, I. RNA folding and unfolding. Curr. Opin. Struct. Biol. 2004, 14, 374−379. (17) Henkels, C. H.; Oas, T. G. Ligation-State Hydrogen Exchange: Coupled Binding and Folding Equilibria in Ribonuclease P Protein. J. Am. Chem. Soc. 2006, 128, 7772−7781. (18) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; et al. AMBER 11; University of California: San Francisco, CA, 2010. (19) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; et al. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926. (20) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089. (21) Hornak, V.; Abel, R.; Okur, A.; et al. Comparison of multiple amber force fields and development of improved protein backbone parameters. Proteins: Struct., Funct., Bioinf. 2006, 65, 712−725. (22) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; et al. Improved sidechain torsion potentials for the Amber ff99SB protein force field. Proteins: Struct., Funct., Bioinf. 2010, NA−NA. (23) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327−341. (24) Pearson, K. On Lines and Planes of Closest Fit to Systems of Points in Space. Philos. Mag. 1901, 2, 559−572. (25) Amadei, A.; Linssen, A. B. M.; Berendsen, H. J. C. Essential dynamics of proteins. Proteins: Struct., Funct., Genet. 1993, 17, 412− 425. (26) Hayward, S.; Groot, B. L. d. Normal Modes and Essential Dynamics. Methods Mol. Biol. 2008, 443, 89−106. (27) Kitao, A.; Go, N. Investigating protein dynamics in collective coordinate space. Curr. Opin. Struct. Biol. 1999, 9, 164−169. (28) Berendsen, H. J.; Hayward, S. Collective protein dynamics in relation to function. Curr. Opin. Struct. Biol. 2000, 10, 165−169. (29) Miller, B. R., III; McGee, T. D., Jr.; Swails, J. M.; et al. MMPBSA. py: An efficient program for end-state free energy calculations. J. Chem. Theory Comput. 2012, 8, 3314−3321. (30) Chen, H.; Luo, R. Binding induced folding in p53-MDM2 complex. J. Am. Chem. Soc. 2007, 129, 2930−2937. (31) Yu, Q.; Ye, W.; Wang, W.; et al. Global Conformational Selection and Local Induced Fit for the Recognition between Intrinsic Disordered p53 and CBP. PLoS One 2013, 8, e59627. (32) Wang, W.; Ye, W.; Jiang, C.; et al. New force field on modeling intrinsically disordered proteins. Chem. Biol. Drug Des. 2014, 84, 253− 269. (33) Sharp, K. A.; Honig, B. Electrostatic interactions in macromolecules theory and applications. Annu. Rev. Biophys. Biophys. Chem. 1990, 19, 301−332. (34) Koshland, D. E. Application of a Theory of Enzyme Specificity to Protein Synthesis. Proc. Natl. Acad. Sci. U. S. A. 1958, 44, 98−104. (35) Boehr, D. D.; Wright, P. E. Biochemistry. How do proteins interact? Science 2008, 320, 1429−1430. (36) Csermely, P.; Palotai, R.; Nussinov, R. Induced fit, conformational selection and independent dynamic segments: an extended view of binding events. Trends Biochem. Sci. 2010, 35, 539−546. (37) Kumar, S.; Ma, B.; Tsai, C. J.; et al. Folding and binding cascades: dynamic landscapes and population shifts. Protein Sci. 2000, 9, 10−19. (38) Ma, B.; Kumar, S.; Tsai, C. J.; et al. Folding funnels and binding mechanisms. Protein Eng. 1999, 12, 713−720. (39) Ma, B.; Shatsky, M.; Wolfson, H. J.; et al. Multiple diverse ligands binding at a single protein site: a matter of pre-existing populations. Protein Sci. 2002, 11, 184−197. (40) Tsai, C. J.; Kumar, S.; Ma, B.; et al. Folding funnels, binding funnels, and protein function. Protein Sci. 1999, 8, 1181−1190. 12434
dx.doi.org/10.1021/jp5079289 | J. Phys. Chem. B 2014, 118, 12426−12434