Effects of the PIWI/MID domain of Argonaute protein on the association of miRNAi's seed base with the target.

The small interfering RNAs (siRNA) or microRNAs (miRNA) incorporated into the RNA-induced silencing complex with the Argonaute (Ago) protein associates with target mRNAs through base-pairing, which leads to the cleavage or knockdown of the target mRNA. The seed region of the s(m)iRNA is crucial for target recognition. In this work, a molecular dynamic simulation was utilized to study the thermodynamics and kinetic properties of the third seed base binding to the target in the presence of the PIWI/MID domain of Ago. The results showed that in the presence of the PIWI/MID domain, the entropy and enthalpy changes for the association of the seed base with the target were smaller than those in the absence of protein. The binding affinity was increased due to the reduced entropy penalty, which resulted from the preorganization of the seed base into the A-helix form. In the presence of the protein, the association barrier resulting from the unfavorable entropy loss and the dissociation barrier coming from the destruction of hydrogen bonding and base-stacking interactions were lower than those in the absence of the protein. These results indicate that the seed region is crucial for fast recognition and association with the correct target.

INTRODUCTION

RNA interference (RNAi) is induced by double-stranded RNA (dsRNA) and results in gene silencing of the target RNA. There are two types of RNAi: One is known as small interfering RNAs (siRNA), which is triggered by dsRNA helices having been introduced exogenously into cells; the other one is microRNA (miRNA), which is produced endogenously from small noncoding RNAs (Hammond et al. 2000; Zamore et al. 2000; Liu et al. 2004; Wang et al. 2008b). One strand of the s(m)iRNA duplex is loaded into the Argonaute (Ago) protein, forming the RNA-induced silencing complex (RISC) (Bartel 2004; Song et al. 2004; Tang 2005; Peters and Meister 2007; Hutvagner and Simard 2008; Kawamata and Tomari 2010; Wilson and Doudna 2013); and then the s(m)iRNA within the RISC associates with target mRNAs through base-pairing (Hutvágner and Zamore 2002; Martinez et al. 2002; Saito and Sætrom 2010), resulting in the cleavage or translation repression and destabilization of the target mRNA (Doench and Sharp 2004; Gregory et al. 2005). An RNA silence mechanism, such as the interaction and affinity between miRNA and mRNA in the presence of Ago, has attracted considerable interest and has been studied extensively in recent years (Martinez and Tuschl 2004; Vella et al. 2004; van den Berg et al. 2008; Wang et al. 2010; Chen and Zhang 2012; Dieterich and Stadler 2013; Gan and Gunsalus 2013; Khorshid et al. 2013; Schuck et al. 2013). Biochemical and bioinformatic studies indicated that Ago proteins divide s(m)iRNAs into five distinct domains, and the seed region that resides at the 5′-end of miRNAs and spans from nucleotide positions 2 to 7 was crucial for target recognition (Lewis et al. 2005; Bartel 2009; Pasquinelli 2012; Gorski et al. 2017). The crystal structures (see Fig. 1) of the RNA–protein complex show that the seed region of the guide strand is settled by its phosphodiester backbone to the PIWI/MID domain of Argonaute protein, and the seed region in an A-form helix is important for target recognition (Ma et al. 2005; Parker et al. 2005; Wang et al. 2008a). It has been found that the affinity of the seed-target interaction was greatly increased upon the association of the AfPiwi protein, a model of the PIWI/MID domain of the Argonaute protein (Parker et al. 2009). Single molecule experiments showed that incorporation of the guide RNA into RNP greatly increased the association rate constant (Chandradoss et al. 2015; Jo et al. 2015; Salomon et al. 2015). The RNP finds its targets less efficiently and binds to them less stably without a matched seed sequence.

FIGURE 1.

The structure of the model system (PDB ID: 4W5O). The MID domain (yellow), the PIWI domain (orange), the L2 linker (tan), the seed region of miRNA (blue), and the target RNA (purple).

