Effect of arsenate substitution on phosphate repository of cell: a computational study

The structural analogy with phosphate derives arsenate into various metabolic processes associated with phosphate inside the organisms. But it is difficult to evaluate the effect of arsenate substitution on the stability of individual biological phosphate species, which span from a simpler monoester form like pyrophosphate to a more complex phosphodiester variant like DNA. In this study, we have classified the physiological phosphate esters into three different classes on the basis of their structural differences. This classification has helped us to present a concise theoretical study on the kinetic stability of phosphate analogue species of arsenate against hydrolysis. All the calculations have been carried out using QM/MM methods of our Own N-layer Integrated molecular Orbital molecular Mechanics (ONIOM). For quantum mechanical region, we have used M06-2X density functional with 6-31+G(2d,2p) basis set and for molecular mechanics we have used the AMBER force field. The calculated rate constants for hydrolysis show that none of the phosphate analogue species of arsenate has a reasonable stability against hydrolysis.

1. Introduction

Maintaining negative charge at physiological pH leads phosphate (PO or P_i) ester to dominate the living world as it prevents chemically stored energy and genetic material from escaping the cell [1]. (V) is also capable of retaining negative charge in its phosphate analogue arsenate (AsO or As_i) over a range of physiological pH conditions [2] and it is already demonstrated in bacteria that at a high concentration ratio of As_i/P_i the phosphate transporters allow the transport of As_i on the basis of structural similarity [3]. This structural analogy further leads to the integration of As_i into various metabolic processes associated with P_i inside the organisms [2], whereas the primary concern related to the P_i analogue species of As_i (PASA) is their rapid rate of hydrolysis [1]. However, the scientific world is really surprised by the discovery of arsenate substituted DNA backbone in place of phosphate (As_i-DNA) in a bacterium strain GFAJ-1 of the Halomonadaceae family, which was isolated from Mono Lake, California with a very high As_i/P_i ratio (greater than 103) [4]. The theoretical investigations on the possibility of As_i-DNA also revealed that both arsenodiester and phosphodiester linkages have similar geometric and conformational properties [5,6]. Although these theoretical studies support geometrical similarities between As_i-DNA and P_i-DNA, neither of them claims the stability of As_i-DNA against hydrolysis. Whereas the models described by Fekry et al. and Mládek et al. to estimate the kinetic stability of phosphodiester and arsenodiester linkages in DNA lack the fundamental base stacking interactions, also called π − π stacking [6,7]. In order to evaluate the stability of a molecule like DNA, a more realistic model is required because the work of Wang et al. reveals that the stacking between neighbouring bases in DNA single strands raises the activation energy barrier for hydrolysis of DNA [8]. Along with this, it is found that the discovery of As_i-DNA in GFAJ-1 captured the attention of the scientific community from other PASA. However, it is not clear how As_i could be directly incorporated into P_i-DNA because neither As_i nor P_i are the substrates for the DNA synthesizing enzymes and neither of these are directly incorporated into the diester backbone of the DNA. There are various metabolic steps through which monosaccharides are processed using P_i biomolecules and metabolized into nucleoside triphosphate (NTP’s). NTP’s are the monomeric units of nucleic acid polymers (DNA and RNA), and P_i is able to form the phosphodiester backbone of the DNA through the nucleophilic substitution reaction between deoxy-NTP molecules. Similarly, before reaching DNA, As_i must have a reasonable stability in its dNTP analogues, whereas biosynthesis of As_i analogue of dNTP again depends upon the availability of some other PASA like Ribose-5-As_i, Glucose-1-As_i, Pyro-As_i, etc. [9]. Physiologically, these are dianionic monoesters of P_i and highly sensitive to As_i substitution [10]. Hence, to investigate the effect of As_i substitution on the P_i repository of the cell a comparative study at different biomolecular level is needed, that is not merely confined to analysing the impact of As_i substitution on the stability of a particular P_i biomolecule. However, due to a vast number of P_i species inside organisms, it is difficult to compare the rate of hydrolysis for each and every PASA against their respective P_i analogue. In this study, through a unique classification of P_ibiomolecules, we examined the consequence of As_i substitution on the kinetic stability of P_i species spanning from a simpler monoester form like pyrophosphate to a more complex phosphodiester variant like DNA.

