The complex structure of DNA/UDG was downloaded from the Protein Data Bank (PDB ID: 1EMH (Parikh et al., 2000); Figure 1), and we visualized it by using Chimera (Pettersen et al., 2004). In this original structure, the base at the target location (B5, flips out to UDG) is a pseudosubstrate (P2U), which is to replace uracil so that it binds to UDG stably. Since UDG targets uracil in a real DNA base rather than P2U, we mutated the P2U base to uracil (U) using Chimera. We deleted all the water molecules that are involved in the original structures, as DelPhi (Li et al., 2012), DelPhiForce (Li et al., 2017a), and DelPhiPKa (Wang et al., 2015b) implement an implicit solvent model (Poisson–Boltzmann) in the calculations, which have been proved to be successful in previous studies (Jia et al., 2017; Xie et al., 2020a; Xie et al., 2020b; Guo et al., 2021; Lopez-Hernandez et al., 2021).

As we mentioned previously, in the repairing process of UDG applied on uracil-induced DNA, the incorrect base uracil (U) of DNA (location: B5; Figure 1B) is the target base for UDG to hydrolyze. In order to study the presence of this uracil in a DNA chain and compare it with the original base (before the damage) at this location, which is cytosine (C), we generated a new DNA/UDG structure using cytosine to replace uracil at the B5 position, with the help of Chimera. To better discuss the two structures in the following, we named the DNA with uracil as DNA_RU and the DNA with an original cytosine as DNA_C.

The folding energy pH dependency of the entire UDG enzyme is shown in Figure 2, and the detailed discussion is included in Results and Discussions section. In particular, we studied the binding pocket of UDG, which is referred to as the essential binding area (Parikh et al., 1998). The surface of the binding pocket is colored in magenta as shown in Figure 3, and the residues involved in the UDG pocket are Q144, D145, P146, Y147, H148, F158, S169, S247, H268, P269, S270, P271, L272, and S273.

In order to study the electrostatic features, DelPhi (Li et al., 2012; Li et al., 2013) was utilized to calculate the electrostatic potential of DNA and UDG. In the framework of continuum electrostatics, DelPhi calculates the electrostatic potential ϕ (in systems comprised of biological macromolecules and water in the presence of mobile ions) by solving the Poisson–Boltzmann equation (PBE): ∇ · [ ϵ ( r ) ∇ ϕ ( r ) ] = − 4 π ρ ( r ) + ϵ ( r ) κ 2 ( r ) sinh ( ϕ ( r ) / k B T ) , (1) where   ϕ ( r ) is the electrostatic potential, ϵ ( r ) is the dielectric distribution, ρ ( r ) is the charge density based on the atomic structures, κ is the Debye–Huckel parameter, k B is the Boltzmann constant, and T is the temperature. Due to the irregular shape of macromolecules, DelPhi uses a finite difference (FD) method to solve the PBE.

