Theoretical analyses on a flipping mechanism of UV-induced DNA damage.

As for UV-induced DNA damage, which may induce skin cancer in animals and growth inhibition in plants, there are two types of photoproducts, namely cis-sin cyclobutane pyrimidine dimers (CPD) and pyrimidine-pyrimidone (6-4) photoproducts. When they are to be repaired, base-flipping occurs, and they bind to enzymes. However, this process remains relatively unknown at a molecular level. We analyze conformation and interaction energy changes upon base-flipping using classical molecular dynamics (CMD) simulations and ab initio electronic structure calculations. CMD simulations starting with a CPD in the flipped-in and flipped-out states showed that both states were unchanged for 500 ns, indicating the flipped-in and flipped-out processes do not occur spontaneously (without any help of the enzyme) after photo-damage. To deeply understand the reasons, we investigated interaction energy changes among bases upon structure changes during the flipped-in and flipped-out processes using Parallel Cascade Selection-MD (PaCS-MD) simulations at 400 K, followed by a fragment molecular orbital (FMO) method. The total inter-fragment interaction energy (IFIE) between CPD and other bases at the flipped-in state is estimated to be -60.08 kcal/mol. In particular, four bases strongly interact with CPD with interaction energies being -10.96, -13.70, -21.52, and -14.46 kcal/mol each. On the other hand, the total IFIE at the obtained flipped-out state increased to -10.40 kcal/mol by partly losing hydrogen bonds and π-π stacking interactions, respectively. These results clearly indicate that the base-flipping process of DNA lesions occurs with the help of external forces like interactions with appropriate enzymes such as photolyases.

DNA is damaged by irradiation of solar light in the Ultra-Violet (UV) region (290–320 nm) [1,2]. The UV-induced DNA lesions may cause skin cancer in animals and growth delay in plants [3-5]. There are two major types of photoproducts in DNA. One is a cis–syn cyclobutane pyrimidine dimer (CPD) and the other is a pyrimidine–pyrimidone (6-4) photoproduct ((6-4)PP). These photoproducts originate from covalent bond formations between adjacent pyrimidine bases in DNA via photochemical reactions, and mostly exert their detrimental cytotoxic effects by blocking DNA replication and transcription [6,7].

Photolyases and cryptochromes are blue-light photoreceptors and functionally distinct flavoproteins found in some organisms. The photolyases contain two noncovalently bound flavin adenine dinucleotides (FADH−). FADH− serves as a cofactor that carries out the repair function upon excitation by either direct photon absorption or by resonance energy transfer from a second chromophore (methenyltetrahydrofolate or deazaflavin). The second chromophore is an antenna pigment that harvests sunlight and enhances repair efficiency [8-12]. The DNA repair mechanism is thought to be as follows. Firstly, photolyase recognizes a DNA lesion moiety and binds there. Then, an excited state electron transfers to the protein-bound DNA lesion moiety from FADH−* (excited state FADH−) which is produced by an irradiation of visible light (350–450 nm). After the acceptance of the electron by DNA lesion moiety, the covalent bonds in the dimer are cleaved, resulting in a repair of the DNA lesion. Although cryptochromes are similar to photolyases in sequence and structure, cryptochromes cannot perform detectable DNA repair activity. Alternatively, they regulate growth and development of plants and control the circadian clock in animals with near-UV/blue-light [13-16]. Surprisingly, a newly discovered cryptochrome, cryptochrome Drosophila Arabidopsis Synechocystis Human (CRY-DASH) can repair single strand, but not double strand, DNA lesions in vitro [17]. However, the repair mechanism has not been revealed yet.

