Novel pathway for mutagenic tautomerization of classical А∙Т DNA base pairs via sequential proton transfer through quasi-orthogonal transition states: A QM/QTAIM investigation

In this paper we have theoretically predicted a novel pathway for the mutagenic tautomerization of the classical A∙T DNA base pairs in the free state, the Watson-Crick A·Т(WC), reverse Watson-Crick A·Т(rWC), Hoogsteen A·Т(H) and reverse Hoogsteen A·Т(rH) pairs, via sequential proton transfer accompanied by a significant change in the mutual orientation of the bases. Quantum-mechanical (QM) calculations were performed at the MP2/aug-cc-pVDZ//B3LYP/6-311++G(d,p) level in vacuum phase, along with Bader’s quantum theory of Atoms in Molecules (QTAIM). These processes involve transition states (TSs) with quasi-orthogonal structures (symmetry C1), which are highly polar, tight ion pairs (A-, N6H2-deprotonated)∙(T+, O4/O2-protonated). Gibbs free energies of activation for the A∙T(WC) / A∙T(rWC) ↔ A*∙Т(rwWC) / A*∙Т(wWC) tautomeric transitions (~43.5 kcal∙mol-1) are lower than for the A∙T(H) / A∙T(rH) ↔ A*N7∙Т(rwH) / A*N7∙Т(wH) tautomerisations (~53.0 kcal∙mol-1) (rare tautomers are marked by an asterisk; w—wobble configured tautomerisation products). The (T)N3+H⋯N1-(A), (T)O4+H⋯N1-(A) / (T)N3+H⋯N1-(A) and (T)O2+H⋯N1-(A) H-bonds are found in the transition states TSA-·T+A·T(WC)↔A*·T(rwWC) / TSA-·T+A·T(rWC)↔A*·T(wWC). However, in the transition state TSA-·T+A·Т(H)↔A*N7·T(rwH) / TSA-·T+A·Т(rH)↔A*N7·T(wH), the (T)N3+H⋯N7-(A), (T)O4+H⋯N7-(A) / (T)N3+H⋯N7-(A) and (T)O2+H⋯N7-(A) H-bonds are supplemented by the attractive (T)O4+/O2+⋯N6-(A) van der Waals contacts. It was demonstrated that the products of the tautomerization of the classical A∙T DNA base pairs—A*∙Т(rwWC), A*N7∙Т(rwH) and A*N7∙Т(wH) (symmetry Cs)–further transform via double proton transfer into the energetically favorable wobble A∙T*(rwWC), A∙T*(rwH) and A∙T*O2(wH) base mispairs (symmetry Cs).


Introduction
Investigation of microstructural mechanisms for mutagenic tautomerization of the Watson-Crick DNA base pairs occupies an important place in molecular biophysics and molecular biology, enabling an understanding of the nature of genome instability [1][2][3][4][5]. This follows from the 'rare tautomer hypothesis' proposed by Watson and Crick [1] shortly after they established the spatial architecture of DNA [2]. However, achievements in this area remain rather modest despite its long history [6], encouraging further research in this direction.
Löwdin [3,4] first proposed that the electronic structure of the Watson-Crick (WC) DNA base pairs AÁT(WC) and GÁC(WC) permits their transition into the high-energy tautomerized states A Ã ÁT Ã (L) and G Ã ÁC Ã (L), now called Löwdin (L) base pairs. Here and henceforth, rare (in particular mutagenic) tautomers are marked with an asterisk and differ from each other by the location of a particular proton: in the A Ã rare tautomer proton bonds to N1 nitrogen atom, A Ã N7 -to N7 nitrogen atom; T Ã -to O4 oxygen atom and T Ã O2 -to O2 oxygen atom. Löwdin proposed that the AÁT(WC)$A Ã ÁT Ã (L) and GÁC(WC)$G Ã ÁC Ã (L) transitions occur by double proton transfer (DPT) along neighboring intermolecular hydrogen (H) bonds via proton tunneling. These ideas have been prominent in the field of quantum biology and attracted much theoretical study of the mechanisms of spontaneous transitions and transversions arising during DNA replication [7][8][9][10][11][12][13][14].
Recently, it has become clear that Löwdin's mechanism does not provide the generation of sufficiently long-lived mutagenic tautomers of the DNA bases, which escape from the replicative DNA-polymerase transforming into their canonical tautomeric forms. The root cause of this observation is the absence of the reverse barrier of tautomerization ΔΔG in the AÁT(WC) DNA base pair and its small value in comparison with kT (0.62 kcalÁmol -1 under normal conditions) for the GÁC(WC) DNA base pair [8,9,[15][16][17][18].
In previous papers [19][20][21][22][23][24][25][26][27] we proposed an alternative mechanism for mutagenic tautomerization of the AÁT(WC) and GÁC(WC) base pairs into the corresponding wobble base mispairs and vice versa, which mechanism obviates the above difficulties. The chief difference of our mechanism from the Löwdin mechanism is that, in the process of mutagenic tautomerization through sequential proton transfer, the DNA bases shift laterally relative each other into the DNA minor or major grooves, leading to the wobble configuration which contains the mutagenic tautomers [19]. Moreover, it turned out that a similar mechanism works also for the mutagenic tautomerization of purine-purine [21], pyrimidine-pyrimidine [22,23] and purine-pyrimidine [24][25][26][27] DNA base mispairs, which are active players in the field of spontaneous point mutagenesis.
This allows us to assume that it is the intrapair tautomeric transition of the wobble pairs from the main tautomeric form into the rare one with a WC configuration or close to it, and vice versa, which is the key to understanding the microstructural mechanisms for spontaneous transitions and transversions during DNA biosynthesis [19][20][21][22][23][24][25][26][27]. Theoretical analyses of such mechanisms have been experimentally confirmed in part for the AÁC(w) and GÁT(w) purinepyrimidine pairs [28][29][30][31].

