Dynamics of metal binding and mutation in yybP–ykoY riboswitch of Lactococcus lactis

Riboswitch is a regulatory segment of messenger RNA (mRNA), which by binding to various cellular metabolites regulates the activity of mRNA via modulating transcription, translation, alternative splicing, and stability of the mRNA. yybP–ykoY riboswitch of Lactococcus lactis, which is present upstream of the yoaB gene, functions as a Mn2+-specific genetic ON-switch, and modulates expression of proteins which are significant for Mn2+ homeostasis. The P1.1 switch helix of the aptamer domain of the riboswitch contains an intrinsic transcription terminator structure, which gets stabilized with Mn2+ binding and causes disruption of terminator structure and allows the continuation of transcription. The current research work involved the evaluation of structural and dynamical properties of the yybP-ykoY riboswitch of L. lactis in its Mn2+-free, Mn2+-bound (wild-type), and Mn2+-bound mutant (A41U) states by applying molecular dynamics simulations. Based on the simulations, the effects of Mn2+ absence and A41U mutation were evaluated on the structure and dynamics of the riboswitches followed by the computation of the free energy of metal binding in the wild-type and the mutant riboswitches. The simulation results provided insights into the properties of the riboswitch with the focus on the dynamics of the P1.1 switch helix, and the manganese binding site designated as MB site, as well as the relative stability of the wild-type and the mutant riboswitches, which helped to understand the structural and dynamical role of the metal ion involved in the function of Mn2+-sensing riboswitch.


Introduction
Riboswitch is a regulatory segment of messenger RNA (mRNA) which by binding to various cellular metabolites regulates activity of mRNA by modulating transcription, translation, alternative splicing and stability of mRNA. 1,2 To date, nearly 40 different classes of riboswitches are known which respond to different ligands along with metabolites, transfer RNAs (tRNAs), enzyme co-factors, signaling molecules and metal ions such as Mn 2+ and Mg 2+ . [1][2][3] They are largely present at 5 0 untranslated regions of protein coding genes unlike thiamine pyrophosphate (TPP)-sensing riboswitch which occurs at 3 0 untranslated region or in introns of target genes. 1,2,4 Riboswitches act as genetic switch which modulate gene expression in distinct fungi, bacteria, plants, and archaea. 1,2 Genes that are regulated by riboswitches encode proteins involved in biosynthesis, catabolism, transport and signaling of cellular metabolites that bind to the riboswitches. The regulation of gene expression by riboswitches involves either transcription termination by the formation of rho-independent transcription terminator hairpin or transcription continuation with the formation of anti-terminator structure. Likewise, gene regulation takes place either at the level of translation by sequestering the ribosome binding site to restrict translation initiation or by exposing the ribosome binding site to enable translation. 1,2 Some riboswitches cleave themselves acting as ribozymes in the presence of ample concentration of metabolites, 5 and some riboswitches, in eukaryotes and prokaryotes, regulate splicing of the pre-mRNA. 1,2,4,6,7 Riboswitches contain two domains; one is aptamer (sensory domain) whose sequence is highly conserved that directly binds with metabolites, and other is expression platform whose sequence greatly varies 8 that as a regulatory domain can either repress or activate gene expression by experiencing structural modications in response to the changes in the aptamer structure. 9 yybP-ykoY riboswitch that is mostly present in plants and human bacterial pathogens for instance Lactococcus lactis, Bacillus subtilis, Escherichia coli, Xanthomonas oryzae, etc., has high sensitivity for manganese (Mn 2+ ) ion. It senses charge, geometry and Lewis acid hardness of Mn 2+ by making inner coordination sphere contacts directly with ve phosphoryl oxygen atoms and one nitrogen atom prototyped as OP1 G9 , OP1 U44 , OP2 C45 , OP2 U39 , OP2 C40 , and N7 A41 in L. lactis. 10 It also regulates distinct genes and is involved in Mn 2+ homeostasis. 10 There are two metal binding sites designated as M A and M B in the riboswitch of L. lactis. Magnesium (Mg 2+ ) ion, owing to its more abundance than Mn 2+ at physiological conditions, and study has not been reported yet, comparing the dynamics of Mn 2+ -bound yybP-ykoY riboswitch regarded as the wild-type, which contains Mg 2+ and Mn 2+ at the M A and the M B site, respectively, A41U Mn 2+ -bounda mutant in which an invariable Mn 2+ -sensing adenosine is replaced with uracil, and Mn 2+free with Mg 2+ at M A site only of yybP-ykoY riboswitches of L. lactis. The aim of this study was to gain deep understanding of the structural and dynamical properties of the Mn 2+ -free, Mn 2+bound and A41U Mn 2+ -bound yybP-ykoY riboswitches of L. lactis at the atomistic details.
In one experimental study, single molecule kinetic analysis of RNA transient structure (SiM-KARTS) assay showed that the anti-terminating P1.1 switch helix in yybP-ykoY riboswitch gets stabilized as Mn 2+ binds to M B site via signal transduction through metal ion sensing core, thus restraining the formation of transcription-terminator hairpin. This stabilizing effect is lost upon mutation of the invariable Mn 2+ -sensing adenosine (A41) with uridine. 14 Through molecular dynamics simulations, we elaborated the metal binding and the mutation effect in L. lactis yybP-ykoY riboswitch as well as found how the stability of P1.1 switch helix is affected in the absence of Mn 2+ . To achieve the objective, we ran 1 ms all-atom MD simulations of the Mn 2+free, Mn 2+ -bound and A41U Mn 2+ -bound riboswitches and came up with the detailed description of the structural and dynamical properties with the focus on the dynamic exibility of the P1.1 switch helix.