Recently, the crystal structure of a photolyase (PHRs)-DNA lesion complex has been solved by X-ray measurement [18-20]. In this structure, the DNA lesion deviates from the double strand DNA and interacts with a binding pocket in the PHRs (see Fig. 1 (a)). The structural changes from the regular double strand DNA to that with the DNA lesion exposed is referred to as a forward “base-flipping” process (see Fig. 1 (b)). Since the double strand DNA structure is strongly stabilized both by hydrogen bonding and base stacking interactions among bases, it is expected that the energy required to induce the base-flipping process is quite high. A recent computational study on the base-flipping process by Neil and Wiest reported free energy differences between flipped-in and -out states by using molecular dynamics (MD) simulation [21,22]. The estimated free energy differences for several different sequences ranged from 5.25 to 7.5 kcal/mol, meaning that the free energy required for the base-flipping process was more than 5 kcal/mol. Moreover, the NMR measurements on an imino proton-exchange process have shown that the dissociation constant of the base-flipping of a guanine base in a 5′–GCGC–3′ tetramer is estimated to be 3.3×10−7, which corresponds to a free energy difference of ~9 kcal/mol [23]. These results strongly suggest that the base-flipping process does not occur spontaneously and that the help of PHRs is required. However, the simulation time of the previous MD simulations was too short to track the spontaneous base-flipping process at MD simulation, so the cause of flipping not occurring has not been clarified yet.

Figure 1

(a) Binding structure of damaged DNA that contains a DNA lesion with a photolyase. (b) Conformational change from a flipped-in to -out state.

In this paper, we precisely analyze the base-flipping process of a DNA lesion without PHR in aqueous solution based on both MD simulations and ab initio quantum chemical calculations. In particular, the conformational changes were induced with our recently developed conformational search method, i.e. Parallel Cascade Selection MD (PaCS-MD) [24]. Using the fragment molecular orbital (FMO) method, the inter-fragment interaction energies (IFIEs), which show the interaction energies among a CPD and other bases, between flipped-in and -out structures of a duplex DNA with the CPD were estimated, and important bases for the base-flipping process were detected. Then, the total energy change due to the base-flipping process was calculated by the FMO method to elucidate the chance of the base-flipping of CPD. Finally, we give a brief remark on a necessary element for the base-flipping process found in our analyses.

1. Theoretical Methods

1.1. Classical MD simulations

All MD simulations were performed using the Amber 14 program suite [25]. For the force field parameters of nucleic acid bases, such as adenine (A), thymine (T), guanine (G), and cytosine (C), amber force filed 14 SB is used [26]. We performed MD simulations for a model lesion DNA which includes CPD (PDBID: 1SNH). Note here that in this structure, the phosphate group PO4− of CPD moiety was replaced by CH2O2, and two bases, which form hydrogen bonds to CPD, were replaced by G, but they are normally A in natural DNA. 1SNH is in the flipped-in state and we need to create the structure of the flipped-out state of the same base sequence. Therefore, the DNA structure for the flipped-out state was constructed from 1SNH with winmostar [28]. For the force field parameters of CPD, we calculated the electronic structure of the CPD at B3LYP/6-31G(d) level using Gaussian 09 to obtain atomic charges on the CPD and converted the data using the Antechamber module [27] for the force field parameter. The systems were neutralized using Na+ counter ions and solvated using the TIP3P water model as provided in tLeap. The solvent box is extended by 10 Å from the center of mass of DNA. The final systems consisted of 11,765 and 11,147 atoms for the flipped-in and -out state simulations, respectively (including 19 bases and 3,705 and 3,499 water molecules and 18 Na+ ions).

At the first step, rough structure minimizations were performed for 10,000 steps by imposing restrictions on heavy atoms in DNA (a force constant of 10 kcal/mol·Å2) to relax positions of water and hydrogen atoms. Then, minimizations without any restriction were performed for 10,000 steps. In the initial MD simulation, the systems were heated from 0 K – to 300 K for 100 ps with the constant-pressure, isothermal ensemble (NPT, P=1 bar and T=300 K) by imposing restrictions on the heavy atoms again (the force constant of 10 kcal/mol·Å2). Next, we performed the equilibration run with the NPT ensemble for 10 ns with restrictions on four edges of DNA (the force constant of 10 kcal/mol·Å2). Here, we combined two thymine bases of residue 5 and 6 as CPD. Therefore, a total number of the residues is 19, where the four edges of the sequences are residue 1 (cytosine) and 9 (cytosine) on A chain and 10 (guanine) and 19 (guanine) on B chain in 1SNH (Supplementary Fig. S1). After the equilibration, the production run was performed for 500 ns with the NPT ensemble.

