Effect O6-Guanine Alkylation on DNA Flexibility Studied by Comparative Molecular Dynamics Simulations Mahmut Kara,1 Tomas Drsata,2,3 Filip Lankas,2 Martin Zacharias1 1

Physik-Department T38, Technische Universit€at M€ unchen, James-Franck-Strasse, D-85748 Garching Germany

2

Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, Flemingovo namesti 2, 166 10 Prague, Czech Republic

3

Department of Physical and Macromolecular Chemistry, Faculty of Science, Charles University Prague, Albertov 6, 128 43 Prague, Czech Republic Received 29 May 2014; revised 4 August 2014; accepted 11 August 2014 Published online 16 August 2014 in Wiley Online Library (wileyonlinelibrary.com). DOI 10.1002/bip.22535

ABSTRACT: Alkylation of guanine at the O6 atom is a highly mutagenic DNA lesion because it alters the coding specificity of the base causing G:C to A:T transversion mutations. Specific DNA repair enzymes, e.g. O6-alkylguanin-DNATransferases (AGT), recognize and repair such damage after looping out the damaged base to transfer it into the enzyme active site. The exact mechanism how the repair enzyme identifies a damaged site within a large surplus of undamaged DNA is not fully understood. The O6alkylation of guanine may change the deformability of DNA which may facilitate the initial binding of a repair enzyme at the damaged site. In order to characterize the effect of O6-methyl-guanine (O6-MeG) containing base pairs on the DNA deformability extensive comparative

were performed. The simulations indicate significant differences in the helical deformability due to the presence of O6-MeG compared to regular undamaged DNA. This includes enhanced base pair opening, shear and stagger motions and alterations in the backbone fine structure caused in part by transient rupture of the base pairing at the damaged site and transient insertion of water molecules. It is likely that the increased opening motions of O6-MeG:C or O6-MeG:T base pairs play a decisive role for the induced fit recognition or for the looping out of the C 2014 Wiley Periodidamaged base by repair enzymes. V

cals, Inc. Biopolymers 103: 23–32, 2015. Keywords: DNA damage; DNA alkylation; DNA repair; molecular

simulation;

molecular dynamics

simulation

molecular dynamics (MD) simulations on duplex DNA with central G:C, O6-MeG:C or O6-MeG:T base pairs Additional Supporting Information may be found in the online version of this article. Correspondence to: Martin Zacharias; e-mail: [email protected] Contract grant sponsor: Deutsche Forschungsgemeinschaft (DFG) Contract grant number: SFB749/project C5 Contract grant sponsor: Leibnitz Rechenzentrum, Germany Contract grant number: pr86pu Contract grant sponsor: Academy of Sciences of the Czech Republic Contract grant number: RVO61388963 (T.D. and F.L.) Contract grant sponsor: Grant Agency of the Czech Republic Contract grant number: 14-21893S (T.D. and F.L.) Contract grant sponsor: Grant Agency of the Charles University Contract grant number: 584213 (T.D.) C 2014 Wiley Periodicals, Inc. V

Biopolymers Volume 103 / Number 1

This article was originally published online as an accepted preprint. The “Published Online” date corresponds to the preprint version. You can request a copy of any preprints from the past two calendar years by emailing the Biopolymers editorial office at [email protected].

INTRODUCTION

A

common type of damage in DNA is caused by alkylating agents which can lead to chemical modifications at several different positions, including the heterocyclic bases and the backbone.1–5 If not repaired such damage can result in mutations upon

23

24

Kara et al.

FIGURE 1 Geometry of the central base pairs in the simulation starting structures (upper panels) and illustration of the duplex DNA with the central O6-MeG nucleotide in van der Waals representation (lower panel).

DNA replication. One of the most highly mutagenic modifications on guanine is the attack of the O6 atom in the base (e.g. by N-methyl-N-nitrosourea and N-methyl-N0 -nitro-Nnitrosoguanidine) resulting, for example, in O6methylguanine (O6-MeG). The additional methyl-group changes the tautomerism of the base which allows pairing with a thymine or cytosine. In subsequent replication rounds, this can lead especially to G:C to A:T transversion mutations. DNA containing an O6-MeG can be recognized and repaired by different repair enzymes.3–10 Among these are oxygenases of the AlkB family that directly repair the damage and glycosylases that excise the damaged base followed by a DNA polymerase mediated template-based synthesis of the damaged strand.4–6 In case of specific O6-alkylguanin-DNA-transferases (AGTs), the alkyl group is transferred to a nucleophilic cysteine residue in the repair enzyme active site.4,5,8,9 This transition is irreversible and the enzyme is degraded after the repair by the ubiquitin pathway (also termed “suicide enzyme”). The active O6-Alkylguanine-DNA alkyltransferase has a helix-turn-helix (HTH) motif that interacts with the minor groove of the DNA double strand while the O6-alkylation is located in the major groove.8,10 The interaction introduces a conformational change, so that it increases the size of the minor groove of DNA such that the intercalation of an arginine finger in the helix becomes possible coupled to transition of the alkylated base to an extrahelical state in the enzyme active site.8–10 The structure of O6-MeG:C and O6-MeG:T base pairs has been investigated using nuclear magnetic resonance (NMR) spectroscopy11–13 and X-ray crystallography14–17 as well as

