Introduction {#Sec1}
============

Genetic mutations, alterations in the DNA sequence, have been identified as the cause of genetic disorders, tumor development and drug resistance in pathogenic microorganisms^[@CR1]--[@CR3]^. There are several types of mutations that can occur in DNA, such as copy number variations, duplications, deletions, insertion and single base substitutions, and the latter is the most common type (two-thirds of human genetic diseases are due to single-base alterations)^[@CR4]^. Therefore, research on the mechanism of point mutations and single-base editing gene therapies has become an attractive topic^[@CR4]--[@CR6]^.

In double-helical DNA macromolecules, the complementarity between pyrimidines (cytosine and thymine) and purines (guanine and adenine) are due to hydrogen bonds. Each hydrogen bonding is essentially a proton trapped in an asymmetric double well potential between two electronegative atoms^[@CR7]^. It has been hypothesized that proton transfer within the double well potential is an important mechanism of base substitution^[@CR7],[@CR8]^. Tautomerism in guanine--cytosine base pairs resulting from intermolecular proton transfer has been suggested to be responsible for the universal guanine--cytosine to adenine--thymine mutation frequently observed in bacteria, fungi, plants, and animals^[@CR9]^. The tautomers, indicated by an asterisks (C\* and G\*), resulting from the transfer of cytosine H~4a~ to guanine and the back transfer from guanine H'~1~ to cytosine cannot form hydrogen-bonded complexes with their natural counterparts but can form mismatches with the bases A and T, respectively, during DNA replication (Fig. [1](#Fig1){ref-type="fig"})^[@CR7],[@CR10]^. Under normal circumstances, the relatively high energy barrier within the double well potential prevents the forward double-proton transfer (DPT) reaction, which accounts for the stability of DNA and the fidelity with which the genetic code is preserved and transmitted to daughter cells^[@CR11],[@CR12]^. However, when the quantum nature of the proton is taken into account, the proton has a small yet finite probability of tunneling through the barrier^[@CR12]^. When such double proton tunneling occurs during DNA synthesis, single base substitutions will occur. Fu *et al*. showed that among several possible DNA mutation processes, only the DPT mechanism could fully explain the universal mutation bias^[@CR9]^.Figure 1Structures and atomic numberings of the base pairs generated by GaussView 5.0. Top: Canonical C-G base pair and C\*-G\* base pair resulting from DPT; bottom: tautomers C\* and G\* mismatched to the bases A and T, respectively.

Previous studies have focused on various environmental factors that might facilitate intermolecular proton transfer within base pairs. Zhang, J. D. *et al*. showed that when an additional hydride is placed on the C~6~ and C~4~ positions of cytosine, the anionic complex formed, which facilitates guanine H~1~ transfer to the N~3~ site of cytosine^[@CR13]^. Noguera, M. *et al*. found that protonation at the N~7~ and O~6~ sites of guanine, affording H^+^ G~N7~C and H^+^G~O6~C, strengthens the binding of the base pair and facilitates the N~1~--N~3~ single-proton-transfer reaction^[@CR14]^. The influences of metal cation (Cu^+^, Ca^2+^ and Cu^2+^) coordination to the CG and AT base pairs on intermolecular proton transfer were studied using density functional theory (DFT) methods^[@CR15],[@CR16]^. Alya, A. A. *et al*. investigated the DPT reaction of a GC base pair under the effect of uniform electric fields on the order of 10^8^ to 10^9^ Vm^−1^. They considered that fields applied along the axis of the double proton transfer in the -x (defined in the C to G direction) direction favor the canonical over the rare tautomers^[@CR12]^. Chen, H. Y. *et al*.^[@CR17]^ showed that the substantial effect of GC stacking originates from the electrostatic interactions between the dipoles of the outer GC base pairs and the middle GC^•−^ base-pair radical anion, the extent of the charge delocalization is very small and has little effect on proton transfer in GC^•−^. The effect of the surrounding water molecules on the DPT in GC was investigated using DFT methods, and the results demonstrate that water is crucial to the proton reactions. It does not act as a passive element but actually catalyzes the DPT^[@CR18]^.

Thus far, few studies have examined the effects of neutral cytosine analogues on the DPT reaction in CG base pare. The purpose of this study was to use existing physic-chemical understanding and predictive capability to identify neutral cytosine analogues that can induce DPT reactions with guanine under physiological conditions and provide new strategies for gene therapy.

Results and Discussions {#Sec2}
=======================

Molecular structures of the modeled cytosine analogues {#Sec3}
------------------------------------------------------

To provide minimized structures to mimic the Watson--Crick base pairing in duplex DNA or RNA, the sugar-attachment sites of cytosine and guanine were methylated^[@CR19]^. To build a library of cytosine analogues, multiple atoms on the cytosine were replaced. To improve its ability to donate protons, the amino group at C~4~ was replaced by a hydroxy moiety. To improve the ability of N~3~ to accept protons, the carbonyl oxygen at C~2~ was replaced by a hydrogen, and separately, the C~2~ carbon atom was replaced by a boron atom. In addition, C~5~ was replaced by a nitrogen atom. The carbon atom at C~6~ was replaced with methine, methylene, and carbonyl moieties. A total of 12 cytosine analogues were modeled (Fig. [2](#Fig2){ref-type="fig"}).Figure 2Molecular structures of the modeled cytosine analogues (Ca). The atoms that had being replaced were shown in red. The box showing the conformers of Ca~1~ and Ca~2~.

Energy barriers and rate constants of the intramolecular proton transfer {#Sec4}
------------------------------------------------------------------------

First, the structures of all the monomers, Ca~0-12~, were optimized by DFT M06-2X/def2svp method. The vibration analysis showed no imaginary frequencies in those monomers, indicating that the structures are stable at the current calculation level. Because both the O~4~ atom's ability to provide protons and the N~3~ atom's ability to acquire them are enhanced, the modeled cytosine analogues may undergo intramolecular proton transfer (H~4~ transfer from O~4~ to N~3~). Therefore, the intramolecular single proton transfer (SPT) reactions of all monomers were investigated. The calculated SPT energy barriers (Δ*G*^≠^) within the cytosine analogues ranged from 24.51 to 32.23 kcal/mol (Table [1](#Tab1){ref-type="table"}). The forward reaction rate constants (calculated from Δ*G*^≠^) of the intramolecular proton transfer process range from 3.40 × 10^−5^ to 1.23 × 10^−10^ s^−1^, which indicate that the intramolecular SPT process are not facile. Therefore, the barriers of all cytosine analogues are relatively high for the proton to hop, indicating that these molecules can be used as candidate molecules to induce DPT reactions with guanine.Table 1The values of the barrier of the intramolecular SPT reactions (Δ*G*^≠^), the energy difference between analogues and their conformers (ΔE), the forward and reverse barriers of the intermolecular DPT reactions (Δ*G*^≠^~*f*~ and Δ*G*^≠^~*r*~), the forward and reverse reaction rate constants of intermolecular DPT reactions (k~f~ and k~r~), and the equilibrium constants (K~eq~). All energy barrier values are given in kcal/mol.Ca~*n*~Δ*G*^≠^ΔEΔ*G*^≠^~*f*~Δ*G*^≠^~*r*~k~f~ (s^−1^)k~r~ (s^−1^)K~eq~Ca~0~43.13---8.31/8.80/7.600.22/0.71/−0.494.05 × 10^6^2.04 × 10^12^1.99 × 10^−6^Ca~1~30.71−2.991.071.601.14 × 10^12^4.81 × 10^11^2.37 × 10^0^Ca~2~32.23−0.490.871.161.57 × 10^12^9.83 × 10^11^1.60 × 10^0^Ca~3~28.19−4.45−1.203.794.53 × 10^13^1.38 × 10^10^3.28 × 10^3^Ca~4~31.50−4.59−1.54/−0.49/−1.15−0.08/0.97/0.311.43 × 10^13^1.34 × 10^12^1.07 × 10^1^Ca~5~29.86−2.97−0.624.241.77 × 10^13^6.63 × 10^9^2.67 × 10^3^Ca~6~31.10−0.47−0.353.671.14 × 10^13^1.67 × 10^10^6.83 × 10^2^Ca~7~28.61−4.27−1.336.015.60 × 10^13^3.75 × 10^8^1.49 × 10^5^Ca~8~31.68−4.55−0.484.951.41 × 10^13^2.09 × 10^9^6.75 × 10^3^Ca~9~27.13−3.230.005.076.46 × 10^12^1.72 × 10^9^5.09 × 10^3^Ca~10~28.19−0.52−0.205.228.94 × 10^12^1.35 × 10^9^6.62 × 10^3^Ca~11~24.51−4.86−2.028.061.72 × 10^14^1.34 × 10^7^1.28 × 10^7^Ca~12~28.36−4.89−0.745.622.15 × 10^13^7.06 × 10^8^3.05 × 10^4^

Due to the rotational orientation of the hydroxyl group, cytosine analogues have possible conformers in which the H~4~ atom is oriented toward the C~5~ or N~5~ atom, as shown in Fig. [2](#Fig2){ref-type="fig"}. The conformers were optimized by DFT at the current level and no imaginary frequencies were found. The present calculations show all the analogues have lower energy compared with their conformers (Table [1](#Tab1){ref-type="table"}), indicating that these analogues are more stable than theirs conformers.

The vdW surfaces and ESP extrema of the cytosine analogues {#Sec5}
----------------------------------------------------------

The electrostatic potential (ESP)-mapped van der Waals (vdW) surface has been used extensively for interpreting and predicting reactivity and intermolecular interactions of a wide variety of chemical systems^[@CR20],[@CR21]^. The negative charge of the Ca~0~ molecule is concentrated on the O~2~ atom with an ESP extrema value of −71.53 kcal/mol. The positive electrostatic potential is centered around H~4a~ of the amino group with an ESP extrema value of +38.77 kcal/mol (Fig. [3](#Fig3){ref-type="fig"}. Ca~0~). When the amino group at the C~4~ position of the cytosine was replaced by a hydroxy moiety, the ESP extrema value was increased to +47.47 kcal/mol, which is favorable for nucleophilic attack and release of the proton (Fig. [3](#Fig3){ref-type="fig"}. Ca~1~). Furthermore, when the C~5~ carbon atom was replaced by a nitrogen atom and the carbon atom at C~6~ was a methine, methylene, or carbonyl, the ESP extrema value of H~4~ was increased to +52.06 kcal/mol, +52.86 kcal/mol and +69.67 kcal/mol, respectively, which are all favorable for nucleophilic attack and release of the proton (Fig. [3](#Fig3){ref-type="fig"}. Ca~2~-Ca~4~). When the carbonyl oxygen at C~2~ was replaced by a hydrogen, the negative charges and the ESP extrema shift to the N~3~ atom, which is favorable for electrophilic attack and accepting a proton at the N~3~ site (Fig. [3](#Fig3){ref-type="fig"}. Ca~5~-Ca~8~). Furthermore, when the C~2~ carbon atom was replaced by a boron atom, the ESP extrema values of N~3~ are strengthen to −45.65 kcal/mol, −39.81 kcal/mol, −51.50 kcal/mol and −33.59 kcal/mol, which are all more favorable for electrophilic attack and accepting a proton at the N~3~ site (Fig. [3](#Fig3){ref-type="fig"}. Ca~9~-Ca~12~). So, these atomic substitutions, especially the substitutions of a hydroxy moiety for the amino group, a hydrogen for the carbonyl oxygen at C~2~ and a boron atom for the C~2~ carbon atom, are all favorable to induce DPT reactions with guanine.Figure 3ESP-mapped vdW surfaces and ESP extrema of cytosine analogues generated by Multiwfn 3.6 and VMD 1.9.1 (<http://www.ks.uiuc.edu/Research/vmd/>). The red color indicates positive ESP regions on the vdW surface; the blue regions correspond to negative ESP regions of the vdW surface. The significant local maxima and minima of ESP on the vdW surfaces are represented as orange and cyan spheres, and labelled by values given in kcal/mol.

Energy barriers and equilibrium constants of the DPT reaction {#Sec6}
-------------------------------------------------------------

The structures of all the Ca~n~G complexes were optimized by DFT at the current level. The vibration analysis showed no imaginary frequencies in these complexes, indicating that the structures are stable. Intermolecular DPT, which generates hydrogen-bonded pairs of tautomers, was computationally predicted for all the complexes. DPT can take place through a concerted DPT (Ca~1-3~G and Ca~5-12~G) or via a stepwise mechanism involving two distinct SPT steps (Ca~0~G and Ca~4~G). The vibration analysis of each proton transfer reaction showed one imaginary frequency, and the vibration mode of the imaginary frequency corresponds to the reactants and products assigned to the transition state. An intrinsic reaction coordinate was prepared to confirm the existence of the transition states.

The Gibbs free energies of the fully optimized Ca~0~G complexes were defined as having zero energy. The relative Gibbs free energies of the optimized reactants, the first and second single-proton transition states, the double-proton transferred states, the intermediate product of single proton transfer, and the product of double-proton transfer were calculated.

The Ca~0~G complex, the model of the Watson--Crick guanine--cytosine (GC) base pair, undergoes a stepwise DPT process. The dissociation energy (DE) of the Ca~0~G base pair at the current level of theory is 21.74 kcal/mol, which is in good agreement with the reported experimental value of 21.0 kcal/mol^[@CR22]^. However, the previously calculated DE values were 25.4 kcal/mol^[@CR23]^, 23.8 kcal/mol^[@CR24]^, and 24.4 kcal/mol^[@CR25]^. This result indicates that the present calculations are a better mimic of physiological conditions. The energy of the DPT product of Ca~0~G lies 9.15 kcal/mol above the reactant species, which is similar to previously calculated values of 9.8 kcal/mol^[@CR14]^ and disfavors the DPT reaction. The Ca~0~G--Ca~0~\*G\* equilibrium (with a calculated equilibrium constant of 1.99 × 10^−6^, which is in good agreement with previously estimated values of 2.0 × 10^−6^)^[@CR12]^ largely favors the reactants, so this double-proton-transfer reaction will rarely occur. In fact, we also calculated the dissociation energy of Ca~0~G and the equilibrium constant of the DPT reaction using the DFT B3LYP/6-311 + +G(d,p) method. And the calculated values were 13.25 kcal/mol and 4.53 × 10^−8^ respectively, which were quite different from the experimental values (21.0 kcal/mol) of the dissociation energy and the previously estimated values (2.0 × 10^−6^) of the equilibrium constant. So, these results calculated by the DFT M06-2X/def2svp method are in good agreement with previous findings related to the DPT reaction of GC base pairs, indicated the calculations are more reasonable.

The forward DPT free energy barriers of the 12 analogue complexes (Δ*G*^≠^~*f*~) range from −2.02 to 1.07 kcal/mol, which are significantly lower than the value of Ca~0~G (8.80 kcal/mol). The barriers increased in the order Ca~11~G \< Ca~7~G \< Ca~3~G \< Ca~12~G \< Ca~5~G \< Ca~4~G \< Ca~8~G \< Ca~6~G \< Ca~10~G \< 0 \< Ca~9~G \< Ca~2~G \< Ca~1~G\<\<Ca~0~G. Lower energy barriers favor the DPT reaction, especially for complexes with negative barriers. The forward rate constants (calculated from Δ*G*^≠^~*f*~) of the DPT process range from 1.14 × 10^12^ to 1.72 × 10^14^ s^−1^, which are significantly higher than the value of Ca~0~G (4.05 × 10^6^ s^−1^). The reverse DPT free energy barriers of the 12 analogue complexes range from 0.97 to 8.06 kcal/mol, which are higher than the value of Ca~0~G (0.71 kcal/mol). In addition, the reverse rate constants (calculated from Δ*G*^≠^~*r*~) of the DPT process range from 1.34 × 10^7^ to 1.34 × 10^12^ s^−1^, which are lower than the value of Ca~0~G (2.04 × 10^12^). The DPT equilibrium (k~eq~ = k~f~/k~r~) ranges from 1.60 × 10^0^ to 1.28 × 10^7^, which is significantly higher than the value of Ca~0~G (1.99 × 10^−6^). The DPT equilibrium constants increased in the order Ca~0~G\<\<0 \< Ca~2~G \< Ca~1~G \< Ca~4~G \< Ca~6~G \< Ca~5~G \< Ca~3~G \< Ca~9~G \< Ca~8~G \< Ca~10~G \< Ca~12~G \< Ca~7~G \< Ca~11~G. For all 12 analogue complexes with positive equilibrium constants, the reaction may occur in the DPT direction. In particular, the DPT equilibrium constants of Ca~7~G and Ca~11~G (1.49 × 10^5^ and 1.28 × 10^7^) are over 1.0 × 10^5^, which indicates that these two proton transfer reactions will occur adequately. All the values are shown in Table [1](#Tab1){ref-type="table"}.

Geometries of the DPT reactions of the analogue complexes {#Sec7}
---------------------------------------------------------

In the Ca~0~G complex, the intermonomer N~3~ − N′~1~ distance (as labeled in Fig. [1](#Fig1){ref-type="fig"}) is 2.92 Å. This distance is consistent with the calculated distance (2.95 Å) and the experimental distance (3.09 Å)^[@CR12]^, as is the distances of N~4~ − O′~6~. As the reaction progresses from the reactant to transition state 1 (TS~1~), the intermediate product of the single proton transfer (SPT) and TS~2~, the guanine and cytosine monomers approach each other. The N~3~ − N′~1~ distances in TS~1~, SPT and TS~2~ are 2.65 Å, 2.72 Å and 2.78 Å, respectively. As the reaction proceeds, the N′~1~ − H′~1~ bond stretches faster than the N~4~ − H~4a~ bond and completes the proton transfer first. In the products, the N~3~ − N′~1~ bridge elongates again to 2.86 Å, which is close to their original separation in the reactants.

In the Ca~1-12~G analogue complexes, the intermonomer N~3~ − N′~1~ distance and O~4~-O′~6~ distance range from 2.70 to 2.81 Å and 2.46 to 2.61 Å, respectively, which are all shorter than those in the Ca~0~G complex (2.92 Å and 2.85 Å), indicating that the forward DPT reaction of those complexes is more favorable than that of the Ca~0~G complex. As the reaction progresses from the reactant to TS, most of the N~3~ − N′~1~ distances and O~4~-O′~6~ distances are shorter than the distance in the Ca~0~G complex and favor the DPT reaction. In the products, the N~3~ − N′~1~ and O~4~-O′~6~ bridges elongate again, and more than half of these bridges are longer than the distance in the Ca~0~G DPT product, which indicates that the DPT reverse reaction in those complexes is more difficult than it is in the Ca~0~G DPT product (Table [2](#Tab2){ref-type="table"}).Table 2Distances (in Å) between the electronegative atoms involved in the two H-bonds in Ca~n~G during the DPT process.Ca~n~GElectronegative AtomsReactantTS (TS~1~/SPT/TS~2~)ProductCa~0~GN~3~-N′~1~2.922.65/2.72/2.782.86N~4~-O′~6~2.852.62/2.59/2.482.58Ca~1~GN~3~-N′~1~2.812.612.84O~4~-O′~6~2.542.432.59Ca~2~GN~3~-N′~1~2.792.602.82O~4~-O′~6~2.522.462.61Ca~3~GN~3~-N′~1~2.772.612.87O~4~-O′~6~2.502.432.60Ca~4~GN~3~-N′~1~2.762.71/2.69/2.602.81O~4~-O′~6~2.462.40/2.45/2.502.65Ca~5~GN~3~-N′~1~2.742.572.92O~4~-O′~6~2.612.522.52Ca~6~GN~3~-N′~1~2.802.562.94O~4~-O′~6~2.602.492.53Ca~7~GN~3~-N′~1~2.722.572.93O~4~-O′~6~2.582.492.51Ca~8~GN~3~-N′~1~2.782.562.89O~4~-O′~6~2.522.412.59Ca~9~GN~3~-N′~1~2.782.592.94O~4~-O′~6~2.582.452.54Ca~10~GN~3~-N′~1~2.782.582.91O~4~-O′~6~2.562.442.57Ca~11~GN~3~-N′~1~2.702.592.95O~4~-O′~6~2.532.472.55Ca~12~GN~3~-N′~1~2.772.652.91O~4~-O′~6~2.492.402.62

Effects of cytosine atom substitution on the DPT reaction {#Sec8}
---------------------------------------------------------

When the amino group at the C~4~ position of the cytosine was replaced by a hydroxy moiety and the other conditions were held constant, the energy barrier of the DPT reaction decreased by 7.73 kcal/mol. This result indicates that the hydroxy group significantly favored the DPT reaction. Since the electronegativity of oxygen atom is lower than that of nitrogen atom, the attraction of oxygen to proton is less than that of nitrogen, which is the possible reason that favors the transfer of H~4~ from O~4~ to N′~6~.

When the C~2~ carbonyl oxygen atoms of the cytosine analogues were replaced by hydrogen, the energy barriers of the DPT decreased by an average of 1.01 kcal/mol. Replacing C~2~ with a boron atom further decreased the energy barriers of the DPT by an average of 0.06 kcal/mol, further favoring the DPT reaction. Removing the carbonyl group from C~2~ and replacing C~2~ with a less electronegative boron atom can transfer the negative charge to N~3~, which is conducive to the ability of N~3~ to acquire proton and favors the transfer of H'~1~ from N′~1~ to N~3~. A large number of BN-containing heterocycles have been developed by means of the substitution of a carbon-carbon double bond with an isoelectronic boron-nitrogen bond (BN-substitution)^[@CR26],[@CR27]^. BN-substituted nucleobases are very attractive BN-substituted heterocyclic compounds. Bielawski *et al*.^[@CR28]^ reported the synthesis of B(6)-phenyl-BN-uracil, which was characterized by mass spectrometry, IR spectroscopy and elemental analysis. It is interesting to note that Hiroshi *et al*.^[@CR26]^ also reported the synthesis of B(6)-substituted 5-aza-6-borauracils (BN-substituted uracils) and -thymines (BN-substituted thymines). The structures of these BN-substituted nucleobases are similar to the cytosine analogues in present study. Therefore, these boron-containing cytosine analogues can be synthesized as realistic candidates.

When the C~5~ carbon atom of the cytosine analogues was replaced by a nitrogen atom, the energy barriers of the DPT decreased by an average of 0.03 kcal/mol. When the double bonds between N~5~ and C~6~ of these analogues were hydrogenated to single bonds, the energy barriers of the DPT further decreased by an average of 1.62 kcal/mol. However, when the hydrogen on C~6~ was exchanged for a carbonyl oxygen, the energy barriers of the DPT increased again. Therefore, the second result seemed to indicate that the DPT reaction was favored. Replacing C~5~ with highly electronegative nitrogen atom may transfer the negative charge from O~4~ to N~5~, which will reduce the attraction of O~4~ to proton and favors the transfer of H~4~ from O~4~ to N′~6~.

Conclusions {#Sec9}
===========

The abilities of 12 modified cytosine--guanine complexes to undergo DPT were predicted and compared using theoretical calculations. The DE value, DPT barriers, hydrogen bond lengths and equilibrium constants of the Ca~0~G complex are similar to previously calculated values and experimental values, indicating that the present calculation method is reasonable. Because of the intramolecular SPT process are not facile, the 12 modeled cytosine analogues can be used as candidate molecules to induce DPT reactions with guanine. Eight modified complexes (Ca~3~G, Ca~5~G and Ca~7-12~G) were significantly more prone to undergo DPT reactions than Ca~0~G. In particular, Ca~7~G and Ca~11~G may undergo sufficient DPT reactions under physiological conditions according to the present calculations. Here, we present the part of the study on theoretical calculations and predictions of candidate molecules. Next these analogues are expected to be synthesized for incorporation into single-stranded RNAs and to be validated in *in-vivo* model. By binding such modified RNA to DNA or RNA, DPT reactions of guanine may be induced at a fixed point. If cytosine analogues that can spontaneously induce DPT reactions of guanine under physiological conditions are identified experimentally, targeted pathogenic mutations can be used to restore original functions at the level of DNA or RNA.

Theoretical Methods {#Sec10}
===================

DFT with the M06-2X functional and the def2svp basis set has been the primary research method used to investigate proton transfer reactions between cytosine analogues and guanine. The M06-2X functional is a high nonlocality functional with twice the nonlocal exchange (2X), and it is considered suitable for the study of main-group thermochemistry, thermochemical kinetics, noncovalent interactions, and electronic excitation energies to the valence and Rydberg states. The M06-2X functional also gives the best performance for hydrogen-transfer barrier calculations and the lowest values of balanced mean unsigned error (BMUE), which means it gives the best overall performance for barrier calculations^[@CR29]--[@CR31]^. The D3 dispersion correction was used in this study to improve the accuracy of the energy and structure calculations^[@CR32]--[@CR35]^. To simulate the DNA surroundings in a biological environment, all the calculations were carried out with water solvation at *T* = 310.15 K (37 °C) and *p* = 1 atm. The sophisticated polarizable continuum model^[@CR36]^ has been used to investigate solute--solvent interactions in water using the scrf = (solvent = water, pcm) keyword. The Gaussian 09 program package was used throughout this study^[@CR37]^. The structures of all the monomers and complexes were optimized by DFT at the current level. The vibration analyses were performed to determine whether the molecular structures were stable. The H~4~ and H'~1~ hydrogen atoms of the optimized complexes were placed in the middle points between two electronegative atoms (N~3~-N′~1~, N~4~-O′~6~ or O~4~-O′~6~) and then the structures of the transition states were optimized using the opt = (calcall, ts, noeigen) keyword. The vibration analyses were performed to determine the existence of the transition states. The forward energy barrier, Δ*G*^*≠*^~*f*~, is the difference between the Gibbs free energy of the transferred states and the reactants. The reverse energy barrier, Δ*G*^≠^~*r*~, is the difference between the Gibbs free energy of the transferred states and the products.

The GaussView 5.0^[@CR38]^ was used to create molecular structures of the modeled cytosine analogues. The molecular structure of classic cytosine, at first, was created and optimized by the program package. And then the 12 modeled cytosine analogues were created by replacing the atoms in the classic cytosine and optimized by the program package.

To investigate the ability of N~3~ to accept protons and O~4~ to donate protons, the ESP-mapped vdW surfaces and ESP extrema were rendered using the VMD 1.9.1 program^[@CR39]^ based on the outputs from the Multiwfn 3.6 program^[@CR40],[@CR41]^.

To investigate the equilibrium of the DPT reactions and chemical reactions of the analogue complexes, the forward and reverse rate constants (k~f~ and k~r~) were given by the following equation^[@CR42]^$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$k=2.1\times {10}^{10}T{{\rm{e}}}^{-1000\varDelta {G}^{\ne }/(1.9859\times T)}$$\end{document}$$Where k is the forward or reverse rate constants, *T*=310.15 K (37 °C), Δ*G*^≠^ is the forward or reverse energy barriers. The equilibrium constants (K~eq~ values) of the DPT reactions are the quotients of k~f~ and k~r~ (K~eq~= k~f~/k~r~).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We are grateful to Jiahui Xia, Huanxiao An, Lingqian Wu, Desheng Liang, Jigeng Bai, Yongxin Yang, Yu Wang, Wei Wei and the March of dimes birth defects foundation of Shanxi for their support in this study.

Conceptualization: Jinjie Xue. Methodology: Jinjie Xue, Xingbao Wang and Yafeng Xiao. Investigation: Jinjie Xue and Xingping Guo. Writing Original Draft: Jinjie Xue. All authors reviewed and approved the manuscript.

The authors declare no competing interests.
