Insights into the mutations leading to capreomycin resistance in S-adenosyl-L-methionine binding motif in TlyA from Mycobacterium tuberculosis

Abstract Capreomycin is a second line antibiotic used for the treatment of drug resistant Tuberculosis (TB), primary reason of death from a solo infectious organism, Mycobacterium tuberculosis (M.tb). Capreomycin targets the ribosome of bacteria and is known to bind at the interface where the large and small ribosomal subunits interact in M.tb using an S-Adenosyl Methionine (SAM) dependent methyltransferase, TlyA (Rv1794). Besides the methyltransferase activity, TlyA has also been found to show substantial haemolytic activity. The dual activity of TlyA highlights its crucial role in pathogenesis and virulence of M.tb. In the present study, docking and molecular dynamics (MD) simulations were carried out to explore the impact of mutations in a conserved SAM binding motif, 90GASTG94, on the affinity of TlyA enzyme for SAM. Two already reported mutations, A91E and S92L, and the remaining wild type residues, Gly90, Thr93, Gly94 mutated to alanine were taken into consideration resulting in a total of six systems, wild type + SAM, G90A + SAM, A91E + SAM, S92L + SAM, T93A + SAM and G94A + SAM that were subjected to 100 ns MD simulations. Docking scores and MD simulations analyses revealed that in contrast to wild type, mutants reduced the affinity of SAM for TlyA with most prominent effect observed in case of alanine mutants. Mutations also led to the loss of hydrogen bond and hydrophobic interactions and large-scale movement of atoms evident from the principal component analyses indicating their destabilizing impact on TlyA. The present study gives insights into influence of mutations on binding of SAM to TlyA in M.tb and promoting capreomycin resistance. Communicated by Ramaswamy H. Sarma


Introduction
Mycobacterium tuberculosis (M.tb) is the etiological agent of Tuberculosis (TB), a deadly infectious disease accountable for millions of deaths globally. As stated by the World Health Organization, in 2019, approximately 10 million people fell ill with TB worldwide amongst which 3.2 million were women, 1.2 million children, and 5.6 million were men (WHO, 2020). Altogether 30 countries accounted for 87% of total new cases with two thirds of cases from India (highest number of cases) followed by Indonesia, China, Philippines, Pakistan, Nigeria, Bangladesh, and South Africa. The emergence of drug resistant TB has further deteriorated the situation. A total of 0.20 million people with multidrug resistant TB (MDR-TB) were found and reported in 2019 adding up to a 10% increase from 0.18 million in 2018 (WHO, 2020).
In addition to being resistant to first-line anti-TB drugs i.e. MDR-TB, M.tb strains resistant to any of the second-line injectable drugs, (amikacin, capreomycin or kanamycin) and fluoroquinolone are known as extensively drug-resistant (XDR-TB). Among the second-line anti-TB drugs, capreomycin is an essential second-line drug which belongs to the aminoglycoside family of antibiotics and is used for treatment of MDR-TB and persistent infection (Lin et al., 2014). This drug binds to 70S ribosomal subunit of bacteria using 2 0 -O-methyltransferase encoded by tlyA gene thus preventing protein synthesis. TlyA catalyses the S-Adenosyl Methionine (SAM) dependent methylation of two cytidines at position 1409 in 16S rRNA present in the small ribosomal subunit and at position 1920 in 23S rRNA located in the large ribosomal subunit (Witek et al., 2017). Although, resistance in drugs targeting ribosomes is largely related to addition of methyl group, in case of TlyA, bacterial antibiotic resistance occurs due to the loss of function which puts TlyA in an exceptional group of methyltransferases (Georghiou et al., 2012;Pablos-Mendez et al., 2002). In recent studies, the mycobacterium TlyA also exhibited hemolysis by destabilizing bacterial membrane thus making it a dually active protein, a ribosomal methyltransferase as well as a hemolysin (Rahman et al., 2010). Fascinatingly, TlyA is restricted to pathogenic strains of M.tb, making TlyA a vital protein to be investigated in detail and understanding its function in M.tb virulence (Kumar et al., 2015). Mutations in TlyA affect the SAM dependent methyltransferase activity of TlyA which leads to the capreomycin resistance in M.tb (Feuerriegel et al., 2009;Stanley et al., 2010). Even with the vital role of TlyA in capreomycin resistance, the knowledge of the structure and mechanism of action of TlyA in capreomycin resistance remains unexplored.
Earlier a novel conserved motif 90 GASTG 94 has been shown to be indispensable for SAM-binding (Witek et al., 2017) and thus mutations in this motif were of prime focus in the present study. Extensive computational analyses involving docking, molecular dynamics (MD) simulations and principal component analyses (PCA) was performed to investigate the conformational changes in TlyA as well as its SAM binding affinity as a result of mutations in 90 GASTG 94 motif. Alanine scanning mutagenesis was used during which residues, Gly90, Thr93 and Gly94 were mutated to alanine along with the two reported mutations, A91E and S92L (Engstrom et al., 2011;Maus et al., 2005). Our analysis reveals changes in conformation of TlyA due to the point mutations in SAMbinding motif thus providing detailed understanding of the enzymatic activity of TlyA and its role in capreomycin resistance.

