Cong-Min Ji(祭聪敏), Yusong Tu(涂育松), and Yuan-Yan Wu(吴园燕)

College of Physics Science and Technology,Yangzhou University,Yangzhou 225009,China

Keywords: G-quadruplex DNA,hydration,diffusion,reorientation dynamics

G-quadruplexes (GQs) are guanine-rich, non-canonical nucleic acid structures that play fundamental roles in genomic stability and the regulation of gene expression and also are enriched in promoter sequences of growth regulatory genes.They have been a subject of great interest due to their suggested potential therapeutic applications.[1-4]Hydration of GQ DNA plays an essential role in both preserving the DNA structure and ensuring its proper biochemical function.[5-10]For example,Nagatoishiet al.reported that the enhanced hydration of two thymine sites in loops may underlie the mechanisms of binding between non-canonical DNA structures and proteins that play important roles in biological processes.[11]Milleret al.reported that the hydration layer is a major determinant of GQ stability and conformation of the human telomere 3＇sequence of d(AG3(TTAG3)3).Compact GQs appear to be stabilized by dehydration, which has relevance to intracellular conditions where the absolute water activity must be much lower than 55.5 M.[12]In recent years,the properties of biomolecular hydration shells have been extensively studied,with special attention devoted to their dynamics.Liet al.[13]have surveyed the available folded right-handed quadruplex crystal structures in the Protein Data Bank.Almost all report the presence of water molecules, which is unsurprising given the typical water content of at least 40%for an oligonucleotide crystal.[14]The dynamics of water in a hydration shell have remained highly controversial,and different perspectives have been inferred based on different methods.A new concept called ‘biological water’ has been proposed, which is a more rigid water structure with a pronounced slowing down of structural fluctuations and rotational motion with a concomitant increase in H-bond lifetime.[15,16]In contrast, both simulations at the molecular level and experimental work have challenged this concept and reported a moderate slowing down of water dynamics by the rather modest factors of two to three compared with the bulk.[17-19]

The dynamics of water in the biomolecular hydration shell clearly present a slowdown compared with the bulk,even though the magnitude of this slowdown continues to be debated.However,is this slowdown identical everywhere within the hydration layer or is the impact of some sites more significant than others? Due to the chemically and topologically heterogeneous of exposed biomolecular surfaces, the second option seems more likely.[20]For GQ DNA in particular, the chemical groups,including phosphate group with highly negatively charged in backbone and different types of aromatic rings in bases, are all of special interest as H-bond acceptors and donors; in addition, the shape of the exposed surface is quite complex including pockets inside loops,an inner channel and grooves.Identifying which sites produce greater and smaller perturbations requires a site-resolved mapping of the dynamics, and this is extremely challenging with all experimental techniques.In this work, we employed all-atom molecular dynamics(MD)simulations to study the hydration patterns of GQ DNA.The Drude oscillator model was used in MD simulation as a computationally efficient method for modeling electronic polarization in DNA ion solutions.Hydration structure was analyzed in terms of radial distribution functions and three-dimensional high-density hydration sites.Analysis of hydration dynamics focused on self-diffusion rates and orientation time correlation at different structural regions of GQ DNA.

2.1.MD system construction and force field

In general, G4 sequences are non-randomly distributed but mainly clustered in pivotal genomic regions, namely,telomeres, gene promoters and DNA replication origins, and their topology can be influenced by the nature and length of loops and their location in the sequence.[21,22]To describe the general hydration properties of GQ DNA, three GQ DNAs were analyzed in our study, including a human telomere GQ DNA (PDB ID 2gku),[23]a c-myc promoter GQ DNA (PDB ID 2mgn).[24]and a GQ DNA formed by a synthesized DNA oligonucleotide(PDB ID 1u64).[25]The sequence of the latter GQ DNA (1u64) loop is only thymidine, and the thymidine residues Thy4-Thy7 form a diagonal loop, whereas Thy15-Thy18 form an edge type loop (Fig.1).These are two very typical and concise loop types that are beneficial for exploring the general hydration properties of loops.Therefore,we mainly analyzed the hydration patterns around the last GQ DNA (1u64) in the main text, and discuss the other two molecules(2gku and 2mgn)in Appendix A.