1.1. Hydrolysis pathways

Pathways of the reaction for the hydrolysis of As_i and P_i are well documented in the literature [11–13]. The associative pathway is suggested to be dominant for both the esters [11] and follows the S_N2 mechanism with first-order kinetics. The initial step of the pathway as shown in figure 1 involves the attack of water as a nucleophile on the central atom (As/P) of the ester. It is found that the attack of water has the highest barrier in terms of activation energy and hence determines the rate of hydrolysis of the ester [6,11,12]. In addition to the nature of reactant and attacking nucleophile, the chemistry of the leaving group is also the key factor for the kinetics of nucleophilic substitution reactions. However, the complete mechanism of the reaction for the nucleophilic attack on different types of phosphate ester was investigated in detail by Lopez et al. who concluded that the rate determining step is always that for the attack of nucleophile rather than the departure of the leaving group [12]. In other words, the activation energy barrier for nucleophilic substitution is mainly influenced by attacking nucleophile and steric congestion around a central atom [14,15]. We employed the activation strain model [14] to relate the height of the activation energy barrier with the geometrical deformation and rigidity of reactants (figure 1). Our study aims at comparing the kinetic stability of As and P esters against water, hence we only modelled the reactants and transition state (TS) structures which belong to the rate determining step. In an attempt to understand the TS, it was found that after reaching the TS for monoanionic monoesters and diesters of P form a reaction intermediate, having a pentacoordinated centre with trigonal bipyramidal geometry [12]. However, direct product is observed in the case of dianionic monoester of P [12]. We also observed that the TS is followed by a pentacovalent intermediate, during the hydrolysis of all classes of PASA, irrespective of the type of ester and anionic form. The reaction via pentacovalent intermediate is further proceeded by internal proton transfer that leads to the breaking of the As/P-O bond [6,11,12]. The structures of reactant and TS as shown in figure 1 for all three classes were modelled (available with the electronic supplementary material). We carried out intrinsic reaction coordinate calculations (movie files available with the electronic supplementary material) to confirm coordinates of the reactant and intermediate for the initial step as shown in figure 1.

1.2. Classification

In order to determine the effect of As_i substitution on the P_i repository of the cell, we first need to understand the physiological functions of P_i species. P_i has a variety of roles inside a cell, and is incorporated into biomolecules mainly through kinases. This process is also known as phosphorylation and the reverse of this is de-phosphorylation, removal of P_i by hydrolases. This cycle of addition and removal of P_i governs biological phenomena such as activation or deactivation of proteins, synthesis of DNA and removing falsely paired nucleotide bases from it, breaking down hexoses into trioses, synthesis and utilization of ATP, etc. However, all of these happen at specific sites in a well-controlled way. Phosphorylation/de-phosphorylation on proteins usually takes place on the amino acids that have a hydroxyl group on their side chains, e.g. serine, threonine and tyrosine. Similarly, on hexoses the hydroxyl group at terminal carbon atoms is the site of phosphorylation/de-phosphorylation. Whereas in ADP, the hydroxyl group attached to the terminal beta phosphorus (Pβ) atom is used to synthesise ATP [9]. These specific sites are conserved throughout the organisms, and the rate of phosphorylation/de-phosphorylation at these sites is well synchronized with each biological process. We classified physiological P_i species based on the type of ester and conserved sites through which P_i anchors with biomolecules (table 1). The same sites also behaved as leaving groups during dephosphorylation or hydrolysis of P_i ester.

The modelled compound from each class is shown in figure 2. Our classification broadly covers most of the P_i biomolecules, Class 1A representing the energy currency of the cell, like ATP and ADP, while examples of Class 1B members are sugar or glycerol moiety containing biomolecules like Glucose-1-P_i, Fructose-1,6-Bi-P_i, Phosphoglycerate, etc. Whereas the molecules capable of storing genetic information belong to Class 2, like DNA and RNA. This classification not only enabled us to determine the relative stability of PASAs against hydrolysis with their respective P_i analogue, it also helped us to understand the consequence of As_i substitution at two different levels of biological esters (monoester and diester) and the role of conserved sites for phosphorylation/de-phosphorylation in the kinetics of hydrolysis.

