In silico Analysis of Novel Mutation Ala102Pro targeting pncA Gene of M. tuberculosis Computer Science & Systems Biology

This study reports on the structural and functional basis of pyrazinamide (PZA) resistance conferred by a novel mutation Ala102Pro in pncA gene as sequenced from a PZA resistant Mycobacterium tuberculosis strain. Molecular modeling studies of Wild Type (WT) and Mutant Type (MT) of Pyrazinamidase (PZase) showed the mutation at Ala102Pro does not impact on the conformation of the protein. However, the docking studies infer that MT has a higher inhibitory constant (Ki-990.0m) compared to WT (Ki-822.42m), which is indicative of drug resistance in MT. Furthermore, molecular interaction studies also reveal that WT forms 4 hydrogen bonds involved in PZA-WT inhibitory interactions, whereas, in case of MT, it formed 5 hydrogen bonds with PZA. However, Ala102 in WT was found to be less fluctuating and more stable in molecular dynamics simulation when compared to Pro102 in MT which was highly fluctuating and unstable. This implies that Ala102 shall be a key residue involved in PZA inhibitory interactions. Moreover, MT does not show hydrogen bonding with PZA with Pro102 and also deviating in terms of PZA binding pose in comparison with WT. Hence, the observed deviations in terms of MT-PZA interactions shall be attributed to the drug resistance conferred.


Introduction
Pyrazinamide (PZA) is the first-line drug used in the treatment of tuberculosis along with isoniazid and rifampicin and it inhibits semidormant mycobacteria only at low pH in vitro.
The pncA gene encodes pyraziniamidase (PZase), and mutations in pncA are associated with resistance to PZA [1,2] or loss of PZase activity. PZA acts by targeting the fatty acid synthase/synthetase enzyme, and is responsible for the killing of persistent tubercle bacilli in the initial intensive phase of chemotherapy. It is a prodrug that is converted to its active form namely, pyrazinoic acid (POA) by the catalytic action of PZase enzyme, encoded by the pncA gene in M. tuberculosis. Interestingly, PZA is active only at low pH since acidic environment favours accumulation of POA in the cytoplasm due to an ineffective efflux pump, thereby leading to improper efflux out of the amidase from the cell to the exterior [2,3].
Although a large number of mutations have been described, no mutational hotspots have been identified so far [4]. This can be explained by the fact that mutations occur along the entire length of the pncA gene [5]. However, studies show that mutations that confirm the PZA resistance occur mainly in the putative promoter region of the pncA gene. Scorpio and Zhang in 1996 [6] had identified the PZase gene (pncA) from M. tuberculosis and had shown that pncA mutations are a major mechanism of PZA resistance [6]. The identified pncA mutations are largely missense mutations causing amino acid substitutions, and in some cases nucleotide insertions or deletions and nonsense mutations in the pncA structural gene or in the putative promoter region of pncA [7]. The uniqueness of the mutations of pncA gene is its diversity and scattering along the whole gene though there does appear to be some degree of clustering at three regions of pncA protein (3 to 17, 61 to 85, and 132 to 142). These regions are likely to contain catalytic sites for the Pzase enzyme [8]. The catalytic residues comprise the active site (D8, K96, A134 and C138) and the metal-binding site (D49, H51 and H71) [9,10]. Cys-138, Ala-134, Thr-135, Trp-68, and Asp-8 in the M. tuberculosis PZase could be key residues for hydrolysis of PZA [8]. At the protein level, these regions were found to be well conserved among the amino acid sequences of pncA proteins from different bacterial species.
The crystal structure of the M. tuberculosis pncA protein has been determined, showing significant differences in the substrate binding cavity when compared to the pyrazinamidases from Pyrococcus horikoshii and Acinetobacter baumanii. In M. tuberculosis, this region was found to hold a Fe 2+ ion coordinated by one aspartate and three histidines, the most crucial structural element in this loop appears to be the specific positioning of residue His57 which is directly involved in the coordination of the Fe 2+ ion. The overall architecture of the pyrazinamidase of M. tuberculosis is similar to that reported for the other pyrazinamidases of A. baumanii and P. horikoshii [11].
In the pncA model, the putative catalytic centre would be located in a pocket formed by one α-helix (αE), four β-strands (β1, β2, β3 and β4) consisting of β1 (Asp-8 and Phe-13), β2 (Asp-49), β3 (Lys-96), β4 (Ala-134 and Thr-135) and αE (Cys-138). In this pocket, the conserved active cysteine residue Cys-138 is located close to the conserved residues: Asp-8, Trp-68, Lys-96, Ser-104, Ala-134 and Thr-135. In the pncA model, the side chains of the two residues Asp-8 and Lys-96 are found to point towards Cys-138 of the active-site. The modification of the amino acid residues Asp-8, Lys-96 and Ser-104 in the mutants D8G, K96T and S104R resulted in enzymes showing specific activities drastically impaired (%0.004 unit mg), thus suggesting that these Abstract This study reports on the structural and functional basis of pyrazinamide (PZA) resistance conferred by a novel mutation Ala102Pro in pncA gene as sequenced from a PZA resistant Mycobacterium tuberculosis strain. Molecular modeling studies of Wild Type (WT) and Mutant Type (MT) of Pyrazinamidase (PZase) showed the mutation at Ala102Pro does not impact on the conformation of the protein. However, the docking studies infer that MT has a higher inhibitory constant (Ki-990.0m) compared to WT (Ki-822.42m), which is indicative of drug resistance in MT. Furthermore, molecular interaction studies also reveal that WT forms 4 hydrogen bonds involved in PZA-WT inhibitory interactions, whereas, in case of MT, it formed 5 hydrogen bonds with PZA. However, Ala102 in WT was found to be less fluctuating and more stable in molecular dynamics simulation when compared to Pro102 in MT which was highly fluctuating and unstable. This implies that Ala102 shall be a key residue involved in PZA inhibitory interactions. Moreover, MT does not show hydrogen bonding with PZA with Pro102 and also deviating in terms of PZA binding pose in comparison with WT. Hence, the observed deviations in terms of MT-PZA interactions shall be attributed to the drug resistance conferred.
In silico Analysis of Novel Mutation Ala102Pro targeting pncA Gene of M. tuberculosis residues are essential for the pncA activity. The amino acids found at positions 8, 13, 61, 69, 96, 103, 104 and 146 are functionally and or structurally important in pncA [9]. Hence, in this study we attempt to utilize this empirical structural data to analyze to the impact of Ala102Pro mutation on PZA resistance.