Before the DelPhi calculations, the PQR files of DNA and UDG were generated by the PDB2PQR (Dolinsky et al., 2004) tool. We used the AMBER force field for PDB2PQR calculation and removed the water molecules in the process, using the PDB2PQR web server (https://server.poissonboltzmann.org/pdb2pqr). The PDB2PQR built the new hydrogen atoms in proper distances with existing atoms to avoid clashes which also optimized the hydrogen bonding network.

For the calculation parameters in DelPhi, the grid resolution was set to be 2.0 grids/Å. The dielectric constants were set as 2.0 for protein and 80.0 for the water environment. The probe radius for generating the molecular surface was 1.4 Å. The salt concentration was 0.15 M. The boundary condition for PBE was set as a dipolar boundary condition. After the calculation, the values of electrostatic potential on the surface were visualized with Chimera (Figure 4). In order to visualize the electric field lines between DNA and UDG, the separation distance of the DNA from UDG was set to 20Å with respect the direction of their mass centers connection line. Visual molecular dynamics (VMD) (Humphrey et al., 1996) was implemented based on the electrostatic potential map from DelPhi calculations, and the color scale range was set from −3.0 to 3.0 kT/e. For a better representation of electric field lines, we chose the line size to be 4 and delta value to be 0.25, and set the gradient magnitude value to be 3.64 (which shows the lines within the volumetric). Besides, those selected field lines for display have the minimum line length of 4.12 and maximum line length of 35.31.

The net charges of proteins at the unfolded state were calculated using the following equation: Q u ( pH ) = ∑ i = 1 N 10 − 2.3 y ( i ) ( pH − pKa ( i ) ) 1 + 10 − 2.3 y ( i ) ( pH − pKa ( i ) )   , (2) where the summation is of all the titratable groups, y(i) value is −1 for acidic groups, and +1 for basic groups, respectively. As for the folding free energy, the next equation was implemented: ΔN ( p H folding ) = 2.3 RT ∫ p H i p H f ( Q f ( pH ) − Q u ( pH ) d ( pH ) ) , (3) where Q f ( p H ) and Q u ( p H ) stand for the net charge of folded and unfolded state, respectively. R is the universal gas constant taken as 1.9872 × 10 − 3 k c a l M o l ∗ K and T is the temperature with the value of 300 K.

DelPhiPKa (Wang et al., 2015a; Wang et al., 2015b) was used to calculate the pH dependence of folding energy for UDG, given the pH ranging from 0 to 14 with the pH interval of 0.5. During the calculations, we used the AMBER force field (Wang et al., 2004). Water molecules and HETATM were removed because the implicit solvent model is used in DelPhiPKa. Variance of Gaussian Distribution was set to be 0.7, salt concentration was 0.15 M, reference dielectric was 8.0, and external dielectric was 80.0. Note that this method calculated relative folding energies. The folding energy at pH = 0 was set as reference (0 kcal/mol). Lower folding energy at certain pH value indicates higher stability at that pH.

To compare the strengths and directions of electrostatic forces between DNA and UDG, DelPhiForce (Li et al., 2016b) was implemented to perform the force calculations. During DelPhiForce calculations, grid solution was set to be 2.0 grids/Å and salt concentration was 0.15 M. In order to study the binding process of DNA to UDG, we separated the DNA from UDG in the direction of their mass centers connection line with the distances ranging from 20Å to 40Å with the step size of 4Å.

The electrostatic binding forces calculated by DelphiForce were visualized with VMD and represented by arrows. The arrows in Figure 5 represent the directions of forces between DNA and UDG, as they were normalized to be of the same size. In order to study the base B5 (RU or C) and the UDG pocket in particular, apart from the calculations between DNA and the whole UDG, we also did the same calculations of B5 and UDG, as well as B5 and the UDG pocket (Figures 5C,D).

As discussed in the Introduction section, UDG is able to detect the damaged DNA base pair AU in a double-stranded DNA by identifying the unusual kink of 45° (Satange et al., 2018). The uracil base then flips out to the UDG binding pocket so that UDG hydrolyzes the uracil from DNA. Figure 1 shows the binding state of uracil to the UDG pocket.

To better understand the stability of UDG in different environments, especially at different pH values, we calculated the pH dependence of UDG folding energy by using DelPhiPKa.

The calculation was performed at different pH values ranging from 0 to 14 with an interval of 0.5 (Figure 2). From the trends in Figure 2, we observed that the folding energy decreases from 0 to 5, then it becomes relatively more stable from 5 to 10, and increases from 10 to 14. Figure 2 indicates that UDG is stable at pH ranging from 5 to 10. Therefore, the average value of optimal pH is 7.5, which matches the storage conditions of UDG in the laboratory (Wu et al., 2015). Note that the folding energies in Figure 2 are relative values because we set the reference energy to be 0 kcal/mol when pH is equal to 0. We did not calculate the absolute folding energies (energy difference between folded and unfolded states) since we focused on the pH dependency of the folding energies (energy difference between folded energies at a certain pH and pH 0).

DNA always binds to UDG at the pocket side (colored as magenta in Figure 3; sequence: Q144, D145, P146, Y147, H148, F158, S169, S247, H268, P269, S270, P271, L272, and S273) rather than the other side. This is crucial for the binding process since only the pocket area is able to “cut” uracil in DNA instead of other regions on UDG. But the factors to guide DNA binding with the UDG binding pocket efficiently are not fully understood; here, we illustrated the electrostatic potential on the surface of UDG to demonstrate the binding mechanism.

To study the electrostatic features, DelPhi was utilized to calculate the electrostatic potential on the surfaces of DNA and UDG. The electrostatic potential distribution on DNA is shown in Figure 4B and movie 1 (see the Supplementary Material), which were rendered by Chimera with a color scale from −3.0 to 3.0 kT/e. The charge distribution on UDG is shown in Figures 4DF and movie 2 (see the Supplementary Material), which were rendered by Chimera with a color scale from −3.0 to 3.0 kT/e as well, negatively and positively charged areas are colored in red and blue, respectively.

By comparing the electrostatic potential on the surfaces of the pocket side and the non-pocket side of UDG, it is clear that UDG has a polar charge distribution. The pocket side has an overall positively charged surface, while the non-pocket side has dominantly negative surface. This charge distribution helps to increase the binding efficiency and decrease the binding direction errors, which ensures the DNA binds directly to the pocket, followed by other processes. Similarly, such polar distributions are commonly found in many other protein–protein interactions, such as molecular motors binding with microtubules (Li et al., 2016b), viral capsid binding with each other (Xian et al., 2019), and enzyme binding with inhibitors (Li et al., 2017b).

Next, by analyzing the electrostatic potential on the surfaces of DNA and UDG, it is obvious that DNA has a negatively charged surface (Figure 4B) while UDG (pocket side; Figure 4D) has a positively charged surface. This fact indicates the attractive forces between DNA and UDG (pocket side). In order to investigate more about the attractive interactions, field lines were generated between DNA and UDG, which is discussed in the Electric field lines section.

After looking at the overall potential distribution, we colored charged amino acids in Figures 6C,F to visualize the charge distribution.

In Figure 6, black squares indicate the area of binding interface of UDG when DNA binds to it. By looking at Figure 6B, the binding interface is overall blue, which is positively charged. While by looking at Figure 6D, there is only one negatively charged residue (ASP in pink) inside the pocket, while the surrounding area of the binding interface has several positively charged residues (ARG in blue and LYS in green). This fact indicates that the binding pocket itself does not provide the attractive forces to DNA, since DNA is overall negatively charged. However, the whole binding interface does provide the positively charged environment for DNA to be attracted and bound to UDG.

In order to investigate more about the pocket area in particular, we calculated the forces between DNA and UDG with different partial complex structures that are DNA/UDG, DNA (B5 base only)/UDG, and DNA (B5 base only)/UDG (pocket only). Those results are discussed in the Electronic forces section.

Electric field lines between UDG and DNA were calculated. To better visualize the field lines between the binding interface, DNA was separated from UDG by 20 Å (Figure 7).

The field lines distribution confirmed that UDG and DNA have attractive forces between each other. In the analysis of field lines, the density of the distribution indicates the strength of the electrostatic binding forces, which means that the denser distribution has the stronger interactions. From Figure 7B, we noticed that the flipping out base uracil has a very dense field lines connected to the UDG binding pocket.

Electrostatic forces of DNA and UDG were calculated by DelPhiForce. The calculated results were visualized with arrows (Figure 5) and line graphs (Figure 8). Arrows in Figure 5 represent the net forces between DNA and UDG by shifting the DNA away from UDG by variable distances ranging from 20 Å to 40 Å with the step size of 4 Å. The direction of arrows represents the force directions. To better visualize the direction of the net forces, the magnitudes of the net forces were normalized to be of the same size, which means that the size of the force does not represent the force strength. Force strengths are discussed in the Force strengths section and visualized in Figure 8.