Recently, the CHARMM[26,27]and AMBER[28]force fields (FFs) have been the most commonly used in studies of nucleic acids.The CHARMM additive all-atom FF for nucleic acids treats electrostatics within the framework of the fixed-atomic-charge approximation, where effective charges assigned to particles are independent of a system’s configuration.While this approach is very attractive from the viewpoint of computational efficiency, the underlying Coulomb electrostatic interaction potential does not allow for the system to respond to changes in the polarity of the environment via redistribution of the electron density.For such a complex biological system as polyanionic DNA immersed in an aqueous salt environment, the conformational behavior of which is determined to a significant extent by solvation effects and interactions with the surrounding ionic atmosphere, the omission of polarization effects may preclude a physically correct description of its conformational behavior.In the all-atom CHARMM polarizable FF[29]for DNA based on the classical Drude oscillator model,electronic polarizability to represent electronic induction is explicitly included in the FF using classical Drude oscillators, where an auxiliary(Drude) charged particle is attached to each polarizable atom by a harmonic spring.In MD simulations using the Drude model,charge redistribution as a response to the change in the local electrostatic field is approximated by updating self-consistently the positions of Drude particles.The SWM4-NDP polarizable water model,[30]based on classical Drude oscillators, is re-optimized for negatively charged Drude particles,and is calibrated to reproduce important properties of the neat liquid at room temperature and pressure: vaporization enthalpy,density,static dielectric constant,self-diffusion constant,the liquid shear viscosity and free energy of hydration.In this paper, the simulations were performed in Nanoscale Molecular Dynamics (NAMD) using a Drude polarizable FF and the SWM4-NDP polarizable water model.

Fig.1.Cartoon representation of the NMR structure of GQ DNA(PDB ID 1U64) with bound Na+ ions (gold).The guanine bases of the GQ core are colored by tetrad(1,red; 2,blue; 3,green).End 1: inner region of the loop outside tetrad 1(red);end 2,inner region of the loop outside tetrad 3(green).

The starting structure for GQ DNA was taken from the nuclear magnetic resonance (NMR) ensemble in PDB entry 1U64.[31]The first of the eight deposited models was used.The deposited structure lacks bound Na+expected in the tetrad core, so two Na+ions were added to the GQ structure using the CHARMM program.[32]The positions of the ions were assigned using the average coordinates of guanine base carbonyl oxygen (O6) atoms in consecutive tetrads, yielding symmetric Na+coordination.Guanines 2,8,12,and 20 comprise tetrad 1, guanines 1, 13, 9, and 21 comprise tetrad 2 and guanines 10, 14, 19, and 22 comprise tetrad 3.A cartoon representation of the structure with bound Na+ions is shown in Fig.1.A solvated system was initially prepared with the additive C36m nucleic acid FF.[27]The GQ was centered in a 70 °A×70 °A×70 °A cubic cell, which was filled with CHARMM-modified TIP3P water[33]and 0.15 mol·l-1NaCl, including neutralizing Na+counterions, which is the typical intracellular NaCl concentration in eukaryotes.Standard CHARMM ion parameters were applied to the NaCl.Energy minimization in NAMD was performed to relax the solvated system.

2.2.MD simulations

Equilibration of the systems was performed in NAMD[34]for 1 ns under an NPT ensemble while all non-hydrogen atoms of the GQ and the bound Na+ions were restrained and water molecules and bulk NaCl ions were free.The NPT ensemble was enforced using the Langevin thermostat method at 298 K,pressure was kept constant at 1 atm 1.01325×105Pa) and periodic boundary conditions were applied in all dimensions.The short-range van der Waals forces were smoothly switched to zero from 10 °A to 12 °A, and the particle mesh Ewald method[35]was used to calculate electrostatic interactions with a real-space cutoff of 12 °A and a Fourier grid spacing of 1 °A.All bonds containing hydrogen atoms were constrained using the SHAKE algorithm and water molecules were kept rigid with SETTLE,allowing an integration time step of 1 fs.