crystal structure analysis of O6-MeG containing DNA in complex with DNA polymerases.18 The O6-MeG:C base pair adopts a (regular) anti-glycosidic dihedral state with a “wobble” hydrogen-bonding and with the O6 methyl group solvent exposed (Figure 1). At acidic pH, a different hydrogen bonding pattern is possible due to a different tautomer structure of the O6-MeG base.17 In dodecamer B-DNA crystal structures containing a O6-MeG:T base pair two hydrogen bonds between the two bases can be assigned (Figure 1). In contrast, a O6MeG:T pair with only a single hydrogen bond, was inferred from NMR spectroscopy.11 The overall shape of the O6-MeG:T base pair is similar in both the NMR and crystallographic structures.11,17 The efficient localization of a base lesion and base extrusion is a long-standing general question. There are three different hypotheses how the enzyme locates the base lesion.4,5 One possibility is, that unspecific binding of the enzyme to DNA helps to flip out the base and distinguishes if damaged or not after inspection in the active site or some other recognition region. A second possibility is a stable binding of the enzyme only to an already looped out damaged base. This would imply that the enzyme has no influence on the barrier of looping out and just has to “wait” for the spontaneous looping out of a damaged base. Another third possibility is the specific recognition of some structural property of the DNA damage with the damaged base still in the intrahelical state. Enhanced binding to such motif would increase the probability for looping out during this transient binding and accelerates the formation of a complex with the enzyme and a damaged base. Molecular dynamics (MD) simulations are a useful tool to compare the fine structure and dynamics of O6-MeG containing DNA and undamaged DNA at atomic resolution and under realistic solution conditions including explicit solvent and surrounding ions.19 This could help to distinguish between the above discussed possible mechanisms of damage recognition. Previous simulation studies on alkylated DNA have focused on the sequence specificity19 and the mechanism of looping out the damaged nucleotide in the presence of a repair enzyme.20 In the present study, comparative MD simulations on O6MeG:C, O6-MeG:T, and G:C base pairs located at the center of a dsDNA and in the absence of a repair enzyme have been performed. In addition to classical MD simulation, we applied a Hamiltonian-replica exchange (H-REMD) simulation technique to specifically enhance the sampling of possible substates of DNA molecules [termed biasing potential (BP)REMD].21,22 Significant differences in the helical mobility of the base pairs containing O6-MeG compared to regular base pairs were observed. This included especially enhanced base pair opening and shear motions at the damaged site and transient disruption of base pairing coupled to insertion of water Biopolymers

Effect O6-Guanine Alkylation on DNA Flexibility

25

FIGURE 2 Comparison of root-mean-square deviation (RMSD) vs. simulation time of DNA heavy atoms from B-DNA starting structure during continuous MD simulations (left panels) and in the reference replica (without biasing) of the BP-REMD simulations (right panel).

molecules. The enhanced helical mobility likely plays a decisive role for both initial recognition of the damage by the bulky arginine (Arg) finger residue of the alkyltransferases which may probe local opening flexibility and for the energy barrier of the damage looping out process itself.

MATERIALS AND METHODS All MD simulations were started from regular B-DNA of the sequence d(5’-GCCATG(6OG)CTAGTA) and a complementary strand with either a C or T opposite to the O6-MeG (Figure 1). For comparison, simulations on regular undamaged B-DNA of the same sequence d(5’-GCCATGGCTAGTA) and complementary strand) were also performed. The sequence was chosen because it corresponds to an experimentally known structure (pdb1T83) in complex with the O6alkylguanine-DNA alkyltransferase (AGT).9,10 The Amber12 Molecular Dynamics Package23 was used for the simulations. All simulations were performed in explicit water (TIP3P)24 with a truncated octahedral box and a minimum distance of 10 A˚ between DNA and box boundary. Sodium counter ions were included to neutralize the system. All classical MD simulations were carried out with the pmemd module of the AMBER12 package using the parmbsc0 force field for DNA25 Force field parameters specific for the O6-MeG nucleobase were taken from the modified nucleic acid data base.26 The nucleic acid backbone of the O6-MeG included the pambs0 correction. The simulation systems were first energy minimized (5000 steps) with restraints on DNA followed by 5000 steps of unrestrained energy minimization. The system was heated up in three steps (each step 100 ps)