Since its discovery, many methods have been developed to predict miRNA and its targets through thermodynamic properties (Enright et al. 2003; Lewis et al. 2003; Krek et al. 2005; Hsu et al. 2006; Krüger and Rehmsmeier 2006; Maragkakis et al. 2009; Peterson et al. 2014). However, due to lack of thermodynamic parameters of RNA with the protein, the seed's effects could not efficiently incorporate the methods, and the accuracies are limited. The methods could improve their performance by using more accurate thermodynamic parameters of RNA with the protein. Molecular dynamic simulation as a powerful tool has been used to study the thermodynamic and kinetic properties of RNA (Várnai et al. 2004; Zhang et al. 2008, 2011; Gong et al. 2011; Nguyen et al. 2017). Recently, the simulation has been used to obtain the quantitative thermodynamic and kinetic parameters of the terminal base pair upon opening-closing transition (Wang et al. 2016, 2018). In this paper, we obtained the thermodynamic and kinetic properties of the third seed base binding to the target in the presence of the Argonaute protein through simulating the association-dissociation switch of the base pair. The results showed that the binding affinity was increased due to the reduced entropy penalty, which resulted from the preorganization of the seed base into the A-helix form. The association barrier resulting from the unfavorable entropy loss and the dissociation barrier coming from the disruption of hydrogen bonding and base-stacking interactions in the presence were lower than those in the absence of the protein, leading to the increase of the rate constant with which it binds to the target.

RESULTS AND DISCUSSION

The nucleotide of the miRNA in the seed region is preorganized

As has been shown, the terminal base pair would go through the open–close transition through a transition state when there is no protein (Wang et al. 2016, 2018). The open state, closed state, and transition state can be classified from the root mean square deviation (RMSD) of the two terminal nucleotides relative to their starting structures (see Figs. 2A, 4A), the backbone torsion angles ζ (see Figs. 3A, 4B), and the dihedral angle of the four atoms C3′(n)-O3′(n)-P(n + 1)-O5′(n + 1) (Hershkovitz et al. 2003), where n represents the nth nucleotide in a polynucleotide chain.

FIGURE 2.

(A) The RMSD of guide and target bases during the whole simulation time at the temperatures of 530, 540, 560, and 580 K. (B) The distribution of RMSD for the seed and target bases at 540 K. (C) The distribution of RMSD for the seed base at 540 K. (D) The distribution of RMSD for target base at 540 K.

FIGURE 4.

(A) The RMSD of the seed and target bases during the simulation times of 0–8 nsec at 540 K. (B) The torsion angles ζ of the target base for the simulation times of 0 –8 nsec at 540 K. (C) The RMSD and torsion angles ζ near the transition states. (D) The atomic structures of the association, dissociation, and transition states (ata and dtd).

FIGURE 3.

(A) The torsion angle of the target base during the whole simulation time at temperatures of 530, 540, 560, and 580 K. (B) The distribution of the torsion angle of the seed base at 540 K. (C) The distribution of the torsion angle of the target base at 540 K.

Closed state: The conformations are considered as a closed state when the RMSD of the terminal bases is ∼0.7 Å and the torsion angle ζ is between −50° and −100° (see Fig. 4A,B). At this state, the two terminal bases are only vibrating at their starting locations, and the base pair and base stacking are intact.

Open state: The conformations are classified as open state when the RMSD is larger than 2 Å and the torsion angle ζ is between 50° and 100° (see Fig. 4A,B). At this state, the two bases deviated away from the starting positions and the terminal bases flipped into the solvent.

Transition state: The conformations are denoted as transition state when the RMSD of terminal bases is larger than 2 Å, but the torsion angle ζ remains between −50° and −100° (see Fig. 4C). In this state, the torsion angle was still in the closed state region while the bases were flipped out into the solvent (see Fig. 4D). The transition states can be divided into ctc and oto (Dokholyan et al. 2000; Zhang et al. 2006); ctc is defined as the transition state between the closed state, which transited from the closed state and then back to the closed state; and oto is the transition state between the open state, which transited from the open state and then back to the open state.