DNA extraction
DNA was extracted from the M. tuberculosis isolate by boiling at 80°C for 10 minutes followed by centrifugation at 3,000 rpm. 5 µl of supernatant was used as the template DNA.

PCR targeting pncA gene
The amplification reactions contained 200 μM of each dNTPs [dATP, dTTP, dGTP, dCTP (Bangalore Genei)], 1 μM of outer and inner set of primers, 1×buffer [10 mM Tris -HCl (pH 8.3), 50 mM KCl, 1.5 mM MgCl 2 ], 3 units of Taq DNA polymerase enzyme and 5 μl of template DNA. PCR targeting pncA gene was performed with the following forward and reverse primers pncA1: 5'GGCGTCATGGACCCTATATC3' and pncA2: 5'CAACAGTTCATCCCGGTTC3' [13] with the thermal profile for 30 cycles of denaturation at 95°C for 1 min, annealing at 60°C for 30 sec and extension at 72°C for 1 min which yielded a 670 bp product. For each PCR run, negative control and positive control were included in the study. The PCR results were considered valid only when the negative control was negative without amplicon and the positive control yielded specific band.

Detection of amplified product
10 µl of the amplified product was subjected to electrophoresis on 2% agarose gel incorporated with 0.5 µg/ml ethidium bromide for visualization by UV transilluminator (Vilber Lourmet, France).

Standard Precautions for PCR
Strict and adequate precautions were taken to prevent amplicon getting contaminated. Separate rooms were used for the preparation of DNA, its amplification and analysis of amplified products. PCR cocktail was prepared in a laminar flow workbench; the microfuge tubes and milliQ water used for PCR were double sterilized.

DNA sequencing and Data analysis
DNA Sequencing of amplicons was carried out using an ABI prism 3110 automated DNA sequencer. Cycle sequencing of the amplified products was carried out using ABI Prism BigDye terminator kit (Applied Biosystems, USA) following manufacturer's instruction. The sequences were analyzed by Bio Edit sequence alignment software [14]. The sequences generated were compared with the wild type sequence (Genbank accession nos. for X68081 for katG) by using Multalin software [15] to identify the presence of mutation or polymorphism.

