Abstract

We studied for the first time 16 tautomers/rotamers of diphosphocytosine by four computational methods. Some of these tautomers/rotamers are isoenergetic although they have different structures. High-level electron correlation MP2 and MP4(SDQ) ab initio methods and density functional methods employing a B3LYP and the new M06-2X functional were used to study the structure and relative stability of 16 tautomers/rotamers of diphosphocytosine. The dienol tautomers of diphosphocytosine are shown to be much more stable than the keto-enol and diketo forms. The tautomers/rotamers stability could be ranked as PC3 = PC12 < PC2 = PC11 < PC1 < PC10 < PC8 < PC9 < PC15 < PC16 < PC6~PC7 < PC13 < PC4~PC14 < PC5. This stability order was discussed in the light of stereo and electronic factors. Solvation effect has been modeled in a high dielectric solvent, water using the polarized continuum model (PCM). Consideration of the solvent causes some reordering of the relative stability of diphosphocytosine tautomers: PC3~PC12~PC2~PC11 < PC1 < PC10 < PC8 < PC9 < PC15~PC16 < PC13 < PC6~PC7~PC14 < PC4~PC5.

1. Introduction

DNA is a naturally occurring biological macromolecule, containing thousands of nucleic acid bases, and it is of prime importance in genetic determination [1]. The five nucleic acid bases, cytosine, thymine, uracil, adenine, and guanine, found in DNA and RNA control the replication of DNA, store information required to synthesize proteins, and translate this information to the protein. Tautomerism is a well-known phenomenon occurring in nucleic acid bases [216], in which proton transfer from the heterocyclic ring center to an exocyclic oxo- or imino- group leads to the formation of either an –OH or an –NH2 groups. These processes are known as keto-enol or imino-amino tautomerism, respectively. Tautomerism, a sort of isomerism, plays an important role in organic chemistry, medicinal chemistry, pharmacology, and molecular biology. Tautomerism partially explains the structure of nucleic acids and their mutations [17]. In DNA bases, tautomerization results in altered base pairing configurations or mispairing due to changes in hydrogen-bonding capabilities. DNA mutations are likely to be caused by such alterations. Cytosine is one of the building pyrimidine nucleobases of RNA. A large amount of experimental [1834] and theoretical work [3558] has been carried out in order to elucidate the structure of cytosine and its tautomers. Chemically modified bases have attracted extensive interest due to their numerous pharmacological, biochemical, and biological capabilities [5962]. Replacing nitrogen in pyrimidine and purine DNA bases by the biological element of interest, phosphorus, may produce new synthesizable bases of special interest. Modification of cytosine base to phosphorus-containing cytosine analogues, diphosphocytosine, is a good start point in this topic. Also, their existence in many tautomeric forms, like other nucleoside bases, seems to be crucial in order to explain the mutation occurring during DNA duplication [6366]. Searching the literatures confirms that no previous studies include all numbers of cytosine tautomers. This is the first time to include 16 isomers of cytosine. Similar to cytosine, the diphosphocytosine should also have 16 tautomers. There are many interesting questions on the diphosphocytosine tautomers concerning, for example, their relative stability order, energy criteria, proton transfer processes, electronic structure, and chemical properties. In an attempt to answer some of these questions using computational studies, some DFT methods (B3LYP, M06-2X) and ab initio MP2 and MP4(SDQ) methods are chosen to study the stability order and tautomerization processes of all the possible 16 diphosphocytosine tautomers in the gas phase as well as in water. Compared to the popular B3LYP method, studies of tautomerism using the recent method M06-2X are very rare and it is worthwhile to consider them here in this work to assess their accuracy compared to MP4(SDQ) data.

Among cytosine tautomers, experimental investigations performed in gas phase, solid state, or in solution indicate that the canonic amino-oxo tautomer of cytosine is the most stable one [67]. One of our interesting computational results [67, 68] of diphosphocytosine tautomers is the reverse stability order as compared to the parent analogue, cytosine. The origin of this stability order needs to be interpreted as it is not fully discussed. One of the main objectives of this work is to investigate this behavior on the basis of the unique bonding properties of the phosphorous atom.

2. Computational Methods

All quantum chemical calculations were performed using the Gaussian 09 suite of programs [69]. Geometry optimizations for all tautomers have been performed using ab initio density functional theory (DFT) at the B3LYP functional [7072] and the new M06-2X functional in conjunction with the 6-31G(d) basis set. For each stationary point, we carried out a force constant harmonic frequency calculation at the same levels to characterize their nature as minima or transition states and to correct energies for zero-point energy and thermal contribution. The frequencies are scaled by a factor of 0.98. The vibrational modes were animated using the ChemCraft program [73]. Natural charges were computed within full Natural Bond Orbital Analysis, using NBO implemented in Gaussian 09 [69]. The vibrational modes were examined by using the ChemCraft program [73]. Partial charge distributions were calculated using the natural population analysis (NPA) method [74]. Solvation effect has been modeled in water using the polarized continuum model (PCM) of Tomasi et al. [7578] at the B3LYP/6-31+G(d,p) level for all tautomers/rotamers.