To elucidate the association and dissociation mechanism of the seed base 5′-A and the target base 3′-U, we also calculated the RMSD of the two bases relative to their starting positions and the corresponding backbone torsion angles ζ. Figure 2A shows the time-dependent RMSD relative to the starting structure of the base pair of the seed base 5′-A and the target base 3′-U. Similar to that of the AU base pair in the absence of the PIWI/MID domain (Wang et al. 2016), the values of the RMSD distributed in a wide range (from 0 to 12 Å) but clearly centered on two regions, as shown in Figure 2B. In the first region, the RMSD fluctuated about 0.7 Å, indicating that the atoms of the seed base 5′-A and the target base 3′-U only vibrate around their initial positions and the base pair and base stacking were intact. In the other region, the RMSD was centered around 2 Å and it distributed in a wide range. The time-dependent RMSD of the seed base 5′-A and the target base 3′-U relative to the starting structures showed different distribution behavior, as shown in Figure 2C,D. The RMSD of the target base 3′-U relative to its starting structure was distributed in a wide range and the distributions were clearly centered on 0.7 and 3 Å regions, as shown in Figure 2D. However, the RMSD of the seed base 5′-A was almost only centered at 0.7 Å, as shown in Figure 2C, which denotes that the corresponding conformations of the seed base 5′-A vibrated only at its starting position.

We also calculated the backbone torsion angles ζ (see Figs. 3A, 4B) of the seed base 5′-A and target 3′-U. Similar to the case without the protein, the torsion angle ζ of the target base 3′-U was clearly centered on −75° and 50° regions (Fig. 3C), whereas the torsion angle ζ of the 5′-A was always in the region of −75° (Fig. 3B). The RMSD and the torsion angle ζ showed consistency (Figs. 4A,B): When RMSD was in the first region (0.7 Å), the torsion angle ζ was also in the first region (−75°); and when RMSD was in the second region (2 Å), the torsion angle ζ was also in the second region (50°). It has been found that the torsion angle ζ in the helix and coil state is just distributed about −75° and 50°, respectively (Hershkovitz et al. 2003). So according to the RMSD and the torsion angle, all the conformations of the seed base 5′-A and target 3′-U base pair are classified into three states as that without protein: association state (closed state), dissociation state (open state), and transition state. So even the 5′-A of the seed nucleotide was not paired with the target nucleotide 3′-U; it still kept its position as it paired, denoting a preorganization of the seed region.

The reduced entropy penalty enhances the affinity of the seed-target interaction

At each simulation temperature, the probability of the association, dissociation, and transition states can be obtained as

respectively, where τ is the total simulation time, and Na, Nd, and Nt are the total number of occurrences of the corresponding conformations that are in the association, dissociation, and transition states, respectively. , , and are the period in which the conformations are staying at the ith times in the association, dissociation, and transition states, respectively. Figure 5A showed the occupied probabilities of the association state during the simulation at all the simulation temperatures, in which each point is calculated at the time interval from initial to corresponding simulation time. It can be seen that for all the simulation temperatures, when the simulation time gets to 3600 nsec, the occupied probability of the association state reached a stable value, denoting the association-dissociation equilibration. In order to make sure the simulation reaches the association-dissociation switch equilibrium, each simulation temperature lasted 4000 nsec. The probabilities of the association, dissociation, and transition states at the simulation temperatures are listed in Table 1.

FIGURE 5.

(A) The probability of the association state during the simulation time at the following simulation temperatures: 530, 540, 560, and 580 K. (B) The ln(p/p) at different simulation temperatures. Symbol: simulation results; line: linear fit.

TABLE 1.

The average lifetime τ(nsec), the occupied probability p, and the total number N of occurrences of conformations at association, dissociation, and transition states at all simulation temperatures (K) during the whole 4000 nsec simulation time

Table 1 showed that the transition states only occupied negligible probabilities. So the conformations of the seed-target base pair were first treated as association and dissociation two-states. Therefore, the free energy difference ΔG for the association of the third seed nucleotide with the target can then be obtained through the equilibrate probabilities of the association and dissociation states: where p and p are the probability of the dissociation state and association state at temperature T, respectively; k is the Boltzmann constant.