1.2. PaCS-MD for reproducing the flipped-out of CPD

To reproduce the flipping process of the CPD moiety, the PaCS-MD [24] was applied to the DNA lesion. The PaCS-MD is an enhanced conformational sampling method for generating structural transition pathways between end-point structures under a condition that a set of reactant and product is known a priori. In the PaCS-MD, cycles of restarting short-time MD simulations from different initial structures are repeated to generate the structural transition from the reactant to the product. The PaCS-MD starts from the reactant with distinct initial velocities and continues until some of initial structures reach sufficiently close to the product. To generate the initial structures for the first cycle, preliminary MD simulation is performed starting from the reactant with different initial velocities. Then, generated snapshots are rank-ordered for approaching the product. Herein, appropriate reaction coordinates (RCs) are utilized as a ranking measure. Top N snapshots are selected as the initial structures for the first cycle. Then, the short-time MD simulations starting from the initial structures are independently performed. After that, the generated snapshots are rank-ordered again and highly ranked N snapshots are newly selected as the initial structures for the next cycle. This cycle is subsequently repeated until the highly ranked snapshots reach sufficiently close to the product. As the cycles go, the initials selected at each cycle would be gradually close to the product. For the details of PaCS-MD, see the original reference [24] and see also Supplementary Figure S2 for the flowchart of the PaCS-MD.

In this study, a set of the reactant and the product were preliminary given, i.e. a model structure (flipped-in state) as the reactant and that taken from X-ray structure (PDBID: 1TEZ [18]) as the reference structure (flipped-out state). As an RC for representing the base-flipping process of the CPD moiety, an inner product of characteristic angles between the flipped-in and -out states is adopted. For the definitions of the characteristic angles, see references [21,22] and Supplementary Figure S3. To reproduce structural transitions from the flipping-in to -out states, the PaCS-MD was repeated until the 100th cycle. In each cycle, the snapshots with top 10 higher inner product values were selected as the initial structures, and 100-ps MD simulations under NPT ensemble (P=1 bar and T=400 K) were performed. To keep the shape of DNA at 400 K, a set of positional restraints was imposed for the DNA backbones except for CPD using harmonic potentials with a force constant (10 kcal/mol·Å2). Note here that there exist two different rotation directions, a clockwise direction and an anticlockwise one (Supplementary Fig. S4), and they were treated separately.

1.3. Fragment molecular orbital method

The fragment molecular orbital (FMO) method enables one to perform ab initio electronic structure calculations for large size systems [29-34]. In this method, the DNA is divided into small fragments and performs ab initio electronic structure calculations for all monomer fragments and fragment pairs within a two-body FMO (FMO2) approximation. The total energy is approximated using the total energies of the monomer fragment and the fragment pairs as

where N is the number of the monomer fragments and E is the total energy of the I th monomer fragment and E is that of a fragment pair IJ, respectively. The total energy is rewritten as

In the above expression, E′ (x=I for a monomer and IJ for a pair) is the energy without the environmental electrostatic energy contribution defined as

Vx and Dx are the matrix elements of the environmental potential and the atomic orbital-based density matrix, respectively. Δ is a difference density matrix defined as

with

The environmental potential Vx consists of the nuclear attraction, u, and the Hartree term, J, from the surrounding monomers K≠x given by

Here, ΔE=E′−E′−E′ is defined as the inter-fragment interaction energy (IFIE) [34,35]. To treat solvent effects around the DNA, a polarizable continuum model was adopted. For the calculations of FMOs, we used GAMESS program package [ver Dec 5, 2014] [36,37].

