Sequence-dependent response of DNA to torsional stress: a potential biological regulation mechanism.

Torsional restraints on DNA change in time and space during the life of the cell and are an integral part of processes such as gene expression, DNA repair and packaging. The mechanical behavior of DNA under torsional stress has been studied on a mesoscopic scale, but little is known concerning its response at the level of individual base pairs and the effects of base pair composition. To answer this question, we have developed a geometrical restraint that can accurately control the total twist of a DNA segment during all-atom molecular dynamics simulations. By applying this restraint to four different DNA oligomers, we are able to show that DNA responds to both under- and overtwisting in a very heterogeneous manner. Certain base pair steps, in specific sequence environments, are able to absorb most of the torsional stress, leaving other steps close to their relaxed conformation. This heterogeneity also affects the local torsional modulus of DNA. These findings suggest that modifying torsional stress on DNA could act as a modulator for protein binding via the heterogeneous changes in local DNA structure.

INTRODUCTION

Understanding how the underlying organization of genome contributes to biological regulation is an important question. One important element in this organization is linked to the torsional strain that is imposed on the DNA double helix by many essential biological processes. This strain leads to supercoiling, namely axial bending combined with over- or undertwisting, and potentially to the formation of interwound plectonemic structures. One example of such a process is transcription where immobilized RNA polymerase within transcription factories forces DNA to rotate around its axis as the double helix threads through the transcription machinery (1,2). This leads to undertwisting upstream DNA (negative supercoiling) and overtwisting downstream DNA (positive supercoiling). These changes in supercoiling can then propagate along a chromatin fiber, with ranges and speeds of propagation that are dependent on the underlying nucleotide sequence (3). Modifying supercoiling will impact chromatin structure at many scales, locally modifying DNA conformation (notably, bending and twisting), stabilizing or destabilizing nucleosome core particles (4), and changing higher-order chromatin structures, with a resultant impact on protein–DNA interactions (5).

In vivo studies, using various DNA-mapping probes, show that the overall bulk of genomic DNA in mammalian cells is torsionally relaxed, while supercoiling exists at some loci (6,7). More recent investigations by Naughton et al. showed that genomes are organized into supercoiling domains (3). These domains are delimited by CG/AT boundaries, but are constantly remodeled by cellular machines. The CG-rich domains are, on average, significantly negatively supercoiled, highly transcribed, and enriched in open chromatin fibers, transcription initiation sites, RNA polymerase and topoisomerase I binding sites. In view of these observations, a comprehensive understanding of the sequence-dependent energetics and mechanics of DNA under torsional strain is important.

Despite being a crucial factor for chromosomal structure and function, little is known about the local structure of torsionally constrained DNA. Experimentally, supercoiling is studied using ensemble-average techniques, such as electron microscopy or electrophoretic separation (8,9), or single-molecule manipulation techniques, such as magnetic tweezers (10,11). While the first group of techniques probe overall supercoiling-mediated changes in compaction and overall geometry, the second category can provide dynamic information, for example, by monitoring supercoiling-diffusion, or the behavior of nucleosome arrays (12–14). These techniques provide valuable insights into the mechanical properties of DNA on the mesoscale, but little information on sequence-dependent structural rearrangements that may have important biological consequences, not least for the formation of regulatory protein–DNA complexes.

The mechanics of DNA supercoiling has also been addressed by various theoretical approaches. A number of phenomenological models of different complexity have been developed, including the worm-like-chain and rod-like-chain approaches (15–18) that describe DNA as an isotropic elastic polymer with constant persistence length, diameter and charge density. However, these models are limited to a description of the overall mechanical response of sufficiently long DNA molecules, typically hundreds to thousands of base pairs (bp) under moderate torsional stress (10,19).