Protein preparation
The crystal structure of M.tb TlyA (PDB ID: 5kyg) was retrieved from Protein Data Bank (Berman et al., 2000). The protein structure was pre-processed with Accelrys Viewer Lite (Accelrys Inc., San Diego, CA, USA) to remove the water molecules, heteroatoms, and ligands. Further Protein Preparation Wizard accessible through Schrodinger suite was used to add hydrogen bonds, create disulfide bonds, assign bond orders, cap termini, and remove lousy contacts (Sastry et al., 2013).

Alanine mutant generation and receptor grid formation
In order to identify key SAM binding amino acids, known mutations, A91E and S92L (Engstrom et al., 2011;Maus et al., 2005), in SAM binding motif, 90 GASTG 94 of TlyA were investigated. Further alanine mutants were generated to study the impact of mutation on other residues of binding pocket where Gly90, Thr93 and Gly94 were substituted by alanine. The mutants were generated by introducing point mutations in wild type PDB structure through Schrodinger Maestro interface. The Receptor Grid Generation module of Schrodinger Glide (Friesner et al., 2004) was used to generate a grid of size 30 Â 30 Â 30 Å centered on SAM binding residues of wild type and mutant TlyA.

Ligand preparation
The 3D chemical structures of S-adenosyl-L-methionine (SAM) and the drug, capreomycin were obtained from PubChem Compound database, CIDs: 34755 and 3000502, respectively. The compounds were then corrected chemically using Schrodinger LigPrep which generated accurate conformations of the compounds by modifying Lewis structures and removing other mistakes from the structure. Various conformers of both the compounds were generated at different ionization states with pH range of 7-9 and OPLS2005 force field (Jorgensen & Tirado-Rives, 1988) was used for energy minimization.

Molecular docking
Docking studies of wild type and mutant TlyA with SAM and capreomycin were conducted using extra precision (XP) mode of Schrodinger's Glide Ligand Docking (Friesner et al., 2004) module. Glide XP mode is based on anchor-and-grow sampling approach and uses a stringent scoring function that also sets penalty on ligands which do not fit well into the receptor conformation. Furthermore, validation of docking with SAM was performed using AutoDock and PatchDock. Interaction analyses between SAM and wild type and mutant protein was performed using LigPlus (Laskowski & Swindells, 2011) and the docked complexes were visualized using UCSF Chimera (Pettersen et al., 2004) and PyMOL (Llc, 2010).

Molecular dynamics simulations of unbound and SAMbound TlyA
To understand the time-dependent movement of the atoms and molecules surrounded by salts and solvent, MD simulations of unbound TlyA and bound mutant, G90A þ SAM, A91E þ SAM, S92L þ SAM, T93A þ SAM and G94A þ SAM were completed using GROMACS version 5.0.7 (Abraham et al., 2015) and applying GROMOS3a1 force field. The wild type and mutated unbound and bound proteins were placed inside the cubic box packed with SPC water model. In order to balance the total charge inside the cubic box, Na þ and Clions were added, subsequently energy minimization was carried on until the system reached a threshold force of 1000 kJ/mol/nm. To reach a chemical stasis, the system was equilibrated using NVT and NPT ensembles. Finally, 100 ns production run for all individual systems (unbound and SAM bound mutant proteins) were conducted at 300k and 1 bar pressure.

MD Trajectory analysis
Various GROMACS tools, g_rmsd, g_rmsf, g_gyrate, g_sasa and g_hbond, were used for MD simulations analyses. To analyze the average variation of atoms' displacement in an individual system, the root mean square deviation (RMSD) analysis was performed by using the g_rmsd tool. To determine the variation in mutant structure from the wild type, the root mean square fluctuation (RMSF) was performed using g_rmsf. g_gyrate was used to estimate the radius of gyration (Rg) of the wildtype and mutated protein, which provides the distance between the atoms and center of mass for the individual system thus providing a better understanding of compactness in the structure. Further to find out the area of the protein within reach of the solvent and its impact on protein folding and conformation, g_sasa was used. g_hbond was used to compute hydrogen bonds between SAM and wild type and mutant TlyA.