2. Results and Discussion

2.1. Conformational analyses by conventional MD and PaCS-MD simulations

We first performed classical MD simulations starting from both flipped-in and -out states of the DNA lesion (Fig. 1) for 500 ns. No base-flip (both forward and backward) was observed in either of the MD simulations. The time series of the root-mean-square deviation (RMSD) are shown in Figure 2 (a) and (b). The RMSDs for the flipped-in state kept almost the same values. On the other hand, those for the flipped-out state frequently fluctuated, though the drastic structural transition from the flipped-in to -out did not occur. The averaged RMSDs measured from the optimized flipped-in and -out state structures for the flipped-in (flipped-out) simulation were 1.60±0.18 Å (2.39±0.18 Å) and 3.02±0.12 Å (2.59±0.21 Å), respectively, indicating that the flipped-out state is less stable due to a flexible partial single strand moiety exposed outside the double strand. The time series of two characteristic dihedral angles, i.e. 5′-pseudodihedral and 3′-pseudodihedral angles [21,22], are shown in Figure 2 (c) and (d). In the simulation of the flipped-out state, these characteristic angles remarkably changed during the simulation in comparison with those in the simulation of flipped-in state. These changes correspond to the frequent fluctuations found in Figure 2 (a) and (b). The average dihedral angles in simulation of the flipped-in state were 70.17±6.32 (5′ side) and 7.78±8.18 (3′ side) degrees, and those in the simulation of the flipped-out state 109.90±68.79 (5′ side) and −10.05±19.61 (3′ side), again indicating the high flexibility of the flipped-out state characterized by these angles. Judging from these results, the flipped-in state remained in a relatively long time scale, meaning that the forward base-flipping process of CPD cannot occur spontaneously without any external forces. Although more fluctuations of these angles were observed in the simulation of the flipped-out state, the flipped-out state does not go back to the flipped-in state within this time scale, once the lesion CPD moiety was exposed to the bulk water. These classical MD simulations imply that both the forward and backward base-flipping processes of the DNA lesion should scarcely occur and should be carefully treated with efficient conformational sampling techniques to search structural transition pathways.

Figure 2

On the left side, is the root-mean square deviation (RMSD) in MD simulations for 500 ns. (a) RMSD of trajectories starting from the flipped-in state with the optimized flipped-in (blue line) and flipped-out states (green line). (b) RMSD trajectories starting from the flipped-out state with the optimized flipped-in and optimized flipped-out states. On the right side, is the dihedral angle in MD simulation for 500 ns. Time series of 5′pseudodihedral (blue line) and 3′pseudodihedral (green line) angles of trajectories starting from the flipped-in (c) and flipped-out (d) states.

Next, the PaCS-MD was applied to sample the structural transition pathways between the flipped-in and -out states. Note here that the PaCS-MD simulation at room temperature (T=300 K) did not induce the structural transition with the same simulation time (data not shown). This result suggests that the base-flipping process does not occur at ambient conditions. Instead, we performed the PaCS-MD simulation at a high temperature (T=400 K), which is higher than a typical temperature of DNA denaturation experiments (T~360 K) to promote the base-flipping process. Figure 3 (a) and (b) show trajectories obtained by the PaCS-MD projected onto a subspace spanned by the characteristic pseudodihedral angles for a clockwise and an anticlockwise rotation, respectively. Note here that the clockwise and anticlockwise rotations correspond to the movements, where the CPD moves outward and inward from the DNA double strand, respectively. During the 100 cycles, the PaCS-MD successfully sampled several conformational areas on the conformational subspace. The upper-right area corresponds to the flipped-in state, while the lower-left area corresponds to the flipped-out state, where the dashed lines in these figures indicate boundaries of the two states. After the conformational sampling with the PaCS-MD, the representative snapshots were selected along the structural transitions pathways between the flipped-in and -out states (see points in Fig. 3 with indication flipped-in state, (A), (B), (C), and (D)). These snapshots are depicted in Figure 4 and characteristic angles of them were listed in Table 1.

