Twisting DNA by salt

Abstract The structure and properties of DNA depend on the environment, in particular the ion atmosphere. Here, we investigate how DNA twist -one of the central properties of DNA- changes with concentration and identity of the surrounding ions. To resolve how cations influence the twist, we combine single-molecule magnetic tweezer experiments and extensive all-atom molecular dynamics simulations. Two interconnected trends are observed for monovalent alkali and divalent alkaline earth cations. First, DNA twist increases monotonously with increasing concentration for all ions investigated. Second, for a given salt concentration, DNA twist strongly depends on cation identity. At 100 mM concentration, DNA twist increases as Na+ < K+ < Rb+ < Ba2+ < Li+ ≈ Cs+ < Sr2+ < Mg2+ < Ca2+. Our molecular dynamics simulations reveal that preferential binding of the cations to the DNA backbone or the nucleobases has opposing effects on DNA twist and provides the microscopic explanation of the observed ion specificity. However, the simulations also reveal shortcomings of existing force field parameters for Cs+ and Sr2+. The comprehensive view gained from our combined approach provides a foundation for understanding and predicting cation-induced structural changes both in nature and in DNA nanotechnology.


Molecular dynamics simulations
We performed independent simulations of the dsDNA for each ion type and salt concentration. Simulations included periodic boundary conditions and electrostatic interactions were treated using particle-mesh Ewald summation S1 with tin foil boundary conditions. We used a 2 fs time step with constraints on the hydrogen bond atoms using the LINCS algorithm. S2 The long-range electrostatic interactions were treated with cubic interpolation and a Fourier space grid of 0.12 nm. Lennard Jones interactions and close Coulomb real space interactions were cut-off at 1.2 nm. Errors from the truncation of LJ interactions were accounted for by long-range dispersion correction for energy and pressure.
Simulation protocol: The dsDNA structure was placed in an orthorhombic dodecahedral box, assuring a minimal distance of 1.5 nm to the edge. The simulation box was filled with 71011 TIP3P water molecules. S3 Subsequently, the number of ions required for a given concentration was estimated based on the number of water molecules in the simulation box plus the ions necessary to neutralize the charge of the DNA (see Tables S1 and S2). Water molecules were randomly replaced by ions to obtain a neutral system with the desired salt bulk concentration.
A pre-equilibration protocol, consisting of energy minimization, NVT equilibration and NPT equilibration, was completed prior to production runs. During the pre-equilibration process, the heavy atoms of the nucleic acids were constraint with a soft harmonic potential with a force constant of 1000 kJ/(mol nm 2 ) to allow the equilibration of the solvent around the nucleic acid. Energy minimization used the steepest descent algorithm with a maximum of 50000 steps. Later, we employed 1 ns NVT and 1 ns NPT simulations to further equilibrate the system while keeping the positions restraints.
Finally, unrestrained runs were performed in the NPT ensemble. Trajectories were 3 µs long for monovalent cations and 5 µs long for divalent cations. NPT simulations used S-2 the isotropic Parrinello-Rahman barostat S4 and the velocity rescaling thermostat with a stochastic term. S5 Analysis of trajectories: Helical structural parameters were analyzed with the broadly used software tools 3DNA S6 and do_x3dna, S7 complemented with in-house scripts. Particularly, we obtained the helical-twist, helical-rise, the base-pair (bp-)twist, and bp-rise, the radius r, the sugar pucker P . r was defined as the mean distance of the phosphorous atoms to the helical axis. Finally, the helical crookedness β, previously introduced, S8 is measured from the ratio between the sum of helical rises h and the sum of bp-rises d (cos β = h/d).
For the analysis of the helical parameters, the last three bases at each end were not considered. The reported error values correspond to the standard error of the mean (SEM) after block averaging. The cations distributions were analyzed using the canion tool of the software Curves+, S9 which calculates the concentration profiles in untwisted helicoidal coordinates, allowing a detailed view of the interaction of the ions at specific volumes in the DNA. The final bulk concentration was obtained from canion (see Table S3). Complementary, gromaρs was used to obtain the three-dimensional distribution of the ions, S10 and the results were visualized with pymol. S11

S-3
Definitions of the twist in the simulations  Figure S1: Comparison of ∆Tw obtained from three different definitions: (i) the sum of the bp-twist (white symbols), (ii) the sum of the helical-twist (filled symbols), and (iii) the end-to-end twist definition as described by Kriegel et al. S12 (points). Empty symbols are the experimental results. The dotted lines are the fitting of the experimental results to the square root of the cation concentration. To compare the different methods, the dashed line and the transparent area corresponding to the fitting of the helical-twist to the square root of concentration and its error are plotted.
To compare the structural changes of the DNA with the changes observed in experiments, we calculated the change in twist using three different definitions, as previously discussed: S12 (i) the sum of the bp-twist, that is, the sum of local twist of a dinucleotide step. (ii) the sum of the helical-twist. It corresponds to the sum of the rotation angles about the helical axis that brings successive base pairs into coincidence. (iii) the end-to-end twist definition described by Kriegel et al.. S12 In this definition, a base-fixed coordinate reference frame S13 is assigned to both ends of the helix using quaternions averaging S14 over three base pairs at each end.
The total twist corresponds to the rotation of the X-Y plane of the reference frames. S12 For each definition, we obtained the change in the twist for every recorded snapshot.
The three methods yielded similar relative changes for the change in twist ∆Tw, as shown in Figure S1. In the following, we report the end-to-end twist because it has the property of being invariant with respect to initial constant rotation offsets about the z-axis, S12 and it reproduces the closest the experimental setup where the initial torsional offset of the beads is not known. For a detailed discussion, see Ref. S12

