Structural Studies of the 3′,3′-cGAMP Riboswitch Induced by Cognate and Noncognate Ligands Using Molecular Dynamics Simulation

Riboswtich RNAs can control gene expression through the structural change induced by the corresponding small-molecule ligands. Molecular dynamics simulations and free energy calculations on the aptamer domain of the 3′,3′-cGAMP riboswitch in the ligand-free, cognate-bound and noncognate-bound states were performed to investigate the structural features of the 3′,3′-cGAMP riboswitch induced by the 3′,3′-cGAMP ligand and the specificity of ligand recognition. The results revealed that the aptamer of the 3′,3′-cGAMP riboswitch in the ligand-free state has a smaller binding pocket and a relatively compact structure versus that in the 3′,3′-cGAMP-bound state. The binding of the 3′,3′-cGAMP molecule to the 3′,3′-cGAMP riboswitch induces the rotation of P1 helix through the allosteric communication from the binding sites pocket containing the J1/2, J1/3 and J2/3 junction to the P1 helix. Simultaneously, these simulations also revealed that the preferential binding of the 3′,3′-cGAMP riboswitch to its cognate ligand, 3′,3′-cGAMP, over its noncognate ligand, c-di-GMP and c-di-AMP. The J1/2 junction in the 3′,3′-cGAMP riboswitch contributing to the specificity of ligand recognition have also been found.