Moving toward an atomistic view of DNA deformations, Olson et al. addressed the local plasticity of the double helix by analyzing structural fluctuations in crystal structures of protein–DNA complexes (20) and derived a set of empirical energy functions for deformations of the double helix. More recently, progress in all-atom molecular dynamics (MD), and in nucleic acids force fields, has enabled more detailed studies of torsionally stressed DNA. Using a simple twist restraint and relatively short MD sampling, Kannan et al. were able to generate free energy profiles with respect to twist and found evidence for sequence-dependent response to torsional stress, with pyrimidine(Y)-purine(R) (also denoted YpR) steps showing the highest response to under- and overtwisting (21). Pettitt and co-workers (22) also exposed the double helix to extreme torsional stresses, while restricting bending, and reported spontaneous sequence-dependent base flipping upon underwinding and the formation of Pauling-like DNA structures (23) upon overwinding. Lastly, Harris and co-workers explored the conformational space of supercoiled DNA minicircles and also observed deformations of the double helix, including the formation of single-stranded bubbles and ‘wrinkles’, induced by undertwisting (24). These studies provide valuable insights, but have not addressed sequence-dependent effects in detail. Other authors have used coarse-grain DNA models including sequence effects, but these assume that base pair steps deform harmonically and thus do not explicitly treat transitions between local conformational substates (25,26).

That such substates are important became clear during recent and comprehensive studies of the sequence-dependent conformational dynamics of DNA by the laboratories belonging to the ABC consortium (27,28). These studies have shown that certain dinucleotide steps are bimodal and can spontaneously adopt two different conformations. Helical twist is one of the variables involved, and it was found that YpR steps are particularly prone to twist bimodality (with substates typically separated by roughly 20° in twist). In addition, some purine(R)-purine(R) (denoted RpR) steps also show twist bimodality, but with a somewhat smaller range of twists. However, it was also noted that the sequences flanking such steps can significantly modify the degree of bimodality and also the preference for a low or high twist state. It is clear that such behavior may play an important role in how DNA reacts to torsional strain and the aim of this study is to address this question.

In order to do this, we have developed a new structural constraint, based on the Curves+ conformational analysis methodology (29), which controls the total helical twist between two chosen base pairs, allowing both under- and overtwisting to be imposed. This restraint has been implemented within PLUMED free energy library (30) and can thus be used in conjunction with a variety of standard all-atom MD software packages. While determining the total twist of a chosen segment of DNA, our restraint has no influence on other helical variables and it does not control how a segment will react locally to torsional stress. While a single restraint is typically applicable to a few turns of the double helix, several restraints can be chained together to study twisting in longer DNA fragments. The restraint can be applied to DNA alone and or to DNA complexes with other molecules.

As a preliminary test of how DNA responds to torsional stress in a sequence-dependent manner, we present results on four 17 bp oligomers (see Table 1), whose central segments contain repeat sequences involving the tetranucleotides ACGT, ACGA, CCGA and AGCT. Using the latest Parmbsc1 force field (31) coupled with our new restraint, we now investigate how these oligomers respond to over- and undertwisting. The results show that torsionally stressed DNA is highly heterogeneous. Sequence-dependence plays a significant role and leads to the notion of ‘twist capacitor’ dinucleotides that, in a given flanking sequence-environment can locally store and release torsional stress.

Table 1.

DNA oligomers studied, showing the four-letter code used to refer to them and the presence of potentially bimodal YpR and RpR steps

Oligomer	Sequence (5′→3′)	Restrained YpR steps	Restrained RpR steps
ACGT	GCACGTACGTACGTAGC	RCpGY, RTpAY, (RTpAR)
ACGA	GCACGAACGAACGAAGC	RCpGR	RApAY,(RApAR)
CCGA	GCCCGACCGACCGACGC	YCpGR	YGpGY
AGCT	GCAGCTAGCTAGCTAGC	YTpAR	YApGY

The twist restraint is applied to the bold and underlined segment of each oligomer. YpR and RpR steps in brackets occur at the end of the restrained segment and have different flanking base pairs.

MATERIALS AND METHODS

Restraining the helical twist