Methods
The crystal structure of the aptamer domain of yybP-ykoY riboswitch of L. lactis was retrieved from protein databank (PDB ID: 4Y1I) 11 that was an asymmetric unit containing two identical chains A and B with a total RMS of 0.58Å. 11 Chain A was extracted from the crystal unit, and all atoms except ribonucleotides, 11 Mg 2+ ions, and a Mn 2+ ion, were removed from the chain, which was taken as an initial structure that corresponded to Mn 2+ -bound riboswitch, whereas the Mn 2+ -free riboswitch was constructed by removing Mn 2+ ion from M B site in the metal bound riboswitch and the mutant form of riboswitch was obtained by introducing A41U point mutation in the metal bound riboswitch with the help of Chimera visualization program. 15 Initial conguration of the simulation systems was obtained by using xLEaP module of AMBER20 in which ff99bsc0cOL3 AMBER force eld [16][17][18] was employed for modeling the ribonucleotides, and non-bonded model approach was adopted to treat the metal ions in which the metal center was described as simple Lennard-Jones sphere with an integer charge considering only electrostatic and van der Waals (vdW) terms which were used to describe the interaction between the metal ions and its surroundings. Here, the metal-ligand interactions were described through electrostatic and van der Waals parameters, with s and 3 values of 1.40Å and 0.0168 kcal mol À1 , and 1.36Å and 0.0102 kcal mol À1 , for Mn 2+ and Mg 2+ , respectively. 19 The implementation of the force eld parameters to these three states of the riboswitch was followed by protonation, neutralization with Na + as counter ions and solvation with 12Å buffer of TIP3P water model 20 in a truncated octahedral box with periodic boundary conditions.