As shown in Figure 5B, the ln(p/p) and the reciprocal of temperature (1/T) presented a linear relation, indicating the linear relationship of the free energy difference ΔG and the temperature T. As ΔG = ΔH − TΔS, where ΔH and ΔS are the enthalpy and entropy changes for the transition between the association state and the dissociation states, respectively, the thermodynamic parameters of the seed-target AU base pair are found to be ΔH = −4.95 kcal/mol and ΔS = −6.78 eu, while they are ΔH = −7.3 kcal/mol and ΔS = −18.5 eu for the same terminal base pair without the PIWI/MID domain (Wang et al. 2016). Although the enthalpy change was decreasing in the presence of the protein, the free energy difference ΔG = ΔH − TΔS between the association and dissociation state in the presence/absence of protein are −2.85 kcal/mol and −1.56 kcal/mol at T = 310 K, respectively. So the presence of the PIWI/MID domain enhances the interaction between the seed and target RNA through a reduced entropy penalty.

The PIWI/MID domain of Argonaute protein reduces the barrier for association/dissociation of the seed and target base pair

The average lifetimes of the association state, dissociation state, and transition state can be obtained through , where τ is the average lifetime, N is the total number of occurrences of the 5′-A seed, and 3′-U target base pair in the association state or dissociation state or transition state, and τ is the ith dwelling period of the conformation in the corresponding state. Figure 6 showed the distribution of the lifetime of the association, dissociation, and transition states. The average lifetimes of the association and dissociation states at different temperatures are shown in Figure 7A. The average lifetime of the association state exhibited strong temperature dependence. However, that of the dissociation state was only weakly dependent on temperature, which was in agreement with the experimental finding that the folding and unfolding rates of RNA had different temperature-dependent behaviors (Craig et al. 1971).

FIGURE 6.

Histograms of the lifetime distribution of the association state, dissociation state, transition state (ata and dtd) at temperatures of 530, 540, 560, and 580 K.

FIGURE 7.

(A) The average lifetime of the association state (square) and the dissociation state (circle) at the simulation temperatures. Symbol: simulated results. Lines: fitted with Equation 2. (B) The ratios of (square) and (triangle). Symbol: simulated results. Lines: fitted via Equations 2 and 3.

For the seed-target base pair association-dissociation two-state transition, the association rate k+ and the dissociation rate k− of the seed and target base pair can be obtained from the average lifetime of the dissociation and association states as , , where and are the average lifetimes of the dissociation state and association state, respectively.

The transition rates from the transition state to the association state k and to the dissociation state k can be acquired from the average lifetime of the transition states: and , where and are the average lifetimes of the transition states ata and dtd, respectively.

As shown in Table 1, the average lifetime of the transition states was much shorter than the average lifetimes of the association and dissociation states, and it was also weakly dependent on temperature. Based on the transition-state theory (Berne et al. 1988; Hänggi et al. 1990; Chung and Eaton 2013), the average lifetime of the association state (t) and the dissociation state (t), the average transition path time from the dissociation state to the association state () and from the association state to the dissociation () are given by the following equations: where D* is the diffusion coefficient at the top of the free energy barrier; (ω*)2, (ω)2, and (ω)2 are the curvatures of the free energy surface at the top of the barrier, in the association and dissociation states, respectively; ΔG is the free energy barrier height from the association state to the dissociation state, ΔG is the free energy barrier height from the dissociation state to the association state, β = 1/k T, k is the Boltzmann constant; T is the temperature; and γ is Euler's constant.

It can be seen from Equations 2 and 3 that the ratios of and are independent of the diffusion coefficient D*, but only depend on the free energy barriers ΔG, ω*/ω, and ω*/ω. We found that was nearly temperature independent (Fig. 7B). Considering that the ratio ω*/ω is constant, as observed in the protein folding (Chung and Eaton 2013), then the free energy barrier of association for the seed 5′-A and target 3′-U base pair from the dissociation state should be directly proportional to temperature as ΔG ∝ T. The free energy difference of the association and dissociation states at temperature T was ΔG = ΔG − ΔG = ΔH − TΔS. By fitting the two curves of and with the temperature, the free energy barrier from the association state to the dissociation state was ΔG = −4.95 kcal/mol, which coincides with the enthalpy change ΔH between the association state and dissociation state. The free energy barrier for the association of the seed 5′-A and target 3′-U was ΔG = TΔS, where ΔS is exactly equal to the entropy difference between the dissociation and association states. So the barrier for the association of the seed miRNA and target mRNA was the reduction in entropy, ΔS, and the barrier for the dissociation of the seed miRNA and target mRNA was the increase in enthalpy, ΔH, by destroying the hydrogen bond and the base-stacking interactions. Therefore, the PIWI/MID domain would decrease the barrier for the seed-target association due to the reduced entropy penalty, which is conductive to binding the target quickly, and the decrease of the dissociation barrier could be helpful to rapidly dissociate from the wrong target.