To prepare the Drude systems, first the final equilibrated additive configuration was converted to the Drude polarizable model[29]in CHARMM by adding Drude oscillators to all non-hydrogen atoms in the system.Lone pairs were also added and the TIP3P water molecules were converted to the polarizable SWM4-NDP model.[30]Secondly,the Drude-2017 FF.[29]was applied to the GQs and ions, the Drude oscillators were relaxed using steepest descent minimization and adopted-basis Newton-Raphson energy minimization,[36]followed by 1 ns NPT equilibration at 298 K and 1 atm pressure.Extended Lagrangian integration, implemented in NAMD as Langevin dynamics, was used to perform the equilibration simulation.[37]

The simulation was performed in NAMD using a Drude polarizable FF.Production simulations were carried out three times for 100 ns each.The same NPT ensemble was maintained as described above and temperature was maintained using an Andersen thermostat with a collision frequency of 1 ps-1.

The structure and dynamics of water around GQ DNA were analyzed from MD simulation of a GQ DNA salt solution.In the following, for different regions on the GQ DNA surface including the phosphate group, loop and groove, we will first describe the hydration structure and then discuss dynamics.

Water structure To describe the water structure around GQ DNA,we first analyze radial distribution functions of water at different regions on the DNA surface.Furthermore,we describe water structure through the three-dimensional density distribution.

Fig.2.The radial distribution functions of oxygen atoms with respect to GQ DNA.(a)Around phosphate group: the inset shows a snapshot of the water in the first hydration shell.(b) In the groove: the resultant curve shows the splitting of the first maximum into two maxima at 2.05 °A and 2.95 °A:the inset shows a snapshot of the water in groove.(c)The inside of the first loop(marked as end 1 in Fig.1): the resultant curve shows the splitting of the first maximum into two maxima at 2.25 °A and 2.85 °A.The inset shows a snapshot of a water in this loop.(d)Inside the other loop(marked as end 2 in Fig.1): the resultant curve shows the splitting of the first maximum into two maxima at 1.95 °A and 2.75 °A.The inset shows a snapshot of water in this loop.

Radial distribution functions The radial distribution functions (RDFs) were calculated for the distances between water oxygen atoms and the nearest heavy atoms of different regions of the GQ DNA and were built for every picosecond of the simulation trajectory.To find the most common features of water distribution around different regions of the GQ DNA,the mean RDFs were built by averaging the RDFs over three 100 ns simulation trajectories.As shown in Fig.2(a),the obtained mean RDF is characterized by two maxima revealing the shells of hydration.The first maximum is the most prominent and corresponds to the first hydration shell where oxygen atoms of water molecules are in direct contact with the phosphate groups.The second maximum is about two times lower than the first and corresponds to the second hydration shell of the phosphate group.There is also a third maximum.It is essentially lower than the first two and corresponds to water molecules that are almost bulk water but are still constrained by interaction with the phosphate group.In the case of groove,it is interesting that there is a ‘spine of water’ in the groove(Fig.2(b),inset),just like the minor groove of double-stranded DNA.The resultant RDF curve shows the splitting of the first maximum into two maxima at 2.05 °A and 2.95 °A(Fig.2(b)),and the position of the first maximum of the RDF is much closer than that of the phosphate group,which corresponds to hydration water close to the bottom of the groove.The second maximum corresponds to water in the‘spine of water’.In the case of the two loops above the two ends of the G-tetrad, the two RDF curves also show the splitting of the first maximum into two maxima at 2.25 °A and 2.85 °A for end 1 and 1.95 °A and 2.75 °A for end 2, indicating that there are two hydration patterns:the first peak corresponds to the water close to center of end of G-tetrad that strongly interacts with the O6 atoms of guanines and the second peak corresponds to the water inside the loop that interacts with the base of the loop.The RDFs indicate the different water geometries around the phosphate group, groove and loop and show a heterogeneous hydration pattern due to the rugged surface and heterogeneous electrostatic potential of GQ DNA.

Three-dimensional water densities The threedimensional water densities were calculated to further compare the hydration structure between different regions of GQ DNA (Fig.3).Most of the high-density sites appear to be homogeneous around the G-tetrad, including the phosphate group and groove.However, there are also some regions that are significantly different.One such region is the outside surface of the loop (marked with a dashed circle shown in the left of Fig.3) where the hydration inside the loop appears to be more rigid and involves a different water structure.On the other hand, there is an insular high-density site in the inner channel of the G-tetrad (marked with a dashed circle in the right of Fig.3).These results suggest that the distribution of water in the first hydration shell largely corresponds to the dynamic and surface structure of GQ DNA.