Principal component analysis and free energy landscape computation
The collective movement of atoms in wild type and mutant protein conformations calculated using Principal Component Analysis (PCA). The diagonalized covariance matrix from atomic variations was created using g_covar and eigen vectors of covariance matrix denoting the direction of the atomic movement were analyzed using g_anaeig of GROMACS. Further Gibbs free energy landscape (FEL) plots were generated from the first two principal components (PC1 and PC2) using g_sham package.

Results and discussion
The present study was designed to explore the impact of substitutions in SAM binding pocket of TlyA protein from M.tb. Two reported mutations, A91E (Maus et al., 2005) found in Beijing strain and S92L found in M.tb H37Rv strain (Engstrom et al., 2011), were considered whereas the other residues, Gly90, Thr93 and Gly94, were substituted by alanine to look into the impact of these substitutions on binding of SAM with TlyA. In the present study, MD simulations were performed at two levels. Prior to molecular docking, MD simulations were performed on the wild type and mutated TlyA (without SAM) proteins. Further, the stable protein structures were extracted from MD trajectories of wild type and mutated proteins and these structures were used for conducting docking studies using SAM. Next wild type and mutant TlyA-SAM complexes were subjected to MD simulations.

Molecular dynamics simulation analyses of the wildtype and mutants in unbound TlyA
The amino acids, Gly90, Ala91, Ser92, Thr93 and Gly94 have been reported to form the SAM binding motif in M.tb TlyA protein (Witek et al., 2017). Amongst these, Ala91 and Ser92 have already been reported to confer resistance to capreomycin antibiotic (Maus et al., 2005). At codon position 91 st , non-polar, aliphatic alanine was replaced by acidic and polar glutamic acid whereas polar serine at 92 nd position, was replaced by aliphatic non-polar leucine. The other amino acid, non-polar aliphatic glycine at positions 90 and 94, and polar threonine at 93 rd codon position was substituted by alanine.
RMSD analysis for the unbound wild type TlyA as well as all the mutant forms indicated that the system attained an equilibrium and protein was stable during the course of 100 ns simulation period. The average RMSD of the backbone in case of wild type and mutants was found to be in the range of 0.22 to 0.29 nm ( Figure 1A and B).
Further in RMSF analyses, an overall decreased fluctuation was seen in TlyA mutants as compared to the wild type  ( Figure 1C and D). The mutant G94A had the lowest overall average RMSF value equivalent to 0.16 nm. Moreover, RMSF value for the key residues in SAM binding region was lower for wild type in comparison to the mutants except T93A (Table 1).
The RMSF values in the table indicate that wild type was more rigid with respect to the mutants during MD simulations except in case of T93A, where at 93rd substitution position, RMSF for wild type Thr was 0.35 nm and RMSF reduced to 0.26 nm upon substitution by Ala. As compared to wild type, all the residues, except Thr93 had higher RMSF at the substitution position indicating that the mutants increased the flexibility of the SAM binding site thus decreasing the stability of TlyA. Gly90 in wild type had RMSF of 0.14 nm and when substituted to Ala, RMSF increased to 0.20 nm. In case of A91E, wild type RMSF of Ala at 91st position was 0.12 nm whereas upon substitution to Glu, RMSF increased to 0.16 nm. At position 92nd for Ser in wild type, RMSF was 0.20 and upon substitution to Leu, RMSF increased to 0.26 nm. In G94A mutant, RMSF for both, wild type Gly and mutant Ala was 0.21 nm.

Molecular docking and interaction analyses of wild type and mutant SAM-TlyA complex
Wild type M.tb TlyA protein as well as five, Gly90, Ala91, Ser92, Thr93 and Gly94, mutant protein structures were docked with SAM to study the impact of mutations on binding of SAM with TlyA. Figure 2 describes the chemical structures of SAM and the drug, capreomycin.
Molecular docking results showed strong binding affinity of SAM with wild type TlyA as is evident from high Glide docking score of À4.84 kcal/mol and AutoDock score of À6.18 kcal/mol as well as five hydrogen bonds and four hydrophobic interactions. As compared to wild type, mutants showed poor binding affinity, the Glide docking scores being À4.64, À4.05, À4.20 and À4.51 kcal/mol for G90A, S92L, T93A and G94A respectively ( Table 2). The global binding energies obtained from PatchDock were also less for all TlyA variants as compared to the wild type.
An overall reduction in hydrogen bonding and hydrophobically interacting residues was seen in case of all mutants except A91E. Supplementary Table 1 lists the hydrogen bonding and hydrophobic interactions between SAM and wild type and mutant TlyA. The mutant A91E had high Glide docking score, À5.56 kcal/mol contrary to wild type. This can be due to the tendency of glutamate to pair with positively charged amino acids and form hydrogen bonds thus stabilizing the protein (Mj & Rb, 2003).