3. Results and Discussion

3.1. Geometrical Aspects of Diphosphocytosine Tautomeric Forms

Geometry of the 16 tautomers and rotamers was optimized at four different levels. Optimized structures of the 16 diphosphocytosine tautomers (at MP4(SDQ)/6-31+G(d)) are presented in Figure 1. For the notations of the tautomers and rotamers of diphosphocytosine, we use PC1, PC2, PC3, 16 for the sixteen forms. All structures with their notations and atomic numbering are shown in Figure 1. Structure geometrical parameters of all the 16 tautomers are listed in Table 1. Bond length changes among tautomers show the same trend as expected from formal bond orders. In tautomers PC2, PC3, PC11, and PC12, all CC and CP bond lengths are in the range characteristic for aromatic rings. C1-C2 bond length changes among tautomers depending on the amount of π-bond character the tautomer has. In tautomers PC2, PC3, PC11, and PC12 this bond distance is about 1.40 Å while in the rest of the tautomers this distance increased as shown in Table 1. In all the dienol forms, the P3-C2 and P3-C4 and also P5-C6 and P5-C4 distances are nearly similar while those for the other keto-enol and diketo forms show considerable differences of about 0.147 Å. This reflects the extent of conjugation and the increasing aromatic character of the six-membered ring. An inspection of the dihedral angles in Table 1 proves that all the tautomers/rotamers structures are nonplanar. It is reported that [79] results of electron correlation calculations (MP2/6-31G*) give nonplanar geometries for several nucleotide bases specifically, for cytosine, which showed a significant pyramidalization of the amino group from MP2/6-31G* calculations. Using enormous computer resources, a complete MP2/6-31G** frequency calculation was carried out [80] for cytosine, confirming that the planar structure is a saddle point. Bond angles in the ring show quite dramatic changes between tautomers and rotamers, by up to more than 10°. Owing to the lack of polarization functions on hydrogens, geometry parameters of the latter are expected to show larger differences, but these less important parameters were not listed in.

In the dienol form of diphosphocytosine, there are three double bonds giving it the six π-electrons needed to gain its aromatic character as in case of benzene. The phosphorous atom has an additional lone pair of electrons but not involved in the conjugated system. Once the keto form is formed in the diphosphocytosine, one resonance form is lost as the phosphorus atom is strongly pyramidal in the keto forms and its lone pairs are no longer available for conjugation. The extent of pyramidalization emerges from the sum of bond angles around the phosphorus atom. An extra destabilization occurs upon forming the second keto tautomer. Accordingly, diketo form is accompanied by destruction of conjugated π-electrons system in the ring of the diphosphocytosine. In case of cytosine, on the other hand, formation of keto or diketo forms still keeps planarity at the nitrogen atoms where they still able to resonate with the adjacent double bonds. Aromaticity is lost in both cytosine and diphosphocytosine, but the number of resonance structures decreased in the case of the latter system.

3.2. Relative Energies and Thermodynamic Parameters of Tautomers

Since our major objective of this study was to obtain information on the relative stabilities of diphosphocytosine tautomers and rotamers, we tried to go as high as possible, with the level of theory and use both DFT and high correlation ab initio MP2 and MP4(SDQ) methods. The results on the relative stabilities of diphosphocytosine tautomers and rotamers using different computational levels are listed in Table 2. For a better view, this is also shown in Figure 2. The results of the four computational methods are consistent with slight differences. The 16 diphosphocytosine tautomers and rotamers could be classified, on the basis of total relative energy, into three groups: (1) low-energy tautomers group containing four tautomers namely, PC2, PC3, PC11, and PC12. These tautomers: are those satisfying the aromaticity of the pyrimidine ring (the benzene-like structure), (2) moderate-energy tautomers group containing also four tautomers, PC1, PC8, PC9, and PC10. (3) high energy tautomers group containing the rest eight tautomers, PC4, PC5, PC6, PC7, PC13, PC14, PC15, and PC16. DFT methods predict PC2 to be the most stable lowest-energy tautomer. High correlation ab initio methods, on the other hand, predict the most stable is the hydroxy-amino form PC3. Energy difference is less than 0.50 kcal/mol. The structural difference between the two tautomers/rotamers is the proximity of the positively charged H to the positively charged phosphorus atom. The same argument is applied to PC11 and PC12, where the positively charged H is repulsively interacting with P3 (0.503, Table 3) in PC11 and P5 (0.531, Table 3) in PC12. Previous studies suggest that the canonical amino-oxo structure is the most stable one [81]. PC5 is the most unstable highest-energy tautomer as predicted by all methods. This could be due to the destruction of the aromaticity and also attributed to the repulsive interaction of the positively charged H and the great positive charge accommodated on P3 (0.657, Table 3). As could be seen from Table 2 that including zero-point energy and thermal energy correction has a considerable effect on the molecular energy. The corrected energy is lower than the uncorrected by about 4.0 kcal/mol. Also, considering these corrections causes reordering of the tautomers stability. Figure 3 shows the MP2-corrected for zero-point and thermal energy versus uncorrected one for all the 16 diphosphocytosine tautomers. Based on the corrected energy, both B3LYP and M06-2X have the same order as PC2 < PC3 = PC12 < PC11 < PC1 < PC10 < PC8 < PC9 < PC15 < PC16 < PC7 < PC6 < PC13 < PC5 ~ PC4 ~ PC14. On the other hand, the relative stability order as predicted by ab initio Moller Plesset MP2 and MP4(SDQ) is PC3 = PC12 < PC2 = PC11 < PC1 < PC10 < PC8 < PC9 < PC15 < PC16 < PC6 ~ PC7 < PC13 < PC4 ~ PC14 < PC5. However, MP4(SDQ) method shows slight differences like PC15 ~ PC16 and PC4 < PC14 ~ PC5. Figure 4 presents tautomers’ order by DFT and MP2 methods. Enthalpy and free energy are consistent with the total energy as can be detected from Table 2. One of the electronic factors that affect the stability of the tautomers is the repulsion between the neighboring –OH and –PH or –CH groups, while the other factor is the attraction between the O–H or P–H bonds and the lone pair of the neighboring phosphorous or oxygen atoms. Another factor is the attractive interaction between the O–H or P–H bonds, on one side, and the lone pair(s) of the adjacent P or O atoms, on the other one [82].