Simulation protocol
The initial congurations of all systems (Mn 2+ -free, Mn 2+ -bound and A41U Mn 2+ -bound riboswitches) were initially subjected to energy minimization in a total 1500 steps to eliminate bad intraand inter-molecular contacts between atoms. First, 500 steps minimization was carried out with positional restraint of 25 kcal mol À1ÅÀ2 on all nucleotides including all the ions except water molecules. The positional restraint was then gradually reduced by 5 kcal mol À1ÅÀ2 for every 100 steps followed by 500 steps unrestrained minimization. Aer minimization, each system was gradually heated up to $300 K over 80 ps under NVT (constant particle number, volume and temperature) conditions using Langevin thermostat 21 with constraints applied to bonds involving hydrogen using SHAKE algorithm 22 and the positional restraint with force constant of 25 kcal mol À1 A À2 was applied to all solutes except water molecules. Each system was then equilibrated for 140 ps in NPT ensemble (T ¼ $300 K, P ¼ 1 atm) employing Langevin thermostat and Berendsen barostat. 23 SHAKE algorithm was remained active till the production MD stage to constrain position of hydrogen atoms with force constant of 25 kcal mol À1ÅÀ2 applied on all atoms except water molecules. Aerward, the positional restraint was gradually released in steps with a difference of 5 kcal mol À1ÅÀ2 . Long-range electrostatic interactions were calculated with particle mesh Ewald (PME) 24 algorithm using cut-off distance of 12Å. Finally, each system was subjected to production MD in the NPT ensemble using CUDA version of PMEMD module implemented in AMBER 20 suite of program for 1 ms simulation with time-step of 2 fs at $300 K and 1 atm set for the temperature and pressure, respectively. 25 The relative stability of the wild-type (Mn 2+ -bound) and the mutant (A41U Mn 2+ -bound) riboswitches were estimated in terms of relative binding free energy (DDG bind ) that was calculated as the free energy (DG AB ) difference of the transformation of wild-type form into mutant (A41U) form with Mn 2+ ion at M B site in process 1 and without Mn 2+ ion at M B site in process 2 employing thermodynamic integration (TI) approach (Scheme 1). TI based free energy estimation method computes the free energy difference between the two states A and B, in this case wild-type and mutant riboswitches, by coupling them via a parameter l that serves as an additional, non-spatial coordinate, which allows the free energy difference between the states to be computed as: (1) Scheme 1 : Thermodynamic cycle employed for the computation of relative binding free energy difference between wild type and mutant yybP-ykoY riboswitches.
The Hamiltonian, H that indicates the extent of change came about between the wild-type (state A) and mutant (state B), was calculated as a function of coupling parameter l. TI simulations were performed for each process employing one-step protocol that is based on disappearing one unique residue and appearing other unique residue simultaneously. Atoms of both unique residues were in socore region and were individually specied during calculation, and both van der Waals (vdW) and charge interactions of disappearing or appearing unique atoms with neighboring atoms were illustrated by socore potentials. We used dual topology approach in which both the wild-type and the mutant molecules were considered, simultaneously and atoms of common residues in both states were having same coordinates. To remove redundant bonding term calculations, both the wild-type and the mutant molecules were merged in such a way that the new topology contained only the common atoms and the socore atoms which were atoms of both unique residues. TI simulations were run for each l value and gradient of Hamiltonian H with respect to the coupling parameter, l, was averaged over numerous equilibrated congurations generated during simulations that was then numerically integrated to get free energy difference (DDG bind ) of the two states. In the current study, the free energy simulation protocol employed 6 l points with an array of 0.2 from 0.0 to 1.0, and 1 ns TI simulation was performed for each l point at the isobaric-isothermal ensemble to obtain DDG bind . The simulation systems were minimized at each l point under periodic boundary conditions, employing steepest-descent algorithm prior to TI simulation. Aerward, the systems were heated up to $300 K over 20 ps under NVT conditions using weak-coupling algorithm with constraints applied to bonds involving hydrogen atoms using SHAKE algorithm, and positional restraint with force constant of 5 kcal mol À1ÅÀ2 was also applied to all atoms except water molecules and hydrogen atoms. Finally, each simulation system was simulated in the isobaric-isothermal (NPT) ensemble at each l point for 1 ns with time-step of 2 fs employing CPU version of pmemd module implemented in AMBER 20; 26 the stretching of bonds involving hydrogen atoms was constrained with SHAKE algorithm, temperature at 300 K and pressure at 1 atm were regulated by Langevin thermostat with collision frequency of 2 ps À1 and Berendsen barostat with pressure relaxation time of 2 ps, respectively. Socore potential was employed, and a and b parameters of the socore potential were set to 0.5 and 12Å 2 , respectively. All essential DG values were obtained by numerically integrating vH/vl values collected during the free energy simulations.
The simulation trajectories were primarily analyzed in terms of root mean square deviation (RMSD) of the backbone atoms (P, O5 0 , C5 0 , C4 0 , C3 0 & O3 0 ) and heavy atoms belonging to nucleotides of the M B site employing crystallographic structure as a reference via CppTraj, 27 an analysis tool implemented in AMBERTOOLS 20. Moreover, radius of gyration (R g ), root mean square uctuation (RMSF), and dynamic cross correlation matrix (DCCM) were also evaluated. A correlation plot between R g and RMSD was obtained with the help of pytraj, a python version of CppTraj. Essential dynamics were obtained via principal component analysis (PCA) using the gromacs analysis tools. [28][29][30][31][32] The trajectories were visualized employing VMD program and the images were captured with the help of CHIMERA and PYMOL, 33 along with VMD. 34