Computational methods
The geometries of all the investigated DNA base pairs and transition states (TSs) were optimized using the Gaussian'09 package [50]. The B3LYP/6-311++G(d,p) level of theory [51-55] was used. This level of theory has successfully proved itself for calculations of similar systems [56][57][58][59][60][61][62][63]. The study included harmonic frequency calculations (using a scaling factor of 0.9668 [64][65][66]) and intrinsic reaction coordinate (IRC) analysis in the forward and reverse directions from each TS using a Hessian-based predictor-corrector integration algorithm [67] at the B3LYP/6-311++G(d,p) level of theory successfully applied in the previous studies [16,17,68,69]. Local minima and TSs (localized by the synchronous transit-guided quasi-Newton method [70]) were confirmed as such by the absence or presence, respectively, of one imaginary frequency. Standard TS theory was applied to estimate the activation barriers for the tautomerisation reactions [71]. Single point electronic energy calculations were performed for the B3LYP geometries at the MP2/aug-cc-pVDZ level of theory [72,73]. MP2 has been successfully applied to gain chemical information about similar proton transfer reactions in DNA systems [74][75][76][77][78][79]. The choice of the MP2 level of theory is caused by the insignificant errors in comparison with CCSD(T) method, that was convincingly shown in the benchmark works of Hobza and Šponer [80,81].
We have performed investigations for the isolated H-bonded pairs of nucleotide bases, that adequately reflects the processes occurring in real duplex environment [14,30,31]. At this we relied on experience received in the previous works [82][83][84][85] devoted to related topics and systems, where the negligibly small impact of the stacking and sugar-phosphate backbone on the tautomerisation processes has been shown.
The Gibbs free energy G for all structures was obtained in the following way: where E el = electronic energy, while E corr = thermal correction to Gibbs free energy. Electronic interaction energies ΔE int were calculated at the MP2/6-311++G(2df,pd) level of theory as the difference between the total energy of the base pair and energies of the monomers and corrected for the basis set superposition error (BSSE) [86,87] through the counterpoise procedure [88,89] without consideration of the deformation energies of the monomers due to their relatively small values [90].
The energies of the attractive van der Waals contacts [101,102] in TSs for tautomeric transitions of the base pairs were calculated by the empirical Espinosa-Molins-Lecomte (EML) formula [103,104] based on the electron density distribution at the (3,-1) BCPs of the specific contacts: where V(r) = value of a local potential energy at the (3,-1) BCP. Energies of conventional AHÁÁÁB H-bonds were evaluated by the empirical Iogansen formula [105]: where Δν = magnitude of the stretching frequency shift for the AH H-bonded group involved in the AHÁÁÁB H-bond relative to the unbound group. Partial deuteration was applied in order to avoid the effect of vibrational resonances [106][107][108][109][110][111][112][113][114]. The atom numbering scheme for the DNA bases is as per convention [108]. Conformers of the AÁT base pairs remain plane symmetric structures along the entire IRC of tautomerization. This also holds for base pairs tautomerising via proton transfer along intermolecular H-bonds as per currently known mechanisms for mutagenic tautomerization of WC pairs [16,17,19,49].