Biopolymers

to 300 K in steps of 100 K followed by gradual removal of the positional restraints from 25 to 0.5 kcalmol21A˚22 (in five steps) and a 1ns unrestrained equilibration at 300 K. The equilibrated structures served as starting structures for continuous (c)MD and BP-REMD simulations (see below). Analysis of helical parameters was performed with Curves127 or 3DNA.28 In order to control the effect of the sequence context on the dynamics of base pairs involving O6-MeG bases additional cMD simulations of 100 ns were performed on G:C, O6-MeG:C, and O6-MeG:T base pairs flanked by A:T base pairs. In this case, the DNA sequence was d(5’-GCCATA(6OG)TTAGTG) and a complementary strand with either a C or T opposite to the O6-MeG and sequence d(5’-GCCATAGTTAGTG) in case of the reference with a regular central G:C. The same equilibration and production protocol as described above was used.

Biasing Potential Replica-Exchange Based on Pseudotorsion Angles For the biasing potential replica-exchange MD simulations (BPREMD), a two-dimensional (2D) biasing potential acting on two dihedral angles g and h was employed for a given nucleotide. The approach followed a protocol used previously in a study of oxidized guanine in DNA.22 The pseudo-dihedral angles are defined by the nucleic acid backbone atoms O3’-O5’-C4’-O3’ in case of g (here the first O3’ belongs to the previous nucleotide) and C4’-O3’-O5’-C4’ in case h (O5’ and C4’ belong to the next nucleotide). Changes in these pseudo dihedral angles correlate with the two most important conformational substates compatible with B-DNA: Coupled changes of the a and c dihedral angles (also called crankshaft motions or a/c-flip)29,30

26

Kara et al.

and coupled changes of e and f dihedral angles.22,29,31 The main coupled e/f states are called BI (corresponds to the most common substate of regular B-DNA) and BII.31 The addition of a 2D biasing potential for the two coupled g and h dihedral angles centered at the position of known substates can result in a destabilization of favorable regions and hence accelerates transitions to alternative nucleic acid backbone substates (helps to promote transitions over energy barriers separating favorable substates). A BP-REMD simulation consisted of a set of nine parallel replica MD simulations with different levels of the above 2D biasing potential applied to the g, h pair of a selected nucleotide. The biasing potentials were centered on the known centers for BI and BII states as well as for a/c-flipped states. The width of each potential was 60 . In each BPREMD simulation, one of the replicas was simulated using the original force field without biasing and served as the reference replica used for subsequent analysis. Exchanges between neighboring replicas were attempted every 1000 steps (2 ps) and accepted or rejected according to a Metropolis criterion.22 The magnitude of the applied biasing potential levels in each replica differed by 0.66 kcalmol21. Test simulations indicated that this choice resulted in an acceptance probability for replica exchanges of 50%. The added biasing potentials in the replica runs were applied to the central alkylated nucleotide as well as to the opposing nucleotides.

Analysis of DNA Stiffness In order to systematically compare the conformational flexibility of regular DNA and DNA containing O6-MeG nucleotides the MD trajectories were analyzed in terms of a nonlocal harmonic rigid base model.32–34 The trajectories are used to extract the variances and covariances of all intra- and interhelical parameters associated with the sampled DNA conformations resulting in a stiffness matrix to describe the deformability of the duplex in terms of effective force constants for helical deformations. Details of the approach have been described in previous publications.33–35 The overall stiffness is characterized by conformational entropy per rigid base coordinate, as described earlier.36

RESULTS AND DISCUSSION Conformational Flexibility of O6-Methylated Guanine Containing DNA and Undamaged DNA In order to compare the conformational flexibility of O6alkylated DNA and undamaged DNA, MD simulations in explicit solvent were performed (200 ns) starting from regular B-DNA structures with a central G:C, an O6-MeG:T, or an O6MeG:C base pair (Figure 1, sequence given in Materials and Methods section). The average root mean-square deviation (RMSD) of backbone atoms from B-DNA start structures was stable reaching up to 4 A˚ which can be attributed to global bending fluctuations of the helix and in part to fraying (transient opening) of the terminal base pairs (Figure 2). Both the RMSD distribution with respect to B-DNA as well as to the DNA conformation found in complex with the AGT enzyme (leaving out the central base pair) were similar for all three