Data sets
The protein sequence for Pyrazinamidase (Pzase) pncA (UniProtKB id: Q50575) was obtained from UniProtKb to perform in silico sequence analysis. The 3D atomic coordinates (PDB ID:3PL1) were obtained from Protein Data Bank (http://www.pdb.org/pdb/home/home.do) for structure analysis, mutation modeling (Ala102Pro), molecular docking and molecular dynamics simulation studies.

Analysis of secondary structure using Protein Structure Prediction Server (PSIPRED v3.0)
The PSIPRED is a reliable sequence based Protein Secondary Structure Prediction Server and was used in this study to predict the secondary structure of the mutant sequence so as to assess the impact of the novel mutation at the secondary structure elements. The output gives the details of presence of helix (H), sheets (E) and coils(C) in the protein sequence with graphical representation (http://bioinf.cs.ucl. ac.uk/psipred/).

Detection of conserved amino acids in the human pathogenic organisms containing pyrazinamidase/ nicotinamidase
Consurf analysis can be utilized to reveal the conservation of residues among the group of organisms with common metabolic functions. Here, in this study, the bacteria which exhibit PZase activity were grouped together in terms of sequence identity and were further analyzed for conserved residues through consurf analysis.

Homology modeling of mutant Pzase
Homology modeling helps in predicting the 3D structure of a protein with unknown structure by comparing it with a known structure sharing high sequence similarity. In this study, MODELLER9v7 [16] was used for modeling the novel mutant protein. Since the crystal structure of the pyrazinamidase M. tuberculosis (PDB id: 3PL1) was available [12], the same was used as template to predict the mutant form. Finally, the stereochemcial property of WT and MT of the protein structure was predicted using SAVS server (http://nihserver. mbi.ucla.edu/SAVES/).

Structure optimization and validation
The WT and the modeled MT PZAse structures were subjected to energy minimization by steepest descent algorithm using GROMOS96 force field [17]. The quality of the protein structure modeled was again checked using Q-Mean server [18] and ProQ (http://www.sbc. su.se/~bjornw/ProQ/ProQ.cgi)

Active site analysis of Pzase enzyme
The active site of Pzase from M. tuberculosis was analyzed using Q-SiteFinder, CASTp (http://sts-fw.bioengr.uic.edu/castp/calculation. php) and Pocket Finder server (http://www.modelling.leeds.ac.uk/ pocketfinder/help.html) to identify the binding site for docking studies. The results obtained were correlated with the active site residues from published literature [7,8,19].

Ligand optimization
The initial structure of PZA was obtained from PubChem in 3D SDF format. It was converted to PDB format using Open Babel (http://www.molecularnetworks.com/online_demos/corina_demo. html). Hydrogen bonds were added to the structure and was geometry optimized (Broyden-Fletcher-Golfarb-Shanno line search method set to 1000 steps) using ArgusLab [20], wherein, the energy minimization was carried out using Universal Force Field (UFF) [21].

Docking of PZA-Pzase
PZA was docked into the structures of WT and MT using AutoDock 4.0 [22]. In this docking simulation, semiflexible docking protocols [23] were used, in which the protein structures were kept rigid and the PZA being docked was kept flexible. Blind docking was performed using grid point value (X, Y and Z) of 100Ǻ and spacing between the gird points was 0.375Ǻ. Lamarckian genetic algorithm was selected for ligand conformational searching and default docking parameters were used. The best docked complexes based on the lowest binding energy were further analyzed for hydrogen bonding interactions and the binding energy of WT and MT was compared. Finally, the Proteinligand complexes were analyzed using Discovery Studio Visualizer [24] and PYMOL visualization tool (http://www.pymol.org/).

Docking complex simulation
To study the stability of PZA ligand with wild and mutant forms, molecular dynamics simulation was carried out for the docked Proteinligand complexes. All simulations were performed using GROMACS [19], with GROMOS96 43al force field. The topology file for the ligand was generated using Dundee PRODRG Server [http://davapc1.bioch. dundee.ac.uk/prodrg/].
The system was placed at the center of a cubic periodic box with an area of 37×37×37 Å and was solvated with a simple point charge (SPC213) water model. The distance between solute and edge of the box was 0.9 Å for the system. It included water molecules for all the complexes. In order to neutralize the systems charge to zero, adequate number of Na + or Clions were added. Each protein-ligand complex was subjected 2000 steps of energy minimization using steepest descent algorithm. Equilibration of the ensembles was performed by position restraint of the system for 100 ps. A constant temperature of 300 K and 1 atm pressure was maintained in the system. For long range electrostatics calculation, Partial Mesh Ewald coulomb type was used with a cut off of 1.0 nm. All bonds including H-bonds were constrained using LINCS algorithm. The MD simulation was carried out for 1ns for each complex. After simulation, the trajectories were analyzed by ngmx software package. All the molecular modeling and simulation studies were carried on OPEN Discovery 2 Linux Platform [25].

Results and Discussion
Phenotypic drug susceptibility testing

SNP Analysis of I-Mutant Server
In this study, wild type sequence of pncA at 102 nd position Ala was replaced by Pro to predict protein stability changes through I-mutant server. The results infer loss of stability by the mutant protein with negative Gibbs free energy value of -1.10 at pH 5.5 and 37°C.

Secondary structure analysis for Wild and Mutant Type of PZA (Pdbsum SERVER)
PSIPred results for the MT shows no significant change when compared to WT and MT (Figure1A and figure 1B).

Homology modeling and loop refinement of MT Pzase
The structure of MT Pzase (Ala102Pro) was modeled using Modeler9v7 with WT structure as template. The modeled structure was found to be highly plausible as it had 99% sequence identity with that of the template. Moreover, the Ramachandran plot also showed 83.5% of residues in most favored regions with no residues in disallowed region ( Figure 2).

Structure optimization and Validation
The structures of WT and MT were energy minimized using GROMOS96 force field by Gromacs software package. The potential energy of WT and MT were found to be -2.6391059e+04 Kcal/mol and MT -2.3403371e+04 Kcal/mol, respectively. The optimized structure was validated using Q-mean and ProQ Server (Table 1 and figure 3).

Docking studies of WT &MT type of Pzase with PZA
The docking studies were performed using the Autodock software. The aminoacid residues Asp8, Phe13, Thr61, Pro69, Lys96, Tyr103, and Cys138 were assigned as catalytic region as these residues were proven to be potentially active site [8,9,19].
The WT had binding energy of -4.21Kcal/mol and theoretical inhibitory constant (Ki) of 822.42 µM ( Table 2). The residues Asp8, Lys96, and Ala102 interacted through 4 hydrogen bonds with PZA. The interaction was found to span within the chosen active site [8,9,19].
The MT had binding energy of-4.1 Kcal/mol and theoretical inhibitory (Ki) constant of 990.0 µM ( Table 2). The results of docking study infer that both WT and MT interact with PZA at the same binding pocket region. However, MT showed deviation in terms of interaction with PZA by excluding hydrogen bond formation with Lys96 and Ala l02 and including Ala134 and Ile133 (Figures 5A and 5B).

Dynamics Simulation analysis
To further validate the molecular docking studies, high performance molecular dynamics simulation protein-ligand complex were performed and from the resultant trajectory, root mean square fluctuation (RMSF) and gyration analysis were performed. The RMSF results showed that the MT has high fluctuation when compared to WT due to the acquired mutation ( Figure 6A). Further, the Radius of gyration analysis also showed MT to be slightly compact than WT which shall interfere with the flexibility of the protein ( Figure 6B).

Conclusion
The novel mutation Ala102Pro was found to span in the putative active site (Lys96-Tyr103) region of PZase. The SNP analysis by I-mutant server showed that the stability of the protein was affected with negative Gibbs free energy value-1.10 at pH 5.5 and 37ºC. Molecular modeling and docking analysis showed that MT has a higher inhibitory constant than WT in terms of PZA binding, which indicates that the drug affinity is highly affected in the MT due to structural changes. The molecular dynamics simulation and gyration analysis together infer the change in flexibility at the active site cavity of MT. Hence, this study gives insight on the impact of novel mutation on the activity of this protein which can be attributed to the drug resistance observed. This study also presents a need for the discovery of potential newer drug molecules against both mutant and wild types.