The helical twist of DNA should be measured around the helical axis, which cannot always be assumed to be a straight line. To take this problem into account we follow the general approach adopted by Curves+ (29). This requires defining local helical axes U and U at either end of the segment to be restrained. In order to obtain the U vectors, we begin by creating a local frame for each nucleotide b whose three unit vectors are derived from atoms belonging to each base: (i) the bond vectors N9-C1′ for purines, or N1–C1′ for pyrimidines; (ii) the cross product between this vector and the difference vector C8-C1′ for purines, or C6-C1′ for pyrimidines; (iii) the cross product between the two former vectors. The rotation matrix between two successive base frames can be obtained as Q = b is a 3 × 3 rotation matrix associated with three eigenvectors and three eigenvalues. The first eigenvector defines the U screw axis linking the two nucleotide frames by a rotation θ12 around the axis and a translation λ12 along the axis, where Cos θ12 = [trace(Q) - 1]/2.

The local helical axis U is then obtained by averaging the inter-base screw axes U and U in each strand of the double helix (i.e. four unit vectors in total). U is similarly the average of U and U within each strand. By combining U (and similarly U) with the unit vector along the N1-N9 axis of base pair i (oriented in the direction from the 5′-3′ to the 3′-5′ strand), we are able to create an axis frame B (and similarly B) associated with each end of the restrained DNA fragment. Solving the equation Q = B provides the total rotation between the two axis frames θij in the same manner described for the inter-base frames above. Finally, the total twist (the component of the total rotation corresponding to rotation around the helical axis) between B and B is given by twist = θij(U) where U is now the screw vector relating the B and B frames.

During the MD simulation, we can restrain the twist between base pairs i and j to the desired value twist using a simple quadratic function, k(twist-twist)2. Following initial trials, the force constant k was chosen as 0.06 kcal mol–1 degrees–2, the smallest value which was able to achieve the desired twist. Note that since θij is only defined in the range ±180°, the twist restraint can only be applied to oligomeric fragments of DNA, typically a few turns, where, in addition, the limited curvature implies that the restrained twist will remain close to that measured with the curvilinear axis derived by Curves+. Initial trials showed this to be true to within a few tenths of a degree. If it is necessary to control the total twist of longer DNA fragments, several restraints acting on short segments can be chained together.

The new twist restraint was implemented within version 2.2.1 of the PLUMED free energy library environment (30). This allows it to be used in combination with a variety of MD software packages. The C++ code and a user guide for the twisting constraint can be found here: https://github.com/annareym/PLUMED_DNA-Twist.

Over- and undertwisting DNA oligomers