FIGURE 3 Simulation snapshots of the central O6-MeG:C (upper panel) and O6-MeG:T (middle panel) base pairs (stick representation). The first and the second snapshots correspond to frequently sampled conformations close to the starting geometry. The snapshots including a water molecule between the two central bases represent transient states with durations of a few ps. (Lower panel) Distribution of the sampled opening angle for the central base pair and for the adjacent (previous) base pair.

DNA molecules indicating no selective global preference of alkylated DNA toward the DNA conformation found in complex with the repair enzyme (not shown). In the initial geometry the O6-MeG:T and O6-MeG:C base pairs formed two possible hydrogen bonds (Figure 1a). For the central regular G:C case, the typical Watson–Crick hydrogen bonding pattern remained stable throughout the whole simulation time whereas fluctuations between 0, 1 or 2 hydrogen bonds were observed in base pairs of the DNA with central O6MeG (Figure 3). The dominant base pair state resembled the conformation found in X-ray structures of dsDNA containing O6-MeG opposite to T or C, respectively. Although the central damaged base pair showed increased fluctuations compared to the intact G:C base pair (see next paragraphs) no transitions to Biopolymers

Effect O6-Guanine Alkylation on DNA Flexibility

27

fully looped out structures with completely solvent exposed bases were observed. However, within short intervals of 5–10 ps transient significant deviations from the starting base pair geometry were observed (Figure 3). In these states, a transient rupture of the base pair geometry and partial opening toward the solvent occurred with water molecules diffusing into the space between the central nucleo-bases (Figure 3). This leads to a shift and broader distribution of the base pair opening angle for the O6-MeG opposite to T or C compared to a regular base pair (Figure 3). Interestingly, the shift is frequently for the opening of O6-MeG:T in the opposite direction compared to O6-MeG:C. In case of a O6-MeG:C, the water insertion resulted in a movement of the O6-MeG base toward the major groove whereas the opposing T moves toward the major groove in case of O6-MeG:T base pair (Figure 3). Such transient states may represent intermediates on a transition pathway toward a looped out extrahelical conformation of the damaged base and indicate an increased looping out tendency at the damage site. Alternatively, repair enzymes may recognize such flexible sites by an induced fit mechanism which allows tighter initial binding in case of a base pair with an alkylated base.

Nucleic Acid Backbone Structure at Damaged Nucleotide Complex formation between DNA and proteins often involves transitions in the sugar-phosphate backbone as well as global DNA conformational changes.29–31,36 Beside a direct interaction with the alkylated DNA recognition may in part also be mediated by a specific readout of an altered helical or backbone DNA structure containing alkylated bases.17,37 Unfortunately, due to energy barriers separating alternative backbone states the sampling of such states during conventional MD simulations may be limited. We have recently developed a Hamiltonian replica exchange MD (H-REMD) method that employs a specific biasing potential to promote transitions between nucleic acid substates.22 This technique, termed biasing potential (BP)-REMD, was shown to significantly improve the sampling of relevant conformational backbone substates in case of oxidatively damaged Guanine in DNA.22 Following the protocol used previously the method was applied using 9 parallel simulations for 8 ns with varying 2D biasing potentials and one reference replica without biasing potential that served for the evaluation of substates (see Methods for details). The most important fine structure substates in DNA involve coupled changes in the e/f torsion angles (giving rise to BI/BII substates), the a/c angles (crankshaft motion) and the pucker states of the deoxy ribose sugar. The analysis of characteristic nucleic acid substates indicated no significant differences in e/f substates between a central regular G:C base pair or the Biopolymers

FIGURE 4 Comparison of sampled a/c dihedral angle states observed during cMD simulations (left panels) and during BPREMD (right panel) at the G (upper two panels), O6-MeG opposite to a C (middle panels) or O6-MeG opposite to T (lower panels). Each point corresponds to a sampled dihedral combination of the 200-ns cMD or 8-ns BP-REMD (same point density in each plot).

alkylated modifications, with the BI state dominating in all simulations (Supporting Information Table S1). Note, that also in experimentally determined DNA structures with O6-MeG bases a dominance of the BI state was found (Supporting Information Figure S1). The BP-REMD simulations indicated sampling of two additional alternative a/c states in case of a central O6MeG:C and one additional substate in case of the central O6MeG:T pair (Figure 4). However, the abundance of these substates was