- 1Key Laboratory for Molecular Enzymology and Engineering of Ministry of Education, School of Life Science, Jilin University, Changchun, China
- 2State Key Laboratory of Theoretical and Computational Chemistry, Institute of Theoretical Chemistry, Jilin University, Changchun, China
- 3Department of Electric Engineering and Computer Science, C.S. Bond Life Sciences Center, University of Missouri, Columbia, MO, United States
- 4College of Computer Science and Technology, Jilin University, Changchun, China
Thermus thermophiles Argonaute (TtAgo) is a complex, which is consisted of 5′-phosphorylated guide DNA and a series of target DNA with catalytic activities at high temperatures. To understand why high temperatures are needed for the catalytic activities, three molecular dynamics simulations and binding free energy calculations at 310, 324, and 338K were performed for the TtAgo-DNA complex to explore the conformational changes between 16-mer guide DNA/15-mer target DNA and TtAgo at different temperatures. The simulation results indicate that a collapse of a small β-strand (residues 507–509) at 310 K caused Glu512 to move away from the catalytic residues Asp546 and Asp478, resulting in a decrease in catalytic activity, which was not observed in the simulations at 324 and 338 K. The nucleic acid binding channel became enlarged at 324 and 338K, thereby facilitating the DNA to slide in. Binding free energy calculations and hydrogen bond occupancy indicated that the interaction between TtAgo and the DNA was more stable at 324K and 338K than at 310 K. The DNA binding pocket residues Lys575 and Asn590 became less solvent accessible at 324 and 338K than at 310 K to influence hydrophilic interaction with DNA. Our simulation studies shed some light on the mechanism of TtAgo and explained why a high temperature was needed by TtAgo during gene editing of CRISPR.
Introduction
Argonaute (Ago) proteins preform a critical role in guide target RNA recognition, cleavage, and product release, which depend on key components of a RNA-induced silencing complex (Peters and Meister, 2007; Hutvagner and Simard, 2008; Kawamata and Tomari, 2010; Sheng et al., 2014). These proteins belong to the TNRC6/GW182 protein family, coordinates downstream silencing events, which also named TNRC6A-C and GW proteins in humans. Ago proteins are provided with four domains with distinct functions (Huntzinger and Izaurralde, 2011; Pfaff and Meister, 2013; Chen et al., 2014; Mathys et al., 2014; Hauptmann et al., 2015). The N-terminal domain plays a key role for small RNA binding. The P element-induced wimpy testis (PIWI)–Argonaute–Zwille (PAZ) domain is essential for recruiting the middle domain (MID) and the 5′ end of a small RNA. The PIWI domain has the same function as the RNase H domain, and some Ago proteins can serve as small RNA-guided endonucleases (Hauptmann et al., 2015). The PIWI domain is essential for the cleavage activity of Ago in which the Asp–Asp–Asp/His catalytic triad of Ago processes the identification of double stranded RNAs and the cleaving of their strand (Liu et al., 2004; Song et al., 2004; Ma et al., 2005; Rivas et al., 2005; Yuan et al., 2005; Parker and Roe, 2014). It also acts on guide-target RNA duplexes to cleave their target strand. Mg2+ cation mediates the endonucleolytic cleavage of a target RNA strand (Martinez and Tuschl, 2004; Schwarz et al., 2004; Lingel and Sattler, 2005; Jinek and Doudna, 2009; Parker, 2010), thereby forming 3′-OH and 5′-phosphate ends (Elbashir et al., 2001). A series of experimental structures of Thermus thermophilus Argonaute (TtAgo) binding to various complexes of different nucleic acid strands provide insights into conformational changes in proteins and DNAs (Jung et al., 2013; Zhu et al., 2016a,b; Sheng et al., 2017). Two structures, namely, PDB ID 3F73 and PDB ID 3HM9, whose substrates have different lengths (Wang et al., 2008a,b), suggested that the PAZ domain undergo a series of motions during a catalytic cycle. They also reveal the PAZ domain disassociates from the 3′ end of the guide upon target binding (Xia et al., 2012; Wang et al., 2013; Nam et al., 2014; Swarts et al., 2014; Jiang et al., 2015).
This study focused on TtAgo. Figure 1 shows the crystal structure of TtAgo (PDB ID 4NCB) with bound 16-mer guide DNA and 15-mer target DNA (Sheng et al., 2014). TtAgo is also made up by four functional domains, similar to other Ago proteins, which named N, PAZ, MID and C-terminal PIWI domain. These four domains are connected by two linker domains L1 (Linker 1) and Linker 2 (L2) (Figures 1A–F). MID bind to 5′ terminals of the guide and PAZ bind to 3′ terminals of the guide, both of them define the nucleic acid binding channel (Figure 1C and Figure S1). Driven mainly by positively charged residues of TtAgo and negatively charged DNA backbone atoms (Zander et al., 2014; Swarts et al., 2015; Figure S1A), the nucleotides 2–8 of the guide of 5′ segment, which called seed region, involved the combination of the channel binding (Figures S1B,C). The correct positioning of the seed region is important for the binding between the target nucleic acid and TtAgo and human Ago (Lee et al., 2016; Khin et al., 2017). The distance between Glu512 and active center (Asp478, Asp546, and Asp660) directly affects the activity of TtAgo. When Glu512 is far from the catalytic pocket (Figure S1D; Wang et al., 2008a,b), it turns out to be a cleavage-incompatible conformation. Glu512 positioned on loop L2, close to the catalytic pocket, represents a key residue for a cleavage-compatible conformation (Figure 1G and Figure S1E).
Figure 1. Argonaute protein of TtAgo. (A) Domain composition (left) and nucleic acid binding channel (right) of TtAgo (PDB 4NCB). (B) Sequence and pairing of guide (blue) and target (red) strands in the ternary complex of 4NCB. (C–E) The 3-D structures of TtAgo ternary complex were colored as (A,B). In (C), the TtAgo protein is shown as a surface, whereas the nucleic acid is shown as a cartoon. In (D,E), the TtAgo protein is shown as a cartoon, whereas the nucleic acid is shown as a stick. (G) Catalytic residues of TtAgo. The catalytic residues in a stick representation were highlighted by pink sticks and a cyan mesh.
As CRISPR/Cas9 becomes a popular molecular biology technology for gene editing, there is a continuous interest in searching for new genome-editing methods. Recently, it has been reported that argonaute proteins prossess DNA-guided sequence-specific DNA endonuclease activity in Pyrococcus furiosus, Thermus thermophilus, and Methanocaldococcus jannaschii (Beshnova et al., 2014; Cherstvy and Teif, 2014; Harikrishna and Pradeepkumar, 2017) and it can be as a host defense system to eliminate invading nucleic acids. These argonaute proteins usually needs high temperature to be fully active, which may restrict the application in many aspects. In 2016, the argonaute protein (NgAgo) from Natronobacterium gregoryi was reported to as mesophilic microbe whose enzymes, which is functional at 310K. Compared to Cas9, it have higher genome editing activity in mammalian cells (Chandradoss et al., 2015). Since NgAgo does not require a PAM motif and the authors reported the PAM motif is not essential for NgAgo and high efficiency to targets of high guanine-cytosine content, NgAgo was once considered promising as a useful tool for genome editing. However, as soon as researchers set off on a whirlwind of experiments (Chandradoss et al., 2015; Salomon et al., 2015), questions started to emerge and the hope was quickly faded. Nevertheless, it is unclear why NgAgo cannot perform the function of TtAgo, although the two proteins have significant sequence similarity, and whether it is possible to use or modify argonaute as a target for gene editing.
It has been reported that TtAgo is active at high temperatures, not at low temperatures. But how temperature affects TtAgo activity, it has not been clarified so far. A few molecular dynamics (MD) studies on siRNA–TtAgo and siRNA–hAGO2 complexes have been conducted (Wang et al., 2013; Nam et al., 2014; Hanlun et al., 2015). Jiang et al. constructed from MD simulations and found that MSMs can elucidate the conformational dynamics of AGO at biologically relevant timescales (Hanlun et al., 2015). Jinhyuk Lee pointed out that they revealed different structural conformations of the RNA duplex in Ago2, depending on Mg2+concentration and demonstrated that cation effects on Ago2 structural flexibility are critical to its catalytic/functional activity, with low Mg2+ favoring greater Ago2 flexibility and less miRNA/mRNA duplex stability, thus favoring slicing (Nam et al., 2014). Wang et al. discovered three findings during MD simulations: (1) three important (PAZ, Mid and PIWI) domains existed in Argonaute which defined the global dynamics of the protein; (2) the interdomain correlated movements were so crucial for the interaction of Ago-RNAs that they not only facilitated the relaxation of the interactions between residues surrounding the RNA binding channel but also induced certain conformational changes; and (3) it was just these conformational changes that expand the cavity of the active site and open putative pathways for both the substrate uptake and product release (Wang et al., 2013). However, several fundamental questions remain unanswered: (1) How does the temperature affect the catalytic activity of Ago? (2) What factors cause the conformational changes in the catalytic center at different temperatures?
In the present work, MD simulations and binding free energy calculations (molecular mechanics–Poisson–Boltzmann surface area, MM-PBSA) were performed on 16-mer guide DNA, 15-mer target DNA, and TtAgo at temperatures of 310, 324, and 338K to find possible answers to these questions. Principal component analysis (PCA) was also applied to investigate the interactions between TtAgo and DNA. Our MD results provided some informative clues for the mechanism of conformational changes in the catalytic center at different temperatures.
Materials and Methods
Protein Preparation
In 2014, Sheng et al determined the crystal structures of ternary Thermus thermophilus Argonaute (TtAgo) complexes with 5′-phosphorylated guide DNA and a series of DNA targets (PDB Id 4NCB) (Sheng et al., 2014), in which the complex structures solved at improved resolution up to 2.2 Å. It was very interesting reliable target for MD simulation. And so we chose TtAgo (PDB Id 4NCB) to perform simulations. The initial structure of TtAgo–nucleic acid complex in the MD simulations was obtained from the Protein Data Bank (PDB ID 4NCB at 2.19 A resolution). TtAgo is a 685-amino acid protein, and its binding nucleic acid is composed of 16-mer guide DNA and 15-mer target DNA. The Mg2+ in the 4NCB structure were also used in the simulations.
Conventional MD Simulations and Analysis
The three systems with different temperatures (310, 324, and 338K) were utilized in 200 ns MD simulations with NAMD2.10 b1 by using CHARMM27 all-force field parameters (Best et al., 2012). In TtAgo, for Mg2+ and nucleic acids (guide DNA and a series of DNA targets), generalized CHARMM27 all-force field parameters (Mackerell et al., 2004) were applied. In three system, the complexes were placed in the explicit TIP3P water model (Jorgensen et al., 1998), and extended at 15 Å from water molecules in a cubic periodic box. The systems were added with sodium and chloride ions to obtain a final ion concentration of 0.15 mol/L. To avoid steric clashes or improper geometries, the steepest descent algorithm with 50,000 steps was carried out to was energetically minimized for each system. The Langevin dynamics (Schlick, 2002) was used to constant temperature control with a damping coefficient (gamma) of 1.0 ps. The long-range electrostatic interactions was calculated by Particle Mesh Ewald summation algorithm (Darden et al., 1993). Each simulation was performed for 200 ns for the protein–nucleic acid complexes at different temperatures (310, 324, and 338K) by using a pressure coupling time constant of 2.0 ps and a temperature coupling time constant of 0.1 ps and a pressure coupling time constant of 2.0 ps. The timestep for each simulation was 1ps.
Cluster Analysis
Clusteringis a general data-mining technique (Cormack, 1971; Jain et al., 1999) using different methods (algorithms) that can be applied to judge structure similarity. In this study, clustering was performed on three systems (the 310K, 324K, and 338K) molecular dynamics simulation trajectories using the RMSD as a metric according to hierarchical method (Rokach and Maimon, 2005) into distinct groups. The cut off value of RMSD was set with 2.0Å for each systems. Then the TtAgo structure in the 200 ns simulation trajectories which had the similar conformation were divided into the same cluster. The proportion of each cluster were calculated and representative frame of the two largest clusters were extracted for the further analysis.
Binding Free Energy Calculation
The Equation (1) was used to calculate the binding free energy as follows:
where ΔGb is the total binding free energy; The van der Waals interaction energies and electrostatic between a protein and its ligand are Δ and Δ, respectively; TΔS is the contribution of conformational entropy to the binding; ΔGsol is the solvation energy, and. Δ and Δ were calculated by the same parameter set as in the molecular dynamics simulation, but no cutoff was applied to the calculation.
The binding free energies of nucleic acids and TtAgo in different temperature systems was performed by MM-PBSA calculation. The Amber 14 package was used to perform 2 ns simulation with the Amber ff99 force field parameters (Kollman et al., 2000). In our study, 500 snapshots were extracted from the last 1.5 ns MD trajectory for all of the energy components calculating.
Results and Discussion
Stability of TtAgo in Different Temperature Systems
To analyze the influence of temperature on TtAgo activity, the 200 ns molecular dynamics simulation prepared at three temperatures (310, 324, and 338K) were prepared and the simulation of each temperature was repeated three times. The root mean square deviation (RMSD) was monitored with the three different temperature systems of MD trajectories to measure the stability of each domain of TtAgo with temperature variation and to examine the protein stability during the MD simulations. The RMSD of protein backbone atoms and different domain backbone atoms was conducted (Figure 2 and Figure S2). In Figure S2, the RMSD variations of the different temperature systems for three times at the 200 ns time scale indicated that the complexes had the similar atomic coordinates and the initial structures. The RMSDs of the 310, 324, and 338K systems converged to 1.80, 2.40, and 2.20 Å, respectively, indicating that the three systems were stable. The RMSDs of different domains in the three systems are shown in Figures 2A–F. The Linker 1 domain of TtAgo became more flexible in the 324K system (Figure 2B) than that of the 310 and 338 K systems, whereas the PIWI domain was more flexible in the 338K system than in the two other systems (Figure 2F). The increase in temperature causes the protein residues to fluctuate dramatically. In this study, rising the temperature to 324K and 338K resulted in an increase of flexible for Linker 1 or PIWI domains. This may affect the relative position between the Linker 1 and PIWI and thus regulated the catalytic activity of TtAgo. In Figure S3, the total energy vs. time was also calculated and it can be seen that the total energy for three systems keep steady during 200 ns simulation. In addition, We chose the four structures with lower energy from three simulation trajectories and calculated the Ramachandran Plot. As the Figures S4–S6. The four lower energy structures for each system had the same Ramachandran conformation which indicated that the whole simulation process was stable.
Figure 2. Comparison of the RMSD of the 310K, 324K, and 338K systems at different TtAgo domains. (A–F) RMSD plot of backbone atoms in different domains in three temperature systems.
Then the root-mean-square fluctuation (RMSF) for the backbone atoms of the nucleic acids with TtAgo protein in different temperature systems was analyzed to determine the mobility of the complexes. Figure 3 illustrates that the large punctuations of residues mainly occurred in the 324 and 338K systems. The RMSF fluctuations of Gly38, Asp62, and Gly83 located at the N terminus (Figure 3A), Arg123 and Arg137 located at Linker1 (Figure 3B), Pro218 located at PAZ (Figure 3C), Ile297 and Arg324 located at Linker 2 (Figure 3D), Gln387 located at MID (Figure 3E), and Asp478, Glu507, Gly547, and Gly612 located at PIWI (Figure 3F) were higher in the 324 and 338K systems than in the 310 K system. The results indicated that the protein at 324 and 338K exhibited a higher mobility than at 310 K during MD simulations. Asp478, Glu507, and Gly547 near the active center at the PIWI domain that exhibited a high mobility in the 324 and 338K systems might affect the catalytic activity of TtAgo. To further confirm the RMSF results, we calculated the RMSD values for each residue, and then colored each according to the RMSD value (Figures 3G–I). It can be seen that the intensity of the movement of various residues in the 324 or 338K system was higher than that in the 310K system, especially in the Linker 2 and active pocket domain. In order to further analyze residual atomic flexibility, an isotropic temperature factor (B-factor) has been calculated and the results was shown in Figure S7. As we can see, the change trend of B-factor was consistent with that of RMSF and RMSD of individual residues. Figure S8 shows that the radius of gyration (Rg) in different domains in the three systems during the 200 ns time scale MD was stable.
Figure 3. RMSF and backbone flexibility during 200 ns MD at different TtAgo domains. (A) RMSF values of the N terminus calculated in three temperature systems for 200 ns MD trajectories. (B) RMSF values of Linker 1 calculated in three temperature systems for 200 ns MD trajectories. (C) RMSF values of PAZ calculated in three temperature systems for 200 ns MD trajectories. (D) RMSF values of Linker 2 calculated in three temperature systems for 200 ns MD trajectories. (E) RMSF values of MID calculated in three temperature systems for 200 ns MD trajectories. (F) RMSF values of PIWI calculated in three temperature systems for 200 ns MD trajectories. (G–I) Visualizations of backbone flexibility of TgAgo at 310 K (G), 324K (H), and 338 K (I). The different degrees of motion correspond to tubes with different colors and thicknesses.
The solvent-accessible surface area (SASA) predicts the number of residues present in the outlier regions (surface) of the protein and the number of residues that are in the hydrophobic core (buried). In this study, the time-dependent SASA in the different domains of TtAgo from the simulations was also calculated (Figure S9). The results showed that the ensembles of the different domains in the three temperature systems were similar, which indicated that the change in temperature hardly influenced the SASA value of each domain. This analysis further calculated the SASA values of the residues that interact with the guide DNA and target DNA. As we can see in Figure S10A, the SASA of the hydrophilic residues Lys575 and Asp590 decreased in the 324 and 338K systems compared with that in the 310K system, thereby illustrating increased interaction with the DNA. Asp590 interacted with the G1 of the target DNA (Figure S10B), whereas Lys575 interacted with the C11 and C12 of the target DNA (Figure S10C), forming a tighter interaction with DNA. This result indicated that TtAgo at 324 and 338K might have undergone drifting, thereby leading to its transition into an advantageous conformation and binding with a nucleic acid.
Conformational Changes in Different Temperature Systems
It has been reported that the binary complex of a yeast Ago with a bound 5′-phosphorylated guide RNA strand have found that the distance between Glu512 and catalytic pocket composed of three Asp residues, namely, Asp478, Asp546, and Asp660, determines the activity of Ago (Nakanishi et al., 2012). The activated Ago needs to insert Glu512 into the catalytic pocket. The change in the distance between Glu512 and the catalytic pocket during MD simulations in the three temperature systems was also analyzed to observe the influence of temperature change on enzyme catalytic activity. In Figure 4A, the distance between Glu512 and Asp478 in the 310K system increased from 10 to 21 Å at 80 ns. The distance between Glu512 and Asp546 in the 310K system also increased from 9.50 to 23 Å at 80 ns (Figure 4B). In the 324 and 338K systems, the distance from Glu512 to Asp546 and Asp478 was approximately 10 Å during MD simulations. In addition, the distance between Glu512 and Asp478 or Asp546 for other two repeated simulations were calculated and the results were prepared in Figure S11. It can be seen that the distance between Glu512 and Asp478 or Asp546 showed a remarkably increase at 47 and 90 ns in 310K system respectively. But for the 324 and 338K systems, the distance from Glu512 to Asp546 and Asp478 was approximately 10 Å during two MD simulations. This finding implied that the distance from the Glu512 to the active center in the 310K system was larger than that in the 324 and 328K systems. To confirm these results, we compared the structures between 0 and 150 ns at different temperatures (Figures 4C–H).The distance from Glu512 OE2 to Asp546 OD1 and Asp478 OD1 at 150 ns was 2.5 times than at 0 ns in the 310K system (Figures 4C,D). The distance between 0 and 150 ns changed slightly in the 324 and 338K systems compared with that in the 310 K system (Figures 4E, 5H). To illustrate the dynamic process of distance change, we created animations 1–3 in Supporting information of the distance between the Glu512 OE2 and Asp478 OD1 of the visual dynamic variation at three different temperatures. The results revealed that Glu512 was far from the catalytic pocket at 310K, and the enzyme activity at 324 and 338K was stronger than that at 310K.
Figure 4. Distance between Glu512 and active pocket in 310K, 324K, and 338K systems. (A) Time evolution of distance between Glu512 OE2 and Asp478 OD1 in three temperature systems. (B) Time evolution of distance between Glu512 OE1 and Asp546 OD1 in three temperature systems. (C,D) Distance between Glu512 (blue) and active pocket (Asp546, cyan; Asp478, yellow) in 310K system at 0 and 150 ns. (E,F) Value of distance between Glu512 (blue) and active pocket (Asp546, cyan; Asp478, yellow) in 324K system at 0 and 150 ns. (G,H) Distance between Glu512 (blue) and active pocket (Asp546, cyan; Asp478, yellow) in 338K system at 0 and 150 ns.
Figure 5. Comparison of secondary structure in different temperature systems. (A,B) Structure of TtAgo in 310K system at 1.5 ns (A) and 150 ns (B) during 200 ns MD, and residues Glu507 to Arg513 were highlighted with orange (1.5 ns) and green (150 ns) cartoons. (C,D) Structure of TtAgo in 324K system at 1.5 (C) and 150 ns (D) during 200 ns MD, and residues Glu507 to Arg513 were highlighted with split-pea (1.5 ns) and yellow (150 ns) cartoons. (E,F) Structure of TtAgo in 324K system at 1.5 (E) and 150 ns (F) during 200 ns MD, and residues Glu507 to Arg513 were highlighted with cyan (1.5 ns) and yellow (150 ns) cartoons.
To explain the cause of the change in distance, we investigated the secondary structure modifications during 200 ns MD simulations in the different temperature systems (Figures S12A–C and Figure 5; Table 1). The β-strand of Glu507 to Gln509 near Glu512 disappeared at 150 ns compared with that at1.5 ns in the 310K system (Figure S12A). The β-strand contents of Glu507 to Gln509 in the 310, 324, and 338K systems were 1.50, 37.60, and 46.80%, respectively (Table 1, the β-strand contents of Glu507 to Gln509 in the 310, 324, and 338K systems for other two repeated simulations were calculated and the results were prepared in Tables S1, S2). Secondary structure analysis revealed that the short intermediate β-sheet at 324 and 338K during 200 ns molecular dynamics simulation represented the ordered structure that might help maintain the close distance between Glu512 and the catalytic pocket. Glu512 moved far from the catalytic pocket when β-sheet disappeared.
Changes in Nucleic Acid Binding Channel in Different Temperature Systems
The minimum distances from Linker 1 to PIWI and from PAZ to PIWI have been reported to determine the nucleic acid binding channel in TtAgo. In our study, DCH1 was used to describe the center of mass (c.o.m.) distance from the part of L1 (residues 165–174, highlighted in red) to the long helix on PIWI (residues 654–668, highlighted in blue). DCH2 was the c.o.m. distance from the part of PAZ (residues 192–202, highlighted in pink) to the part of PIWI (residues 572–581, highlighted in green), as shown in Figure 6A. These two distances were used to measure the diameter of the nucleic acid binding channel. To analyze the effect of temperature change on the nucleic acid binding channel, we detected the changes in DCH1 and DCH2 during 200ns MD in the different temperature systems (Figures 6B,C). DCH1 varied at 27.5, 28.5, and 29.5 Å in the 310, 324, and 338K systems, respectively. These findings suggested that the minimum distance from Linker 1 to PIWI was larger in the 324 and 338K systems than that in the 310K system. Similarly, DCH2 reached 31 and 30 Å in the 324 and 338K systems and in the 310K system, respectively.
Figure 6. Analysis of nucleic acid binding channel in three temperature systems. (A) DCH1 and DCH2 were used to describe the diameter of the nucleic acid binding channel. DCH1 measured the c.o.m. distance in residues 165–174 (red cartoon) and in residues 654–668 (blue cartoon). DCH2 was the c.o.m. distance in residues 192–202 (orange cartoon) and in residues 572–581 (green cartoon). The guide-DNA and target-DNA were highlighted by white and marine cartoon respectively. (B) Time series of DCH1 in three temperature systems (left) and the probability distribution of DCH2 calculated from 0 to 200 ns trajectories (right). (C) Time series of the value for DCH2 in three temperature systems (left) and the probability distribution of DCH2 calculated from 0 to 200 ns trajectories (right). (D–F) 2-D(top) and 3-D (bottom) visualization of the channel obtained with CHEXVIS in the 310K (D), 324K (E), and 338K systems (F) for 150 ns structure. For 2-D channel profile visualization, bulkiness and average flexibility are shown in split visualization over the radius profile as a color map. The bottleneck radius was highlighted by rectangles. For 3-D channel profile visualization, the channel was highlighted with blue, and TtAgo was represented by a cartoon.
CHEXVIS is a tool used for molecular channel calculation and visualization based on alpha complex representation (Masood et al., 2015). In addition to predictive channels, channel residues and their flexibility can be predicted by CHEXVIS. In our study, this tool was utilized to analyze the distinction of the nucleic acid binding channels in the different temperature systems at 150 ns for 200 ns MD (Figures 6D–F; Table 2). The length of the channel in TtAgo in the 310K system (60.1Å) was longer than those in the 324 (50.37Å) and 338K (50.76Å) systems. Studies on bottlenecks have recommended promising areas for the modification of channel properties because their substitutions can affect the informative tunnel geometry when protein dynamics is considered. The results showed that the bottleneck of TtAgo in the 310K system was 3.452 Å and changed to 4.375 and 4.362 Å in the 324 and 338K systems, respectively. The analysis of bulkiness and average flexibility of channel residues, which were presented in a 2D profile, suggested that the flexibility of the channel residues was higher in the 324 and 338K systems than in the 310K system. With such a high flexibility, the channels easily become wide. These results indicated that the diameter of nucleic acid binding channels were larger and the binding of nucleic acid and protein may occur more easily in the 324 and 338K systems than in the 310K system as temperature increased. Our theoretical finding is consistent with the computational study made by Wang et al. (2013).
Cluster Analysis
In addition, clusteranalysis was performed on each trajectory of the three systems to reveal a clear-cut structural difference in the nucleic acid and TtAgo complex in the different temperature systems. As shown in Figure 7A, conformations for 310K systems during the simulations were divided into six clusters and the percentage population of clusters was 43.1, 34.3, 9.2, 6.7, 4.3, 8, 2.4, 5, and 4% respectively. Then the representative structures of most populated cluster in cluster1 and cluster2 were selected to analyze the conformation and channel. It can be seen that the residues Glu507 to Gln509 in both the representative structures of cluster1 and cluster2 at 310K formed a loop. The DCH1 and DCH2 in cluster1 representative structures were 27.43 and 29.37Å while they were 27.28 and 29.69Å in cluster2 representative structures, respectively.
Figure 7. Analysis on clusters in the 310K (A), 324K (B), and 338K systems (C). Two representative structures of most populated clusters in three systems are represented as cartoons. The DCH1 and DCH2 residues were highlighted with mesh. The beta strands at residues Glu507 to Arg513 were colored in red, while the loops at residues Glu507 to Arg513 were colored by orange.
Figures 7B,C exhibited the cluster analysis for 324K and 338K systems, respectively. The percentage population of cluster1 and cluster2 for 324K and 338K systems were 46.1, 27.3, 53.2, and 21.4%, respectively. The sum of cluster1 and cluster2 both in those two systems were over 70%. Two representative structures based clustering approach at 324K (cluster1 and cluster2) were extracted the results were shown in Figure 7B. In those conformations, the residues Glu507 to Gln509 remained β-sheet and the DCH1 were over 29.56Å in two representative structures while the and DCH2 were over 30 Å. As the same as 324K systems, residues Glu507 to Gln509 in cluster1 and cluster2 of 338K systems for representative structures also presented β-sheet. The DCH1 and DCH2 at 338K were more than 29 Å (Figure 7C). These findings revealed that the β-sheet (Glu507 to Gln509) at 310K formed a loop but remained in its original state in the 324 and 338K systems. These results were consistent with those illustrated in Figure 6. The different nucleic acid binding channels in the three temperature systems in the cluster1 and cluster2 of TtAgo illustrated that DCH1 and DCH2 in the 310K system were lower than those in the 324 and 338K systems.
Hydrogen Bond Network Analysis
The hydrogen bond occupancy between nucleic acids (target DNA and guide DNA) and TtAgo during the 200 ns MD is listed in Tables S3, S4 to obtain insights into the molecular interactions between TtAgo and nucleic acid at different temperatures. The MD simulation results showed that the number of hydrogen bonds between nucleic acids and TtAgo in the 310K system was more than that in the 324 and 338K systems. For example, the hydrogen bonds, R59 (HN2)—DC5 (O2), D590 (OD1)—DG19 (N1) that connected the protein and the target DNA, and R615 (NE)—DT6 (O5′) and Y457 (NZ)—DA3 (O2P) that linked TtAgo and guide DNA disappeared in the 324 and 338K systems. The decrease in hydrogen bond occupancy was due to the increase in the distance between hydrogen bond donors and receptors. The distances of the donors and receptors were examined during 200 ns MD (Figure 8). In Figure 8, the distances from the donors to the receptors were larger in the 324 and 338K systems than in the 310K system. The rise in temperature increased the entropy, thereby promoting the relative motion of proteins and DNA. Such an increase reduced the stability of hydrogen bonds. Hydrophobic reduction greatly increased the solubility of TtAgo and nucleic acids in solvents, resulting in TtAgo was easier to complete the cutting action of nucleic acid. Overall, the binding of nucleic acid and TtAgo became tighter at higher temperature and TtAgo has high activity at high temperature.
Figure 8. Distance between hydrogen bond donors and receptors at different temperatures. (A) Distance from R59 (HN2) to DC5 (O2) in the three systems and the probability exhibited in (A1). (B) Distance from D590 (OD1) to DG19 (N1) in the three systems and the probability exhibited in (B1).(C) Distance from R615 (NE) to DT6 (O5) in the three systems and the probability exhibited in (C1). (D) Distance from Y457 (NZ) to DA3 (O2P) in the three systems and the probability exhibited in (D1).
MM-PBSA Calculations
It was well known that the arguably most popular end point method to calculate the free energy of binding between the enzyme and the ligands is the topic of MM/PBSA molecular mechanics [MM] with–Boltzmann [PB] and surface area solvation (Guimarães, 2011). The method has been used in a range of settings, including protein design by free energy calculations. The accuracy is better than for docking and scoring methods (Rifai et al., 2018), comparable to other end point methods, such as LIE methods. Furthermore, it is very fast and thence it was chosen to used for calculated the free energy of binding. The binding free energy of the TtAgo–nucleic acid complex was evaluated using MM–GBSA to confirm the results. The non-bonded van der Waals (ΔEvdW), non-bonded electrostatic (ΔEelec) interactions, and binding free energy (ΔGbin) of the TtAgo–nucleic complex in the three temperature systems are shown in Table 3. The binding free energies were driven primarily by electrostatic interactions and polar solvation, whereas the vdW interaction, polar solvation, and entropy basically remained the same. No significance was found in the variations in the binding affinity between nucleic acid and TtAgo in the different temperature systems. Table 3 shows that the binding free energy of the TtAgo–nucleic complex at 310K (−419.3667 kcal/mol) was higher than those at 324K (−461.1173 kcal/mol) and 338K (−434.8954 kcal/mol). (Seen from Tables S5, S6, it listed the other two repeated simulations results for MM-GBSA Calculations. The results were similar to the results listed in Table 3). These data indicated that the interaction between TtAgo and nucleic acid was more stable at 324 and 338K than at 310K.
To analyze the effects of individual amino acids on TtAgo and DNA binding, the per-residue free energy decomposition was performed to characterize and identify the crucial residues that mainly contributed to the binding free energy (Figure 9). Six amino acid residues, namely, Tyr171, Arg172, Ile173, Arg174, Thr174, and Met413 had relatively low energies in the three temperature systems compared with other residues and might play a pivotal role in combining with guide DNA (Figure 9A). It can be seen that Asp546, Lys575, and Asp590 had lower energy which instructed the important role in the process of the combination of target DNA.
Figure 9. Free energy decomposition analysis at three different temperatures (310 K, 324 K, and 338 K) for TtAgo–nucleic complex. (A) The total energy for residues which involved in the interaction with guide DNA. (B) The total energy for residues which involved in the interaction with target DNA. The residues that contributed significantly were highlighted.
Conclusion
TtAgo is active at high temperatures. The mechanism by which the temperature affects the activity of TtAgo remains poorly understood. In the present study, MD simulations revealed that the disappearance of β-sheet from Glu507 to Gln509 at 310K could cause Glu512 to move far from the catalytic pocket, thereby reducing the activity of TtAgo. The binding channel of the nucleic acid widened at 324 and 338K. Consequently, the binding of nucleic acids and proteins occurred more easily at 324 and 338 K than at 310 K. Binding free energy (MM-PBSA) calculations and hydrogen bond occupancy showed that the binding between TtAgo and nucleic acid was more stable at 324and 338K than at 310 K. These results suggested that the temperature influenced the activity of TtAgo by regulating the binding channel width of TtAgo and the distance from Glu512 to the active center. The activity of TtAgo was higher at 324 and 338K than that at 310K. Our results would be useful as a basis for further studies on the influence of temperature on TtAgo and design of CRISPR methods.
Author Contributions
YL wrote and revised this paper. ZY prepared the tables. JZ made the Supplementary Materials. SW provided some revision advice. DX and WH provided the ideas and modified the papers.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
This work was supported by the Major Scientific Research Projects of Jilin Province (20140203025NY), the Natural Science Foundation of Jilin Province (3B812C201465), and the US National Institutes of Health (GM100701). Our simulations were conducted at the High Performance Computing Center of Jilin University.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2018.00223/full#supplementary-material
References
Beshnova, D. A., Cherstvy, A. G., Vainshtein, Y., and Teif, V. B. (2014). Regulation of the nucleosome repeat length in vivo by the DNA sequence, protein concentrations and long-range interactions. PLoS Comput. Biol. 10:e1003698. doi: 10.1371/journal.pcbi.1003698
Best, R. B., Zhu, X., Shim, J., Lopes, P. E., Mittal, J., Feig, M., et al. (2012). Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ1 and χ2 dihedral angles. J. Chem. Theory Comput. 8, 3257–3273. doi: 10.1021/ct300400x
Chandradoss, S. D., Schirle, N. T., Szczepaniak, M., Macrae, I. J., and Joo, C. (2015). A dynamic search process underlies MicroRNA targeting. Cell 162, 96–107. doi: 10.1016/j.cell.2015.06.032
Chen, Y., Boland, A., Kuzuoglu-Öztürk, D., Bawankar, P., Loh, B., Chang, C. T., et al. (2014). A DDX6-CNOT1 complex and W-binding pockets in CNOT9 reveal direct links between miRNA target recognition and silencing. Mol. Cell 54, 737–750. doi: 10.1016/j.molcel.2014.03.034
Cherstvy, A. G., and Teif, V. B. (2014). Electrostatic effect of H1-histone protein binding on nucleosome repeat length. Phys. Biol. 11:044001. doi: 10.1088/1478-3975/11/4/044001
Cormack, R. M. (1971). A review of classification. J. R. Stat. Soc. 134, 321–367. doi: 10.2307/2344237
Darden, T., York, D., and Pedersen, L. (1993). Particle mesh Ewald: an N·log(N) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. doi: 10.1063/1.464397
Elbashir, S. M., Harborth, J., Lendeckel, W., Yalcin, A., Weber, K., and Tuschl, T. (2001). Duplexes of 21-nucleotide RNAs mediate RNA interference in cultured mammalian cells. Nature 411, 494–498. doi: 10.1038/35078107
Guimarães, C. R. (2011). A direct comparison of the MM-GB/SA scoring procedure and free-energy perturbation calculations using carbonic anhydrase as a test case: strengths and pitfalls of each approach. J. Chem. Theory Comput. 7, 2296–2306. doi: 10.1021/ct200244p
Hanlun, J., Fu Kit, S., Lizhe, Z., Xin, G., Julie, B., and Xuhui, H. (2015). Mutations in the binding pocket for miRNA(5') show different effects on hAgo2-miRNA binding. PLos Comput. Biol. 11:e1004404.
Harikrishna, S., and Pradeepkumar, P. I. (2017). Probing the binding interactions between chemically modified siRNAs and human Argonaute 2 using microsecond MD simulations. J. Chem. Inf. Model. 57, 883. doi: 10.1021/acs.jcim.6b00773
Hauptmann, J., Schraivogel, D., Bruckmann, A., Manickavel, S., Jakob, L., Eichner, N., et al. (2015). Biochemical isolation of Argonaute protein complexes by Ago-APP. Proc. Natl. Acad. Sci. U.S.A. 112, 11841. doi: 10.1073/pnas.1506116112
Huntzinger, E., and Izaurralde, E. (2011). Huntzinger, E., Izaurralde EGene silencing by microRNAs: contributions of translational repression and mRNA decay. Nat. Rev. Genet. 12, 99–110. doi: 10.1038/nrg2936
Hutvagner, G., and Simard, M. J. (2008). Argonaute proteins: key players in RNA silencing. Nat. Rev. Mol. Cell Biol. 9, 22–32. doi: 10.1038/nrm2321
Jain, A. K., Murty, M. N., and Flynn, P. J. (1999). Data clustering: a review. ACM Comput. Surv. 31, 264–323. doi: 10.1145/331499.331504
Jiang, H., Zhu, L., Gao, X., Bernauer, J., and Huang, X. (2015). Markov state models reveal a two-step mechanism of miRNA loading into the human Argonaute protein: selective binding followed by structural re-arrangement. PLoS Comput. Biol. 11:e1004404. doi: 10.1371/journal.pcbi.1004404
Jinek, M., and Doudna, J. A. (2009). A three-dimensional view of the molecular machinery of RNA interference. Nature 457, 405–412. doi: 10.1038/nature07755
Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1998). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. doi: 10.1063/1.445869
Jung, S. R., Kim, E., Hwang, W., Shin, S., Song, J. J., and Hohng, S. (2013). Dynamic anchoring of the 3'-end of the guide strand controls the target dissociation of Argonaute-guide complex. J. Am. Chem. Soc. 135, 16865–16871. doi: 10.1021/ja403138d
Kawamata, T., and Tomari, Y. (2010). Making RISC. Trends Biochem. Sci. 35, 368–376. doi: 10.1016/j.tibs.2010.03.009
Khin, N. C., Lowe, J. L., Jensen, L. M., and Burgio, G. (2017). No evidence for genome editing in mouse zygotes and HEK293T human cell line using the DNA-guided Natronobacterium gregoryi Argonaute (NgAgo). PLoS ONE 12:e0178768. doi: 10.1371/journal.pone.0178768
Kollman, P. A., Massova, I., Reyes, C., Kuhn, B., Huo, S., Chong, L., et al. (2000). Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 33, 889–897. doi: 10.1021/ar000033j
Lee, S. H., Turchiano, G., Ata, H., Nowsheen, S., Romito, M., Lou, Z., et al. (2016). Failure to detect DNA-guided genome editing using Natronobacterium gregoryi Argonaute. Nat. Biotechnol. 35, 17. doi: 10.1038/nbt.3753
Lingel, A., and Sattler, M. (2005). Novel modes of protein-RNA recognition in the RNAi pathway. Curr. Opin. Struct. Biol. 15, 107–115. doi: 10.1016/j.sbi.2005.01.010
Liu, J., Carmell, M. A., Rivas, F. V., Marsden, C. G., Thomson, J. M., Song, J. J., et al. (2004). Argonaute2 is the catalytic engine of mammalian RNAi. Science 305, 1437–1441. doi: 10.1126/science.1102513
Ma, J. B., Yuan, Y. R., Meister, G., Pei, Y., Tuschl, T., and Patel, D. J. (2005). Structural basis for 5'-end-specific recognition of guide RNA by the A. fulgidus Piwi protein. Nature 434, 666–670. doi: 10.1038/nature03514
Mackerell, A. D. Jr., Feig, M., and Rd, B. C. (2004). Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 25, 1400–1415. doi: 10.1002/jcc.20065
Martinez, J., and Tuschl, T. (2004). RISC is a 5′ phosphomonoester-producing RNA endonuclease. Genes Dev. 18, 975–980. doi: 10.1101/gad.1187904
Masood, T. B., Sandhya, S., Chandra, N., and Natarajan, V. (2015). CHEXVIS: a tool for molecular channel extraction and visualization. BMC Bioinformatics 16:119. doi: 10.1186/s12859-015-0545-9
Mathys, H., Basquin, J., Ozgur, S., Czarnocki-Cieciura, M., Bonneau, F., Aartse, A., et al. (2014). Structural and biochemical insights to the role of the CCR4-NOT complex and DDX6 ATPase in microRNA repression. Mol. Cell 54, 751–765. doi: 10.1016/j.molcel.2014.03.036
Nakanishi, K., Weinberg, D. E., Bartel, D. P., and Patel, D. J. (2012). Structure of yeast Argonaute with guide RNA. Nature 486, 368–374. doi: 10.1038/nature11211
Nam, S., Ryu, H., Son, W. J., Kim, Y. H., Kim, K. T., Balch, C., et al. (2014). Mg2+ effect on Argonaute and RNA duplex by molecular dynamics and bioinformatics implications. PLoS ONE 9:e109745. doi: 10.1371/journal.pone.0109745
Parker, J. S. (2010). How to slice: snapshots of Argonaute in action. Silence 1:3. doi: 10.1186/1758-907X-1-3
Parker, J. S., and Roe, S. M. (2014). Crystal structure of a PIWI protein suggests mechanisms for siRNA recognition and slicer activity. Embo J. 23, 4727–4737. doi: 10.1038/sj.emboj.7600488
Peters, L., and Meister, G. (2007). Argonaute proteins: mediators of RNA silencing. Mol. Cell 26, 611. doi: 10.1016/j.molcel.2007.05.001
Pfaff, J., and Meister, G. (2013). Argonaute and GW182 proteins: an effective alliance in gene silencing. Biochem. Soc. Trans. 41, 855–860. doi: 10.1042/BST20130047
Rifai, E. A., Vermeulen, N. P. E., and Geerke, D. P. (2018). Binding free energy predictions of farnesoid X receptor (FXR) agonists using a linear interaction energy (LIE) approach with reliability estimation: application to the D3R Grand Challenge 2. J. Comput. Aided Mol. Des. 32, 239–249. doi: 10.1007/s10822-017-0055-0
Rivas, F. V., Tolia, N. H., Song, J. J., Aragon, J. P., Liu, J., Hannon, G. J., et al. (2005). Purified Argonaute2 and an siRNA form recombinant human RISC. Nat. Struct. Mol. Biol. 12, 340–349. doi: 10.1038/nsmb918
Rokach, L., and Maimon, O. (2005). Clustering methods. Data Mining Knowl. Discov. Handb. 3, 321–352. doi: 10.1007/0-387-25465-X_15
Salomon, W. E., Jolly, S. M., Moore, M. J., Zamore, P. D., and Serebrov, V. (2015). Single-molecule imaging reveals that Argonaute reshapes the binding properties of its nucleic acid guides. Cell 162, 84–95. doi: 10.1016/j.cell.2015.06.029
Schlick, T. (2002). Molecular modeling and simulation. Interdisc. Appl. Math. 26, 1–31. doi: 10.1007/978-0-387-22464-0
Schwarz, D. S., Tomari, Y., and Zamore, D. P. (2004). The RNA-induced silencing complex is a Mg-dependent endonuclease. Curr. Biol. 14, 787–791. doi: 10.1016/j.cub.2004.03.008
Sheng, G., Gogakos, T., Wang, J., Zhao, H., Serganov, A., Juranek, S., et al. (2017). Structure/cleavage-based insights into helical perturbations at bulge sites within T. thermophilus Argonaute silencing complexes. Nucleic Acids Res. 45, 9149–9163. doi: 10.1093/nar/gkx547
Sheng, G., Zhao, H., Wang, J., Rao, Y., Tian, W., Swarts, D. C., et al. (2014). Structure-based cleavage mechanism of Thermus thermophilus Argonaute DNA guide strand-mediated DNA target cleavage. Proc. Natl. Acad. Sci. U.S.A. 111, 652–657. doi: 10.1073/pnas.1321032111
Song, J. J., Smith, S. K., Hannon, G. J., and Joshua-Tor, L. (2004). Crystal structure of Argonaute and its implications for RISC slicer activity. Science 305, 1434–1437. doi: 10.1126/science.1102514
Swarts, D. C., Hegge, J. W., Hinojo, I., Shiimori, M., Ellis, M. A., Dumrongkulraksa, J., et al. (2015). Argonaute of the archaeon Pyrococcus furiosus is a DNA-guided nuclease that targets cognate DNA. Nucleic Acids Res. 43, 5120–5129. doi: 10.1093/nar/gkv415
Swarts, D. C., Jore, M. M., Westra, E. R., Zhu, Y., Janssen, J. H., Snijders, A. P., et al. (2014). DNA-guided DNA interference by a prokaryotic Argonaute. Nature 507, 258–261. doi: 10.1038/nature12971
Wang, Y., Juranek, S., Li, H., Sheng, G., Tuschl, T., and Patel, D. J. (2008a). Structure of an argonaute silencing complex with a seed-containing guide DNA and target RNA duplex. Nature 456, 921–926. doi: 10.1038/nature07666.
Wang, Y., Li, Y., Ma, Z., Yang, W., and Ai, C. (2013). Mechanism of microRNA-target interaction: molecular dynamics simulations and thermodynamics analysis. PLoS Comput. Biol. 6:e1000866. doi: 10.1371/journal.pcbi.1000866
Wang, Y., Sheng, G., Juranek, S., Tuschl, T., and Patel, D. J. (2008b). Structure of the guide-strand-containing argonaute silencing complex. Nature 456, 209–213. doi: 10.1038/nature07315
Xia, Z., Clark, P., Huynh, T., Loher, P., Zhao, Y., Chen, H. W., et al. (2012). ERRATUM: molecular dynamics simulations of Ago silencing complexes reveal a large repertoire of admissible ‘seed-less’ targets. Sci. Rep. 2, 569. doi: 10.1038/srep00909
Yuan, Y. R., Pei, Y., Ma, J. B., Kuryavyi, V., Zhadina, M., Meister, G., et al. (2005). Crystal structure of A. aeolicus argonaute, a site-specific DNA-guided endoribonuclease, provides insights into RISC-mediated mRNA cleavage. Mol. Cell 19, 405–419. doi: 10.1016/j.molcel.2005.07.011
Zander, A., Holzmeister, P., Klose, D., Tinnefeld, P., and Grohmann, D. (2014). Single-molecule FRET supports the two-state model of Argonaute action. RNA Biol. 11, 45–56. doi: 10.4161/rna.27446
Zhu, L., Jiang, H., Cui, X., Wang, Y., Xin, G., et al. (2016a). Understanding the core of RNA interference: the dynamic aspects of Argonaute-mediated processes. Prog. Biophys. Mol. Biol. 128, 39–46 doi: 10.1016/j.pbiomolbio.2016.09.008
Keywords: Thermus thermophilus Argonaute, gene editing, CRISPR, molecular dynamics simulation, MM-PBSA, conformational change
Citation: Liu Y, Yu Z, Zhu J, Wang S, Xu D and Han W (2018) Why Is a High Temperature Needed by Thermus thermophilus Argonaute During mRNA Silencing: A Theoretical Study. Front. Chem. 6:223. doi: 10.3389/fchem.2018.00223
Received: 13 February 2018; Accepted: 28 May 2018;
Published: 14 June 2018.
Edited by:
Nino Russo, Università della Calabria, ItalyReviewed by:
Tugba Kucukkal, Gallaudet University, United StatesAndras Szilagyi, Hungarian Academy of Sciences (MTA), Hungary
Copyright © 2018 Liu, Yu, Zhu, Wang, Xu and Han. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Dong Xu, xudong@missouri.edu
Weiwei Han, weiweihan@ jlu.edu.cn