Effect of mutations on binding of ligands to guanine riboswitch probed by free energy perturbation and molecular dynamics simulations

Abstract Riboswitches can regulate gene expression by direct and specific interactions with ligands and have recently attracted interest as potential drug targets for antibacterial. In this work, molecular dynamics (MD) simulations, free energy perturbation (FEP) and molecular mechanics generalized Born surface area (MM-GBSA) methods were integrated to probe the effect of mutations on the binding of ligands to guanine riboswitch (GR). The results not only show that binding free energies predicted by FEP and MM-GBSA obtain an excellent correlation, but also indicate that mutations involved in the current study can strengthen the binding affinity of ligands GR. Residue-based free energy decomposition was applied to compute ligand-nucleotide interactions and the results suggest that mutations highly affect interactions of ligands with key nucleotides U22, U51 and C74. Dynamics analyses based on MD trajectories indicate that mutations not only regulate the structural flexibility but also change the internal motion modes of GR, especially for the structures J12, J23 and J31, which implies that the aptamer domain activity of GR is extremely plastic and thus readily tunable by nucleotide mutations. This study is expected to provide useful molecular basis and dynamics information for the understanding of the function of GR and possibility as potential drug targets for antibacterial.


INTRODUCTION
Riboswitches have been regarded as a new class of genetic control elements founded in the 5 -leader region in multiple bacterial messenger RNAs (mRNA) and play signif-icant roles in modulation of gene expression in bacteria, plants and fungi (1,2). Riboswithches consist of two mutually interacting domains, namely the aptamer domain and the expression platform. The aptamer can directly bind to the metabolite molecules and is responsive to intracellular ligand concentration. The expression platform ensures the structural transformation in response to the changes in the aptamer so as to modulate either ribosome binding or transcription antitermination. Presumably, riboswitches modulate gene expression by an allosteric rearrangement due to binding of small metabolite molecules (3)(4)(5)(6)(7)(8)(9). Riboswitches have currently attracted interests as potential drug targets for antibacterial (10).
By now, the metabolite molecules recognized by riboswitches mainly include amino acids (11,12), nucleotides (13,14), vitamin cofactors (15,16) and metal ions (17). All riboswitches can not only bind their respective targets with high affinity and selectivity, but also discriminate even against very closely related compounds (18,19). The purine riboswitches, including the guanine-specific riboswitch (GR) and the adenine-sensing riboswitch (AR), were found in 2003 (13,14). The crystallographic structures suggest that these two classes of purine riboswitches consist of three helices P1, P2 and P3 that surround a three-way junction J12, J23 and J31 connecting them, with phylogenetically conserved loops L2 and L3 capping P2 and P3 (20), respectively ( Figure 1A and B). The ARs and GRs share highly conserved primary and secondary structure (14), it is observed that a cytosine 74 (C74) is conserved in all GRs and a uridine 74 (U74) is conserved in all ARs. Mutation of C74 into U74 makes the GRs lose the binding ability to guanine and inversely have a strong affinity to adenine (14). Thus, it is of significance to probe molecular mechanism involving the conformational change and ligand binding mechanism of riboswitches for further understanding the role of riboswitches as potential drug targets for antibacterial. Based on significant roles played by purine riboswitches in the modulation of gene expression, many experimental studies have focused on its structure and function. The previous structural study suggested that the mRNA rearrangement upon ligand binding will build a pocket which can interact with all functional groups of the purine and form a Watson−Crick base pair with C74/U74 of AR/GR (14). This information was further supported by the crystal structures of the ligand-binding domain of the purine riboswitches (20,21). Significantly, the work from Gilbert et al. contributed multiple crystal structures of the Bacillus subtilis xpt-pbuX GR and its mutant C74U (GR C74U) associated with different ligands and the related binding data (22). Moreover, several groups around the world have performed different experiments to probe the conformational changes of purine riboswitches induced by ligand bindings (23,24). For example, Brenner et al. applied single molecule fluorescence resonance energy transfer spectroscopy to explore the conformational changes of the xpt GR aptamer (23). Stoddard et al. adopted X-ray crystallography and chemical probing to study the effect of modest sequence alterations on the activity of the purine riboswitches, and their results suggested that the introduction of non-natural compositions induces instability to regulate gene expression in vivo and the aptamer domain activity is highly plastic and thus readily tunable to meet cellular needs (25). Batey et al. performed a detailed mutagenic survey on the xptpbux GR from B. subtilis and the results indicated that the conserved nucleotides, mainly arising from the formation of the ligand-GR complex, promote an open state of the apo RNA and take part in the secondary structural switch with the expression platform (19). Therefore, for deeply understanding the role of riboswitches as potential antibacterial drug targets, it is of significance to further investigate the conformational changes of riboswitches and the energy basis of ligand binding to riboswitches.
Currently, MD simulation (26)(27)(28)(29)(30)(31)(32) and binding free energy prediction (33)(34)(35) have been extensively adopted to explore the interaction mechanism of biomolecules with small organic ligands, such as the ligand-induced conformational change and the free energy basis of the binding. Moreover, MD simulation has been applied to successfully probe the conformational dynamics of RNA influenced by ligand bindings or mutations in sequence (18,(36)(37)(38)(39)(40). Recently, multiple computational methods with various levels of computational accuracy and expense, including the free energy perturbation (FEP) (41)(42)(43)(44), molecular mechanics Poisson−Boltzmann surface area (MM-PBSA) (45)(46)(47)(48)(49), and the thermodynamic integration (TI) (50-52) methods, were employed to study the binding affinity of ligands to RNAs. Hu et al. combined the TI and MM-PBSA methods to study the binding of four guanine analogues to the GR and GR-C74U, and their results indicated that the loss in binding free energies upon mutation is mainly driven by electrostatic interactions (53). Sund et al. integrated the FEP with linear interaction energy (LIE) method to calculate the binding affinities of 14 different purine analogues with GR and AR, and their results demonstrated that the AR specifically uses the electrostatic preorganization to discriminate against guanine by preventing the formation of a G-U wobble base pair (54).
In this study, the influence of sequence mutation on the binding mechanisms of the xpt-pbuX GR with ligands are explored by the combination of the MD simulation, FEP, and MM-GBSA method. Two ligands were selected, including hypoxanthine (HPA) and 9H-purine-2,6-diaminopurine (6AP), and their structures were shown in Figure 1C and D. Meanwhile, several mutations of xpt-pbuX GR, including A24U, U25A/A46G, A24U/U25A/A46G, AU25A/A46G/C74U, A24U/U25A/A46G/C74U were considered to comparatively probe changes in binding affinities and conformation of the xpt-pbuX GR induced by mutations. This study is expected to provide significant energy basis and dynamic information of GRs relating with mutations so as to better understand the role of riboswitches as potential drug targets for antibacterial.

Construction of molecular systems
Atomic coordinates of the wild-type (WT) GR complexed with HPA were obtained from protein data bank (PDB ID: 4FE5) (20). The X-ray structures of HPA complexed with several mutants of GR (A24U, U25A/A46G and A24U/U25A/A46G) can also be obtained from protein data bank (PDB ID: 4FEJ, 4FEL and 4FEN) (25). Atomic coordinates of 6AP complexed with mutated GR (U25A/A46G/C74U and A24U/U25A/A46G/C74U) were taken from protein data bank (PDB ID: 4FEO and 4FEP) (25). The initio structure of 6AP-WT GR complex can be constructed based on the complex structure of HPA with WT GR. The missing hydrogen atoms were connected to the corresponding heavy atoms by using the Leap module in Amber (55,56). The Amber ff99SB force field (57) was employed to produce the topology files of the WT and mutant GRs. The TIP3P model was selected for water molecules. The structures of HPA and 6AP were optimized with the Gaussian 09 program (58) at the HF/6-31G* level, and then the restrained electrostatic potential (RESP) charges and the general Amber force field (GAFF) (59) were assigned to them. The force field parameters of these two ligands are provided as the supplemental data (Files S1-S4). The initialized structures of each complex system was solvated in a truncated octahedral box of TIP3P water molecules with a 14.0Å buffer along each dimension (60). The appropriate number of sodium counterions were added to neutralize each system. Details of the system preparation were listed in Supplementary Table S1.

Conventional MD simulations
To remove high energy contacts formed during the initialization of the simulated systems, each system was optimized by energy minimization using a combination of the steepest descent and conjugated gradient methods. Subsequently, each system was gently heated from 0 to 300 K in 1 ns at constant volume and equilibrated at 300 K for another 1 ns. Finally, a 1-s MD simulation was conducted without any restrictions at constant pressure and 300 K to fully relax each system. In the current work, all MD simulations were ran by using the PMEMD program (61,62) in Amber. During MD simulations, all chemical bonds contain hydrogen atoms were restricted using the SHAKE algorithm (63) and the time step was set to 2 fs. The Langevin thermostat with a collision frequency of 2.0 ps −1 was adopted to regulate the temperature of the systems (64). The particle mesh Ewald (PME) (65,66) method was used to compute the long-range electrostatic interactions with the cutoff distances of 12.0Å, which is also suitable for the calculations of van der Waals interactions. The cross-correlation map and principal component analysis (PCA) (67)(68)(69)(70)(71) were carried out using the CPPTRAJ module (72) in Amber and the details involving these analyses have been presented in our previous studies (73,74). Meanwhile, the cluster analysis was also performed using the iMWK-Means method proposed by Salsbury, Jr et al. to reveal conformational changes of GRs induced by mutations (67,75,76). The PyMOL (http: //www.pymol.org.) and VMD (77) softwares were utilized to check MD trajectories and draw pictures.

Free energy perturbation calculations
The free energy perturbation (FEP) method proposed by Shirts et al. (78) was adopted to calculate the absolute binding free energies ( G bi nd ) of ligands to the WT and mutated GRs. The thermodynamics cycle for the FEP was shown in Figure 2 and G bi nd can be determined by the following equation where 1 and 0 indicate the states of the real ligands and the dummy ones without any molecular properties, respectively. describes the alchemical coupling parameters transforming from 0 to 1 and . . . λ denotes an ensemble average in states . ∂ V(λ)/∂λis the derivative of the hydrid potential function on . As G 2 = 0, G bi nd can be represented as in which G 1 and G 3 represent the change of free energies associate with the disappearance of ligand from GR and bulk water, respectively. In the current study, all FEP calculations were performed by using Gromacs 5.0.4 (79). The last snapshot taken from the first 200 ns of the conventional MD simulation of each system was used as the starting conformation of FEP calculation. The van der Waals and electrostatic interactions between ligands and GRs were annihilated by using a linear alchemical pathway with a -value interval of 0.05 for van der Waals interactions, and a -value interval of 0.1 for electrostatic interactions (42). In addition, 12 nonuniformly distributed values (0.0, 0.01, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.5, 0.75, 1.0) were employed to add restrictions on ligands. In this study, a total of 42 and 31 windows were used for simulations of the ligand-GR complexes and ligands in the bulk water, respectively. For each window, energy minimization of 3000 steps was conducted to remove high-energy contacts between atoms by using the steepest descent algorithm. Subsequently, each system was relaxed for 1ns in the canonical ensemble and the non-hydrogen atoms of the solute were restricted with a force constant of 1000 kJ·mol −1 ·nm −2 . The temperature in the simulation was controlled at 300K by using a V-rescale thermostat with a response time of 1.0 ps (80). Meanwhile, the pressure of the simulated system was regulated by using the Parrinello-Rahman coupling algorithm with a target pressure of 1 atm (81). For each window, the Berendsen weak coupling algorithm was applied to deeply relax the system for 1 ns in the isothermalisobaric ensemble in each window (82). Finally, 5 ns produc-tion simulation was performed without any restrictions for each window to collect data using Hamiltonian-exchange Langevin dynamics in the NPT ensemble. 1 million swaps between the state pair in an interval of 4 ps were carried out by using the Gibbs sampling scheme proposed by Shirts et al. (78,83). During simulations, the LINS method (84) was adopted to constrain the chemical bond that contain hydrogen atoms, and the time step of simulation was set to 2 fs. The electrostatic interactions were treated by using the Particle Mesh Ewald (PME) method with a cut-off distance of 12Å, a spline order of 6, a relative tolerance of 10 −6 and a Fourier spacing of 1.0Å. The van der Waals interactions were computed using the soft-core potential and a switch function between 10 and 12Å (85). The above scheme treating non-bonded interactions aim to relieve singularities, numerical instabilities, and additional minima in the potential energy for all combinations of non-bonded interactions in all intermediate alchemical states (86). Based on the success of FEP method in predicting binding ability of ligands to receptors (42,43,52,(87)(88)(89), thus FEP method was employed to evaluate the effect of several mutations on binding ability of HPA and 6AP to GRs.

MM-GBSA calculations
For in-depth comparison, the molecular mechanics generalized Born surface area (MM-GBSA) method (49,90,91) has also been applied to compute the binding free energies of HPA and 6AP to GRs: in which the first three terms represent binding free energy components in gas phase and the last two terms denote solvation free energies. E vdW and E ele are van der Waals and electrostatic interactions between ligands and GRs and these two components were calculated by using molecular mechanics and ff99SB force field. −T S indicates the contributions of the entropy changes to ligand associations with GR, which is computed by using the mmpbsa py nabnmode program (92). G egb is polar solvation energy and can be calculated with the GB model developed by Onufriev et al. (93). The last term G esur f is non-polar solvation free energy, which is determined using the following equation where the parameter γ and SASA respectively denote the surface tension and the difference in the solvent accessible surface areas caused by ligand associations. The parameter β reflects a regression offset of the linear relationship. For our current work, γ and β were assigned as 0.0072 kcal·mol·Å −2 and 0 kcal·mol −1 , respectively (94). Due to the high cost in calculation of the entropy, only 50 snapshots extracted from the equilibrated trajectories were adopted for the entropy calculation, while for calculations of the other free energy components, the total of 200 conformations taken from the equilibrated MD trajectories were used.

Effect of mutations on binding ability of ligands to GR
To comparatively evaluate the effect of mutations on the binding ability of ligands to GRs, FEP and MM-GBSA methods were applied to compute the binding affinity of ligands 6AP and HPA to the WT and mutated GRs. The calculated results from FEP and MM-GBSA methods are independently listed in Table 1 and 2. Due to the lack of experimental free energy information involving mutations on the ligand-GR bindings, the correlation relationship between binding free energies computed by FEP method and MM-GBSA method is estimated to evaluate the reliability and the results are depicted in Supplementary Figure S1. It is observed that the correlation coefficient (R 2 ) reaches 0.83, which suggests that our current free energy calculations on predicting the binding affinity of HPA and 6AP to the WT and mutated GRs are reliable. In FEP calculations, 5-ns production simulations were performed on each window to obtain fully statistical samplings and the last 3-ns trajectories of each window simulation were utilized to estimate the binding free energies by using the g bar tool in Gromacs 5.0.4 (Table 1). In Table 1, G 3 and G 1 represent the changes of free energies in decoupling ligands from water and GRs, respectively. As shown in Table 1, the change of free energy caused by decoupling HPA from water is close to that caused by decoupling 6AP from water and the difference between them is only 0.23 kcal/mol. However, the free energy change induced by decoupling 6AP from the binding site of the WT GR is 5.68 kcal/mol stronger than that caused by disappearing HPA from the WT GR, which provides significant contribution to the difference in the binding ability of HPA and 6AP to the WT GR. According to Table 1, the binding affinity of 6AP to the WT GR is strengthened by 5.45 kcal/mol relative to that of HPA to the WT GR, suggesting that binding strength of 6AP to the WT GR is much stronger than HPA and also implying that the WT GR has a binding preference and selectivity for 6AP over HPA. Comparing to the WT GR, the mutations A24U, U25A/A46G and A24U/U25A/A46G separately lead to the strengthening of 0.77, 0.53 and 1.16 kcal/mol in binding abilities of HPA to GR. Similarly, the mutations U25A/A46G/C74U and A24U/U25A/A46G/C74U respectively result in the enhancing of 3.09 and 1.14 kcal/mol in the binding affinity of 6AP to the mutated GR relative to the WT one. Based on the above results, it is concluded that the GR mutants involved in the current work obviously strengthen the binding ability of ligands HPA and 6AP to the mutated GRs.
In MM-GBSA calculations, 1-s MD simulations were performed on the seven studied systems. Root-mean-square deviations (RMSD) of atoms P, O3 , O5 C3 , C4 and C5 in GR were computed relative to the first structure through the entire MD simulations to evaluate the reliability of system equilibrium and the results were displayed in the supporting information (Supplementary Figure S2). It is observed that all systems basically reach the equilibrium after 200 ns of MD simulations. To obtain a good conformational sampling, a total of 200 conformations extracted from the last 800 ns trajectories of MD simulations with a time interval of 4 ns were used to compute the binding free energies. Binding free energies and separate free energy components were listed in Table 2.
According to Table 2, binding free energies are decomposed into van der Waals interactions ( E vdW ), electrostatics interactions ( E ele ), entropy contributions (−T S), polar solvation free energies ( G egb ) and non-polar solvation free energies ( G esur f ). It is found that the binding free energy of 6AP to the WT GR is strengthened by 5.22 kcal/mol comparing to that of HPA to the WT GR, indicating that the binding ability of 6AP to the WT GR is stronger than HPA. Meanwhile, comparing to the WT GR, binding free energies of HPA to the mutated GRs are separately enhanced by 1.16, 1.65 and 5.44 kcal/mol due to the mutations A24U, U25A/A46G and A24U/U25A/A46G, while that of 6AP to the mutated GRs are improved by 4.61 and 3.89 kcal/mol because of the mutations U25A/A46G/C74U and A24U/U25A/A46G/C74U, suggesting that these mutations strengthen the binding of two ligands to the mutated GRs. These results agree well with the previous FEP calculations. As shown in Table 2, although A24U and A24U/U25A/A46G lead to the decrease of 2.40 and 0.73 kcal/mol in van der Waals interactions of HPA with GR, the polar interactions of HPA with GR ( G ele+egb ), an unfavorable factor generated by the sum of electrostatic interactions and polar solvation free energies, are decreased by 3.35 and 6.79 kcal/mol in the two mutated systems, which completely compensate the decrease in van der Waals interactions. Thus, the decrease in polar interactions of HPA with GR induced by A24U and A24U/U25A/A46G is a main force of enhancing the binding ability of HPA to the mutated GRs. According to Table 2, comparing to the WT GR, the increase of 2.60 kcal/mol in the polar interaction of HPA with the U25A/A46G mutated GR weakens the binding of HPA to GR, but the van der Waals interaction of HPA with the U25A/A46G mutated GR is improved by 2.40 kcal/mol, besides the entropy change (−T S) due to the mutations U25A/A46G is also decreased by 1.69 kcal/mol relative to the WT one. Thus, it totally strengthens the binding ability of HPA to the U25A/A46G mutated GR. In the case of 6AP, mutations U25A/A46G/C74U generate a decrease of 0.81 kcal/mol in van der Waals interactions of 6AP with the mutated GR relative to the WT one, which slightly impairs the association of 6AP with GR. However, the mutations U25A/A46G/C74U respectively lead to the decrease of 4.00 kcal/mol and 1.41 kcal/mol in polar interactions and entropy changes compared to the WT GR, which not only completely compensate the decrease in van der Waals interactions of 6AP with the U25A/A46G/C74U mutated GR, but also totally strengthen the binding of 6AP to the mutated GR. As seen in Table 2, the decrease of 0.72 kcal/mol in the entropy changes induced by the mutations A24U/U25A/A46G/C74U softly strengthens the association of 6AP with the mutated GR. However A24U/U25A/A46G/C74U produces the decrease of 1.56 kcal/mol in van der Waals interaction and 1.54 kcal/mol in polar interaction of 6AP with GR, which totally enhances binding of 6AP to the A24U/U25A/A46G/C74U mutated GR. Based on the above analysis of the binding free energies from FEP and MM-GBSA calculations, the mutations mentioned in this study enhance the bindings of HPA and 6AP to GRs, implying that these mutations may exert a dramatic impact on the activity of GR. The work from Stoddard et al. also suggested that the aptamer domain activity of GR is highly plastic and thus readily tunable to meet cellular needs, which is in good agreement with our current conclusions (25).

Impact of mutations on ligand-GR interaction networks
To reveal the molecular mechanism under the strengthened binding ability of HPA and 6AP to mutant GRs, residuebased free energy decomposition method was utilized to compute the interaction spectrum of 6AP and HPA with the separate nucleotides of GRs (Figure 3 and Supplementary Figure S3). Geometrical positions of 6AP and HPA relative to the key nucleotides in GRs were depicted in Figures 4, 5 and Supplementary Figure S4 by using the averaged structures taken from the last 10 ns of MD trajectories. Hydrogen bonding interactions of 6AP and HPA with nucleotides were analyzed with the program CPPTRAJ in Amber and the information was listed Table 3 and Supplementary Table  S2.
In addition, it is worth noting that several water molecules were observed in the binding pocket of all averaged structures taken from the last 10 ns trajectories of MD simulations, which indicates that water molecules may be of significance for the binding of ligand to GRs and function of GRs. Consistent with our current results, Sund et al. also suggested that several water molecules appear in the averaged structure of GR from MD simulations (54). Meanwhile, Batey et al. also found that 5% of the surface in the center of a three-way helical junction is still solvent accessible and a conserved water molecule is observed in the highest-resolution crystal structures to form the bridge interaction between the carbonyl groups of U/C74 and U51 (20,95), which may show an induced fit mechanism of binding. Except for the crystal structure, the molecular-docked structure performed by Daldrop et al. also detected the existence of water molecules in ligand-GR complexes (96). All of these works indicated that water molecules may be requisite for the binding of ligands to GR. Based on the above analyses, it is concluded that the changes in the ligandnucleotide interactions induced by multiple mutations enhanced the binding ability of 6AP and HPA to GRs, implying that the aptamer domain activity of GR is extremely plastic and thus readily tunable by nucleotide mutations so as to meet cellular needs, which basically agrees with the experimental study of Stoddard et al. (25).

Influence of mutations on dynamics of GR
To evaluate the influence of mutations on the flexibility of GRs, root-mean-square fluctuations (RMSFs) were computed by using C1 atoms of nucleotides in GR based on the equilibrated MD trajectories ( Figure 6). In the 6AP-GR complexes, the flexibility of nucleotides in GR is heavily affected by mutations ( Figure 6A). Mutations U25A/A46G/C74U and A24U/U25A/A46G/C74U extremely strengthen the flexibility of J23 in GR and slightly enhance the flexibility of P1 in GRs. Mutations U25A/A46G/C74U and A24U/U25A/A46G/C74U generate opposite impact on P2 and L3 of GR. The former one strengthens the fluctuations of P2 and L3, while the latter one highly weakens the flexibility of P2 and L3 in GR. As a whole, the flexibility of the rest in GR is also restrained by mutations U25A/A46G/C74U and A24U/U25A/A46G/C74U, including J12 and J31. Interestingly, J12, J23 and J31 consist of binding pocket of 6AP to GR, thus the changes in the flexibility of J12, J23 and J31 caused by mutations certainly affect the binding ability of 6AP to GR. Similar to the 6AP-GR complexes, mutations A24U, U25A/A46G and A24U/U25A/A46G also give rise to the important effect on the flexibility of GR complexed with HPA ( Figure 6B). Apart from L2 and L3, the flexibility of the rest part of GR is extremely constrained by mutations, especially for J12, J23 and J31 around the binding site of GR, which may regulate the binding ability of HPA to GR. The above analyses suggest that modest sequence alterations in nucleotides have a dramatic impact on flexibility of GR.
To obtain the details involving the influence of mutations on motion modes of GR, the correlation maps of C1 atoms in nucleotides of GR were computed using the program CPPTRAJ in Amber and the results were depicted in Figure 7 and 8 in color-coded modes. The red and yellow reflect strongly correlated motions between specific nucleotides of GR, while the blue and dark blue are indicative of strongly anti-correlated motions between nucleotides. The diagonal regions represent the motion of a nucleotide to itself. It is found that mutations induce obvious changes on the motion modes of GR.
For the HPA-WT GR complex, the structural components J12 and J23 of the WT GR produce strongly anticorrelated motions relative to nucleotides 15-20, while the components P1, P3 and J31 separately form strongly correlated motions relative to nucleotides 20-25, 20-25 and 15-20 ( Figure 8A). It is also observed that strongly correlated movements occur in the regions R1 and R2, and the region R1 involves the structural components J23 and P3, while the region R2 mainly includes the structural component J31. The region R1 reflects the correlated motions of J23 and P3 relative to themselves and the region R2 mostly describes the correlated movements between J31 and nucleotides 20-25 ( Figure 8A). Compared to the WT GR with HPA, mutation A24U not only heavily weakens the correlated motions of the regions R1 and R2, but also the anti-  Figure 8B). Although mutation A24U strengthens the correlated motions between J31 and nucleotides 20-25, A24U also weakens the correlated motions between P1 and nucleotides 20-25 ( Figure 8B). AS shown in Figure 8C, mutations U25A/A46G produce obvious effect on motion modes of GR with HPA. Mutations U25A/A46G not only highly weaken the correlated motions of the regions R1 and R2, but also the correlated motions of P3, J31 and P1 separately relative to nucleotides 20-25, 20-25 and 15-20 compared to the WT GR. In addition, mutations U25A/A46G also obviously weaken the anti-correlated motions of J12 and J23 relative to nucleotides 15-20 ( Figure 8C). According to Figure 8D, the correlated motions in the regions R1 and R2 and that of P3, J31 and P1 separately relative to nucleotides 20-25, 20-25 and 15-20 are extremely decreased by mutations A24U/U25A/A46G compared to the WT GR. Meanwhile, the anti-correlated movements of J12 and J23 relative to nucleotides 15-20 are also heavily reduced due to mutations A24U/U25A/A46G relative to the WT GR.
Based on the above analyses of RMSF, several mutations involved in this study produce obvious effect on the flex-ibility of GR, especially for J12, J23 and J31 around the binding site of ligands. The calculated correlation maps also suggest that mutations highly affect motion modes of GR and change the correlated extent of motions between nucleotides. These results show that modest sequence alterations in nucleotides can exert a dramatic influence on flexibility and dynamic behavior of GR, implying the aptamer domain activity of GR and binding ability can be easily regulated by nucleotide mutations, which basically agrees with the study of Stoddard et al. (25).

DISCUSSIONS
In this work, 1-s molecular dynamics simulations were performed on seven systems involving the WT and mutated GRs complexed with two ligands 6AP and HPA to probe the effect of mutations on dynamics behavior of GRs. In a summary, RMSD of atoms P, O3 , O5 C3 , C4 and C5 in GRs calculated using the equilibrated MD trajectories indicate that all simulated systems reach the equilibrium after 200 ns of MD simulations. The averaged RMSD of seven systems are lower than 4.21Å and the fluctuation range of RMSDs is also <1.13Å, suggesting that the equilibrated MD trajectories can be reliably used to predict the binding affinities, selectivity of two ligands toward the WT and mutated GRs as well as dynamics changes of GRs due to mutations.
The cross-correlation analyses based on the equilibrated MD trajectories show that mutations change motion modes and internal dynamics of GR. Meanwhile, the calculated RMSFs of C1 atoms in nucleotides of GR indicate that mutations highly affect the structural flexibility of GR. To better discuss the impact of mutations on conformations of GRs, cluster analysis was conducted on the last 600 ns of MD trajectories by using the iMWK-Means method, depicted in Supplementary Figures S5 and S6. The results show that mutations induce changes in cluster redistributions (including C1, C2 and C3) of GRs. The superimpositions of structures extracted from main clusters indicate that 6AP and HPA produce obvious slide and rotation in the binding pocket of the mutated GRs relative to the WT GR, which may affect the binding activity of GR. Meanwhile, PC analysis was also run on the equilibrated MD trajectories. Supplementary Figure S7 displays the function of eigenvalues as index of eigenvectors. In the complexes of the mutated GRs with HPA, mutations A24U, U25A/A46G and A24U/U25A/A46G weaken the motion strength of GR compared to the WT GR, while for that of the mutated GRs with 6AP, mutations U25A/A46G/C74U and A24U/U25A/A46G/C74U strengthen the motion of GR relative to the WT GR. The porcupine plots captured by the first eigenvector arising from PC analysis (Supplementary Figures S8 and S9) suggest that the mutations in this work exert significant influence on dynamics behavior of the regions P1, P2, J23 and J31 in GR. Free energy landscapes were built by using projections of MD trajectories on the eigenvectors of the first two principal components PC1 and PC2 and structures corresponding to different free energy basins E1 and E2 were superimposed to understand conformations of GR due mutations (Supplementary Figures  S10 and S11). The results indicate that mutations not only induce the redistribution of free energy basins, but also lead to obvious slide and rotation of two ligands in the binding pocket of the mutated GRs relative to the WT GR, which may tune binding ability of ligands to GR. Therefore, the results of cluster and PC analyses further support the conclusion that the aptamer domain activity of GR is extremely plastic and thus readily tunable by nucleotide mutations.
FEP and MM-GBSA methods were applied to compute binding free energies of 6AP and HPA to the WT and mutated GRs so as to comparatively evaluate influence of mutations on binding ability of ligands to GR. The results not only suggest that binding free energies predicted by two different methods exist encouraging correlations, but also clarify that mutations in the current works strengthen bindings of ligands to GR. Meanwhile, mutations generate important influence on enthalpy changes and entropy changes. As shown in Table 2, the enthalpy change ( H) and entropy change (−T S) of HPA-GR binding is respectively increased by 0.91 kcal/mol and decreased by 0.15 kcal/mol due to A24U compared to the WT GR, totally leading to stronger association of HPA with the mutated GR. By comparison with the WT GR, despite almost no change in binding enthalpy of HPA to GR caused by U25A/A46G, a decrease of 1.69 kcal/mol in binding entropy of HPA to GR induced by U25A/A46G provides main contribution for stronger bindings of HPA to the U25A/A46G mutated GR. Although the entropy change of HPA-GR binding is increased by 0.8 kcal/mol due to A24U/U25A/A46G relative to the WT GR, this unfavorable change is completely compensated by an increase of 6.24 kcal/mol in enthalpy change. By comparison with the 6AP-WT GR, mutations U25A/A46G/C74U and A24U/U25A/A46G/C74U not only increase the enthalpy of 6AP-GR binding, but also decrease the entropy change, thus two mutations highly strengthen binding of 6AP to the mutated GRs. It is observed from Table 2 that binding strength of 6AP to the WT GR is much stronger than HPA, indicating a more favorable binding selectivity of 6AP toward the WT GR over HPA. The comparison of separate free energy components computed by MM-GBSA method shows that the difference in van der Waals interactions and nonpolar solvation energy is small, while the unfavorable polar interaction of 6AP with the WT GR is reduced by 10.38 kcal/mol compared to that of HPA with the WT GR, which is mostly responsible for driving contributions to more favorable selectivity of 6AP toward the WT GR than HPA.
The ligand-nucleotide interaction spectrum calculated by using residue-based free energy decomposition method demonstrate that mutations highly affect interactions of ligands with key nucleotides U22, U51 and C74 in GR, which provides significant contributions for the changes of binding ability of ligands induced by mutations. More importantly, the mutations involved in the current study induce the formation of new additional hydrogen bonding interactions of HPA and 6AP with GR, moreover multiple published works also observed this phenomenon (19,22,25). These changes of ligand-nucleotide interactions should be responsible for certain contributions to conformational changes of GRs due to mutations and play significant role in the plasticity and tuning of the aptamer domain activity of GR. We expect that this study contribute molecular basis and dynamics information for understanding function of GR and possibility as potential drug targets for antibacterial.