Results and discussion
The structural and dynamical properties of the Mn 2+ -free, Mn 2+bound (wild-type) and A41U Mn 2+ -bound (mutant) riboswitches were primarily assessed via time evolution root mean square deviation (RMSD) for the backbone atoms shown in Fig. 2a. For the case of Mn 2+ -free riboswitch, RMSD was stabilized aer 20 ns, whereas stabilization of the RMSD occurred aer 100 and 40 ns for the case of wild-type and mutant riboswitches, respectively. The averaged RMSD was computed to be approximately same for Mn 2+ -free (6.28 AE 0.78Å) and the wild-type riboswitch (6.27 AE 0.81Å), which were also deduced from the probability plots shown in Fig. 2b. On the other hand, the RMSD value for the mutant riboswitch was found to be lower (4.39 AE 0.60Å) than those for the case of other two riboswitches. The very close RMSD values for the case of Mn 2+ -free and wild-type riboswitches indicated that the magnesium ions in the solution diffused to occupy the M B site, thus stabilizing the RMSD prole of Mn 2+ -free riboswitch as observed in case of the wild-type. The coordination of the magnesium ion to the Mn 2+ -free riboswitch was also visualized in the simulation trajectory via VMD. The similar structure behaviour of Mn 2+ and Mg 2+ toward the riboswitch was inferred as both ions have same charge, Lewis acid hardness and octahedral coordination geometry. 35 The RMSD prole of the mutant riboswitch was signicantly lower than that of the wild-type riboswitch thus indicating higher stability of the mutant riboswitch compared to the latter one as well as reecting that the mutation induced global structural changes to the riboswitch. Fig. 2c illustrated the RMSD proles of the M B site nucleotides in Mn 2+ -free, wild-type and mutant riboswitches. In case of Mn 2+ -free riboswitch, the RMSD prole exhibited more uctuation compared to the wild-type case, even though the averaged RMSD of the Mn 2+ -free riboswitch (0.59 AE 0.09Å) was close to the wild-type case (0.56 AE 0.06Å) but lower than that of the case of mutant riboswitch (0.78 AE 0.08Å). The variation in RMSDs pointed toward a contrast dynamical behaviour of the two ions, however, both were reportedly having very similar structural properties. 35 Moreover, the variation in the backbone RMSD was lower in both cases (Mn 2+ -free and wild-type riboswitches) compared to that of the mutant case thus suggesting that the overall dynamics of the M B site in Mn 2+ -free riboswitch was lower than that of the mutant riboswitch, however it was conversely observed for the M B site RMSD as the mutant riboswitch case demonstrated high uctuations than that of the wild-type thus pointing toward the lower stability of the metalbinding site due to mutation. The contrasting behaviour was attributed to the mutation in the riboswitch that increased dynamic exibility of the metal-binding site but caused the overall structure of the riboswitch less exible. The RMSD data thus reected that the dynamics of yybP-ykoY riboswitch was altered with the mutation and absence of manganese ion, which therefore led to further evaluate the effect of mutation and Mn 2+ absence.
Time evolution radius of gyration (R g ) of the backbone atoms of the three riboswitches was computed to obtain information of the overall dynamics in the absence and presence of the metal ion, and due to the mutation shown in Fig. 3a. R g gives estimation of a system's compactness as the higher R g value indicates lower compactness. 36 The average R g value for the Mn 2+ -free riboswitch (24.82 AE 0.35Å) was very close to that of the wild-type (24.64 AE 0.35Å) but higher than the mutant riboswitch (23.94 AE 0.27Å), which indicated that the volume of the Mn 2+ -free and the wild-type riboswitches was indifferent demonstrating that Mn 2+ absence has little effect on the overall structure of the riboswitch, however, their volume was slightly higher compared to the mutant riboswitch. Fig. 3b displayed the R g vs. RMSD correlation plots; in case of Mn 2+ -free riboswitch, distribution of the sampled conformers was located at RMSD from 4.05 to 8.05Å corresponding to R g value between 24.0 and 25.7Å, whereas for the wild-type case, the corresponding distribution showed similar RMSD from 4.01 to 8.29Å and R g from 25.69 to 23.8Å. In the case of the mutant riboswitch, RMSD appeared between 2.99 and 5.60Å which corresponded to R g value from 23.29 to 24.63Å. The correlation plots exhibited that the dynamics of the Mn 2+ -free and wild-type riboswitches yielded to be similar, whereas the mutant case demonstrated distinct dynamics than the other two riboswitches. 37 Effect of Mn 2+ absence and A41U mutation on the dynamics of the riboswitch, in particular the P1.1 switch helix, was analyzed by root mean square uctuation (RMSF), which facilitated to understand the exibility of different nucleotide regions during simulations. 38 Fig. 4a showed the RMSF proles of the Mn 2+ -free, the wild-type and the mutant riboswitch; in the case of Mn 2+ -free riboswitch, residual uctuations were higher than that of the wild-type riboswitch except for U66, G67, U68 and U86 nucleotides, and the M B site nucleotides (G9, U39, C40, A41, U44 and C45) exhibited high uctuations depicted in Fig. 4d-f. The higher exibility of different regions in Mn 2+ -free riboswitch compared to that in the wild-type was attributed to the absence of the manganese ion which caused the riboswitch more exible. A similar uctuation trend was observed for the P1.1 switch helix comprising of A5, G6, G7, G8, C96, C97, U98, U99, U100 and C101 nucleotides in the Mn 2+ -free riboswitch demonstrating high exibility as deduced from the high intensity uctuations in the RMSF proles shown in Fig. 4b and c. This was interpreted as manganese and magnesium ions reportedly have distinct dynamical properties therefore, the binding of Mg 2+ ion to the Mn 2+ -free riboswitch was assumed altered thereby destabilizing the M B site and the P1.1 switch helix. The stabilization of these regions can only be achieved only if the suitable metal ion like Mn 2+ essentially binds to the metal-binding site as in the wild-type riboswitch, since according to an experimental study, the binding of Mn 2+ ion disrupts transcription terminator hairpin structure and thus stabilizes anti-terminating P1.1 switch helix leading to transcription continuance. 11,20 The RMSF prole of the mutant riboswitch showed mixed residual uctuations in comparison to the wild-type case as shown in Fig. 4a, and the M B site residues in particular G9 and U41 nucleotides showed higher uctuations, and U39 and C40 nucleotides showed lower uctuations, whereas U44 and C45 nucleotides possessed similar uctuations as observed in case of wild-type riboswitch (Fig. 4df). Nucleotides of the P1.1 switch helix, in the mutant riboswitch showed distinct uctuation pattern compared to the wild-type, since A3 and A4 nucleotides produced low uctuations compared to G6, G7, U100, C101 and C102 nucleotides having high uctuations, whereas A5, C96 C97 U98 and U99 nucleotides possessed similar uctuations like those of the wild-type riboswitch ( Fig. 4b and c). This implied that the mutation had mixed effect on the overall structure exibility of the riboswitch, including the M B site and the P1.1 helix.  It was therefore observed that the mutation perturbed the M B site despite an overall structural stabilization occurred in the mutant riboswitch. To evaluate the contrast in the structure and dynamics of the wild-type riboswitch in comparison to other two riboswitches, comparative dynamic cross-correlation matrix (DCCM) analysis was carried out for Mn 2+ -free, wild-type and mutant riboswitch to gure out the correlation or anticorrelation between nucleotides of the M B site and other regions, in particular the P1.1 switch helix (Fig. 5). DCCM plots depicted correlated movements between residues belonging to different regions of the riboswitches which were observed throughout 1 ms trajectories of each riboswitch. In DCCM plot, correlation range was from +1 to À1 representing with variation in the color intensity; red-colored contours correspond to positively correlated movements, white-colored regions represent zero correlation or uncorrelated movements and bluecolored contours signify negatively or anti-correlated movements. The characteristic features in the DCCM plots were depicted as boxed regions labelled as 1 to 16 showing effects of Mn 2+ absence and A41U mutation on the structure and dynamics of the riboswitch. In the case of Mn 2+ -free riboswitch, region 8, and 3 and 5 corresponded to positively correlated and zero correlated motions, respectively (Fig. 5a), whereas the presence of Mn 2+ ion in the wild-type riboswitch caused the riboswitch to attain mostly positively correlated motions (Fig. 5b). The A41U mutation was demonstrated to affect the atomic motions causing the mutant riboswitch to adopt distinct motions for instance regions 12, 13, 15 and 16 corresponded to zero correlation (Fig. 5c). Moreover, a contrast was observed between Mn 2+ -free and wild-type cases, as regions 1 and 7 in the Mn 2+ -free case depicted positively correlated motions compared to the corresponding region with zero correlated motions. Similarily, a constrast was observed for the case of mutant riboswitch when comparing with the wild-type case; regions 7, and 2, 9, 14 and 11 in the mutant case corresponded to positively and negatively correlated motions, respectively, whereas Mn 2+ -free case had mixed correlations consisting of both positive and negative correlations for the atomic motions corresponding to respective regions 1 and 7, and 2, 4, 6, 9 and 10. Nevertheless, the difference in the DCCM plots was too small to observe the obvious effect of Mn 2+ -absence and A41U mutation on the correlation dynamics among different regions, particularly, between M B site nucleotides (G9, U44, A41, C45, U39 and C40) and nucleotides of P1.1 switch helix (A3, A4, A5, G6, G7,  C96, C97, U98, U99, U100, C101 and C102) of the riboswitches, which provoked to perform essential dynamics (ED) or principal component analysis (PCA) in order to evaluate distinct dynamics of the riboswitches.
PCA reduces the dimensionality of huge data set and helps to identify signicant motions which are incumbent for the biological activity of a molecule. 39 As the essential motions are extracted from sampled conformations, the method of applying PCA to a simulated trajectory is called essential dynamics (ED). [40][41][42] With the removal of roto-translation motions of all atoms as of the nucleic acid, conformations are collected throughout the MD trajectory. Consequently, covariance matrix is set up that upon diagonalization of system's atomic coordinates targets the signicant motions of conformation which are accountable for the mass variance in atomic positions as eigenvalues. Eigenvectors (essential modes) are determined by ranking eigenvalues of principal components (PCs) represented by scree plot shown in Fig. 6. Aerward, essential conformational space is determined by projecting the original conguration of the system on the PCs. 39 Essential subspace comprises a few degree of freedom, in which positional uctuations occur an-harmonically, therefore, essential degree of freedom points out the motions relevant to the function of biomolecules. 42 Here, eigenvalues from covariance matrix were constructed from all-atom coordinates of the Mn 2+ -free, the wild-type and the mutant riboswitches which are shown in descending order. Eigenvalues of the Mn 2+ -free case showed steeper decrease than those of the wild-type and the mutant cases, and the wild-type case also showed steeper decrease than the mutant case. Trace of the covariance matrix of the Mn 2+ -free case happened to be 332.71 nm 2 whereas for the wild-type and mutant cases, the values yielded to be 203.66 and 205.39 nm 2 , respectively. Since the high values of trace of the covariance matrix generally pointed toward more exible system, 39 therefore, high value of trace of covariance in case of Mn 2+ -free riboswitch demonstrated high dynamic exibility due to Mn 2+ absence in comparison to other two riboswitches, and the mutant case with slightly high trace of covariance matrix was assessed to have slightly more exible structure compared to the wild-type riboswitch. Moreover, large eigenvalues for the case of Mn 2+free riboswitch suggested to be less compact structure compared to the wild-type and the mutant riboswitches, thus giving a more clear description of the volume of Mn 2+ -free and wild-type riboswitches as described earlier by R g calculations (see Fig. 3a). The rst two eigenvectors were selected for the evaluation of essential dynamics of the riboswitches since they corresponded to the highest eigenvalues demonstrated by visible kink in the scree plot (Fig. 6).
Conformational sampling of Mn 2+ -free, wild-type and mutant riboswitches in the essential subspace was depicted in Fig. 7, illustrating the conformations along the two eigenvectors (EV1 and EV2) projected by all-atoms. It was observed that Mn 2+free riboswitch occupied larger and distinct conformational subspace than wild-type riboswitch, however, atomic positions in the Mn 2+ -free case was also overlapping with that of the wildtype case showing some differences in the sampling of Mn 2+free and wild-type riboswitches, thus pointing toward the change in the conformational space of the riboswitch due to Mn 2+ ion (Fig. 7a). On the other hand, the mutant riboswitch occupied smaller conformational subspace than the wild-type riboswitch and atomic positions nearly overlapped with that of the wild-type riboswitch demonstrating that the mutation stabilized the overall structure of riboswitch and also decreased conformational subspace occupied by the mutant riboswitch (Fig. 7b). To further ascertain the mutation effect on the riboswitch dynamics in particular the M B site and the P1.1 switch helix, atomic uctuations of nucleotides, and dynamic correlation between the M B site and the P1.1 switch helix, accountable for the essential subspace dynamics were analyzed via RMSF calculations along EV1 during PCA, since it reduced the functionally irrelevant motions and unmasked functionally signicant motions (Fig. 8). The PCA-based RMSF helped to gain further insight into the essential dynamics of the riboswitch that was not achieved earlier by simple RMSF calculations and DCCM analysis (see Fig. 4 and 5). Fig. 8a showed that the mutation decreased the exibility of the overall structure of the riboswitch, as low uctuations were observed in the RMSF prole of the mutant riboswitch thus revealing that the mutation stabilized the overall structure which also correlated the backbone RMSD of the mutant riboswitch showing enhanced structural stability (see Fig. 2a).
PCA based RMSF analysis also demonstrated low uctuations of the M B site nucleotides in the mutant riboswitch such as G9, U39, U44 and C45 except C40 and U41 depicted in Fig. 8d-f. The RMSF prole also revealed that the mutation affected the overall dynamics of the riboswitch but not the exibility of the M B site nucleotides rather than make the metal site more rigid. On the other hand, all nucleotides (A4, A5, C96, C97, U98, U99, U100, C101 and C102) except A3 in the P1.1 switch helix of the mutant riboswitch showed relatively high uctuations depicted in Fig. 8b and c, demonstrating that the mutation increased the helix exibility. This nding complemented an experimental study reporting Mn 2+ binding at the M B site that stabilized P1.1 switch helix in yybP-ykoY riboswitch, and this stabilizing effect was lost upon Mn 2+ sensing adenosine mutation. 14 The rst two principal components, PC1 and PC2 of each Mn 2+ -free, wild-type and mutant riboswitch were employed to obtain information of the global minimum-energy (stable) conformations responsible for the essential dynamics, which were depicted in the form of Gibbs free energy landscapes (FEL) shown in Fig. 9. The lowest energy states or stable conformations represented by blue colored maps whereas red colored maps corresponded to high-energy conformations. Each riboswitch produced different FEL pattern; the Mn 2+ -free case showed two free energy basins (valleys) with pointed and at ends, which were less broad than the single energy basin observed in the wild-type case ( Fig. 9a and b). The minima with at ends corresponded to clusters of stable conformations, the minima with pointed end represented a stable conformation, and energy barriers were large between basins unlike between the minima in basins. Since narrow energy valleys indicated small number of conformations, the Mn 2+ -free case yielded less number of stable conformations in comparison to the wild-type riboswitch. FEL of the mutant riboswitch depicted broader minima with little energy barrier than that of the wild-type riboswitch; minima with at ends indicated clustering of the stable or low energy conformations, and a canonical end suggested the presence of a stable conformation. In the case of mutant riboswitch, large number of stable conformations were produced than the wild-type riboswitch as reected by FEL pattern (Fig. 9c). The stable conformations indicated by the minima were responsible for essential dynamics in the riboswitch, and as the number of stable conformations decreases exibility of the system increases and vice versa. Consequently, FEL pattern served to gure out the effect of the mutation and Mn 2+ absence on the sampled essential subspace.
The coordination sphere of the metal binding M B site was also analyzed by evaluating dynamics of its geometric conguration in terms of bond lengths and bond angles. Fig. 10 illustrated time evolution coordination distances between the metal ion and the coordinating atoms of the nucleotides in the Mn 2+ -free, the wild-type and the mutant riboswitches. Based on the distance plots, the averaged coordination distances were computed for the metal-nucleotide atomic pairs in each riboswitch and listed in Table 1. In case of Mn 2+ -free riboswitch, distances between Mg 2+ and the M B site residues were averaged to be less than 2Å throughout the simulation time except Mg-N7 A41 distance (2.31 AE 0.11Å) that was same as the Mn-N7 A41 distance (2.30 AE 0.09Å) in the wild-type riboswitch, thus demonstrating the similar structure behavior of these two metal ions. Nevertheless, based on the distance plots, the average coordination distances were computed between the metal ions and nucleotides of the M B site in each riboswitch. In case of Mn 2+ -free riboswitch, distances between Mg 2+ ion and the M B site residues were averaged to be less than 2Å throughout the simulation time except the Mg-N7 A41 distance (2.31 AE 0.11Å) that was same as the Mn-N7 A41 distance (2.30 AE 0.09Å) in the wild-type riboswitch. However, the amplitude of these two distance plots had a characteristic constrast as the distance in Mn 2+ -free riboswitch was slightly larger than the wild-type riboswitch (Fig. 10a), reecting less stability of the M B site in the absence of essential Mn 2+ ion. Likewise, a relatively large amplitude was observed for Mg-OP1 G9 distance (1.96 AE 0.05Å) compared to Mn-OP1 G9 distance (2.02 AE 0.05Å) indicating less stability of the metal coordination site even though the distance in the Mn 2+ -free riboswitch was low than the wild-type case (Fig. 10b). Other metal-nucleotide distances in Mn 2+ -free riboswitch were also smaller than those of the wild-type riboswitch but the Mg-nucleotide distances had similar amplitudes (Fig. 10c), suggesting its rigidity than those in the wild-type riboswitch. Based on the distance plots, it was evaluated as the Mg 2+ ion unlike Mn 2+ ion showed constant variations while making interactions with the M B site nucleotides, whereas Mn 2+ being an essential constituent of the riboswitch as of the wild-type, forms dynamically stable metal coordination site which could be attributed to the function of the Mn 2+ -sensing riboswitches. The rigidity of the metalnucleotide distances in Mn 2+ -free riboswitch also corroborated the RMSF data, demonstrating structural exibility of the riboswitch (see Fig. 4d-f), however, irrespective of very similar RMSD values for the M B site nucleotides in Mn 2+ -free and wild-type riboswitches. The average distances between Mn 2+ and M B site residues in the wild-type riboswitch were very close as that of the mutant riboswitch as listed in Table 1, except Mn-N7 A41 distance (2.30 AE 0.09Å) in the wild-type that was larger than Mn-O4 U41 distance (2.12 AE 0.09Å) in the mutant form shown in Fig. 10d, in which U41 nucleotide was not having N7 atom. This indicated that the mutation led to Fig. 10 (a) Distance between N7 A41 and Mg 2+ , and N7 A41 and Mn 2+ of Mn 2+ -free and wild-type riboswitches, respectively, (b) distance between OP1 G9 and Mg 2+ , and OP1 G9 and Mn 2+ of Mn 2+ -free and wild-type riboswitches, respectively, (c) distance between Mg 2+ and OP1 U44 , OP2 U39 , OP2 C45 , OP2 C40 of Mn 2+ -free, and Mn 2+ and OP1 U44 , OP2 U39 , OP2 C45 , OP2 C40 of wild-type riboswitch, (d) distance between O4 U41 and Mn 2+ of mutant riboswitch, and N7 A41 and Mn 2+ of wild-type riboswitch, (e) distance between OP1 U44 and Mn 2+ of mutant riboswitch, and OP1 U44 and Mn 2+ of wild-type riboswitch, (f) M B site of mutant and wild-type riboswitches. shorter and rigid contacts between Mn 2+ and M B site residues which also correlated the PCA results signifying that the exibility of M B site decreases with mutation (see Fig. 8d-f). Fig. 10d depicted a fragile coordination distance of Mn 2+ with O4 atom of uracil in the mutant riboswitch that was formed at the beginning of the simulation time as deduced from the distance plot for the rst 5 ns, which upon equilibration at $1 ns got decreased thus causing the coordination shell to be stabilized, and this nding was found to be correlating with an experimental study 11 (Fig. 10f). According to the experimental study, the A41U mutation led Mn 2+ coordination scheme to adopt certain changes in such a way that the coordination scheme of M A site nucleotides and G10, U39, C40, and C45 of M B site remained unaffected, whereas the backbone phosphate atoms of U44 was shied away and replaced by a water molecule as resolved by 2.2Å electron density. 11 Conversely, the simulation study showed that the contact between OP1 atom of U44 nucleotide and Mn 2+ remained stable and was not substituted by water as demonstrated by the distance plots illustrated in Fig. 10e and f. Aer having analyzed the effects of Mn 2+ absence and the mutation in terms of structure and dynamics of yybP-ykoY riboswitch, and in particular, the exibility (stability) of the P1.1 switch helix and the M B site, the A41U mutation effect was quantied as the binding free energy of the mutant riboswitch. The free energy was estimated with respect to the Mn 2+ ion bound to the wild-type and the mutant riboswitches through relative binding free energy difference (DDG bind ) via free energy simulations based on TI approach, 43 in order to evaluate the relative stability of the two Mn 2+ -bound riboswitches. Fig. 11 illustrated the plot of relative free energy for each interval of l that is between adjacent Hamiltonian for both riboswitches, where process 1 illustrated binding free energies in the presence of Mn 2+ , computed as a function of coupling parameter, l, whereas binding energy calculation (without Mn 2+ ) was carried out in process 2 to obtain useful information employing non-physical thermodynamic cycle (Scheme 1) for the relative binding free energy difference based on the following relation, Table 2 listed the free energies computed during process 1 (DGÃ B 1 ) and process 2 (DGÃ B 2 ), and the relative binding free energy difference (DDG bind ) revealed that the mutant riboswitch was more stable than the wild-type riboswitch with relative binding energy difference of 23.75 kcal mol À1 . The increase in the binding free energy for the mutant riboswitch was attributed to the more rigid contact between O4 U41 and Mn 2+ compared to the Mn-N7 A41 coordination bond in the wild-type riboswitch (Table 1 and Fig. 10). Based on the free energy data, it was obvious that the A41U mutation stabilizes A41U yybP-ykoY riboswitch of L. lactis. Fig. 11 Relative binding free-energy between mutant and wild-type riboswitches as a function of l (coupling parameter) employing thermodynamic integration approach. Process 1 (and 2) outlines free-energy difference for the conversion of wild-type system into mutant system with Mn 2+ (and without Mn 2+ ) at binding site M B .