For each of the 17-mer oligomers studied, the twisting restraint was applied between the 3rd and 15th bp (i.e. 12 bp steps, see Table 1). This avoids potential problems linked to deformations within the weaker base pairs close to the oligomer termini. The total twist between these base pairs was increased or decreased in 6° steps (corresponding to an average change in twist of 0.5° per base pair step. Following MD simulation of each step (see below), convergence was improved by using the final structure under each twist restraint as the starting structure for the following step, working outward in the over- and undertwisting directions from the relaxed oligomer. Over- and undertwisting of each oligomer was limited to a maximum of ±5° with respect to 34.9°, the relaxed twist per base pair step averaged over the four oligomers studied. Note that this range is approximately double that sampled by torsionally unrestrained DNA at room temperature.

Molecular dynamics umbrella sampling simulations

The DNA oligomers studied were initially constructed in standard B-DNA conformations using the nucleic acid modeling program JUMNA (32). Unrestrained MD simulations followed by umbrella sampling simulations were both performed using the GROMACS MD software package, version 5.1 (33). Simulations were carried out using the latest AMBER all-atom nucleic acid force field Parmbsc1 (31). Note that this parameter set avoids artefacts that, after a few tens of nanoseconds of simulation, caused a build-up of unusual backbone conformations involving the α (P-O5′) and γ (C5′-C4′) dihedrals and resulted in a steady loss of helical twist. The parmbsc1 force field has been the subject of extensive testing and has been shown to accurately reproduce a wide range of experimental data on DNA (31).

Each oligomer was first neutralized with 32 K+ counterions and solvated with TIP3P water molecules (34) corresponding to a solvent layer of 10 Å, contained within a cubic cell under periodic boundary conditions. Additional K+ and Cl− ions were then added to achieve a physiological salt concentration of 150 mM. The conformation of each oligomer was initially energy minimized with 5000 steps of steepest descent, followed by a 200 ps simulation at constant volume, while raising the temperature to 300 K. Simulations were then carried out at constant pressure and temperature (1 atm., 300 K) using a weak-coupling thermostat (35) with a 0.2 ps coupling constant and an isotropic Parrinello-Rahman barostat (36) with a 2 ps coupling constant. Simulations used a 2 fs time step. Bonds involving hydrogen atoms were constrained with the LINCS algorithm (37) and the non-bonded pair list was updated every 20 fs with the group scheme (38). Electrostatic forces were evaluated with particle-mesh Ewald (39) with a real-space cutoff of 10 Å. The van der Waals forces were truncated at 10 Å and long-range corrections were added. Center of mass movement was removed every 0.2 ps to avoid the building up translational kinetic energy (40).

Since restraints are necessary to study a sufficiently large range of helical twists, it is not possible to directly obtain free energy curves using the inverse Boltzmann procedure that determines the energy of a given state on the basis of its probability of occurrence. It is first necessary to correct the probabilities for the impact of the restraining potential. This can be done in an iterative manner using the so-called Weighted Histogram Analysis Method (WHAM) method (41). We used the version of this approach implemented in PLUMED (30) to calculate the corresponding potential of mean force (PMF) with respect to DNA twisting. Following a 300 ns simulation on each unrestrained oligomer, umbrella sampling was carried out with 100 ns of sampling per window. The initial 20 ns of the sampling time was discarded as equilibration and WHAM analysis was applied to two equal 40 ns blocks of data to test convergence. Overall the deviations between the PMF profiles representing the two sampling blocks were negligible (with maximum deviation of 0.02 kcal mol−1), with the exception of the AGCT oligomer for which sampling was extended to 150 ns per window to reach similar level of convergence. The total simulation time therefore was 2.4 μs for ACGT, ACGA and CCGA and 3.4 μs for AGCT.

Conformational analysis

The conformational analysis of the recorded trajectories was performed in several stages, beginning with pre-processing using the cpptraj program (42) of the AMBERTOOLS 14 package. For each MD snapshot, extracted at 1 ps intervals, DNA helical parameters, backbone torsional angles, and groove geometry parameters were analyzed using Curves+ and Canal (29) (available at http://curvesplus.bsc.es), providing complete, time-dependent information on the response of the DNA oligomers to the imposed torsional stress.

RESULTS AND DISCUSSION

The four oligomers studied in this work are listed in Table 1. These oligomers were chosen so as to include pyrimidine-p-purine (YpR) and purine-p-purine (RpR) steps that have already been shown to be able to exist in high and low twist states (28,43,44). Table 1 lists the presence of these steps within the segment of each oligomer that will be subjected to twist restraints (namely, base pairs 3–15). We also indicate the purine/pyrimidine (R/Y) nature of the flanking base pairs. Note that the restrained segments also contain 7 out of the 10 distinct dinucleotide steps, including two out of the three YpR steps and all four RpR steps.

The first three oligomers contain the YpR step CpG. The first and the fourth oligomers contain another YpR step, TpA. The last three oligomers contain the RpR steps: ApA in ACGA, GpG in CCGA and ApG in AGCT. Based on microsecond-scale ABC and subsequent studies using the AMBER Parmbsc0 force field (28,45) the CpG step within ACGT would mainly be expected to prefer a high twist, that within ACGA should have an intermediate twist and that within CCGA should prefer a low twist. The TpA step in the AGCT oligomer (i.e. CTAG) should also prefer a low twist, while that in the ACGT oligomer (i.e. GTAC) should prefer a high twist. However, all these steps should be able to shift between lower and higher twist substates. The AA step in ACGA (i.e. GAAC), the GG step in CCGA (i.e. CGGT) and the AG step in AGCT (i.e. TAGC) should all prefer high twists. These observations are confirmed for the Parmbsc1 simulations on our oligomers, as shown by 300 ns unrestrained MD twist distributions shown in Figure 1 for the base pair steps belonging to the central tetranucleotides. the single exception is the TpA step in AGCT, which, although it has a clear low twist population, spends more time in the high twist state. It should however be pointed out that twist conformer distributions can be very slow to fully stabilize.

Figure 1.

Twist distributions of dinucleotide steps derived from 300 ns unrestrained MD trajectories of the four oligomers.

Sequence-dependent torsional modulus of DNA

Following 300 ns of unrestrained MD, the new twist restraint was used to modify the total twist of the 12 central base pair steps by ±5° with respect to the average twist per base pair of the restrained segments (34.9°). This corresponds to a change in superhelical density σ of approximately ±0.15, estimating this value for a straight DNA fragment on the basis of the average relaxed total twist. The base pairs of the restrained segments remained intact during both over- and undertwisting, although occasional opening angles beyond 30° was observed in the most undertwisted state (generally for less than 1% of the trajectory and only reaching 4 to 6% for 3 bp). No significant bending was observed, with the bending probability distributions showing maxima of at most 5°, with most values below 2°.

PMF plots, that describe the change in free energy as a function of the imposed twist change (shown in Figure 2), were obtained using the umbrella sampling protocol described above with 100 ns of sampling (150 ns for the AGCT oligomer) for each 6° change in total twist (corresponding to a 0.5° average twist change per base pair step). Although DNA is usually assumed to show a harmonic behavior with respect to small deformations, Figure 2 shows that the PMF curves are not truly symmetric about the minimum. We thus tried fitting the curves with both quadratic and cubic functions. In the case of ACGT and AGCT, the root mean square error (RMSE) was almost identical for quadratic and cubic fits (0.008 versus 0.007 kcal mol−1 for ACGT and 0.010 versus 0.008 kcal mol−1 for AGCT). However, ACGA and CCGA show somewhat more asymmetry around the minimum as evidenced by RMSE’s for quadratic and cubic fits of 0.015 versus 0.004 kcal mol−1 and 0.015 versus 0.007 kcal mol−1 respectively. In all cases, the curves show small irregularities around the minima where exchanges between several conformational substates occur frequently (see below). We have consequently derived the optimal average twist angles and the corresponding force constants from the fits to the overall PMF curves. In order to compare the results more easily, the curves are plotted as changes in the average twist per base pair step with the respect to the corresponding optimal twist value, which varies from 34.4° for CCGA to 35.2° for AGCT (see Table 2). The energies at the optimal twist values have also been shifted to zero. The results in Figure 2 show only a moderate sequence effect. On average, modifying a single base pair step by ±4° (which corresponds to change in superhelical density σ of approximately ±0.12) implies an energy cost of 0.45 kcal mol−1. We can also see from the figure that the asymmetry associated with ACGA and CCGA is due to a slightly higher resistance to overtwisting compared to undertwisting for these oligomers.

Figure 2.

PMF with respect to the average change of twist per base pair step with respect to the corresponding value for the relaxed oligomer.

Table 2.

Calculated average relaxed twists, torsional constants, torsional moduli and torsional persistence lengths for the oligomers studied.

Oligomer	<Relaxed Twist> (deg)	K (kcal mol−1/deg2)	C (pN nm2)	P (nm)
ACGT	34.7	0.054	420	102
ACGA	35.0	0.061	475	115
CCGA	34.4	0.064	495	120
AGCT	35.2	0.057	443	108

The force constants derived from the fitted curves (K) can be used obtain the torsional modulus C, using the isotropic rod equations, T = KΔθ = C Δθ/L, where T is the torque resulting from a change in twist Δθ over a length L. Note that we use a standard value of L = 0.34 nm in order to compare with other published results, although L is weakly sequence dependent and can also vary with the imposed torsion (see later). The torsional modulus can be used to derive the torsional persistence length P using the equation P = C/k (taking k at room temperature as 4.114 pN nm). The results are shown in Table 2. The resulting order of torsional stiffness is CCGA > ACGA > AGCT > ACGT. The calculated values are in the range of sequence-averaged experimental measurements using single molecule techniques (applying torsion while gently stretching the DNA preventing the formation of plectonemes) which yielded = 410 ± 30 pN nm2 (implying