S-4
Influence of duplex length on the calculated change in twist Furthermore, recent experimental S15 and computational evidence S16,S17 have shown that the sequence strongly influences the equilibrium conformation of DNA. In our magnetic tweezers experiments, one would expect this effect to average out as measurements included a mixed sequence with 7900 bps. In contrast, for simulations of 33 bp it is not clear, a priori, how the sequence and length of the simulated helix may affect the relative changes in twist. To assess this effect, we computed ∆Tw for different duplex lengths. Hereby, we subsequently removed the last bp at each end of the simulated helix ( Figure S2). Our results show similar changes up to 15 bp (see Figure S2), which indicates that sequence effects average out with the 33 bp sequence used in the simulations. compete" for occupancy in the ion atmosphere around the nucleic acid S18,S19 and where the resulting concentrations of ions in close proximity to the DNA is intermediate between the corresponding single ion conditions. Comparison to ion counting experiments S18 suggests that the effect on DNA twist is at least approximately proportional to fraction of charge neutralized by the respective ion species in the ion atmosphere. The DNA twist for multiple ions in solution was determined similar to the measurements for single ion species reported in the main text. Rotation curves where first recorded for 100 mM KCl as the reference conditions. Subsequently, rotation curves were recorded in the mixed buffer conditions and the change in DNA twist was determined from the centers of the rotation curves. Bars and error bars are the mean and standard deviation from measurements with 18-20 independent molecules for the mixed ion conditions.

S-8
Supplementary simulation results

Changes in major and minor grooves
As shown in Figure S6, the width of the minor groove changes with the ion type and concentration. Li + compacts the minor groove the most, followed by K + , Na + , and Cs + . In contrast, the major groove width is essentially unaffected by Li + and Na + , but it is compressed by K + and Cs + ( Figure S6), reflecting the accumulation of the ions in the grooves as discussed in the main text.  Figure S6: (a,b) Changes in the minor groove and the major groove width depending on ion type and concentration. Error bars correspond to the SEM. (c) Snapshot of the simulated DNA structure indicating the binding volumes. (d) Snapshot of s segment of DNA with a detail of Cs + bound at the major groove (atoms N7 and O6) pulling two consecutive bases together.

S-9
Helical parameters as function of ion type and concentration  Figure S7: Changes in DNA helical properties as function of monovalent ion type and concentration. (a-d) Changes in the helical rise ∆h, the bp-rise ∆d, the helical radius ∆r, and crookedness ∆β, respectively, with respect to 100 mM of KCl. Dashed lines correspond to fittings proportional to the square root of the cation concentration.

S-10
Binding for K + and Na + as a function of DNA sequence 0 20 40 Occupancy(%) backbone Figure S8: Time averaged occupancy of the main binding sites at the backbone (A) and the nucleobase (atoms N7) of the DNA as a function of the sequence (B). The comparison shows 100 mM of KCl and 100 mM of NaCl. A binding site was considered occupied when a cation was within a cutoff of 3.5Å of the binding site. In agreement with X-ray experiments, S20 our simulations show a preferential binding of K + to the G-tracks of the helix.

S-11
DNA backbone sub-states in the simulations As a possible source of change in twist, we evaluate the DNA backbone sub-states, which have been related to changes in the local twist. S17 In particular, we investigated the substates BI and BII defined from the torsion angles ε and ζ as BI (ε−ζ < 0) and BII (ε−ζ > 0).

S-12
Changes in twist for divalent cations

S-13
Helical properties for divalent cations  Figure S11: Changes of the characteristic DNA properties for different divalent cations: helical rise h, bp-rise d, radius r, and crookedness β for Mg 2+ , Ca 2+ , and Ba 2+ . All results are relative to 100 mM KCl. For Ca 2+ , the data are for the affinity optimized parameter set from Ref. S22 In all cases, the lines are a fit proportional to the natural logarithm of the concentration.

S-14
Change in twist using TIP3P and SPC/E water models We performed a set of four control simulations using the SPC/E S23 water model to observe the effect of the water model on the change of twist. Namely: 100, 500, and 1000 mM of KCl and 1000 mM of NaCl. The SPC/E simulations started from the last frame of the corresponding TIP3P simulation and each calculation was 200 ns long. All other parameters were identical to the TIP3P simulations. For SPC/E simulations, we performed the twisting analysis and obtained the relative change in a twist ( Figure S12).
The SPC/E water model reproduces the increase of twist with concentration for K + (Figure S12A,B) and the relative larger twist of k B T relative to Na + at 1000 mM qualitatively.
However, the change in twist predicted by SPC/E is larger than the change predicted by TIP3P. Even though SPC/E water reproduces the physical properties of water better compared to TIP3P, the agreement with the experimental DNA twist is worse for SPC/E. Note, that the less convincing agreement for DNA twist in the SPC/E simulations is likely caused by the fact that the ion force fields were optimized in combination with TIP3P water (and not SPC/E). To improve the performance in SPC/E, a more systematic comparison with several ion force fields would be required which is beyond the scope of our current work.
S-15  Figure S13: Correlations between helical properties obtained for TIP3P and SPC/E water models. The two-dimensional histograms were calculated for the central 27 bp and for the total simulation time, i.e., 3000 ns for TIP3P and 200 ns for SPC/E. Plots in blue and orange are for K + and Na + , respectively. The SPC/E water model exhibits similar trends for the correlation between helical properties as the TIP3P model.    Table S4: Force field parameters for the metal cations used in this work. The charge and the Lennard-Jones parameters (σ and ) are listed. Values were take from Refs. S21 λ iR σ is the optimized scaling factor for the ion-phosphate oxygen interaction as described in Ref. S22

Ion
Charge (