Figure 3

Projections of trajectories obtained by the PaCS-MD during the 100 cycles, for (a) clockwise and (b) anticlockwise rotations. All the trajectories were projected onto the conformational subspace spanned by the characteristic dihedral angles (5′ and 3′ pseudodihedral angles). The dashed lines indicate the boundaries for the flipped-in and flipped-out states of CPD.

Figure 4

The rotation direction of base-flipping for CPD. On the left side, is each structure (A), (B), (C), and (D) obtained PaCS-MD for clockwise rotation. On the right side, is each structure (A), (B), (C), and (D) obtained PaCS-MD for anticlockwise rotation. While the flipped-in state is obtained by PDBID 1SNH and the flipped-out state is modeled using winmostar software.

Table 1

The dihedral angles for the flipped-in state, flipped-out state, (A), (B), (C), and (D) of clockwise rotation and anticlockwise rotation

Structure	Clockwise rotation [deg.]		Anticlockwise rotation [deg.]
	5′pseudodihedral	3′pseudodiheral	5′pseudodihedral	3′pseudodihedral
Flipped-in	78.15	1.32	78.15	1.32
(A)	91.54	−52.63	74.77	1.74
(B)	95.59	−62.27	7.52	53.37
(C)	146.70	−124.04	−75.88	84.23
(D)	−62.88	75.16	−79.84	84.76
Flipped-out	144.67	−21.80	144.67	−21.80

2.2. Inter-fragment interaction energy (IFIE) and total energy analyses

To investigate how the CPD moiety interacts with other bases during the base-flipping process, we calculated the IFIEs for the NMR flipped-in structures (PDBID: 1SNH), the model flipped-out structure, the structures (A), (B), (C), and (D) obtained by PaCS-MD using the FMO method at PCM-FMO2-MP2/6-31G(d) level. For the FMO calculations, we replaced PO3 with H except for CPD, and H was added to the remaining O (Supplementary Fig. S5). We divided the model DNA lesion into 19 fragments. We hereby show IFIEs for the optimized flipped-in and -out states, (A), (B), (C), and (D) in Figure 5. According to these results, the CPD strongly interacts with the neighbors (DA4 and DA6) and the face-to-face bases (DA14 and DA15) at the initial flipped-in state. The initial flipped-in state has two hydrogen bonds between CPD and DA14 and between CPD and DA15, and the corresponding IFIEs were evaluated to be −21.52 and −14.46 kcal/mol, respectively. Moreover, the IFIEs of the CPD and DA4 or DA6 exhibit strong attractive interactions, −10.96 or −13.70 kcal/mol due to the π-π stacking interaction. A partial sum of the IFIEs between CPD and DNA bases except for DA4, DA6, DA14, and DA15 for the flipped-in (flipped-out) state was 0.563 (−2.355) kcal/mol. Both in the clockwise and anticlockwise rotations, absolute values of IFIEs were gradually reduced as the dihedral angles change toward the flipped-out state.

Figure 5

Inter-fragment interaction energy for each structure in clockwise and anticlockwise rotations. The interaction energy decreases with transitioning from flipped-in state to flipped-out state both clockwise and anticlockwise.

Another partial sum of IFIEs between CPD and the rest of the DNA bases, i.e. DA4, DA6, DA14, and DA15, was −60.08 kcal/mol for the flipped-in state and −10.40 kcal/mol for the flipped-out state, respectively. These results obviously suggest that the base-flipping of DNA cannot occur by thermal fluctuation. The total energy for flipped-out of each structure is shown in Table 2. The energy differences measured from the initial flipped-in state for structures (A), (B), (C), and (D) in the clockwise (anticlockwise) rotation were 339.69 (335.81), 339.50 (356.16), 342.31 (399.90), and 26.87 (50.83) kcal/mol, respectively. These results indicate that the energy barrier from the flipped-in to –out is quite high. In order to reveal the reason, we performed the FMO calculation for the systems without the CPD moiety. The energy differences for structures (A), (B), (C), and (D) in the clockwise (anticlockwise) rotation were 8.71 (−2.71), −4.47 (−2.68), −12.66 (20.82) and 1.80 (0.64) kcal/mol, respectively. This result clearly indicates that the high energy barrier comes from destabilization in a distortion of structures and solvation effects due to the CPD moiety during the structural transition.