Molecular dynamics simulations analyses of TlyA-SAM complexes
Detailed MD simulations were performed to gain insights into conformation of SAM and stability of wild type (wild type þ SAM) and mutant (G90A þ SAM, A91E þ SAM, S92L þ SAM, T93A þ SAM and G94A þ SAM) complexes. RMSD for the wild type and mutant TlyA-SAM complexes was found to be 0.25 nm. The wild type and mutant TlyA-  SAM complexes except G94A had average RMSD value between 0.20-0.25 nm and remained stable throughout the 100 ns simulation period indicating that the systems could be considered for further analyses ( Figure 3A and B). The G94A þ SAM complex had lowest RMSD of 0.16 nm. In particular, the average backbone RMSD values for wild type þ SAM, G90A þ SAM, A91E þ SAM, S92L þ SAM and T93A þ SAM complexes were 0.23 nm, 0.25 nm, 0.20 nm, 0.25 nm, and 0.22 nm, respectively.
Further RMSF which is an indicator of residue wise fluctuation was analyzed. Low level of fluctuation points towards rigidity and protein stability whereas high level of fluctuation means high flexibility of protein during simulations. The overall similar RMSF values were observed in case of wild type and other mutants except G90A ( Figure 3C and D). Huge peaks in the region covering residues 184-198 denoting highest RMSF value, 0.20 nm, were seen in case of G90A mutant complex. In case of SAM binding residues, higher  average RMSF values and greater movement and instability than wild type (0.16 nm) was observed for G90A (0.39 nm) and T93A (0.15 nm) complexes. However low level of flexibility was exhibited in case of A91E, S92L and G94A complexes with average RMSF values corresponding to 0.08 nm, 0.09 nm and 0.14 nm.
The average Rg values for the wild type was 1.65 nm indicating towards the compactness of protein structure through the entire simulation period. In case of mutant SAM-TlyA complexes, the average Rg values were 1.64 for G90A, 1.60 nm for A91E, 1.63 nm for S92L, 1.65 nm for T93A and G94A complexes ( Figure 4A and B). Low Rg values for mutants, G90A, A91E and S92L, showed that these mutant TlyA-SAM complexes were more tightly packed than wild type.
The area of the protein exposed to the solvent was measured using SASA. It was observed that protein was buried deeply in case of mutants, A91E and S92L with average SASA values each of 103.45 nm 2 and 102.87 nm 2 as compared to wild type which had high SASA value of 106.65 nm 2 ( Figure 4C and D). However, the protein was slightly highly exposed to solvent in case of G90A mutant TlyA-SAM complex with average SASA value of 107.55 nm 2 . The average SASA value for mutant G94A-SAM complex was 106.63 nm 2 .
Hydrogen bonds formed between SAM and wild type and mutant TlyA proteins during MD simulations were also measured. Figure 5 shows the number of hydrogen bonds formed between SAM and wild type and mutant TlyA. The average number of hydrogen bonds between SAM and TlyA protein was highest, 2.21 in case of wild type and lowest in case of T93A (0.57) mutant model followed by G94A (0.75). Lower number of hydrogen bonds was also observed in mutants, G90A (0.97), A91E (1.43) and S92L (1.17).
In conclusion, RMSF, Rg, and SASA and hydrogen bond analyses revealed that the mutants had significant impact on overall conformation of TlyA. Furthermore, all the mutations reduced the binding affinity of SAM towards TlyA contrary to wild type as is evident from the hydrogen bond analyses.

Interaction pattern analysis after MD simulations
The representative structures of wild type and mutant TlyA-SAM complexes were extracted from the MD simulations time frame during which complexes were most stable. A change in the orientation of SAM in the binding pocket of wild type TlyA was observed during MD simulations and the complex was stabilized by four hydrogen bonds (Gly190, Val192, Gly195 and His199) and seven hydrophobic interactions (Val187, Gln191, Gly193, Pro194, Val197, Val198 and Leu203) ( Figure 6A). Post MD-simulations, a huge reduction   Table 2. In case of G90A-SAM complex, SAM formed two hydrogen bonds with Gly94 and Gln117 and six hydrophobic interactions with Tyr115, Gln117, Leu139, Phe157 and Phe158 ( Figure 6B).