Conclusions

In conclusion, the thermodynamic and kinetic parameters of the third position in the seed region of the miRNA binding with the target were acquired by simulating the association–dissociation transition of this base pair in the presence of the PIWI/MID domain of the Argonaute protein. The results showed that in the presence of the protein, the entropy and enthalpy changes upon the seed base of the s(m)iRNA binding to the target were less than those in the absence of the protein. The binding affinity was increased due to the reduction of the entropy penalty, which resulted from the preorganization of the seed base into the A-helix form. The presence of the PIWI/MID domain would lower the association barrier resulting from the unfavorable entropy loss and the dissociation barrier coming from the disruption of hydrogen bonding and base-stacking interactions. These results indicate that the seed region is important for fast binding to the target, and would quickly dissociate when it was binding with the wrong target. As more structures of the RNA-protein complexes are given (Schirle et al. 2015; Sheu-Gruttadauria and MacRae 2018), the thermodynamic and kinetic properties of s(m)iRNA and target interactions over and outside the seed regions as well as the sequence and position dependence should be further studied in the future.

MATERIALS AND METHODS

The Amber10 all-atom force field (Zgarbová et al. 2011) and Gromacs 4.5.6 (Pronk et al. 2013) were used to perform all the molecular dynamic simulations. The model system consists of an A-form duplex, which mimics nucleotides 2 to 6 over the seed region of guide RNA base-pairing with the messenger RNA, and the PIWI/MID domain of the Argonaute protein. The RNA sequence was r(AAUUU)·r(UUAAA). The starting structure of the RNA and protein complex (see Fig. 1) was obtained from the crystal structure (PDB ID: 4W5O) (Schirle et al. 2014). The protein contains N terminus, Linker 1, PAZ, Linker 2, MID, and PIWI domains, but the seed region was enveloped by Linker 2, PIWI, and MID domains. We retained the Linker 2 and PIWI/MID domains to build the system. The NaCl concentration was set to 0.5 M by adding Na+ and Cl− ions (Joung and Cheatham 2008) at random initial positions. The TIP3P (Jorgensen et al. 1983; Mahoney and Jorgensen 2000) water model and the SETTLE algorithm was used to keep water molecules completely rigid (Miyamoto and Kollman 1992). The simulating box is set as a 9.9 × 9.1 × 9.0 nm triclinic box and periodic boundary conditions were used, ensuring three or four solvation layers in each direction. The final system had about 68,440 atoms. The simulations were performed at constant temperature and pressure by using velocity rescaling (Bussi et al. 2007) and the Parrinello–Rahman barostat algorithm (Parrinello and Rahman 1981), respectively. The Particle-Mesh Ewald method (Darden et al. 1993; Essmann et al. 1995) was used to treat the electrostatic interactions with a 10 Å direct space cutoff. Lennard-Jones interactions (Lennard-Jones 1931) were truncated at 10 Å. Bond lengths of the solute were constrained by the LINCS algorithm (Hess et al. 1997). The neighboring grid search method (Hess et al. 2008) was used and was updated every 10 steps. The equations of motion were integrated via the Verlet algorithm. The time step was 2 fsec and the coordinates were saved every 2 psec. The system was relaxed by energy minimization to equilibrium over 20 nsec at 290 K. The starting structure for the simulations at high temperatures was selected from the equilibrium structures. To obtain the thermodynamic and kinetic properties of the association and dissociation of the guide 5′-A and target 3′-U base pair in the presence of the PIWI/MID domain of Argonaute protein, all other atoms of the complex except those of the guide 5′-A and target 3′-U nucleotides were fixed with position restraint by using an additional force (force constant: 1000 kJ/mol nm−2). To effectively sample the conformation space of the associated and dissociated states, the simulation temperatures were 530, 540, 560, and 580 K; each simulation lasted 4000 nsec.