Table 2

The total energy for with CPD and without CPD using FMO2-MP2/6-31G(d)

Structure	Clockwise rotation			Anticlockwise rotation
	Total Energy [kcal/mol]			Total Energy [kcal/mol]
	with CPD	without CPD	Energy difference	with CPD	without CPD	Energy difference
Flipped-in	0.00	0.00	0.00	0.00	0.00	0.00
(A)	339.70	8.71	330.99	335.81	−2.71	338.52
(B)	339.50	−4.47	343.97	356.16	−2.68	358.84
(C)	342.31	−12.66	354.97	399.90	20.82	379.08
(D)	26.87	1.80	25.07	50.83	0.64	50.19
Flipped-out	41.79	13.02	28.77	41.79	13.02	28.77

For all cases, the energy without CPD is little less than that of the flipped-in state.

Since the base-flipping of CPD does not occur by itself as shown in this work, it is natural to consider that driving forces exist to promote the base-flipping process and that the differences in the photorepair activity on DNA among the same family of enzymes may originate from the difference in driving forces. In order to characterize them, we hereby show the charges on the surface of PHR and CRY-DASH by a Poisson-Boltzmann method as shown in Figure 6 [38-43]. In comparing PHR with CRY-DASH, the former exhibit more positive charge distribution around the binding pocket than the latter. This result suggests that PHR has stronger interactions with DNA lesions than CRY-DASH does. Therefore, it is expected that PHR can promote the base-flipping process of the DNA lesion and CRY-DASH cannot promote it. In future work, we will calculate the interaction energy between protein and UV-induced DNA and the excited-state electron transfer efficiency of both cryptochromes by means of first-principles method to fully understand the origin of the repair activity of PHR.

Figure 6

Surface charges for Photolyase and Cryptochrome DASH. Blue is the positive charge and red is the negative charge. The positive charge of photolyase distributes around the binding pocket more than cryptochrome DASH.

3. Conclusion

It was experimentally confirmed that the base-flipping process of the duplex DNA with UV-lesions occurs when it binds DNA photolyase. However, its molecular detail has remained unknown. We performed MD simulation for 500 ns for both the flipped-in and -out states and found that both the forward and backward base-flipping processes do not occur during the simulation time. These results indicate that the base-flipping process of the DNA lesion cannot occur without other interactions. We have numerically investigated the reason why the base-flipping process of the DNA legion cannot occur by itself based on the analyses on the conformation sampling by PaCS-MD followed by the inter-fragment interaction energy (IFIE) by FMO. Judging from these analyses, the CPD moiety of the DNA lesion strongly interacts with neighboring bases. The total IFIEs were −60.08 kcal/mol for the flipped-in state and −10.40 kcal/mol for the flipped-out state. This energy difference is larger than the thermal energy. Therefore, we suggest that DNA lesions require other driving forces, such electrostatic interaction, to overcome the energy barrier for the transformation from flipped-in to -out state.

The repair of UV-induced DNA is one of the critical reactions in life. A dysfunction of the DNA repair will cause severe diseases such as skin cancer in animals and growth inhibitor in plants. Although this reaction has been extensively investigated, even molecular mechanism of the base-flipping process, which is an early process of DNA repair, has not been revealed, yet. In this work, we analyzed the base-flipping process of a DNA lesion alone by using molecular dynamics simulations and quantum mechanics calculations. Our results clearly indicate that the base-flipping process does not occur without a help of photolyase.