The amino N accommodates the highest negative charge (up to −0.904); then the O atom comes in the second place (up to −0.778). Positive charges are accommodated on the P and H atoms. Very little negative charge on H bonded to P atoms as appeared in PC4, PC6, and PC8. The values of the NPA charges at hydrogen atoms of –OH groups are larger than the values at the hydrogen atom of the –CH group (see Table 3). As a result, the repulsion between the neighboring –OH and –CH groups is stronger than the repulsion of the neighboring –PH and –CH groups.

3.3. Solvation Study

One of the most crucial factors determining the tautomer distribution in the biological media is the environment. Solvent preferentially stabilizes structures with localized charges. Solvation effect has been explored in water (ε = 78.54) using the polarized continuum model (PCM) at the B3LYP/6-31+G(d) level for all the 16 tautomers. Table 4 collects the relative energies for different tautomers/rotamers (kcal/mol) at B3LYP/6-31+G(d) in water. In water, the most stable diphosphocytosine tautomers are PC3 and PC12 and the most unstable tautomers/rotamers are PC4 and PC5. As seen from Table 4, the relative stabilities in water go in the order PC3 ~ PC12 ~ PC2 ~ PC11 < PC1 < PC10 < PC8 < PC9 < PC15 ~ PC16 < PC13 < PC6 ~ PC7 ~ PC14 < PC4 ~ PC5 at B3LYP/6-31+G(d) (PCM/water). This order is slightly different from that predicted at the same level in gas phase (ε = 1). For example, PC3 ~ PC12 ~ PC2 ~ PC11. This result is easily understandable on the basis of the calculated dipole moment of these tautomers. Polar tautomers/rotamers are favored in polar solvents.

The NPA charges computed for tautomers/rotamers geometry optimized in water are given in Table 5.

The NPA charges computed for tautomers/rotamers geometry optimized in water are generally smaller than those computed for tautomers/rotamers optimized in gas phase. This is understandable as a role of the high dielectric constant solvent, water. The amino N accommodates the highest negative charge (up to −0.846), then the O atom comes in the second place (up to  −0.707). Positive charges are accommodated on the P and H atoms. Very little negative charge on H bonded to P atoms as appeared in PC4, PC6, and PC8. The values of the NPA charges at hydrogen atoms of –OH groups are larger than the values at the hydrogen atom of the –CH group (see Table 5). As a result, the repulsion between the neighboring –OH and –CH groups is stronger than the repulsion of the neighboring –PH and –CH groups.

4. Summary

On the basis of the present theoretical study, we may confirm the following conclusions.(i)We studied for the first time 16 tautomers/rotamers of diphosphocytosine by four computational methods. Some of these tautomers/rotamers are isoenergetic although they have different structure.(ii)From the total energy point of view, the 16 tautomers/rotamers are categorized into three main groups. These are low-energy difference group, moderate-energy difference group, and high-energy difference group.(iii)PC3 and PC12 are the most stable tautomers/rotamers.(iv)Our systematic calculations show that the energy difference between components of a rotamer pair is very stable, independent of the level of theory. This leads to the important conclusion that besides the dominant PC3 form, its pair PC2 may lie not higher than ~0.4 kcal mol−1; therefore, any experimental study should consider this possibility when evaluating spectroscopic experimental studies.(v)Benzene-like structures are the most stable tautomers/rotamers. (vi)Charge is distributed on the ring atoms; −ve charges are on N and O, while +ve charges on P and H.(vii)Solvation does not affect the order significantly; PC3 ~ PC12 ~ PC2 ~ PC11 < PC1 < PC10 < PC8 < PC9 < PC15 ~ PC16 < PC13 < PC6 ~ PC7 ~ PC14 < PC4 ~ PC5.

Acknowledgment

This work was done at the Vienna Supercomputer Centre of the University of Vienna, Austria.