Fig.3.Three-dimensional high-density water sites around GQ DNA(pink).The densities were averaged over 100 ns from sampling of water around DNA after superposition to a common reference structure.Density contours are shown at a level of iso=4.0.The left and right figures show front and top views,respectively.

Water dynamics We will now turn to an analysis of dynamic features,in particular water diffusion and the reorientation slowdown factor.

Diffusion coefficients The diffusion coefficients[38]were calculated from the mean square displacement(MSD)of water oxygen atoms around different region of GQ DNA:

The diffusion coefficient is calculated from the MSD according to

The time evolution of the MSD indicates a slowdown in diffusion for all hydration water of GQ DNA with respect to bulk water (shown in Fig.4, left).Diffusion rates calculated from the long-time slope of MSD versus time show a sequence of the degree of diffusion slowdown: end 2＆gt;end 1＆gt;groove＆gt;phosphate group.The results seem to depend on the depth of hydration water buried in the structure of GQ DNA.

Fig.4.Mean square displacement (MSD) as a function of time for water around different regions of GQ DNA:phosphate group(black),groove(red),end 1(green)and end 2(blue).The resultant curves were built by averaging the MSD over all 100 ns simulation trajectories.

In normal diffusion, MSD is a strictly linear function of time.However, previous studies have suggested that anomalous diffusion may occur in DNA hydration.In anomalous diffusion,MSD obeys the following law:

wheredwis the anomalous diffusion exponent.Ifdw= 2,normal diffusion is recovered.To judge whether anomalous diffusion is occurring,plots of log[MSD(t)/t]versus logtare convenient (shown in Fig.4, right).For the loop region including end 1 and end 2, anomalous diffusion is present for all time scales(100 ps)where log[MSD(t)/t]varies as a function of logt.It appears that the degree of anomalous diffusion increases,but after～30 ps diffusion becomes normal for water around the phosphate group and groove.For water inside the loop, the anomalous diffusion is likely due to a cage effect where water molecules can move freely in a very short region but are more restricted from making long-range displacements.For water around the phosphate and groove,the change from anomalous to normal diffusion essentially depends on the time required to overcome the electrostatic attraction and structural limitations.

Reorientation dynamics We also characterized the reorientation dynamics of water at the GQ DNA interface.The average orientational relaxation of water initially in the DNA hydration shell exhibits a much slower decay than in the bulk(Fig.5(a)).Its pronounced non-exponential character reveals the presence of strong dynamical heterogeneity within the hydration shell.By calculating the reorientation next to each GQ DNA part separately, we found that the reorientation dynamics of water in the groove and loop exhibit similar characteristics,but a much slower decay is found for water around the phosphate group.These essentially correspond to electrostatic interaction between hydration water and the polarized groups at different parts of DNA, where phosphate groups have the strongest electrostatic potential.For a quantitative comparative analysis,the average slowdown factor of hydration water at different parts of DNA were calculated(Fig.5(b)).Hydration water in the groove and loop experience a similar slowdown,with an average slowdown factor of～3.5,and the water around the phosphate group experienced a much slower decay with an average slowdown factor of～7.5.Furthermore,we found that the ‘spine water’ in the groove experienced a slightly slower decay than other hydration water in the groove.

Fig.5.(a)Water orientation time correlation functionC(t)=〈d(t)·d(τ+t)〉τ for water dipoles d(t) as a function of time t from averages over 100 ns for a GQ DNA solution system in different surface regions of DNA.Water molecules initially within the phosphate group(black), groove(red), the end of the G4 stem hydration layer (green) and within the bulk (blue), respectively.(b) The reorientation slowdown factor of water within the first hydration shell of DNA relative to bulk water.

In this study, we investigated the physical properties of water at different regions of the GQ DNA surface including around the phosphate group,loop and groove.We found that the hydration properties were significantly dependent on the structure, dynamics and electrostatic potential distribution of DNA and showed highly heterogeneous hydration patterns.Through extensive calculations and detailed analyses of the physical properties of hydration at different regions of the GQ DNA surface,the following major conclusions were obtained:

