A Molecular Dynamics Simulation of the (dC), (dC), Minihelix Including Counterions and Water THERESA JULIA ZlELlNSKl’ and MASAYUKI SHIBATA2 ‘Department of Chemistry, Niagara University, Niagara University, New York 14109, and 2Biophysics Department, Roswell Park Memorial Institute, Buffalo, New York 14263
SYNOPSIS
The results of a 60 ps molecular dynamics (MD) simulation of (dG), . (dC), including 10 Na+ counterions and 292 water molecules are presented. All backbone angles and helix parameters for the hexamer are reported in this paper along with trajectory plots of selected angles. Hydrogen bonding between the bases along the helical axis was observed to fluctuate with time, showing the dynamic nature of the base-pairing interaction. These fluctuations gave rise to unusual hydrogen-bonding patterns. Good intrastrand base stacking and no interstrand base stacking were also observed. The hexamer minihelix retains an essentially B-DNA conformation throughout the entire simulation even though some helix parameters and backbone angles do not have strict B-DNA values. The most striking feature obtained from the simulation was a high propeller twist, which resulted in a narrow minor groove for the minihelix. It is proposed that (dG), . (dC), sequences are resistant to DNAase I because of this narrow minor groove in dilute aqueous solution.
INTRODUCTI0N Nucleotide sequencing techniques have shown that strings of 10 or more contiguous purine or pyrimidine residues are found four times more frequently than expected for a random distribution.’ In the human P-globin region, 12 runs of 10 or more As or Ts, including one run of 12 Gs, have been found out of a set of 269 strings. Another homopurine . homopyrimidine tract, (dG)16 . (dC),,, found upstream of the chicken P*-globin gene,2>3has hypersensitivity to S1 nuclease activity related to transcription of the gene. Homopurine homopyrimidine sequences have been linked to DNAase I hypersensitive sites (nucleosome-free regions that have been correlated with gene e x p r e ~ i o n ) .More ~ recently, an endonuclease that selectively cleaves DNA a t (dG), . (dC), tracts has been identified.5 This endonuclease, called endonuclease G, was 0 1990 John Wiley & Sons, Inc.
CCC 0006-3525/90/051027-18 $04.00 Biopolymers, Vol. 29, 1027-1044 (1990)
found to be present in cells from different origins. Since the homopolymer tracts appear next to regions where rearrangement, translocation, deletion, and hypervariability of DNA occurs, i.e., recombinational “hot spots”, endonuclease G activity and homopolymer sequences could be involved in recombination. These observations have increased the interest in homopurine . homopyrimidine, mixed purine mixed pyrimidine, and in particular, (dG), * (dC), nucleic acid base sequences. The relationship between the helix geometry and biochemical behavior of these sequences becomes an important issue that can lead to greater insight into their biological function. The x-ray crystallographic techniques have been applied to GC-rich oligomers. The longest run of ~ - the ~ basis of guanines in such a study is f o ~ r . On their results, McCall et al. developed a model for poly(dG) p01y(dC).~This model had a wide shallow minor groove and a deep cavernous major groove. The mean roll, slide, and twist values were 5”, 0.19 nm, and 32.1”, respectively. There were
-
1027
1028
ZIELINSKI AND SHIBATA
11.2 base pairs per turn with a mean rise of 0.288 nm per residue and a tilt of 12". The base pairs were displaced by 0.49 nm from the helix axis, placing the model firmly in the A-DNA family. The P-P distance across the major groove was 1.39 nm. By analogy with the d(G,C,) structure, it was suggested that the major groove of the model would be filled with water and that the ion/water bridges would link phosphates on opposite strands. In the model the poly(dG) strand showed base stacking involving the five-membered ring of one base stacking over the six-membered ring of its neighbor. The poly(dC) strand showed little or no stacking. It has been proposed that DNAase I is sensitive to DNA backbone geometry as determined by base sequence." DNAase I cuts most poorly at A . T or G C homopurine . homopyrimidine stretches. Using the available crystallographic data it was suggested that the minor groove width is too small in (dA), . (dT), and too large in (dG), . (dC), to fit the enzyme. The enzyme was used with the synthetic oligomer (dA),,(dG),, containing flanking extremes of conformational preference, and with a naturally occurring sequence G,,AGA, from downstream of the gene for the 5sRNA of sea urchin. While the results were not conclusive, the cutting patterns indicate an important role for flanking base sequences in determining groove width and consequent enzyme recognition." The importance of groove width is also supported by the effect of antibiotics such as actinomycin and distamycin on DNAase I activity. Distamycin induces a structural change-narrowing of the minor groove-in neighboring G . C-rich sequences, which enhances the cutting rate of the enzyme.12 The Drew and Travers model for the action of DNAase I depended on the structural details obtained from single crystal x-ray diffraction studies. These structural details may not be applicable when the oligomers are found in dilute aqueous solution. This point is raised as a consequence of the Raman spectral studies of poly(dG) . poly(dC) in fiber and s01ution.l~In the concentrated solution the polymer behaves as predominantly A-DNA, but displays a predominantly B-DNA form in dilute solution. The fiber polymer, interestingly, is an intermediate case with a non-A structure, which departs significantly from typical B-DNA, and may involve variants of the C,-endo/anti G including the C,-exo/unti G conformation. In another experiment it was found that the isotopic exchange 1 H-3H for the C8H of guanine of poly(dG) . poly(dC) in aqueous dilute solution with 0.15M NaCl is retarded compared to the guanine in
-
Escherichia coli DNA and dGMP.', This corresponds to a decrease in the accessibility of the C8H in the homopolymer, which can be interpreted in two ways. Either the structure of poly(dG) . poly(dC) in this solution is intermediate between A- and B-DNA forms, or there is an equilibrium between the A- and B-DNA forms with approximately equal concentrations of both. Lesnik et al.', prefer the latter interpretation. In higher 3M NaCl solutions, poly(dG) . poly(dC) forms intermolecular aggregates with a corresponding shift to the ADNA form even with low polymer concentrations. It has been proposed that the G * C sequence is important for stabilizing the A-form of DNA.l5,I6 In light of the structural variability of poly(dG) . poly(dC) found by spectroscopic studies of aqueous solutions and the difference between solution and solid state results, it would be timely to examine this system by a completely independent method. Molecular dynamics (MD) has been shown to be useful for elucidating the structural features of oligonucleotides in vacuum and in ~ o l u t i o n . ' ~ - ~ ~ In this paper we report a model for poly(dG) . poly(dC) obtained by MD simulation of the hexamer (dG), . (dC), with explicitly included Na+ counterions and solvent.
METHOD The initial hexamer structure was constructed from standard B-DNA coordinates with a base-pair step height of 0.34 nm and a helix twist of 36.0°.21The 5' and 3' terminals ended in hydroxy group?. The united atom approximation was used for all carbon hydrogens while all other hydrogens were treated explicitly. Electrical neutrality was achieved by placing octahedrally hydrated Na ions, coordination number of 6,22 a t each of the phosphate positions such that each of the free phosphate oxygens formed a bifurcated H-bonding structure with one of the waters of hydration of a sodium ion. The resulting complex (solute, sodium ions, and their waters of hydration) was then submerged in a box of SPC water.23 The GROMOS (Groningen Molecular Simulation System)24 program PROBOX generated 1581 extra waters around the solute complex. The total number of waters around the minihelix plus sodium ion system was then reduced to 292 by discarding any waters beyond 0.5 nm of a solute atom. The net effect is that of a glove of water molecules over the surface of the solute. These 292 (the 60 Na+ waters are a subset) water +
MD SIMULATIONS OF (dG), . (dC),
molecules should give a reasonable description of the first hydration shell since approximately 300 water molecules hydrate one (B-DNA and counterions) turn of the helix with about 240 of these in the first hydration In the dynamics simulation the minihelix, 274 atoms, and the 10 associated Na' were grouped together as the solute. All water molecules were grouped in the solvent set. In order to eliminate the possibility of fraying of the helix during the simulation the center hydrogen-bonded pair of heavy atoms for the base pairs at the top and bottom of the minihelix were restrained to their initial interatomic distances using SHAKE. The potential energy function of the GROMOS programs was used without m~dification.~~ GROMOS partial charges were also used with a 1.0 dielectric constant. The cutoff range for nonbonded pair interactions was set at 0.8 nm and evaluated every MD step. The electrostatic interaction was calculated using a 1.8 nm cutoff and a pair list that was updated every 10 time steps. The complete solute-solvent system was relaxed with conjugate gradient energy minimization until the difference in energy between two successive steps was less than The initial velocities for the dynamics run were taken from a Maxwellian distribution. Constant temperature was maintained by weakly coupling the system to a 300 K thermal bath. No pressure effects were included. The MD run extended for 60 ps with a step size of 0.001 ps for the first 20 ps in order to equilibrate the system to the chosen temperature. A step size of 0.002 ps was used thereafter. No restrictions were placed on bond lengths,
SYNOPSIS
The results of a 60 ps molecular dynamics (MD) simulation of (dG), . (dC), including 10 Na+ counterions and 292 water molecules are presented. All backbone angles and helix parameters for the hexamer are reported in this paper along with trajectory plots of selected angles. Hydrogen bonding between the bases along the helical axis was observed to fluctuate with time, showing the dynamic nature of the base-pairing interaction. These fluctuations gave rise to unusual hydrogen-bonding patterns. Good intrastrand base stacking and no interstrand base stacking were also observed. The hexamer minihelix retains an essentially B-DNA conformation throughout the entire simulation even though some helix parameters and backbone angles do not have strict B-DNA values. The most striking feature obtained from the simulation was a high propeller twist, which resulted in a narrow minor groove for the minihelix. It is proposed that (dG), . (dC), sequences are resistant to DNAase I because of this narrow minor groove in dilute aqueous solution.
INTRODUCTI0N Nucleotide sequencing techniques have shown that strings of 10 or more contiguous purine or pyrimidine residues are found four times more frequently than expected for a random distribution.’ In the human P-globin region, 12 runs of 10 or more As or Ts, including one run of 12 Gs, have been found out of a set of 269 strings. Another homopurine . homopyrimidine tract, (dG)16 . (dC),,, found upstream of the chicken P*-globin gene,2>3has hypersensitivity to S1 nuclease activity related to transcription of the gene. Homopurine homopyrimidine sequences have been linked to DNAase I hypersensitive sites (nucleosome-free regions that have been correlated with gene e x p r e ~ i o n ) .More ~ recently, an endonuclease that selectively cleaves DNA a t (dG), . (dC), tracts has been identified.5 This endonuclease, called endonuclease G, was 0 1990 John Wiley & Sons, Inc.
CCC 0006-3525/90/051027-18 $04.00 Biopolymers, Vol. 29, 1027-1044 (1990)
found to be present in cells from different origins. Since the homopolymer tracts appear next to regions where rearrangement, translocation, deletion, and hypervariability of DNA occurs, i.e., recombinational “hot spots”, endonuclease G activity and homopolymer sequences could be involved in recombination. These observations have increased the interest in homopurine . homopyrimidine, mixed purine mixed pyrimidine, and in particular, (dG), * (dC), nucleic acid base sequences. The relationship between the helix geometry and biochemical behavior of these sequences becomes an important issue that can lead to greater insight into their biological function. The x-ray crystallographic techniques have been applied to GC-rich oligomers. The longest run of ~ - the ~ basis of guanines in such a study is f o ~ r . On their results, McCall et al. developed a model for poly(dG) p01y(dC).~This model had a wide shallow minor groove and a deep cavernous major groove. The mean roll, slide, and twist values were 5”, 0.19 nm, and 32.1”, respectively. There were
-
1027
1028
ZIELINSKI AND SHIBATA
11.2 base pairs per turn with a mean rise of 0.288 nm per residue and a tilt of 12". The base pairs were displaced by 0.49 nm from the helix axis, placing the model firmly in the A-DNA family. The P-P distance across the major groove was 1.39 nm. By analogy with the d(G,C,) structure, it was suggested that the major groove of the model would be filled with water and that the ion/water bridges would link phosphates on opposite strands. In the model the poly(dG) strand showed base stacking involving the five-membered ring of one base stacking over the six-membered ring of its neighbor. The poly(dC) strand showed little or no stacking. It has been proposed that DNAase I is sensitive to DNA backbone geometry as determined by base sequence." DNAase I cuts most poorly at A . T or G C homopurine . homopyrimidine stretches. Using the available crystallographic data it was suggested that the minor groove width is too small in (dA), . (dT), and too large in (dG), . (dC), to fit the enzyme. The enzyme was used with the synthetic oligomer (dA),,(dG),, containing flanking extremes of conformational preference, and with a naturally occurring sequence G,,AGA, from downstream of the gene for the 5sRNA of sea urchin. While the results were not conclusive, the cutting patterns indicate an important role for flanking base sequences in determining groove width and consequent enzyme recognition." The importance of groove width is also supported by the effect of antibiotics such as actinomycin and distamycin on DNAase I activity. Distamycin induces a structural change-narrowing of the minor groove-in neighboring G . C-rich sequences, which enhances the cutting rate of the enzyme.12 The Drew and Travers model for the action of DNAase I depended on the structural details obtained from single crystal x-ray diffraction studies. These structural details may not be applicable when the oligomers are found in dilute aqueous solution. This point is raised as a consequence of the Raman spectral studies of poly(dG) . poly(dC) in fiber and s01ution.l~In the concentrated solution the polymer behaves as predominantly A-DNA, but displays a predominantly B-DNA form in dilute solution. The fiber polymer, interestingly, is an intermediate case with a non-A structure, which departs significantly from typical B-DNA, and may involve variants of the C,-endo/anti G including the C,-exo/unti G conformation. In another experiment it was found that the isotopic exchange 1 H-3H for the C8H of guanine of poly(dG) . poly(dC) in aqueous dilute solution with 0.15M NaCl is retarded compared to the guanine in
-
Escherichia coli DNA and dGMP.', This corresponds to a decrease in the accessibility of the C8H in the homopolymer, which can be interpreted in two ways. Either the structure of poly(dG) . poly(dC) in this solution is intermediate between A- and B-DNA forms, or there is an equilibrium between the A- and B-DNA forms with approximately equal concentrations of both. Lesnik et al.', prefer the latter interpretation. In higher 3M NaCl solutions, poly(dG) . poly(dC) forms intermolecular aggregates with a corresponding shift to the ADNA form even with low polymer concentrations. It has been proposed that the G * C sequence is important for stabilizing the A-form of DNA.l5,I6 In light of the structural variability of poly(dG) . poly(dC) found by spectroscopic studies of aqueous solutions and the difference between solution and solid state results, it would be timely to examine this system by a completely independent method. Molecular dynamics (MD) has been shown to be useful for elucidating the structural features of oligonucleotides in vacuum and in ~ o l u t i o n . ' ~ - ~ ~ In this paper we report a model for poly(dG) . poly(dC) obtained by MD simulation of the hexamer (dG), . (dC), with explicitly included Na+ counterions and solvent.
METHOD The initial hexamer structure was constructed from standard B-DNA coordinates with a base-pair step height of 0.34 nm and a helix twist of 36.0°.21The 5' and 3' terminals ended in hydroxy group?. The united atom approximation was used for all carbon hydrogens while all other hydrogens were treated explicitly. Electrical neutrality was achieved by placing octahedrally hydrated Na ions, coordination number of 6,22 a t each of the phosphate positions such that each of the free phosphate oxygens formed a bifurcated H-bonding structure with one of the waters of hydration of a sodium ion. The resulting complex (solute, sodium ions, and their waters of hydration) was then submerged in a box of SPC water.23 The GROMOS (Groningen Molecular Simulation System)24 program PROBOX generated 1581 extra waters around the solute complex. The total number of waters around the minihelix plus sodium ion system was then reduced to 292 by discarding any waters beyond 0.5 nm of a solute atom. The net effect is that of a glove of water molecules over the surface of the solute. These 292 (the 60 Na+ waters are a subset) water +
MD SIMULATIONS OF (dG), . (dC),
molecules should give a reasonable description of the first hydration shell since approximately 300 water molecules hydrate one (B-DNA and counterions) turn of the helix with about 240 of these in the first hydration In the dynamics simulation the minihelix, 274 atoms, and the 10 associated Na' were grouped together as the solute. All water molecules were grouped in the solvent set. In order to eliminate the possibility of fraying of the helix during the simulation the center hydrogen-bonded pair of heavy atoms for the base pairs at the top and bottom of the minihelix were restrained to their initial interatomic distances using SHAKE. The potential energy function of the GROMOS programs was used without m~dification.~~ GROMOS partial charges were also used with a 1.0 dielectric constant. The cutoff range for nonbonded pair interactions was set at 0.8 nm and evaluated every MD step. The electrostatic interaction was calculated using a 1.8 nm cutoff and a pair list that was updated every 10 time steps. The complete solute-solvent system was relaxed with conjugate gradient energy minimization until the difference in energy between two successive steps was less than The initial velocities for the dynamics run were taken from a Maxwellian distribution. Constant temperature was maintained by weakly coupling the system to a 300 K thermal bath. No pressure effects were included. The MD run extended for 60 ps with a step size of 0.001 ps for the first 20 ps in order to equilibrate the system to the chosen temperature. A step size of 0.002 ps was used thereafter. No restrictions were placed on bond lengths,
