Elucidating the influence of linker histone variants on chromatosome dynamics and energetics

Abstract Linker histones are epigenetic regulators that bind to nucleosomes and alter chromatin structures and dynamics. Biophysical studies have revealed two binding modes in the linker histone/nucleosome complex, the chromatosome, where the linker histone is either centered on or askew from the dyad axis. Each has been posited to have distinct effects on chromatin, however the molecular and thermodynamic mechanisms that drive them and their dependence on linker histone compositions remain poorly understood. We present molecular dynamics simulations of chromatosomes with the globular domain of two linker histone variants, generic H1 (genGH1) and H1.0 (GH1.0), to determine how their differences influence chromatosome structures, energetics and dynamics. Results show that both unbound linker histones adopt a single compact conformation. Upon binding, DNA flexibility is reduced, resulting in increased chromatosome compaction. While both variants enthalpically favor on-dyad binding, energetic benefits are significantly higher for GH1.0, suggesting that GH1.0 is more capable than genGH1 of overcoming the large entropic reduction required for on-dyad binding which helps rationalize experiments that have consistently demonstrated GH1.0 in on-dyad states but that show genGH1 in both locations. These simulations highlight the thermodynamic basis for different linker histone binding motifs, and details their physical and chemical effects on chromatosomes.


INTRODUCTION
In eukaryotes, chromosomes serve as the primary storage medium of genomic information within an organism and consist predominantly of organized, long condensed fibers of DNA and structural proteins (1). These fibers are made of compacted repeating arrays of DNA-protein complexes collectively known as chromatin (2,3). Despite being tightly condensed, chromatin still allows for enzyme induced replication, repair and transcription (4)(5)(6). The basic building block of chromatin fibers is the nucleosome core particle (NCP) which is comprised of 147 base pairs of DNA wrapped around an octameric core of histone proteins that are built from duplicates of four histones: H2A, H2B, H3 and H4 (1,(7)(8). These histones bind to one another to form H2A-H2B and H3-H4 dimers, while the H3-H4 dimers associate into a tetramer. This tetramer then combines with the H2A-H2B dimers to form the octameric core (9,10).
Linker histones primarily bind to the nucleosome in two states. In the 'on-dyad' location (33)(34)(35) the linker histone is centered on the dyad axis (Figure 1 (top)), whereas in the "off-dyad" configuration the histone binds in a DNA groove off the dyad axis (36)(37)(38) (Figure 1 (bottom)). The model studied here has the on-dyad linker histone interacting with the DNA minor groove, while in the off-dyad state it is adjacent to the DNA major groove at the +0.5 superhelical location, roughly three to seven base pairs from the dyad axis. Variations in the linker histone binding mode may result in differences in the mechanical stability and overall packaging within the greater chromatin architecture, which would naturally affect the accessibility of DNA in nuclear processes (39)(40)(41). This effect was demonstrated in recent coarse-grained simulations by Perišić et al. where they found the off-dyad binding mode, as observed in the experimentally generated cryo-EM image of polynucleosomal arrays (37), is a better chromatin condenser than other binding modes (41). Additionally, it was suggested that on-dyad binding is relatively poor for compaction, even compared to systems with hybrid binding modes, suggesting it plays a role in chromatin transcriptional accessibility and dynamic architecture. In a more recent cryo-EM and crystallography study, Garcia-Saez et al. showed that the on-dyad state can create alternative compact chromatin conformations, and that shifts in the ionic conditions can induce untwisting to reveal a ladder-like structures (42). Altogether, these works suggest that fluctuating linker histone binding modalities can lead to different levels and structures of chromatin compaction in cooperation with varying linker DNA lengths between nucleosomes (43)(44)(45).
Several factors contribute to the thermodynamic preference for on-versus off-dyad binding in linker histones.
Recent work by Zhou et al. examined the binding modes of the wild-type linker histone globular domains of Drosohila melanogaster generic H1 (genGH1), Gallus gallus H1.0 (GH1.0) and an H1.0 pentamutant (36). Paramagnetic relaxation enhancement (PRE) experiments showed that GH1.0 binds on-dyad and that genGH1 binds off-dyad, but also that a small number of mutations are able to shift the equilibrium of the GH1.0 binding state from on-to off-dyad (36). This suggests that the thermodynamic balance between these states is finely tuned by specific linker histone/nucleosome contacts. These studies are supported by cryo-EM experiments of condensed nucleosome arrays which suggested an off-dyad generic H1 binding mode (37). In contrast, in cryo-EM experiments of single nucleosomes, generic H1 has been observed in the on-dyad binding state (34), while there is evidence to suggest that off-dyad binding in the cryo-EM map of the 30-nm fiber may be a result of cross-linker effects (46). Taken together, these experiments paint the picture that linker histones likely bind in an ensemble of on-and off-dyad states, and that the balance of these two conformations is dictated by several factors including the linker histone primary sequence, the chromatosome's stereochemical environment and greater chromatin Nucleic Acids Research, 2020, Vol. 48, No. 7 3593 architecture (47). Note that the linker histone naming convention used here is consistent with that found in the Histone Database and introduced by Talbert et al. and may differ from that found in previous studies (48,49).
The contrasting influence of linker histone variants on the chromatosome structure and energetics, and the extent to which they affect greater chromatin dynamics, remains unclear (50)(51)(52). Using Brownian dynamic docking simulations,Öztürk et al. found that GH1.0 displays a range of conformational flexibility and affects the overall chromatosome dynamics, including the linker DNA (53). With similar techniques they later showed that even slightly varied linker histone sequences, including point mutations and posttranslational modifications, can significantly affect the chromatosome structure (47,54). Moreover, with accelerated molecular dynamics (MD) simulations they found that the GH1.0 ␤-sheet loop (␤-loop), which has both an open and closed-state in the crystal structure (55), favors the closed-state in solution, although the open-state may still be populated. This is in line with the closed-state being the only conformation observed in chromatosome crystal structures (33)(34)46). However, there is evidence suggesting that linker histones may exist in alternative conformations (56,57) and binding orientations (56,(58)(59)(60).
Despite these and many more (61) excellent experimental and computational studies, several questions remain concerning linker histones and their nucleosome binding. For example, to what extent does linker histone plasticity affect its function? What are the effects of on-and offdyad binding on chromatosome dynamics? How do specific thermodynamic forces influence the on-vs off-dyad binding equilibrium? How are these properties influenced by different linker histone variants? To address these questions, we have performed a series of conventional and free energy MD simulations of chromatosomes containing the D. melanogaster generic globular domain of H1 (genGH1) and G. gallus globular domain of H1.0 (GH1.0, which has previously been referred to as H5) bound in both onand off-dyad states. Building off work from the Bai and Wade groups, this work was focused on the globular domain of each linker histone. (36,38,47,(53)(54)56). Results suggest that both genGH1 and GH1.0 readily adopt a single compact configuration. Furthermore, in the off-dyad state linker histones display increased localized sampling while modestly altering the linker DNA dynamics, while linker histones have highly stable binding in the on-dyad state while significantly restricting DNA motions. Energetic analyses shows that the equilibrium between on-and offdyad binding is the result of a balance between Van der Waals and electrostatic interactions that is dictated by the linker histone variant type. Together, these results suggest that, regardless of the variant, on-dyad binding is enthalpically stabilized whereas off-dyad binding is relatively more entropically stabilized. Furthermore, when in on-and offdyad conformations, different linker histone variants have similar effects on chromatosome structures and dynamics, and that the role of linker histone modifications is likely to shift the relative populations between these binding states. The in vitro ensemble of binding modes, and therefore the greater structure of linker histone containing chromatin fibers, is therefore dictated by competing thermodynamic forces which are likely influenced by a myriad of structural and environmental factors in vivo.

System construction
Core histones were modelled based on the 1KX5 crystal structure (resolution 1.94Å (62)). The asymmetric Widom 601 DNA (63,64) was taken from the 4QLC crystal structure, which has a lower resolution (3.50Å) but both DNA and GH1.0 in an on-dyad conformation (33). Missing residues and nucleotides were added using Modeller via the Chimera graphical user interface (65,66). Linker histone coordinates from the 4QLC structure were used for GH1.0 ondyad simulations, whereas for genGH1 on-dyad simulations the GH1.0 primary sequence was mutated to the genGH1 sequence ( Figure 1). The completed on-dyad GH1 had an RMSD of 0.33Å relative to the recently published crystallographic genGH1 structure (PDB ID: 5NL0, resolution: 5.4Å) (34). For simulations of the nucleosome, the linker histone was deleted.
Off-dyad binding models were based on a combination of manual placement, rigid docking, and flexible docking. First, the exit DNA was manually adjusted to allow space for the linker histone to be placed in an off-dyad binding mode. Rigid docking of genGH1 was then performed with the 12Å cryo-EM map as a guide using the Colores module of Situs (37,(67)(68), which was followed by flexible docking using internal coordinates normal mode analysis (iMOD) (69). To validate the linker histone placement, theoretical PRE intensity ratios were calculated and compared to experimental data (see 'Analyses Methods' section below). A model of off-dyad GH1.0 binding was constructed by superimposing and replacing genGH1 coordinates with GH1.0 coordinates.

Molecular dynamics simulations
All systems were prepared with tleap from the Amber-Tools16 (70) software package. Each system was solvated in a TIP3P water box extending at least 10Å from the solute (71,72). Using Joung-Cheatham ions (73,74), the solvent contained 150 mM NaCl, sodium cations to neutralize negative charges, and magnesium ions that replaced the manganese ions in the 1KX5 crystal structure. Only magnesium ions in the DNA grooves were included, whereas those located close to the the linker histones binding locations were excluded so as to not interfere with LH-DNA interactions. The AMBER14SB and BSC1 force fields were used for protein and DNA interactions (75,76). All simulations were performed using NAMD version 2.12 (77). A cutoff distance of 10.0Å with a switching function beginning at 8.0Å was used for non-bonded interactions, and long range electrostatics were treated with particle mesh Ewald calculations (78). For constant pressure calculations a modified NAMD version of the Nosé-Hoover barostat was used with a target pressure of 1.01325 bar while the Langevin thermostat was with a target temperature of 300K (79).
Systems were minimized twice for 5000 steps, first with a 10 kcal·mol -1 ·Å -2 harmonic restraint applied to the solute and then followed by no restraints. Using the Langevin thermostat, systems were then heated in the NVT ensemble from a temperature of 10K to 300K in 1K increments every 4 ps with a 10 kcal/mol kcal·mol -1 ·Å -2 solute restraint. This restraint was then reduced by 0.001 kcal/mol every 60 fs in the NPT ensemble. Equilibration runs with no restraints and a temperature of 300K were then performed. Simulations were conducted on seven systems: four chromatosomes (each containing genGH1 or GH1.0 in either the onand or off-dyad binding mode), one nucleosome and two isolated linker histones (genGH1 and GH1.0). Each simulation was run in triplicate for 250 ns using resources provided by the Extreme Science and Engineering Discovery Environment (XSEDE) (80).
Umbrella sampling. The 2 reaction coordinate was divided into 41 windows spaced every 2 • , which covered a range of 40-120 • . Seed structures for each window were selected from steered MD runs where 2 was adjusted at a rate of 5 deg/ns using a harmonic force constant of 0.1 kcal/mol from the closed-state in the 1HST structure (55). The angle used here, along with 1 , is detailed in Supplementary Figure S2 of Section S4 and Supplementary Table S3 of Section S3. Windows were run for 20 ns each with a force constant of 0.1 kcal/mol/deg 2 , totaling 1.64 s of simulation data, where the first 5 ns of simulation was removed from analysis for equilibration. MD simulation parameters in NAMD were the same as described above for simulations of linker histones in solution, with the umbrella potential provided by the colvars module (81). Trajectories were analyzed using the weighted histogram analysis method (WHAM) (82) with code from the Grossfield group (83). For the Monte Carlo bootstrapping error analysis, 100 trials were run for each distribution and the statistical inefficiency was calculated for each window to calculate the number of statistically independent data points in each window. Convergence of each PMF can be found in Supplementary Figure S3.

Analyses methods
Unless otherwise stated, all analyses were performed using cpptraj with the first 50 ns of the trajectory excluded for equilibration (70). Figures were created using Visual Molecular Dynamics (84). Protein-DNA contacts were defined between the heavy atoms of residues within 4.0Å of one another.
Estimated binding affinity. Binding affinities were estimated with an MM/GBSA analysis (molecular mechanics--generalized Born surface area) using the MMPBSA.py script from the AMBER16 software suite (85). The three trajectory approach was implemented by using the trajectories from chromatosome simulations as the complex, the nucleosome simulations as the receptor, and simulations of linker histones in solution as the ligands. Explicit solvent trajectories were stripped of all solvent molecules while using trajectory frames every 4 ps. The implicit solvent model, GBneck2 with the mbondi3 radii parameters were used as they have been shown to have good agreement with more expensive Poisson-Boltzmann calculations for protein/nucleic acid complexes (86). The salt concentration was set to 150 mM.
Clustering. Binding modes of genGH1 and GH1.0 were compared by calculating RMSD values of the helical ␣carbons in the linker histone with respect to the helical ␣carbons in the core histones. The RMSD analysis was limited to the helical ␣-carbons to reduce noise from the loops and intrinsically disordered tails. Based on the RMSD results, a cutoff of 2.0Å was chosen for the subsequent clustering analysis (Supplementary Figure S4). Clustering was performed using the hierarchical agglomerative approach implemented in the cluster module of cpptraj from the the AmberTools16 software package (87).
Linker DNA dynamics. The linker DNA in-and out-of nucleosomal plane motions were quantified to describe the linker DNA motions. To define the plane, the nucleosomal DNA was divided into four quadrants and the center of mass of the C1' atoms within the two quadrants located distal from the linker DNA were used for two points, while the third point was defined as the C1' center of mass of bases 83 and 250 which are located approximately on the dyad axis (see Supplementary Section S2 for details). The linker DNA vectors were defined as the C1' center of mass of the base pairs at the origin of the linker DNA (bases 20-315 and 148-187) and terminal base pairs (bases 1-334 and 167-168), respectively. The ␣-angles were defined as in-plane and the ␤-angles were defined as out-of-plane motions of this vector. Positive ␣-angles were defined as inward motions towards the dyad axis while positive ␤-angles were defined as outward motions away from the nucleosomal-plane. For reference, the angles shown in Figure 6 are positive.
The normalized mutual information (NMI) between angles was calculated to determine the correlations between each pair of angles (88). The NMI has the advantageous property that all values are in the range of 0-1 and are therefore easier to interpret than standard mutual information (MI). For detail discussion of the NMI see Supporting Data. The change in sampling of bound linker histones from the nucleosome were computed by the Kullback-Leibler (KL) divergence utilizing the nucleosome sampling distributions as the reference set (89): where Q(x) is the normalized reference distribution (nucleosomal linker DNA angles) and P(x) is the normalized data set (chromatosomal linker DNA angles).
Paramagnetic relaxation enhancement (PRE) intensity ratios. Theoretical PRE cross-peak intensity ratios were estimated based on distances between experimentally labeled core histone (H2A T119 and H3 K37) and linker histone methyl-terminated residues (36,38). Since our simulations did not include the MTSL probe, wild type residues were used for distance calculations, similar to Piana et al. (90). For the K37 probe, distances were measured from the lysine terminal nitrogen to the terminal methyl carbon atom of each respective residue, while the T119 distance was measured from the threonine terminal methyl carbon. Note, that each of these neutral residues have two distinct methyl groups and were referred to as methyl-a and methyl-b. Using these distances, the predicted intensity ratios where calculated using known equations and the average of the methyl-a and methyl-b values were used to compared to experimental values (36,38,62). The relationship between between the paramagnetic relaxation cross-peak intensity ratios and the interatomic distances is given by: where, R 2 is the intrinsic relaxation rate (inverse of the transverse time constant (T 2 )), t is the total evolution time of the transverse proton magnetization, R sp is the contribution to the relaxation caused by the paramagnetic probe, and I ox and I red are the peak intensities for the oxidized and reduced states, respectively. This last variable is what connects the cross peak equation to the interatomic distances, r. Distances are defined by the following (91,92): where, K is a constant that describes the spin properties of the MTSL label (1.23 × 10 −32 cm 6 s −2 ), h is the Larmour frequency and C is the apparent correlation time which is estimated from the molecular weight of the protein.

Unbound linker histones favor the closed-state
To quantify the dynamics of linker histones in solution, we measured two angles over our unbound genGH1 and GH1.0 simulations, 1  Therefore, 2 appears to be a stronger metric of ␤-loop structures, whereas 1 is subject to increased mobility throughout the ␤-sheet. To quantify the dynamics of the ␤-loop, umbrella sampling simulations were performed to determine the relative stability of the open-state in solution. Based on the results reported above, we determined that 2 was a more direct metric of the ␤-loop dynamics than 1 . Consequently, here we computed potentials of mean force (PMF) along the one-dimensional 2 coordinate space for both genGH1 and GH1.0 ( Figure 3). The PMFs shows that both genGH1 and GH1.0 contain a single broad free-energy well spanning about 25 • with minima at ∼50 • and ∼62 • . Based on the crystallographic structure, these states are consistent with closed ␤-loop. Free energies corresponding to the open ␤-loop (∼100 • ) suggest that this state is sparsely populated in solution for GH1.0, with a free energy of 3.27 ± 0.48 kcal/mol and virtually nonexistent for genGH1 which has an open-state free energy of 6.97 ± 0.56 kcal/mol. For both systems, angles below 45 • led to unphysical steric clashes, which results in increased free energies for these angles.

Bound linker histone dynamics
Off-dyad linker histones sampled multiple states in the DNA binding pocket. There is no reported high-resolution crystal structure for linker histones bound in an off-dyad state. Therefore, a model for off-dyad genGH1 was constructed based on both the low-resolution cryo-EM map of a polynucleosomal array by Song et al. and NMR PRE experiments by Zhou et al. (see 'Materials and Methods' section) (36)(37)(38). Due to the availability of of the higher resolution PRE data, we used the same genH1 as Zhou et al. instead of the H1.4 represented in the array cryo-EM map and the H1.5 from the 5NL0 structure, and we assumed that offdyad states are similar between linker histones, as has been observed for on-dyad states. The GH1.0 off-dyad model was constructed by mutating the genGH1 model to the appropriate primary sequence. This model was validated by comparing to PRE experiments conducted by Zhou et al. To make a direct comparison to the experimental data, we estimated the observed PRE intensity ratios (I ox /I red ) using distances between the core histone labeled probed residues H2A T119 and H3 K37 and genGH1 methyl-terminated residues (Supplementary Figure S5). In general, we observed suitable agreement for both sites, with mean signed errors of 0.02 and 0.04 for the T119 and K37 probe sites, respectively. Comparing our off-dyad model of the chromatosome to that described by Zhou et al., the linker histones exhibit a slight rotational variation in the DNA pocket (38). We attribute this difference largely to our use of data released after the work cited above, such as the cryo-EM map of the nucleosome array that was used to adjust the linker DNA arms (37). Furthermore, our model exhibits more stabilizing LH-DNA contacts (see Supplementary Figures S6-9) leading to an overall lower energy structure. Based on these results, we found the off-dyad structures of genGH1 and GH1.0 sufficient to begin simulations. For reference, we have provided the contacts of genGH1 and GH1.0 with the DNA prior to simulations (Supplementary Figures S8  and 9).
For both genGH1 and GH1.0, there was considerably more sampling of the linker histone position in the offdyad location relative to on-dyad (Supplementary Figure  S10). Specifically, the linker histone was localized within the DNA minor groove where it rocked in and out of the nucleosomal plane, although in our simulation it did not slide along the DNA between the on-and off-dyad states. This is quantified by the increased root-mean-square deviation (RMSD) values in off-dyad simulations (Supplementary Figure S4). This contrast in RMSD values between binding modes is reflected in a clustering analysis which demonstrates that off-dyad systems have more clusters that are less populated than on-dyad systems (Figure 4). More specifically, when clusters were separated by 2Å from one another, six clusters were found for GH1.0 in an on-dyad binding mode, whereas 11 were found for off-dyad GH1.0. Similarly, six and eight cluster were found for genGH1 ondyad and off-dyad. However, some of these clusters had a low population, and when the clusters representing only the top 90% of frames were analyzed this resulted in three clusters for each on-dyad and five clusters for each off-dyad linker histone (Figure 4). These suggest that off-dyad linker histones are more fluid in the DNA pocket. In contrast, the increased number of linker histone-DNA contacts in the on-dyad pocket (discussed below) leads to a more confined and rigid complex, hence decreased sampling.
This increased sampling in off-dyad systems contributes to increased uncertainty in the genGH1 observed timeaveraged PRE ratios (Supplementary Figure S11 and 12 for GH1.0). These values had a higher discrepancy to experiments with mean errors of 0.11 and 0.16, although for most residues the experimental values were within the 80% confidence interval of the simulation derived ratios. Some of these differences are likely due to localized fluctuations of the linker histones that occur throughout the simulations and the fact that distance changes on the order of 2-3Å in methyl/probe distances can result in difference on the order of 0.10-0.15 in the PRE intensity ratio.
In all bound linker histones simulations the ␤-loop remained in the closed conformation. This is exemplified by the 2 angle distributions, which sampled states that had a free energy in solution below 1 kcal/mol with values of 49.7 • ± 4.2 • and 55.0 • ± 11.5 • for genGH1 in on-and off-dyad states, and 52.4 • ± 5.6 • and 63.7 • ± 7.9 • for GH1.0 on-and off-dyad systems ( Figure 2). As previously noted, 1 is a relatively poor metric to distinguish between the open-and closed ␤-loop states.
On-dyad binding restricts DNA motions. Linker histones interact with both the nucleosomal and linker DNA, which has a direct effect on their motions within the DNA binding pocket. Above, we have emphasized the importance of this interaction by showing how the linker histone binding pose affects experimental results. To further probe these dynamics, we plotted the in-and out-of-nucleosomal-plane motions of both linker DNA arms ( Figure 5 and Supplementary Figure S13). Generally, in on-dyad binding modes the interactions of the linker histone with both linker DNA arms restricted both of their motions compared to the nucleosome. However, when the linker histone was bound offdyad, the entry-linker DNA sampling was similar to the nucleosome, whereas the exit-DNA was shifted out of the nucleosomal plane. Due to the asymmetric nature of the Widom 601 DNA sequence, DNA motions were not symmetric in the nucleosome simulations between the entry and exit DNA segments. Beyond Figure 5, we further quantified the in-and out-of-plane linker DNA motions which are termed here as the ␣ and ␤ angles, respectively, for the entry and exit DNA, as inspired by Bednar et al. (34) (see Figure 6 and Supplementary Figure S1 for definitions). The ␣ angles correspond largely to DNA breathing motions, and over all Linker histone binding had a significant effect on both the equilibrium distributions of these motions as well as their correlations to one another. In addition to the 2D-histograms of the ␣ and ␤ angles (Figure 6; Supplementary Figures S15 and 16), their correlations were analyzed by computing the NMI for each angle pair (Table 1), and the changes induced by linker histone binding were quantified with the KL divergence for each probability distribution relative to the nucleosome (Table 2 and Supplementary Table  S2). Of particular note are the ␣-entry/␣-exit angle distributions which highlight the breathing motions of both DNA ends (Figure 6, middle). In this space, the nucleosome samples a wide range of angles in both the entry and exit DNA with little correlation between the two as shown by the NMI value of 0.04. In off-dyad binding there is a similar range of motions for genGH1, which has a KL divergence of only 1.62 relative to the nucleosome, while GH1.0 has a similar range of sampled states with a shifted mean, resulting in a higher KL value of 3.95. For both cases the correlations between the ␣ angles are relatively low, with a modest increase in the NMI values for genGH1 and a slight decrease for GH1.0. In contrast, on-dyad binding shows a significant reduction in the ␣-entry/␣-exit conformational space, as the range of motion of both the entry and exit DNAs are substantially restricted regardless of the linker histone variant. The KL divergences for these states are high, 4.35 and 5.82 for genGH1 and GH1.0, and the correlations for each state are approximately three times higher than for the nucleosome.
Another distribution of interest is the ␣-exit/␤-exit phase space, which describes the in-and out-of-plane motions of the exit DNA (Figure 6, bottom). In the canonical nucleosome the average ␣-exit and ␤-exit angles are 30.5 • and −2.7 • , respectively, indicating that this DNA arm fluctuates about states that are slightly pointing inward when viewed from the side. These motions are largely uncorrelated, with an NMI value of 0.03. Binding in the off-dyad location has a dramatic effect, forcing the DNA outward and shifting the ␤-exit angles to fluctuate around 9.5 • and 7.1 • for the genGH1 and GH1.0 systems. This alters the sampling distributions, with KL values of these angles of 6.33 and 9.48 for genGH1 and GH1.0, and increasing the correlations between the motions by four to 6-fold. In contrast, on-dyad binding decreases the average ␤-exit value to −12.1 • and −6.1 • for genGH1 and GH1.0 as this binding mode pulls the DNA inward from the side view. This results in more modest changes to the ␣-exit/␤-exit phase space and lower correlation increases over the nucleosome.
In contrast to the exit DNA, the entry DNA is largely unaffected by off-dyad binding with KL values below 1.8 relative to the nucleosome. On-dyad binding creates a larger perturbation, with KL values above 5.1 as ranges of motion of both angles are restricted relative to the nucleosome. The average ␤-exit angles have a slight increase from −5.1 • in the nucleosome to 1.8 • and 0.6 • in genGH1 and GH1.0, which are both significantly greater than the values in the exit DNA. Together, this creates an asymmetry in the profile view of DNA structures, as illustrated in Supplementary Figure S15.

Energetic contributions
On-Dyad binding is enthalpically favored. Although both linker histones have similar physical effects on the nucleo-Nucleic Acids Research, 2020, Vol. 48, No. 7 3599  some when bound in on-or off-dyad locations, the difference between their binding energetics determines their in vitro binding mode preference. To estimate binding affinities, we have used an MM/GBSA analysis to decompose the overall binding affinity into the energetic components that drive this interactions (93). Regardless of the variant, on-dyad binding was found to be significantly more enthalpically favorable than off-dyad binding ( Table 3). GH1.0 favored the on-dyad state over off-dyad by −163.3 ± 36.3 kcal/mol, while genGH1 only favored on-dyad by −89.9 ± 39.0 kcal/mol. This difference in E total values was largely driven by Van der Waals (VdW) interaction energies, with GH1.0 favoring the on-dyad state by −91.0 ± 24.5 kcal/mol and genGH1 favoring off-dyad state by 32.1 ± 25.6 kcal/mol. In contrast, genGH1 favors the on-dyad binding mode more than GH1.0 in the electrostatic interaction energies. However, a large overlap in the errors between variants suggests that both systems have a relatively similar electrostatic on-dyad binding preference. We emphasize that MM/GBSA approach includes a number of approximations and do not include important thermodynamic quantities such as conformational entropy or explicit solvent thermodynamics, therefore these values should be taken as qualitative estimates of binding affinities (94,95).

Van der Waals interactions drive binding mode selectivity.
To further investigate the thermodynamic driving forces between linker histone variants and binding locations, we calculated the energetic strain between each binding species in their complexed and isolated states, as defined by: The results in Table 4 suggest that a combination of electrostatic and VdW interactions from both the linker histones and nucleosomes drive the system conformations. The E tot of genGH1 for on-dyad (17.9 ± 5.2 kcal/mol) and off-dyad (28.2 ± 5.6 kcal/mol) is ∼2.5-fold greater than GH1.0 on-dyad (7.4 ± 5.0 kcal/mol) and off-dyad (10.6 ± 5.3 kcal/mol). These variations are largely defined by differences in the VdW interactions where the E VdW of GH1.0 is 14.9 kcal/mol and 17.2 kcal/mol more favorable Table 2. Kullback-Leibler divergence values for 2D probability distributions of DNA Cells are colored on a scale from blue (lower) to white to red (higher). than genGH1 for on-and off-dyad systems, respectively. Additionally, it is worth noting that electrostatic interactions ( E ele ) for both linker histones actually favor the off-dyad complex by 11.6 kcal/mol and 19.9 kcal/mol for genGH1 and GH1.0, respectively.
Off-dyad nucleosomes showed an increased stability when bound to genGH1 over GH1.0 with a E tot of −40.3 ± 29.2 kcal/mol and 5.7 ± 24.8 kcal/mol for genGH1 and GH1.0 systems, respectively. This contrast in E tot can be also be attributed to the VdW energies with E VdW values of −19.6 ± 16.4 kcal/mol and 67.5 ± 15.5 kcal/mol for genGH1 and GH1.0 systems, respectively.
Contacts between linker histones and the DNA (Supplementary Figure S17) show that most on-dyad interactions come from the ␤-sheet (Figure 7), with 66.8 ± 12.9 and 71.2 ± 11.8 contacts for genGH1 and GH1.0. However, the contrast between variants becomes more evident in the off-dyad binding mode with 30.2 ± 14.8 and 2.2 ± 2.6 ␤-sheet-DNA contacts for genGH1 and GH1.0, respectively. A similar relationship is also observed in the off-dyad ␣3 helix and N'-tail. The genGH1 ␣3 helix in the off-dyad binding mode has 16.4 ± 10.8 more contacts than GH1.0 while the N'-tail has 9.1 ± 18.0 more contacts. Combined, these additional contacts between genGH1 and DNA correlate with its VdW-driven preference for off-dyad binding over GH1.0. Taken together, these results suggest that while on-dyad binding is energetically preferred in both variants, genGH1 has a higher propensity for sampling the off-dyad state than GH1.0, which is largely due to the difference in contacts between the ␤-sheet and the DNA.

DISCUSSION
Here, we have used a combination of conventional and free energy MD simulations to probe the effects of linker histone binding on chromatosome structures and dynamics. Umbrella sampling PMFs show that for both genGH1 and GH1.0 the closed-state is thermodynamically favorable, which is consistent with the GH1.0 accelerated MD results ofÖztürk et al. (53). This suggests that the open-state in the 1HST crystal structure is stabilized by crystal packing forces, as exemplified by the fact that the ␤-loop is inserted into a hydrophobic pocket in the neighboring unit. These results are in general agreement with crystal structures of linker histones bound in the on-dyad state which have the ␤-loop in the closed-state (33)(34)46).
Given that there is no high resolution crystal structure of off-dyad binding, the precise binding structure and orientation of linker histones in this pocket remains inconclusive (36,38). We therefore constructed a model of the off-dyad state based on manual placement and docking into the 30nm cryo-EM structure by Song et al. which we found had  A negative value (−) indicates a favorability of the binding species in complex, as oppose to an isolated state, positive (+). Full binding energies for each state are given in Supplementary Table S1. generally good agreement with the PRE data from Zhou et al. (36,37). Based on our umbrella sampling simulations we used the closed-states of the ␤-loops in these structures.
As highlighted by the clustering results in Figure 4, our simulations show that both linker histones are significantly more fluid in the off-dyad DNA binding pocket relative to on-dyad. This is in line with the results of Brownian dynamics docking studies from the Wade group in which they found that genGH1, GH1.0 and assorted mutants can bind in a variety of sequence dependent orientations (54). Furthermore, these series of simulation results would suggest transitions between on-and off-dyad states might be facilitated by multiple stable binding orientations along the DNA and encouraged by additional linker histone conformational freedom in the binding pocket. One of the central mechanisms by which linker histones inhibit transcription and promote the compaction of chromatin fibers is by altering linker DNA dynamics. Our simulations have shown that one of the primary differences in onand off-dyad binding is that on-dyad binding drastically restricts both the entry and exit DNA segments, whereas offdyad binding has a distinct influence on the exit DNA dynamics with little change to the entry DNA. The latter is in line with previous Brownian and MD simulations from the Wade group which have shown that GH1.0 modifies linker DNA motions in off-dyad binding modes (53,56). Furthermore, our results indicate that these effects are largely independent of linker histone variant type. These differences in DNA dynamics have broad implications for greater chromatin structures. For example, Mishra and Hayes have highlighted how the stoichiometric binding of H1 to nucleosome arrays can severely limit linker DNA accessibility to trans-acting factors (31). Meanwhile, cryo-EM structures of polynucleosome arrays have revealed linker histones in both on-and off-dyad locations, with distinct effects on the greater structure of chromatin arrays (37,42). Furthermore, Perišić et al. recently showed with a highly coarse-grained model that linker histone binding position influences their tail positions, which directly impacts greater chromatin structures, with off-dyad linker histones creating more condensed chromatin fibers (41). At last, numerous studies have shown chromatin fiber pliability to be highly dependent on the linker DNA length and, by extension, the nucleosome repeat length (43)(44)(45). In light of these results, it is likely that in vitro on-dyad mono-chromatosomes are a result of additional linker histone DNA contacts which may not be as prevalent in condensed nucleosome arrays due to limited linker DNA conformational freedom (34).
Given that genGH1 and GH1.0 have similar effects on the structure and dynamics of chromatin when bound in onand off-dyad locations, what is the role of various linker histone variants and modifications in vitro? Results from our energetic and contact analyses show that while they have similar structures, the genGH1 and GH1.0 variants have drastically different energetic preferences for the on-and off-dyad states. Indeed, while both have an enthalpic preference for on-dyad binding, the preference is significantly reduced for genGH1. To further quantify the thermodynamics of binding would require explicit calculations of entropic contributions differences in each system. Unfortunately, entropic calculations on systems of this size are notoriously difficult to converge (96). Therefore, we have extrapolated their influences based on changes in the linker DNA and linker histone dynamics. The results suggest that while the enthalpic reward of on-dyad binding in the GH1.0 chromatosome is enough to overcome the entropic penalty from the drastic reduction in linker DNA and linker histone dynamics, that in the case of genGH1 there is not enough of an enthalpic difference between on-and off-dyad states to compensate for these entropic losses, which is why genH1 has typically been observed in the off-dyad binding mode in vitro. This is in line with experimental results that consistently show GH1.0 in the on-dyad state, but that have provided evidence for genGH1 in both the on-and off-dyad states (34,37,46).
In an excellent example of the finely-tuned nature of the on-and off-dyad binding thermodynamics, Zhou et. al used PRE NMR experiments to show that mutation of five GH1.0 residues can shift its binding equilibrium from on-to off-dyad locations (36). We performed additional simulations of this GH1.0 pentamutant, and observed that dynamics in on-and off-dyad states are similar to those for the GH1.0 and gen1.0 systems examined above (Supplementary Section S5: S5.1-4; Figures S19-23). However, these mutations did have a significant effect on the binding energies. Specifically, they lowered the favorability of the on-dyad state, while at the same time had little effect on the off-dyad binding energies (Supplementary Table S4). This change in energy brings the on-dyad binding energy to within error of the off-dyad state. Taken with the increased dynamics in off-dyad states, likely resulting in an entropic penalty of on-dyad binding, this suggests the binding free energy favors the off-dyad state, as observed by Zhou et al., and demonstrates how even relatively minor changes to linker histone structures can have significant influences on their binding thermodynamics.
More specifically, we observe that the change in energetic preference between the two states is driven by increased ␤-sheet, ␣3-helix and N-terminal tail contacts in genGH1 off-dyad systems. Surprisingly, these changes in contacts contributes mostly to differences in the Van der Waals interactions, highlighting their often overlooked influence in protein-DNA binding thermodynamics. Although we emphasize that electrostatics are vital in the binding process, our results show them to be relatively consistent between the on-and off-dyad binding modes, independent of linker histone variant. In the off-dyad binding orientation of the linker histone, the mostly conserved ␤-sheet is more solvent exposed. Given that each linker histone has the same initial placement, the loss of contacts from the on-dyad state to the off-dyad should be similar, specifically between the ␤-sheet and DNA. However, the linker histones facets that remain in contact with the DNA exhibit a less conserved sequence which would express a greater difference in VdW energies due to sidechain variations. We further investigated this by examining the propensity of linker histone residues to be in contact with particular parts of each nucleotides, such as the backbone, sugar and base. Generally, we found that contacts between the DNA backbone and linker hi-stone residues were more dominant, except in the case of off-dyad genGH1 (Supplementary Figure S18). In this system, residues in the ␤-sheet are consistently in contact with the hydrophobic bases of each nucleotide, overshadowing contacts with the backbone. Together, these results suggest that chromatosome systems are driven by a subtle equilibrium wherein multiple binding states may simultaneously be populated in solution. In vitro, there are likely transitions between the on-and off-dyad binding states, with linker histones diffusing along the DNA while guiding chromatin fiber flexibility. While the on-dyad enthalpic reward seen in the GH1.0 chromatosome may be strong enough to overcome the entropic penalty observed in the linker DNA and linker histone, this is not the case for genGH1 which is likely why it has more often been observed in the off-dyad binding mode.
The initial structures used for MD simulations come from a combination of X-ray diffraction (XRD) and cryo-EM structures. In general, cryo-EM-based structures are often more representative of in vitro systems than XRD structures. Despite this, the DNA and linker histone of 4QLC(33) (XRD) are very similar to the cryo-EM 5NL0(34) structures. Moreover, the core histones of our models come from 1KX5 (XRD) which includes models of the tails and loops which are missing from the 4QLC structure. In this structure the tails are in an extended state away from the complex DNA; however in our experience they typically collapse onto the nucleosome DNA in an ensemble of states, which is in line with experimental results that show the tails are not freely exposed in solution (97). For the purpose our study, the tails provide stability to the overall complex and do little to interact with linker histone. At last, experiments have shown the C-terminal tails of H2A and linker histones may have direct interactions with one another that facilitate binding (98)(99)(100). However, because we are modeling only the globular portion of the linker histones that lack the C-terminal tails, we did not observe these interactions which may have an additional effect in shifting the binding equilibria between on-and off-dyad states. At last, Bednar et al. used cryo-EM and X-ray crystallography to show that the long C-terminal domain (∼100 residues) is oriented on the dyad and localized on a single linker DNA arm. In our off-dyad model (Figure 8), this tail would be aligned distal to the dyad and interacting along the outside of a single linker DNA. Here, the C-terminal linker histone tail could compete for binding with the long H3 tail along the same DNA arm potentially inducing a shift from the off-to on-dyad binding mode. The oppositely oriented N-terminal tail (∼20-100 residues) is much less studied and its role in off-dyad binding remains unclear. In our model, the N-terminal domain finds itself along the inside of the same DNA arm where it could interact with either linker DNA segment.
Overall, our study builds on notions of an ensemble of linker histone binding states within a highly dynamic chromatin fiber while emphasizing the contrasting influence of their variants on those structures and dynamics. Currently, the relative populations of these states within chromatin are still debatable. We subscribe to the hypothesis that the on-and off-dyad binding modes exist as an ensemble of states within chromatin fibers (47). The relative populations of these states in vivo are likely due to a balance of not only the linker histone/nucleosome interactions examined here, but also factors outside the scope of this study such as DNA sequence, nucleosome repeat length, and the greater chromatin architecture. Potentially, coarse-grained models may be more adept at sampling the populations of binding states in mono-and poly-chromatosomes arrays. However, care should be taken in these models as the estimated binding energies calculated here demonstrate the importance of Van der Waals interactions within the chromatosome in ad-dition to the more commonly considered electrostatic energies. Therefore, we emphasize that any model which attempts to recapitulate the physics underlying linker histone binding must carefully balance their electrostatic and Van der Waals components.