Principal component analysis and free energy analysis
In order to enhance our understanding of changes in conformation of TlyA protein, collective atomic movements were analyzed in case of wild type and mutant TlyA-SAM complexes using MD simulations trajectories projected on the first two principal components (PC1 and PC2). Large area covered by dots is indicative of huge collective motions and substantial conformational changes and more flexibility. The dots were highly scattered in case of all the mutants in contrast to wild type suggesting protein stability was decreased in all mutant TlyA-SAM complexes. Figure 7 shows the projection of atomic movements using the first two principal components (PC1 and PC2) for wild type and mutant TlyA-SAM complexes. The trace values computed from covariance matrices for the wild type, G90A, A91E, S92L, T93A and G94A were 87.73 nm 2 , 128.04 nm 2 , 73.73 nm 2 , 82.82 nm 2 , 78.09 nm 2 and 71.48 nm 2 .
Higher trace values in case of G90A with respect to wild type indicate large magnitude of flexibility upon ligand binding.
To gain further insights into structural transitions of wild type and mutant TlyA-SAM complexes, FEL were plotted using the coordinates from the first two principal components, PC1 and PC2, which showed that DG value ranged from 0 to 18.9, 20.5, 18.7, 19.4, 20.9 and 21.1 for wild type, G90A, A91E, S92L, T93A and G94A, respectively (Figure 8).
The shape and size of minimal energy area (blue color region) indicates the stability of the protein where a small and concentrated blue area denotes more stable and low energy conformations. Wild type TlyA traversed a large conformational space suggesting the presence of large number of low energy conformations ( Figure 8A). A single concentrated energy was seen in case of T93A ( Figure 8E) mutant indicating minimum structural changes. On the contrary, several minima were observed in case of wild type and other mutants, G90A ( Figure  8B), A91E ( Figure 8C), S92L ( Figure 8D) and G94A ( Figure 8F).

Secondary structure analysis
Time dependent changes in secondary structure of SAM bound wild type and mutant TlyA were computed (Table 3).
An increase in a-helix was observed in case of all the mutants as compared to wild type. The overall b-sheet content was reduced in all the mutants except A91E in which it was same as to in WT. Further among flexible secondary structure elements, turns and coils, reduced percentage of turns was seen in G90A, A91E and T93A, a slight increase was seen in G94A and equal to wild type turn content was found in case of S92L mutation. Decreased coil percentage was observed in case of all the mutants with the lowest in G90A. The total variations in secondary structural elements in case of mutants with respect to wild type during 100 ns simulation period indicate that mutations led to significant changes in conformation of TlyA protein.

Molecular docking of wild type and mutant TlyA with capreomycin
Docking studies were performed to examine the impact of wild type and mutant (G90A, A91E, S92L and T93A) TlyA on the binding of Capreomycin. Glide XP results showed that all the mutants in SAM-binding motif drastically reduced the binding affinity of Capreomycin for TlyA (Table 4). The two reported mutations, A91E (docking score of À0.27 kcal/mol and S92L À2.19 kcal/mol) had the most significant impact on drug binding as compared to wild type (docking score of À5.66 kcal/mol).

Conclusion
The increase in virulence as well as capreomycin resistant strains emerged as a result of inactivation of tlyA has made the treatment of TB difficult. Earlier reports have linked capreomycin resistance in M.tb to mutations in TlyA thus prompting towards detailed studies to understand the mechanism of resistance (Maus et al., 2005). In the present study, we have explored the interaction between M.tb TlyA and SAM, a compulsory co-substrate for methyltransferase activity, along with the impact of mutations in TlyA on SAM binding. Witek et al. (Witek et al., 2017) had revealed a conserved SAM binding motif, 90 GASTG 94 , present in carboxy-terminal domain of TlyA. Mutations, A91E and S92L, are already known to responsible for capreomycin resistance (Engstrom et al., 2011;Maus et al., 2005). Alanine mutants were generated for the other amino acids, G90, T93 and G94, to elucidate the role of these residues along with A91 and S92 in SAM binding and capreomycin resistance using MD simulations. Compared to wild type TlyA, the mutants, G90A, S92L, T93A and G94A, showed reduced binding affinities for SAM and Capreomycin as demonstrated from the docking scores and interaction profiles. MD simulations analyses also established that mutants, especially alanine mutants, induced large conformational changes and decreased the affinity of SAM for TlyA thus reducing the overall stability of TlyA-SAM complexes. The study reveals that residue, 90 GASTG 94 , are critical for SAM binding and functioning of TlyA and the results obtained in the present study could be used for designing inhibitors that could defy resistance in MDR and XDR TB cases.