(i) The RDFs indicating the different water geometries around the phosphate group,groove and loop show heterogeneous hydration patterns due to the rugged surface and heterogeneous electrostatic potential of GQ DNA.

(ii) Most of the high-density sites appear to be homogeneous around the G-tetrad,including the phosphate group and groove.However, there are also some regions that are significantly different,such as the‘spine of water’in the groove and the insular high-density site in the inner channel of the G-tetrad.

(iii) For the loop region including end 1 and end 2,anomalous diffusion is present for all time scales (100 ps).However, after～30 ps diffusion becomes normal for water around the phosphate group and groove.

(iv)The hydration shell around GQ DNA exhibits a much slower decay than in the bulk.The reorientation dynamics of water in the groove and loop exhibits similar characteristics,but there is a much slower decay for water around the phosphate group due to the strong electrostatic potential.Furthermore,‘spine water’in the groove experience a slightly slower decay than other hydration water in the groove.

Overall,the results show highly heterogeneous hydration patterns in both structure and dynamics.In this study,we only addressed the stability of GQs in 0.15 mol·l-1NaCl solution by MD simulations.Actually, the folding topologies of GQ structures critically depend on the sequences and the nature as well as the concentrations of metal ions present in solution,all of which deserve further study.

Appendix A

The hydration pattern around two other G-quadruplex(GQ)DNAs(PDB IDs 2gku and 2mgn)are discussed below.

Radial distribution functionsThe RDFs were calculated for the distances between water oxygen atoms and the nearest heavy atoms of different regions of the GQ DNA for every picosecond of the simulation trajectory to find the most common features of water distribution around different regions of the GQ DNA.As shown in Fig.A2, the obtained mean RDF is characterized by two maxima, revealing the hydration shells for the phosphate group.The first maximum is the most prominent, and corresponds to the first hydration shell

Fig.A2.The radial distribution functions of the oxygen atoms with respect to different regions of GQ DNA:phosphate group(black),groove(red),loop(green).Left,2gku;right,2mgn GQ DNA.

where oxygen atoms of water molecules are in direct contact with the phosphate groups.The second maximum is about two times lower than the first and corresponds to the second hydration shell of the phosphate group.There is also a third maximum.It is essentially lower than the first two and corresponds to water molecules that are almost bulk water but are still constrained by interaction with the phosphate group.The results for the two DNAs are similar.In the case of the groove,for 2mgn DNA the RDF curve shows the splitting of the first maximum into two maxima,and the position of the first maximum of the RDF is much closer than that of phosphate group,which corresponds to the hydration water close to the bottom of the groove; the second maximum corresponds to water in the‘spine of water’like 1u64 DNA in the main text.For 2gku DNA the second maximum is essentially lower than the first maximum.In the case of the loop and the end of G-tetrad,the RDF curves of the two DNAs show a large difference,indicating that the different hydration patterns correspond to different loop types.The RDFs indicate the different water geometries around phosphate group, groove and loop, and show a heterogeneous hydration pattern due to the rugged surface and heterogeneous electrostatic potential of GQ DNAs.

Diffusion coefficientsThe diffusion coefficients were calculated from mean square displacements (MSDs) of water oxygen atoms around different region of GQ DNA.The MSD time evolution indicates slowdown in diffusion for all hydration waters of the two GQ DNAs with respect to bulk water (shown in Fig.A3, left).Diffusion rates calculated from the long-time slope of MSD versus time show a sequence of the degree of diffusion slowdown: loop＆gt;groove＆gt;phosphate group.The results for 2gku and 2mgn DNA show a slight difference.To judge whether anomalous diffusion occurs, plots of log[MSD(t)/t] versus logtare convenient(shown in Fig.A3,right).Anomalous diffusion is present before 30 ps, where log[MSD(t)/t] decreases as a function of logtfor all three regions of 2gku DNA;however,after～30 ps diffusion becomes normal.For the loop region of 2mgn DNA,anomalous diffusion is present for all time scales(100 ps)due to a cage effect inside the loop.

Fig.A3.Mean square displacement(MSD)as a function of time for water around different regions of GQ DNA:phosphate group(black),groove(red),loop(green).(a)2gku GQ DNA,(b)2mgn GQ DNA.

Acknowledgment

Project supported by the National Natural Science Foundation of China(Grant Nos.11705160 and 11647074).