Results and discussion
tautomerisation reactions occur via the initial migration of proton localized at the N6 atom of the N6H 2 amino group, leading to the formation of the A + ÁTion pair and significant change of the mutual orientation of the bases within the pair, i.e. mutual transformation of the cys / trans$trans / cys-orientation of the N1H and N9H bonds relative to each other (Fig 1). Our new mechanism is controlled by the TSs having quasi-orthogonal structures (symmetry C 1 ). Further proton transfers to the N1/N7 nitrogen atom causing the rotation of the base and formation of the terminal wobble base mispair. Each of these tautomeric conversions is followed by the asynchronous DPT along the intermolecular H-bonds in the wobble base mispairs (Fig 2).
Interestingly, the energy of the intermolecular specific contacts (H-bonds and attractive van der Waals contacts) constitute only a minor part of the electronic energy of monomeric interactions for all these H-bonded structures (~14-0.87%) (see Figs 1 and 2). This agrees with previous results for other H-bonded base pairs [112].
All tautomeric transitions in this work are dipole-active, being accompanied by significant changes in dipole moment of the tautomerizing structures along the IRC (0.38-12.65 D), achieving maximum values for each tautomeric transition at its TS (7.38, 6.83, 12.65 and 10.636 D, accordingly) ( Table 2). The Gibbs free energy of activation for the AÁT(WC)/AÁT (Figs 1 and 2).
Note that only one case of mutagenic tautomerization, the AÁT(WC)$A Ã ÁT(rw WC ) reaction, occurs by participation of the dynamically unstable intermediate A Ã ÁT Ã (L) (a Löwdin's base pair [3,4]). The other three AÁT DNA base pairs-AÁT(rWC), AÁT(H) and AÁT(rH)-do not tautomerise via the Löwdin's mechanism. For these three pairs, the local minima corresponding to the tautomerized A Ã ÁT Ã O2 , A Ã N7 ÁT Ã and A Ã N7 ÁT Ã O2 base pairs are absent on the energy hypersurface. This observation is independent of the level of QM theory used.
It should be noted that three out of four tautomerization processes of the AÁT base pairs do not complete with formation of the A Ã ÁT(rw WC ), A Ã N7 ÁT(rw H ) and A Ã N7 ÁT(w H ) mispairs (Fig 2  and Table 1). These plane-symmetric wobble pairs (symmetry C s ) tautomerise further via the DPT mechanism along neighboring intermolecular H-bonds into the energetically-favorable plane-symmetric AÁT Ã (rw WC ), AÁT Ã (rw H ) and AÁT Ã O2 (w H ) DNA base mispairs, respectively (Fig  2, Tables 1 and 2). These processes occur via a concerted asynchronous mechanism with proton transfer along the intermolecular (T)N3HÁ Á ÁN6(A) H-bonds, which, in fact, is a rate-limiting stage. It is noteworthy that the A Ã N7 ÁT(rw H )!AÁT Ã (rw H ) and A Ã N7 ÁT(w H )!AÁT Ã O2 (w H ) tautomerisations are barrier-less (ΔΔG TS = -1.98 and -1.97 kcalÁmol -1 ) ( It is thus possible to say that the tautomerization processes described here terminate with the mutagenic tautomerization of both T and A DNA bases with further formation of the classical mutagenic tautomers T Ã , T Ã O2 [16,19,27,64,102,106] and A Ã [16, 19-22, 25, 26, 100, 114], respectively. In this case, the A Ã N7 ÁT(rw H )$AÁT Ã (rw H ) and A Ã N7 ÁT(w H )$AÁT Ã O2 (w H ) tautomeric equilibria are completely shifted to the right. For the two other cases, the following conformational transitions (MP2/aug-cc-pVDZ//B3LYP/6-311++G(d,p) level of theory in the continuum with ε = 1 at T = 298.15 K) are presented below complexes in brackets. Dotted lines indicate AHÁÁÁB H-bonds and attractive AÁÁÁB van der Waals contacts-their lengths HÁÁÁB and AÁÁÁB are presented in angstroms (for their more detailed physico-chemical characteristics see Table 2); carbon atoms are in light-blue, nitrogen-in dark-blue, hydrogen-in grey and oxygen-in red.  Notably, the methyl group of the T DNA base does not change its orientation during all these tautomerisation processes without exception. Moreover, the heterocycles of the DNA bases remain planar, despite their ability for out-of-plane bending [118][119][120]. A

Conclusions and perspectives
Novel pathways for mutagenic tautomerization of four classical AÁT DNA base pairs, followed by the significant changes of base orientation within the pair, have been predicted by these QM results. The transition states with quasi-orthogonal structure (symmetry C 1 ) are highly polar Table 1 Table 2 tight ion pairs (A -, N6H 2 -deprotonated)Á(T + , O4/O2-protonated). The tautomerization products-the A Ã ÁT(rw WC ), A Ã N7 ÁT(rw H ) and A Ã N7 ÁT(w H ) pairs-further transform via concerted asynchronous double proton transfer into the energetically favorable wobble AÁT Ã (rw WC ), AÁT Ã (rw H ) and AÁT Ã O2 (w H ) mispairs (symmetry C s ), respectively. Moreover, it was established in our recent papers, that wobble A Ã ÁT(rw WC ) base mispair can also be formed from the reverse AÁT(rWC) base pair [20], A Ã ÁT(wWC) base mispair-from the canonical AÁT(rWC) base pair [19] and A Ã N7 ÁT(w H ) base mispair-from the Hoogsteen AÁT(H) base pair [20]. We are currently engaged in elaborating this topic in order to discover biologically important H-bonded nucleobase pairs, for which the mechanism of mutagenic tautomerization plays a key role. Moreover, we suggested that novel mechanism of mutagenic tautomerization presented in this study could lead to the conversion of an anti-parallel DNA helix to a parallel DNA helix. We also consider investigation of these tautomerisation mechanism by the participation of the modified AÁT base pairs [121][122][123][124] as a task for the future.   [105], Espinose-Molins-Lecomte [103,104] (marked with an asterisk) or Nikolaienko-Bulavin-Hovorun [109] (marked with double asterisk) formulas, kcalÁmol -1 .