2. Methods

For quantum chemical calculations, we need a suitable method to describe the thermochemical kinetics accurately for hydrolysis of inorganic systems like pyro-P_i/As_i and for bigger biomolecular systems like DNA that have non-covalent interactions. Zhao and Truhlars’ work recommended Minnesota density functionals for thermochemistry, thermochemical kinetics and non-covalent interactions [16]. The presence of a polar solvent (ɛ = 78.4, water) was mimicked with the polarizable continuum model using the integral equation formalism variant [17]. We carried out computations using M06-2X Minnesota density functional with a 6-31+G(2d,2p) basis set for the hydrolysis of pyroarsenate, and calculated the rate of the hydrolysis at 298 K using equation (2.1)

2.1

where c° = 1 and represents the free energy of activation. The calculated rate constant value of 0.08 s⁻¹ is in good agreement with the available experimental value of 0.05 s⁻¹ for the non-enzymatic hydrolysis of pyroarsenate. [18]

2.1. Modelling of P_i-DNA and As_i-DNA

The crystal structure of B-DNA dodecamer on PDB (id: 1bna.pdb) was used to construct a single stranded P_i-DNA model system consisting of d(TpTpApA) nucleotide bases. For modelling the As_i-DNA d(TAsTAsApA) structure, we replaced the respective two phosphorus atoms of P_i-DNA d(TpTpApA) with arsenic. It is difficult to deal with the whole system quantum mechanically (QM) due to the size complexity of P_i-DNA and As_i-DNA. Hence, QM was only applied on the portion which is involved in the reaction and the remaining part of the structure was dealt with using molecular mechanics (MM). We used two layers (our Own N-layer Integrated molecular Orbital molecular Mechanics) ONIOM [19], integrated into a Gaussian 09 [20] package, for QM/MM calculations. To avoid the problem of cancellation in ONIOM, a non-polar bond between two sp³hybridized C atoms, which is more than four bonded atoms away from the centre of the reaction, was selected to divide the structure into model layer (QM region) and low layer (MM region) as shown in figure 3.

The MM region was treated with AMBER [21] and frozen to maintain the structural constraints for diester backbone. For the QM region we chose M06-2X with 6-31+G(2d,2p) basis set which provides a good description of base stacking interaction for single-stranded DNA [22]. Our modelled P_i-DNA with M06-2X 6-31+G(2d,2p)/AMBER optimization of the geometry shows a reasonable geometrical similarity with the identical portion of B-DNA crystal structure with a RMSD value of 1.54 Å as shown in figure 4.

3. Results and discussion

The oxygen atoms bonded to central atom (As/P) of the esters are labelled (figure 2) to discuss and compare the parameters of the geometry between modelled reactants and TS compounds. Comparison of the main geometrical data between computed arsenate and phosphate models belong to Class 1A, 1B and 2 is presented in tables 2–4, respectively. Although the tabular data shows differences among the selected bonded parameters of the reactant and TS structures, these parameters were rarely used to relate the two different ester species. However, using the activation strain model we represented these geometrical changes in terms of the following energy terms ( energy, energy, energy), were basically derived from single point energy calculations [14]. The calculated values for distortion, interaction and activation energies for each species of three classes were shown in figure 5. Throughout the classification, irrespective of the anionic form, geometry of P_i esters were more deformed in TS compared to the As_i esters, which is reflected in the values of distortion energy (green bars) in figure 5. The higher deformity in TS for P_i esters might be due to the smaller size of the central P atom, compared to the bigger size of the As atom which allows more space for the incoming water nucleophile. Similarly, we also found that the geometry of dianionic P_i esters was more deformed than its monoionic form because two negatively charged oxygen atoms on the smaller central P atom possess a higher steric barrier for incoming nucleophile. However, the effect of higher deformity in the geometry of dianionic P_i monoester was not reflected in terms of the activation energy barrier due to the higher contribution from the interaction energy. In fact, it is the interaction energy which lowers the activation energy barrier for the hydrolysis of all the dianionic As_i as well as P_imonoesters. Interestingly, in the case of diester of class 2, the activation energy barrier achieved a comparable height equal to the interaction energy and this signifies the importance of diester linkages inside a living system. Compared to As_i-DNA a lot of energy is required to deform the structure of P_i-DNA during hydrolysis and this causes the biology to rely upon the smaller P atom for more rigid and stable geometry compared to the larger As atom. Apart from the comparison of geometrical deformity, the kinetics of hydrolysis for the esters is also discussed further in detail.