Introduction
Riboswitch, located at the 5 untranslated region of mRNAs, are the gene regulatory elements, which can control gene expression induced by small-molecule ligand binding [1][2][3][4][5][6]. So far, more than 20 riboswitch families have been found, due to different types of small-molecule ligands [5,. Riboswitches have been proposed as antibiotic drug targets [30], which have attracted structural biologists and biophysicists great attention.
To get a better understanding of the conformational change for the 3 ,3 -cGAMP riboswitch induced by the cognate and noncognate ligands binding, molecular dynamics simulations were employed on the apo 3 ,3 -cGAMP riboswitch, the 3 ,3 -cGAMP riboswitch bound by the 3 ,3 -cGAMP molecule and by the c-di-GMP and c-di-AMP molecules. Our main goals were (1) to address the apo 3 ,3 -cGAMP riboswitch structural feature at the nucleotide level; (2) to explain specificity of ligand recognition of the 3 ,3 -cGAMP riboswitch.

Results
We monitored the root-mean-square deviation (RMSD) values with respect to the average structure for the apo riboswitch, riboswitch + 3 ,3 -cGAMP, riboswitch + c-di-GMP, riboswitch + c-di-AMP and G20A-GEMM-I + 3 ,3 -cGAMP systems ( Figure S1). The small RMSD fluctuations of one simulation indicate that the system has attained equilibrium. It can be seen that the five systems had the small fluctuations of RMSD values after 20 ns, which indicate that all the studied systems have reached equilibrium after 20 ns, and the energies of these systems were also stable during the equilibrated simulation. The nucleotides that have atoms within 5 Å of 3 ,3 -cGAMP (or c-di-GMP) in the riboswitch + 3 ,3 -cGAMP model (or the riboswitch + c-di-GMP model) were defined as the binding pocket. The RMSD values, solvent accessible surface area (SASA) and radius of gyration (Rg) of the binding pocket for the apo riboswitch, riboswitch + 3 ,3 -cGAMP, riboswitch + c-di-GMP and riboswitch + c-di-AMP systems are measured (Figure 2A-C). The RMSD values of the binding pocket in the riboswitch + 3 ,3 -cGAMP model are relatively stable throughout the trajectory. The increase of the RMSD values of the binding pocket in the riboswitch + c-di-GMP, riboswitch + c-di-AMP and apo riboswitch systems predict the binding pocket has a conformational change occurring. The SASA and Rg values of binding pockets in these four systems verify the occurrence of the conformational changes, which indicate that the binding pocket in the apo riboswitch has a smaller binding pocket compare to that in the ligand-bound state. Computational techniques, and especially molecular dynamics (MD) simulations, play an important role in investigating structure and dynamics of macromelecules at the atomistic level [41][42][43][44]. MD simulations have been carried on the adenine riboswitch [45,46], the SAM-I, SAM-II riboswitches [47][48][49][50][51][52][53], the Glms riboswitch [54,55], the preQ1 riboswitch [56][57][58][59][60][61][62], the guanine riboswitch [63,64], and the c-di-GMP class I GEMM riboswitch studied in our team [65], which are useful for providing interpretations for existing experimental results and producing predictions for novel experiments on these riboswitches [66]. However, the computational investigation on the 3′,3′-cGAMP riboswitch are very limited so far.

Results
We monitored the root-mean-square deviation (RMSD) values with respect to the average structure for the apo riboswitch, riboswitch + 3′,3′-cGAMP, riboswitch + c-di-GMP, riboswitch + c-di-AMP and G20A-GEMM-I + 3′,3′-cGAMP systems ( Figure S1). The small RMSD fluctuations of one simulation indicate that the system has attained equilibrium. It can be seen that the five systems had the small fluctuations of RMSD values after 20 ns, which indicate that all the studied systems have reached equilibrium after 20 ns, and the energies of these systems were also stable during the equilibrated simulation. The nucleotides that have atoms within 5 Å of 3′,3′-cGAMP (or c-di-GMP) in the riboswitch + 3′,3′-cGAMP model (or the riboswitch + c-di-GMP model) were defined as the binding pocket. The RMSD values, solvent accessible surface area (SASA) and radius of gyration (Rg) of the binding pocket for the apo riboswitch, riboswitch + 3′,3′-cGAMP, riboswitch + c-di-GMP and riboswitch + c-di-AMP systems are measured (Figure 2A-C). The RMSD values of the binding pocket in the riboswitch + 3′,3′-cGAMP model are relatively stable throughout the trajectory. The increase of the RMSD values of the binding pocket in the riboswitch + c-di-GMP, riboswitch + c-di-AMP and apo riboswitch systems predict the binding pocket has a conformational change occurring. The SASA and Rg values of binding pockets in these four systems verify the occurrence of the conformational changes, which indicate that the binding pocket in the apo riboswitch has a smaller binding pocket compare to that in the ligand-bound state.

Energy Calculations Revealed the Preferential Binding of the 3 ,3 -cGAMP Riboswitch to 3 ,3 -cGAMP over c-di-GMP and c-di-AMP
To compare the competitive binding characteristics of the 3 ,3 -cGAMP riboswitch to the cognate ligand, the 3 ,3 -cGAMP molecule, and to the noncognate ligand, the c-di-GMP and c-di-AMP molecules, MM-GBSA methods were used to calculate the binding free energies for the riboswitch + 3 ,3 -cGAMP, riboswitch + c-di-GMP and riboswitch + c-di-AMP systems ( Table 1). The binding of the 3 ,3 -cGAMP riboswitch to the 3 ,3 -cGAMP molecule with the binding free energy of −11.62 kcal/mol in the riboswitch + 3 ,3 -cGAMP system is more energetically favorable over the binding of the 3 ,3 -cGAMP riboswitch to the c-di-GMP molecule with that of −8.74 kcal/mol, and the binding of the 3 ,3 -cGAMP riboswitch to the c-di-AMP molecule with the binding free energy of 13.31 kcal/mol in the riboswitch + c-di-AMP system indicate the c-di-AMP molecule can't bind to the 3 ,3 -cGAMP riboswitch, which agrees with the experiment result [39]. The ranking of binding energy predicts that the 3 ,3 -cGAMP riboswitch binds preferentially to its cognate ligand, the 3 ,3 -cGAMP molecule, over its noncognate ligand, the c-di-GMP and c-di-AMP molecules. In addition, the binding of the 3 ,3 -cGAMP molecule to the class I GEMM riboswitch with G20A mutation with the binding free energy of −7.64 kcal/mol, which agrees with experiment result analyzed by the dissociation constants K d [39] (Table S1). Table 1. MM-GBSA free energy (kcal·mol −1 ) components for the riboswitch + 3 ,3 -cGAMP and riboswitch + c-di-GMP models.

Interaction Analysis between the 3 ,3 -cGAMP Riboswitch and Ligands
To address the interactions of the 3 ,3 -cGAMP riboswitch with the cognate ligand (i.e., the 3 ,3 -cGAMP molecule) and with the noncognate liand (i.e., the c-di-GMP molecule), the percentages of all possible hydrogen bonds between the 3 ,3 -cGAMP riboswitch and the 3 ,3 -cGAMP molecule for the riboswitch + 3 ,3 -cGAMP system, and between the 3 ,3 -cGAMP riboswitch and the c-di-GMP molecule for the riboswitch + c-di-GMP system were analyzed from the equilibrated MD simulations, with the error bar analyzed by the block averaging [67,68] (Table 2). The hydrogen bond interaction was defined as follow: The distance between donor and acceptor was shorter than 3.5 Å, and the angle among donor, proton and acceptor was greater than 120 • [69,70]. The are some varied hydrogen bonds between the 3 ,3 -cGAMP riboswitch and the 3 ,3 -cGAMP molecule for the riboswitch + 3 ,3 -cGAMP model, and between the 3 ,3 -cGAMP riboswitch and the c-di-GMP molecule for the riboswitch + c-di-GMP model, namely, the nucleotides G8-C15 around the J1/2 junction and G40-A42 in the J2/3 junction ( Table 2). The hydrogen bond occupancies between nucleotides G8-C15 around the J1/2 junction of the 3 ,3 -cGAMP riboswitch and the 3 ,3 -cGAMP molecule increased 316.58%, compared to those between the 3 ,3 -cGAMP riboswitch and the c-di-GMP molecule in the riboswitch + c-di-GMP model, which might cause the stronger interaction of the 3 ,3 -cGAMP riboswitch with the 3 ,3 -cGAMP molecule. At the same time, the occupancies of hydrogen bond between nucleotide G40-A42 in the J2/3 junction of the 3 ,3 -cGAMP riboswitch and the 3 ,3 -cGAMP molecule in the riboswitch + 3 ,3 -cGAMP model were reduced by 79.08% compared to those between the 3 ,3 -cGAMP riboswitch and the c-di-GMP molecule. Overall, however, the hydrogen bond occupancies have an increase in the riboswitch + 3 ,3 -cGAMP model, which suggest the stronger interaction of the 3 ,3 -cGAMP riboswitch with the 3 ,3 -cGAMP molecule compare with the c-di-GMP molecule, which is in good agreement with the binding energy analysis discussed above. Table 2. The occupancies (%) of hydrogen bonds with the error bars between the 3 ,3 -cGAMP riboswitch and the 3 ,3 -cGAMP molecule for the riboswitch + 3 ,3 -cGAMP model, and between the 3 ,3 -cGAMP riboswitch and the c-di-GMP molecule for the riboswitch + c-di-GMP model.
In addition to hydrogen bond interaction, the base-stacking interactions between A α , G β of the 3 ,3 -cGAMP molecule and the 3 ,3 -cGAMP riboswitch in the riboswitch + 3 ,3 -cGAMP model, and between G α , G β of the c-di-GMP molecule and the 3 ,3 -cGAMP riboswitch in the riboswitch + c-di-GMP model have been detected (Figure 4 and Figure S2). The criteria for the base-stacking was defined as follow: The base-base distance ≤4.0 Å, the base-base angle ≤ 30 • and overlap with each other [71]. As shown in Figure 4 and Figure S2, the base of G40, A α , base of A41, G β , base of G8 in the riboswitch + 3 ,3 -cGAMP model and the base of G40, G α , base of A41, G β , base of G8 in the riboswitch + c-di-GMP model are plane-parallel approximation and overlap between the bases, which meet the angle criteria for base stacking. The percentage of distance ≤4.0 Å between base of G8 and G β of 3 ,3 -cGAMP, base of A41 and A α of 3 ,3 -cGAMP increased 45.54%, 30.58% compared to that between base of G8 and G β of c-di-GMP, A41 and G α of c-di-GMP, which also suggest the stronger interaction of the 3 ,3 -cGAMP riboswitch with the 3 ,3 -cGAMP molecule compare with the c-di-GMP molecule. In addition to hydrogen bond interaction, the base-stacking interactions between Aα, Gβ of the 3′,3′-cGAMP molecule and the 3′,3′-cGAMP riboswitch in the riboswitch + 3′,3′-cGAMP model, and between Gα, Gβ of the c-di-GMP molecule and the 3′,3′-cGAMP riboswitch in the riboswitch + c-di-GMP model have been detected (Figures 4 and S2). The criteria for the base-stacking was defined as follow: The base-base distance ≤4.0 Å, the base-base angle ≤ 30° and overlap with each other [71]. As shown in Figures 4 and S2, the base of G40, Aα, base of A41, Gβ, base of G8 in the riboswitch + 3′,3′-cGAMP model and the base of G40, Gα, base of A41, Gβ, base of G8 in the riboswitch + c-di-GMP model are plane-parallel approximation and overlap between the bases, which meet the angle criteria for base stacking. The percentage of distance ≤4.0 Å between base of G8 and Gβ of 3′,3′-cGAMP, base of A41 and Aα of 3′,3′-cGAMP increased 45.54%, 30.58% compared to that between base of G8 and Gβ of c-di-GMP, A41 and Gα of c-di-GMP, which also suggest the stronger interaction of the 3′,3′-cGAMP riboswitch with the 3′,3′-cGAMP molecule compare with the c-di-GMP molecule. Additionally, all possible hydrogen bond interaction sites of the 3′,3′-cGAMP riboswitch in the riboswitch + 3′,3′-cGAMP system occurred at the nucleotides G8, A9, A11, A12, A14, C15 around the J1/2 junction, G40, A41, A42 around the J2/3 junction and C75 of the J1/3 junction (Table 2), which are in good agreement with energy decomposition discussed above [39]. Such hydrogen bond and base-stacking interactions bridge the P1, P2 and P3 through the J1/2, J2/3 and J1/3 junctions, which may be necessary for the orientation of the P1, P2 and P3 helices.

Structural Characteristics of the 3 ,3 -cGAMP Riboswitch in the 3 ,3 -cGAMP-Bound and c-di-GMP-Bound State
To address the structural change of the 3 ,3 -cGAMP riboswitch in the 3 ,3 -cGAMP-bound and c-di-GMP-bound state, the average structures from the equilibrated simulations of the riboswitch + To investigate the differences of interaction of the 3′,3′-cGAMP riboswitch in the 3′,3′-cGAMPbound and c-di-GMP-bound state, the distances between the J1/2 junction and 3′,3′-cGAMP or c-di-GMP, between the J2/3 junction and 3′,3′-cGAMP or c-di-GMP for the riboswitch + 3′,3′-cGAMP and the riboswitch + c-di-GMP systems have been measured ( Figure 6). The average distance between the J1/2 junction and the 3′,3′-cGAMP molecule for the riboswitch + 3′,3′-cGAMP model is 5.17 Å, while the average distance between the J1/2 junction and the c-di-GMP molecule for the riboswitch + c-di-GMP model is 6.08 Å ( Figure 6). This distance with a decrease in the riboswitch + 3′,3′-cGAMP system indicates that the J1/2 junction move toward the 3′,3′-cGAMP molecule, which result in the increase of the hydrogen bond occupancies of the J1/2 junction discussed above. As shown in Figure 6, the average distance between the J2/3 junction and the 3′,3′-cGAMP molecule for the riboswitch + 3′,3′-cGAMP model is 5.93 Å, while the average distance between the J2/3 junction and the c-di-GMP molecule for the riboswitch + c-di-GMP model is 4.84 Å. This distance with an increase in the riboswitch + 3′,3′-cGAMP system indicates that the J2/3 move away to the 3′,3′-cGAMP molecule, which result in the decrease of the hydrogen bond occupancies of the J2/3 junction and have no influence on the ranking of binding energy discussed above.
In addition, the distances between the J1/2 junction and c-di-AMP, between the J2/3 junction and c-di-AMP, between the J1/3 junction and c-di-AMP for the riboswitch + c-di-AMP system, and between the J1/3 junction and 3′,3′-cGAMP in the riboswitch + 3′,3′-cGAMP system have also been measured ( Figure S4). It can be seen that the distances between J1/2, J2/3 and J1/3 junctions and c-di-AMP in the riboswitch + c-di-AMP system are larger than that in the riboswitch + 3′,3′-cGAMP system, which indicate that the J1/2, J2/3 and J1/3 junctions move away to the ligand in the riboswitch + c-di-AMP system. The movements of the J1/2, J2/3 and J1/3 junctions may result in the unbinding of the c-di-AMP molecule to the 3′,3′-cGAMP riboswitch, which is consistent with the binding energy discussed above. To investigate the differences of interaction of the 3 ,3 -cGAMP riboswitch in the 3 ,3 -cGAMP-bound and c-di-GMP-bound state, the distances between the J1/2 junction and 3 ,3 -cGAMP or c-di-GMP, between the J2/3 junction and 3 ,3 -cGAMP or c-di-GMP for the riboswitch + 3 ,3 -cGAMP and the riboswitch + c-di-GMP systems have been measured ( Figure 6). The average distance between the J1/2 junction and the 3 ,3 -cGAMP molecule for the riboswitch + 3 ,3 -cGAMP model is 5.17 Å, while the average distance between the J1/2 junction and the c-di-GMP molecule for the riboswitch + c-di-GMP model is 6.08 Å ( Figure 6). This distance with a decrease in the riboswitch + 3 ,3 -cGAMP system indicates that the J1/2 junction move toward the 3 ,3 -cGAMP molecule, which result in the increase of the hydrogen bond occupancies of the J1/2 junction discussed above. As shown in Figure 6, the average distance between the J2/3 junction and the 3 ,3 -cGAMP molecule for the riboswitch + 3 ,3 -cGAMP model is 5.93 Å, while the average distance between the J2/3 junction and the c-di-GMP molecule for the riboswitch + c-di-GMP model is 4.84 Å. This distance with an increase in the riboswitch + 3 ,3 -cGAMP system indicates that the J2/3 move away to the 3 ,3 -cGAMP molecule, which result in the decrease of the hydrogen bond occupancies of the J2/3 junction and have no influence on the ranking of binding energy discussed above.
In addition, the distances between the J1/2 junction and c-di-AMP, between the J2/3 junction and c-di-AMP, between the J1/3 junction and c-di-AMP for the riboswitch + c-di-AMP system, and between the J1/3 junction and 3 ,3 -cGAMP in the riboswitch + 3 ,3 -cGAMP system have also been measured ( Figure S4). It can be seen that the distances between J1/2, J2/3 and J1/3 junctions and c-di-AMP in the riboswitch + c-di-AMP system are larger than that in the riboswitch + 3 ,3 -cGAMP system, which indicate that the J1/2, J2/3 and J1/3 junctions move away to the ligand in the riboswitch + c-di-AMP system. The movements of the J1/2, J2/3 and J1/3 junctions may result in the unbinding of the c-di-AMP molecule to the 3 ,3 -cGAMP riboswitch, which is consistent with the binding energy discussed above. Figure 6. Time-dependences of distances between the J1/2 junction and ligand, between the J2/3 junction and ligand of the 3′,3′-cGAMP riboswitch for the riboswitch + 3′,3′-cGAMP (magentas) and riboswitch + c-di-GMP (blue) models. The bottom panel highlights these distance differences between the J1/2 junction and ligand, between the J2/3 junction and ligand for the riboswitch + 3′,3′-cGAMP (magentas dotted line) and riboswitch + c-di-GMP (blue dotted line) models in the two models. The yellow circles indicate the centroid of J1/2, J2/3 and ligand for the two models.

Structural Characteristics of the 3′,3′-cGAMP Riboswitch in the Ligand-Free and 3′,3′-cGAMP-Bound State
To address the structural change of the 3′,3′-cGAMP riboswitch in the ligand-free and 3′,3′-cGAMPbound state, the angle changes among the centroid of P2, P3 and P1 of 3′,3′-cGAMP riboswitch of the apo riboswitch and riboswitch + 3′,3′-cGAMP systems have been analyzed and are shown in Figure 7A, and the average structures are shown in Figure 7B. It can be seen that the average angle among the centriod of P2, P3 and P1 decrease from 92.0° in the riboswitch + 3′,3′-cGAMP system to 80.8° in the apo riboswitch system, which indicates that the 3′,3′-cGAMP riboswitch leads to the rotation of P1 helix with 11.2° induced by the 3′,3′-cGAMP molecule unbinding ( Figure 7A,B). Such rotation results in the tighter packing of riboswitch in its free state that is validated by the radius of gyration (Rg) analysis of the 3′,3′-cGAMP riboswitch in the apo riboswitch and riboswitch + 3′,3′-cGAMP systems (Figure 8). The Rg of the apo riboswitch system was lower than that of the riboswitch + 3′,3′-cGAMP system (Figure 8). The P1 helix rotation in the apo riboswitch system may be because of the destruction of hydrogen bonds and base stacking interactions compared with the riboswitch + 3′,3′-cGAMP system. The destruction of these interactions also disturbs the 3′,3′-cGAMP binding pocket conformation. Figure 6. Time-dependences of distances between the J1/2 junction and ligand, between the J2/3 junction and ligand of the 3 ,3 -cGAMP riboswitch for the riboswitch + 3 ,3 -cGAMP (magentas) and riboswitch + c-di-GMP (blue) models. The bottom panel highlights these distance differences between the J1/2 junction and ligand, between the J2/3 junction and ligand for the riboswitch + 3 ,3 -cGAMP (magentas dotted line) and riboswitch + c-di-GMP (blue dotted line) models in the two models. The yellow circles indicate the centroid of J1/2, J2/3 and ligand for the two models.

Structural Characteristics of the 3 ,3 -cGAMP Riboswitch in the Ligand-Free and 3 ,3 -cGAMP-Bound State
To address the structural change of the 3 ,3 -cGAMP riboswitch in the ligand-free and 3 ,3 -cGAMP-bound state, the angle changes among the centroid of P2, P3 and P1 of 3 ,3 -cGAMP riboswitch of the apo riboswitch and riboswitch + 3 ,3 -cGAMP systems have been analyzed and are shown in Figure 7A, and the average structures are shown in Figure 7B. It can be seen that the average angle among the centriod of P2, P3 and P1 decrease from 92.0 • in the riboswitch + 3 ,3 -cGAMP system to 80.8 • in the apo riboswitch system, which indicates that the 3 ,3 -cGAMP riboswitch leads to the rotation of P1 helix with 11.2 • induced by the 3 ,3 -cGAMP molecule unbinding ( Figure 7A,B). Such rotation results in the tighter packing of riboswitch in its free state that is validated by the radius of gyration (Rg) analysis of the 3 ,3 -cGAMP riboswitch in the apo riboswitch and riboswitch + 3 ,3 -cGAMP systems (Figure 8). The Rg of the apo riboswitch system was lower than that of the riboswitch + 3 ,3 -cGAMP system (Figure 8). The P1 helix rotation in the apo riboswitch system may be because of the destruction of hydrogen bonds and base stacking interactions compared with the riboswitch + 3 ,3 -cGAMP system. The destruction of these interactions also disturbs the 3 ,3 -cGAMP binding pocket conformation.   To investigate the binding pocket conformational change of the 3′,3′-cGAMP riboswitch in the apo riboswitch and riboswitch + 3′,3′-cGAMP systems, the distances between J1/2 and J2/3, and between J1/2 and J1/3 in the riboswitch over the simulation times for the apo riboswitch and riboswitch + 3′,3′-cGAMP models has been also analyzed ( Figure 9A,B). The average distance between the J1/2 junction and the J2/3 junction decreases from 9.95 Å in the riboswitch + 3′,3′-cGAMP system to 6.86 Å in the apo riboswitch system ( Figure 9A). At the same time, the average distance between the J1/2 junction and the J1/3 junction is nearly same in two models ( Figure 9B). This distance with a decrease indicates that the J1/2 junction and the J2/3 junction of the 3′,3′-cGAMP riboswitch move close to each other, which may favor the smaller pocket in the apo riboswitch system and rotation of P1 helix. To investigate the binding pocket conformational change of the 3 ,3 -cGAMP riboswitch in the apo riboswitch and riboswitch + 3 ,3 -cGAMP systems, the distances between J1/2 and J2/3, and between J1/2 and J1/3 in the riboswitch over the simulation times for the apo riboswitch and riboswitch + 3 ,3 -cGAMP models has been also analyzed ( Figure 9A,B). The average distance between the J1/2 junction and the J2/3 junction decreases from 9.95 Å in the riboswitch + 3 ,3 -cGAMP system to 6.86 Å in the apo riboswitch system ( Figure 9A). At the same time, the average distance between the J1/2 junction and the J1/3 junction is nearly same in two models ( Figure 9B). This distance with a decrease indicates that the J1/2 junction and the J2/3 junction of the 3 ,3 -cGAMP riboswitch move close to each other, which may favor the smaller pocket in the apo riboswitch system and rotation of P1 helix. To investigate the binding pocket conformational change of the 3′,3′-cGAMP riboswitch in the apo riboswitch and riboswitch + 3′,3′-cGAMP systems, the distances between J1/2 and J2/3, and between J1/2 and J1/3 in the riboswitch over the simulation times for the apo riboswitch and riboswitch + 3′,3′-cGAMP models has been also analyzed ( Figure 9A,B). The average distance between the J1/2 junction and the J2/3 junction decreases from 9.95 Å in the riboswitch + 3′,3′-cGAMP system to 6.86 Å in the apo riboswitch system ( Figure 9A). At the same time, the average distance between the J1/2 junction and the J1/3 junction is nearly same in two models ( Figure 9B). This distance with a decrease indicates that the J1/2 junction and the J2/3 junction of the 3′,3′-cGAMP riboswitch move close to each other, which may favor the smaller pocket in the apo riboswitch system and rotation of P1 helix. (A)

Principal Component Analysis of the 3′,3′-cGAMP Riboswitch in the Ligand-Free and 3′,3′-cGAMP-Bound State
Principal component analysis was used to analyze the trajectories from the corresponding simulations to examine the dominant 3′,3′-cGAMP riboswitch dynamic motions in the apo riboswitch and riboswitch + 3′,3′-cGAMP models ( Table 3). The first 3 principal components, PC1-PC3, account for about 90% of all the motion modes. In the riboswitch + 3′,3′-cGAMP system, the first three principal modes clearly observed a bending motion between P3 and P1 helices, which may cause the closing and opening of P1 helix along the vertical plane and have no obvious rotation of P1 helix on the horizontal plane (Movie S1, S2 and S3). However, in the apo riboswitch system, the first three principal modes obviously showed a twist motion between P3 and P1 helix (Movie S4, S5 and S6), which may induce the rotation of P1 helix on the horizontal plane. Average structure analyses of the 3′,3′-cGAMP riboswitch support the PCA dominant motions (Figure 7). Table 3. Percentages of occupancy times of the first three principal components during the simulations of the apo riboswitch and riboswitch + 3′,3′-cGAMP models.

Principal Component Analysis of the 3 ,3 -cGAMP Riboswitch in the Ligand-Free and 3 ,3 -cGAMP-Bound State
Principal component analysis was used to analyze the trajectories from the corresponding simulations to examine the dominant 3 ,3 -cGAMP riboswitch dynamic motions in the apo riboswitch and riboswitch + 3 ,3 -cGAMP models ( Table 3). The first 3 principal components, PC1-PC3, account for about 90% of all the motion modes. In the riboswitch + 3 ,3 -cGAMP system, the first three principal modes clearly observed a bending motion between P3 and P1 helices, which may cause the closing and opening of P1 helix along the vertical plane and have no obvious rotation of P1 helix on the horizontal plane (Movie S1, S2 and S3). However, in the apo riboswitch system, the first three principal modes obviously showed a twist motion between P3 and P1 helix (Movie S4, S5 and S6), which may induce the rotation of P1 helix on the horizontal plane. Average structure analyses of the 3 ,3 -cGAMP riboswitch support the PCA dominant motions (Figure 7). Table 3. Percentages of occupancy times of the first three principal components during the simulations of the apo riboswitch and riboswitch + 3 ,3 -cGAMP models.

Allosteric Communication of the 3 ,3 -cGAMP Riboswitch from Binding Sites Pocket to the P1 Helix
To explore the allosteric communications, dynamic cross-correlation map (DCCM) of the allosteric process of 3 ,3 -cGAMP riboswitch induced by the 3 ,3 -cGAMP molecule have been analyzed ( Figure 10). The correlation values vary from −1 (high anticorrelation, blue) to 1 (high correlations, red). Because the 3 ,3 -cGAMP molecule binding sites mainly occurred at the nucleotides G8, A9, A11, A12, A14, C15 around the J1/2 junction, G40, A41, A42 around the J2/3 junction, C75 of the J1/3 junction, the DCCM also appeared in those regions of the 3 ,3 -cGAMP riboswitch. It can be found from Figure 10   Simultaneously, the DCCMs of the 3′,3′-cGAMP riboswitch for equilibrated simulation of the riboswitch + 3′,3′-cGAMP, riboswitch + c-di-GMP and apo riboswitch systems are shown in Figure  S5A-C. It can be seen from Figures 10 and S5A,B that the DCCMs of the 3′,3′-cGAMP riboswitch extracted from the equilibrated simulation of the riboswitch + 3′,3′-cGAMP, riboswitch + c-di-GMP systems show very similar features with that from the allosteric process induced by the 3′,3′-cGAMP molecule, which indicate that the similar allosteric communication networks discussed above also exist in the equilibrated riboswitch + 3′,3′-cGAMP and riboswitch + c-di-GMP systems. Comparison of the DCCM of the equilibrated simulations of the apo riboswitch ( Figure S5C) system with that of the riboswitch + 3′,3′-cGAMP ( Figure S5A) system show distinct features. Namely, in the case of riboswitch + 3′,3′-cGAMP simulation ( Figure S5A), a large correlated motion is observed, in contrast, in the apo riboswitch system, the correlated motions have diminished, which indicate that the allosteric communication network from the 3′,3′-cGAMP binding sites to P1 helix discussed above was reduced or abolished in the apo riboswitch system ( Figure S5C). Overall, the allosteric communication between the 3′,3′-cGAMP binding pocket and the P1 helix can obviously occur induced by the ligand binding, such as the 3′,3′-cGAMP molecule or the c-di-GMP molecule.

Allosteric Communication Network from 3′,3′-cGAMP Binding Sites to P1 Helix
As discussed, the indirect and direct allosteric communication occurred at the 3′,3′-cGAMP riboswitch induced by the 3′,3′-cGAMP molecule binding, which relates to the structural changes/ interactions around the J1/2, J2/3, J1/3 junctions and P1 helix. Namely, the 3′,3′-cGAMP molecule binding to the 3′,3′-cGAMP riboswitch causes the movement of the J1/2, J2/3 and J1/3 junctions and the larger binding pocket, which favors the combination of the 3′,3′-cGAMP riboswitch with the Simultaneously, the DCCMs of the 3 ,3 -cGAMP riboswitch for equilibrated simulation of the riboswitch + 3 ,3 -cGAMP, riboswitch + c-di-GMP and apo riboswitch systems are shown in Figure  S5A-C. It can be seen from Figure 10 and Figure S5A,B that the DCCMs of the 3 ,3 -cGAMP riboswitch extracted from the equilibrated simulation of the riboswitch + 3 ,3 -cGAMP, riboswitch + c-di-GMP systems show very similar features with that from the allosteric process induced by the 3 ,3 -cGAMP molecule, which indicate that the similar allosteric communication networks discussed above also exist in the equilibrated riboswitch + 3 ,3 -cGAMP and riboswitch + c-di-GMP systems. Comparison of the DCCM of the equilibrated simulations of the apo riboswitch ( Figure S5C) system with that of the riboswitch + 3 ,3 -cGAMP ( Figure S5A) system show distinct features. Namely, in the case of riboswitch + 3 ,3 -cGAMP simulation ( Figure S5A), a large correlated motion is observed, in contrast, in the apo riboswitch system, the correlated motions have diminished, which indicate that the allosteric communication network from the 3 ,3 -cGAMP binding sites to P1 helix discussed above was reduced or abolished in the apo riboswitch system ( Figure S5C). Overall, the allosteric communication between the 3 ,3 -cGAMP binding pocket and the P1 helix can obviously occur induced by the ligand binding, such as the 3 ,3 -cGAMP molecule or the c-di-GMP molecule.
When the binding of the 3 ,3 -cGAMP riboswitch to the noncognate ligand, the c-di-GMP molecule, the same indirect allosteric and direct allosteric communication networks have been detected versus to the cognate ligand, the 3 ,3 -cGAMP molecule. However, the structural change and interaction differences of binding pocket have been found between the riboswitch + 3 ,3 -cGAMP and riboswitch + c-di-GMP systems. When the c-di-GMP molecule binding to the 3 ,3 -cGAMP riboswitch causes the movement of the J1/2 junction away from the c-di-GMP molecule and the decrease of hydrogen bonds at nucleotides G8, A9, A11, A12, A14 and C15 of J1/2 versus to the 3 ,3 -cGAMP molecule, which leads to the binding energy decrease of the J1/2 junction. Though the binding energy of the J2/3 junction in the 3 ,3 -cGAMP riboswitch has a slight increase, due to the structural change of the J2/3 junction when the binding of the 3 ,3 -cGAMP riboswitch to the c-di-GMP molecule versus to the 3 ,3 -cGAMP molecule, it does not affect the binding ranking. Finally, the results revealed that the combination of 3 ,3 -cGAMP riboswitch with the 3 ,3 -cGAMP molecule is energetically favorable for the combination of the 3 ,3 -cGAMP riboswitch with the c-di-GMP molecule. The 3 ,3 -cGAMP riboswitch binding to the c-di-GMP molecule causes the difference of structure and interactions of the binding pocket (i.e., J1/2, J2/3, J1/3) versus to the 3 ,3 -cGAMP molecule. Furthermore, the correlations of the J1/2, J2/3 and J1/3 junctions to the P1 helix induced by the c-di-GMP molecule binding causes the rotation of P1 helix with 5.8 • versus that by the 3 ,3 -cGAMP binding.

Initial Structures
From previous studies, the initial structure of the 3 ,3 -cGAMP riboswitch in the 3 ,3 -cGAMP-bound state (assigned as riboswitch + 3 ,3 -cGAMP model) was taken from the crystal structure published by Ren et al. (PDB code: 4YAZ) [39] (Figure 1). To obtain the starting structure of the apo 3 ,3 -cGAMP riboswitch, the riboswitch + 3 ,3 -cGAMP model was edited to strip off the 3 ,3 -cGAMP molecule and was assigned as the apo riboswitch model. In order to investigate the specificity of ligand recognition, the structure of the 3 ,3 -cGAMP riboswitch in the c-di-GMP-bound state (assigned as riboswitch + c-di-GMP) was constructed from the crystal structure (PDB code: 4yb0) [39]. The initial structure of the 3 ,3 -cGAMP riboswitch in the c-di-AMP-bound state (assigned as riboswitch + c-di-AMP) was constructed based on the structure of the riboswitch + c-di-GMP model. The coordinates of the c-di-AMP molecule extracted from the crystal structure (PDB code: 3MUV) [39] were superimposed onto the riboswitch + c-di-GMP model as the substitute for the c-di-GMP molecule. To elucidate the 3 ,3 -cGAMP molecule binding to the class I GEMM riboswitch, the initial structure of the class I GEMM riboswitch with G20A mutation in the 3 ,3 -cGAMP-bound state (assigned as G20A-GEMM-I + 3 ,3 -cGAMP model) was taken from the crystal structure (PDB code: 4YB1) [39]. In the five models, we retained all ions and water coordinates. Na + were added to the negative position of RNA in each model to achieve electroneutrality. The systems were solvated in a rectangular box of TIP3P waters, with 10 Å padding of the solvent shell in all directions. Three MD simulations were repeated for each model to check robustness of the results although the starting structure was constructed from the X-ray crystal result, which yield a non-distinguishable result, e.g., the similar binding energy and average structure results (Table S2 and Figure S6A-C), when simulations reach equilibrium.

Ligand Force Field Parameter
The coordinates of the 3 ,3 -cGAMP and c-di-GMP ligand were extracted from riboswitch + 3 ,3 -cGAMP and riboswitch + c-di-GMP models, respectively, and the ligand structure was optimized and RESP [72] charges were calculated at the HF/6-31G* level of theory using the Gaussian09 program [73]. Generalized Amber force field parameters [74] of the 3 ,3 -cGAMP ligand were generated using the Antechamber program in the Amber.

Molecular Dynamics Simulation Protocols
All energy minimization and MD simulations were performed using the sander module in AMBER16 package [75] with the parm99 [76,77] and parmbsc0 force field [78] for RNA nucleotide. MD simulation in the NPT ensemble at the constant pressure and temperature of 1 bar and 300 K was carried out for each system. For each simulation, a 2 fs integration step was used. The MD protocols details are available in the Supplementary Materials and our previous studies [65,[79][80][81].

Principal Component Analysis
Principal component analysis (PCA) can be used to probe the most prominent characteristic of the dynamics of the studied system. The PCA method details are available in the Supplementary Materials and our previous studies [65].
where ∆G binding is the binding free energy, G complex , G ligand and G receptor are the free energy of complex, ligand and receptor, respectively. Computational details can be found in the Supplementary Materials and our previous studies [65,79,80].

Trajectory Analyses
The root-mean-square deviation (RMSD), the hydrogen interaction, the base-stacking interaction and the correlation coefficients [75,86] were calculated by Cpptraj and Ptraj in the AMBER16 program [75] in this work. Details of the correlation coefficients analyses are available in our previous studies [65,[79][80][81].

Conclusions
Molecular dynamics simulations have been performed for the aptamer domain of the 3 ,3 -cGAMP riboswitch in the ligand-free, the cognate ligand-bound (i.e., 3 ,3 -cGAMP-bound) and the noncognate ligand-bound (i.e., c-di-GMP-bound and c-di-AMP-bound) states. The results demonstrated that in the absence of ligand, the distance between J1/2 and J2/3 with a decrease of 3.09 Å causes a smaller binding pocket containing the J1/2, J2/3 and J1/3 junctions. When the binding of the 3 ,3 -cGAMP riboswitch to the 3 ,3 -cGAMP molecule, the hydrogen bonds and stacking interactions between the nucleotides G8, A9, A11, A12, A14, C15 around the J1/2 junction, G40-A42 around the J2/3 junction, C75 of the J1/3 junction and the 3 ,3 -cGAMP molecule increase the size of the 3 ,3 -cGAMP binding site pocket. Furthermore, the structural change of the binding pocket, due to the 3 ,3 -cGAMP binding causes the P1 helix rotation with 11.2 • versus to the apo 3 ,3 -cGAMP structure. The correlation and interaction analysis presents the direct and indirect allosteric communication networks in the