Conclusion
Based on the simulation data, the presence of Mn 2+ ion was signied as the ion was essentially required for the proper functioning of yybP-ykoY riboswitch of L. lactis. For this purposes, MD simulations were applied to the metal-sensing riboswitches employing the most widely used ff99bsc0cOL3 AMBER force eld to model and describe RNA topology despite having several reports that the force eld was unable to reproduce some tetranucleotides and tetraloops for the conformational features and NMR experimental data, and therefore, a number renements and optimizations of the force eld were suggested. 44 Nevertheless, the data obtained from the simulation was found to be in good agreement with the experimental data reported so far. Furthermore, the application of nonbonded parameter for the manganese and magnesium ions was also found to be compatible as the structural and dynamical properties of the metal-binding sites reproduced well when compared with experimental data based on X-ray crystallography. 11 The absence of Mn 2+ ion at the M B site of the Mn 2+ -free riboswitch was evaluated as the structurally similar Mg 2+ ion from the solution attempted to mimic the functional role of the essentially important Mn 2+ ion in the riboswitch but the dynamical properties of the two ions were reported to be distinct which were suggested to be associated to their different functional role in the biological molecules. The simulation results thus yielded that the presence of Mn 2+ ion was essential for the Mn 2+ -sensing riboswitch which in the absence of Mn 2+ ion undergo unfavourable structural and dynamical modications. Moreover, the effect of the A41U mutation in the riboswitch was examined and it was found that the mutation affected not only the overall riboswitch structure but also the metalbinding site and the P1.1 switch helix, which undergo structural and conformational modications due to mutation. The evaluation of the structure and dynamics related to the metalbinding site and the P1.1 switch helix obtained from the simulations was expected to serve further in the future studies with focus on deciphering the functional role of the Mn 2+sensing riboswitch which in turn could be helpful in unravelling how the metal-sensing riboswitch could be considered for therapeutic purposes.

Conflicts of interest
There are no conicts to declare.