We calculated the rate constant for the hydrolysis of different classes of As/P-esters as shown in table 5. Both Class 1A and 1B dianionic ester of P hydrolyzed much faster than the monoanionic form as reported earlier [12]. We observed the same trend in the case of As esters of Class 1A and Class 1B. Higher stability of monoanionic ester over dianionic ester has biological significance. For instance, a stable molecule, e.g. DNA, with a monoanionic phosphodiester backbone, preserves the genetic information. On the other hand, comparatively facile hydrolyzable Class 1A and 1B esters are available physiologically in dianionic form to provide energetic requirements and assist in several other anabolic processes for an organism [1]. In comparison to the P_i counterpart of Class 1A and 1B, we found that none of the PASA have reasonable stability against hydrolysis. We calculated the rate for the hydrolysis of As_i-DNA to be 1.23 × 10⁻⁶ s⁻¹ which is slower than other PASA like pyroarsente and ribose-1-arsenate, which have first-order rate constant of 0.08 s⁻¹ and 0.07 s⁻¹, respectively, in dianionic form and 0.0003 s⁻¹ and 0.0007 s⁻¹, respectively, in monoanionic form. Higher kinetic stability of As_i-DNA against hydrolysis in comparison to other monoanionic esters of As_i belonging to Class 1A and 1B, was expected because of steric hindrance and base stacking interactions, [8] but it is still a long way from competing with the stability of P_i-DNA against hydrolysis which has a first-order rate constant value of 5.74 × 10⁻²⁰ s⁻¹. Another aspect of selecting P_i over As_i is that P_i is 10⁶ times more stable in monoanionic diester form than in dianionic monoester variants of Class 1B, whereas this difference was found to be narrower in the case of As_i. This difference might favour the integration of dianionic monoesters of P_i into a relatively more stable biomolecule like DNA. We also found that there is no significant difference among the rate constants of monoanionic pyro-As_i and ribose-1-As_i, monoionic pyro-P_i and ribose-1-P_i, similar to the dianionic species of both the esters, which signifies that the conserved sites through which these two esters bonded covalently have nothing to do with the kinetics of the hydrolysis. However, these conserved sites have well-defined regulatory roles in protein–protein interactions [23].

We reported that the different classes of PASA have much less stability against hydrolysis as compared to their P_i analogue species. Now, the rapid hydrolysis of PASA results in free As_i which is again going to compete with P_i. In order to avoid the structural confusion between As_i and P_i, the living organisms require metabolic machineries for the methylation of As_i that results into organoarsenical compounds which have one or more As-C bonds instead of As-O bonds. Methylation not only resolves the structural ambiguity between As_i and P_i, it might also prevent conversion of As_i into a highly reactive form which is known as Arsenite(III). On the other hand, organoarsenicals are highly stable against oxidation as well as against hydrolysis [24], as a result a number of marine organisms accept arsenolipids as a component of their cell membrane which are functionally, but not structurally, similar to phospholipids [25–27]. Arsenolipids primarily consist of As-C bonds which indicates the significant biological stability of As-C bond-based compounds as compared to kinetically unstable As-O bonds containing PASA.

4. Conclusion

Our study concluded that all classes of PASA have a higher hydrolysis rate compared to their respective P_i-analogues. In biological systems, where the rate of hydrolysis of P_i biomolecules is synchronously coupled with another biological process, any alteration in the kinetics of hydrolysis could have a severe effect on the metabolic activity of a cell. Although As_i-DNA has a higher resistance against hydrolysis compared to PASA of Class 1A and 1B, it is still a kinetically unstable molecule, and organisms would face severe consequences by allowing such a highly unstable molecule to substitute P_i-DNA for its genetic material. The rapid rate of the hydrolysis of dianionic PASAs raised another concern regarding their potential to reach and integrate into DNA. However, the mechanism behind the stability of organoarsenicals needs to be investigated for a further understanding of the metabolic process involving methylation of As_i.