DNA deformation energies

Next, we analyse DNA deformation energy, that is the energy of DNA thermal fluctuations, in all eight MD trajectories using a multivariate Ising model, parameterised for all tetranucleotides (Liebl and Zacharias, Reference Liebl and Zacharias2021), derived to include the coupling between all possible sequence-specific conformational substates of DNA. First, we calculate the deformation energy for the region covering both TF REs and the first adjacent flanking nucleotides, namely, 6–35 b.p. (Table 1 and Supplementary Fig. S6A). If a protein binding reduces the deformation energy relative to naked DNA, it indicates that the bound conformation can be readily sampled by naked DNA under physiological conditions, pointing towards a conformational selection mechanism. However, if a protein binding leads to a large increase in deformation energy relative to naked DNA, it indicates that DNA needs an external energy source to adopt the bound conformation, pointing towards an induced fit mechanism. As discussed by Battistini et al. (Reference Battistini, Hospital, Buitrago, Gallego, Dans, Gelpi and Orozco2019), a double increase or more in DNA deformation energy upon a protein binding indicates an induced fit mechanism; whereas a value between the one of naked DNA and the double increase indicates that both mechanisms co-exist. In our case, using the multivariate Isling model instead of a simple harmonic model approximation, the naked DNA thermal energy fluctuation is ~4 kcal/mol*b.p. We observe that in the WT case, DNA deformation energy slightly increases upon binding of the CEBPB dimer to naked DNA (by 0.18 kcal/mol*b.p., Table 1), while the E2F1-DP1–DNA association has no obvious impact. As stated above this is not a large increase, thus the DNA conformation for the CEBPB binding should be readily accessible within the thermal fluctuations of DNA. Upon DNA methylation, the E2F1-DP1–DNA binding further increases the deformation energy, irrespective if CEBPB is present or not (by 0.44–0.48 kcal/mol*b.p., Table 1). Although the increase is relatively low, it is a noteworthy change, which suggests that the binding of E2F1-DP1 to ME-DNA could depend on induced fit changes, which may lower the dimer–DNA affinity.

Table 1. DNA deformation energies (kcal/mol) for the 6–35 b.p. region, calculated using a multivariate Ising model (Liebl and Zacharias, Reference Liebl and Zacharias2021)

Note: Values in the parentheses correspond to the differences in the deformation energy with respect to naked WT-DNA for the entire region and per b.p.

Abbreviations: ME, methylated; WT, wild type.

Second, we divide the DNA sequence into three regions: the E2F1-DP1 binding site, the linker, and the CEBPB binding site, and recalculate the deformation energy for each region (Fig. 2b,c and Supplementary Fig. S6B–D). This segmentation provides further insights. The deformation energy is not equally distributed over the DNA sequence. Although the standard deviations are high which make the calculated deformation energies overlap, the large sampling sets (1 million snapshots) make the standard errors approach zero. Thus, we will focus on comparing the mean values (Fig. 2c). In the WT case, the CEBPB–DNA binding increases the deformation energy of the linker region (by 0.17 kcal/mol*b.p.), suggesting the induced conformational changes – analyses of DNA helical parameters confirm that (see below). This increase also slightly reduces the deformation energy within the E2F1-DP1 binding site (0.06 kcal/mol*b.p.), suggesting it becomes more favourable for E2F1-DP1 to bind to DNA. Nevertheless, the differences in the calculated deformation energies between all systems (Fig. 2c and Supplementary Fig. S6B–D) are relatively small, the only noteworthy change is the increase in deformation energy of the E2F1-DP1-binding site caused by the binding E2F1-DP1 to ME-DNA (0.6–0.8 kcal/mol*b.p.). This could, as suggested, make the binding of E2F1-DP1 to ME-DNA more dependent on induced fit changes, which in turn could make the binding site less recognizable. Taken together, the analyses hint at the cooperative nature of the TF dimers–DNA association, which can be perturbed by the hypermethylation.

Protein–DNA contacts

We continue with the analyses of protein–DNA contacts to further understand the mechanistic impact of CpG methylation on the E2F1-DP1–CEBPB–DNA enhanceosome formation. For the analyses, we employ a dynamic contacts map approach we derived earlier (Hörberg and Reymer, Reference Hörberg and Reymer2020; Hörberg et al., Reference Hörberg, Moreau, Tamás and Reymer2021), where for each MD frame, we calculate a contact strength for every pair of protein–DNA residues that interact specifically and nonspecifically. We follow the time-evolution of the contacts strengths (Supplementary Figs S7–S10) and calculate the average strength value for each protein–DNA contact (Fig. 3 and Supplementary Figs S11, S13 and S14). The analyses of the dynamic contacts maps for the WT versus ME cases and solo- versus enhanceosome systems show that the CEBPB dimer exploits nearly the same network of contacts in all systems (Supplementary Figs S11 and S13). This suggests that the CEBPB–DNA binding is independent of CpG methylation or the binding of the E2F1-DP1 dimer. The CEBPB dimer is a BZIP factor, which recognises its REs through specific interactions formed by a conserved five residues motif (NxxAVxxSR) (Fujii et al., Reference Fujii, Shimizu, Toda, Yanagida and Hakoshima2000), residues 281–289, of each monomer. The studied CEBPB-RE has a non-canonical CEBPB half-site (TAGTGTAA), which explains the observed differences in the specific contacts exploited by the two monomers (Supplementary Figs S11 and S12). Of the five residues motif of monomer one, later CEBPB1, Ala284, Val285 and Ser288 interact hydrophobically with the flanking 5′-GpT, first TpA steps, and the adenine of the TpA step on the Crick strand (GTAGTGTAA), respectively, whereas Arg289 is involved in hydrogen bonding with the central CpA-step (TTACACTA) and the guanine of the first GpT step of the CEBPB1 half-site (TAGTGTAA). Though, the Arg289 contacts are somewhat weaker in the ME-systems. Specific interactions of the five residues motif of monomer two, later CEBP2, include cross-bridging hydrogen bonds by Asn281 to the first two TA b.p. (TTACACTA/TAGTGTAA), hydrophobic contacts by Ala284, Val285 and Ser288 with the flanking 5′-GpT step and the TpT step of the CEBPB2 half-site (GTTACACTA), and hydrogen bonding by Arg289 with the central TpG step (TAGTGTAA). Some of the abovementioned contacts exploited by Ser288 and Arg289 of both CEBPB monomers only appear in some of the simulated systems (Supplementary Figs S11 and S12). However, we believe that these contacts are not system specific, rather they are coupled to the flickering power of long-chained residues, which can fluctuate between different conformational substates.

Fig. 3. Specific contacts between the E2F1-DP1 dimer and DNA. The plots show the strength of specific contacts formed by DP1 [marked with (1)] and E2F1 [marked with (2)] monomers in wild-type (WT) and methylated (ME) systems when bound alone to DNA (solo) (panel a) and together with CEBPB (complex) (panel b). Panels c and d show the strength of specific contacts by DP1and E2F1 monomers in solo- versus complex-systems in WT and ME cases, respectively. For the definition of a contact strength, see Supplementary Material.

In contrast, the DNA contacts exploited by the E2F1-DP1 dimer show a great variation among the simulated systems (Figs 3 and 4a, and Supplementary Fig. S14). We observe a rearrangement of the E2F1-DP1–DNA contacts when comparing the WT- and ME-E2F1-DP1–DNA systems, and a strong impact of TF-cooperativity for the E2F1-DP1–DNA binding. The E2F1-DP1 dimer belongs to the E2F-DP family of TFs, the winged-helix folded proteins that recognise their REs with a conserved four residues motif (RRxYD), residues 167–171 and 165–169 for the DP1 and E2F1 monomers, respectively (Zheng et al., Reference Zheng, Fraenkel, Pabo and Pavletich1999. For DP1, Arg167 is involved in hydrogen bonding with the three guanines of DNA Crick strand (TTTCCCGC/GCGGGAAA) in the WT systems, in contrast the residue forms predominantly hydrophobic contacts with the methyl group of the C_W13 base in the ME-systems. Tyr170 forms hydrophobic contacts with the TpC step (TTTCCCGC) in the WT-CEBPB-E2F1-DP1–DNA system and only with T_W10 base (TTTCCCGC) in the ME-E2F1-DP1–DNA system (Fig. 4b). For E2F1, Arg165 forms hydrophobic contacts with the flanking GpC step (GCGCGGGAAA). For both monomers, Arg168/Arg166 interact nonspecifically with DNA backbone, while Asp171/Asp169 show no direct contact with DNA, instead stabilising the orientation of the Arg167–168/Arg165–166 residues. It should, however, be mentioned that Arg166 of E2F1 shows a reduction in the contact strength of nonspecific interactions with ME-DNA; instead, it forms mainly a salt bridge with Asp171 of DP1. This could affect the affinity of the heterodimer binding to ME-DNA. Furthermore, E2F1 has an N-terminal random-coil-tail that twists around DNA forming mostly nonspecific contacts, except for Arg127 which also specifically interacts with several b.p. of the RE (TTTCCCGC/CGGGAAA). Also, for the WT- and ME-E2F1-DP1–CEBPB–DNA systems there is an additional specific contact by Lys117 of the N-terminal tail. However, as random-coil-tails are highly flexible, exhibiting a large conformational space that is hard to sample completely, the observed variations could be simulation-specific rather than coupled to the impact of methylation.

Fig. 4. (a) Schematic representation of specific protein–DNA contacts exploited by the four residues motif (RRxYD) of E2F1-DP1 in different systems. (b) Cartoon representation of the binding orientations and DNA contacts of the E2F1-DP1 dimer for WT (blue) and ME (orange) systems, when bound alone to DNA (solo) and together with CEBPB (complex).

Overall, we observe a reduction in the strength and number of specific and nonspecific DNA–protein contacts for the WT- and ME-E2F1-DP1–DNA systems (Fig. 3 and Supplementary Fig. S14). In the WT-E2F1-DP1–DNA system, the DP1 recognition helix moves out of the major groove, whereas in the ME-case, the recognition helix changes its orientation and slides downwards (Fig. 4b). The two binding modes of the DP1 monomer in the WT- and ME-E2F1-DP1–DNA systems lead to the loss of specificity of E2F1-DP1–DNA binding. For the WT system, the largest changes include the loss of contacts exploited by DP1 Arg167 and Tyr70, which puts the monomer in a position resembling a dissociation state. E2F1 also shows a reduction in contacts, mostly those exploited by the N-terminal tail. This could indicate that DNA is more globally flexible which makes it harder for the tail to settle. For the ME-system, DP1 shows a reduction in contacts strengths exploited by Arg167 and Tyr70; however, the effects are smaller than for the WT system. The DP1 monomer slides downwards, where Arg167 instead interacts with one of the ME cytosines, which give rise to new contacts exploited by another residue, Lys163. As for the WT system, E2F1 shows a reduction in contacts exploited by the N-terminal tail. The presence of CEBPB, secures a stable E2F1-DP1–DNA binding, with higher DNA sequence-specificity in the WT case, which implies the importance of TF-cooperativity for the E2F1-DP1–CEBPB–DNA enhanceosome formation and suggests that the CEBPB–DNA binding should occur first. Nevertheless, for the ME-systems the impact of the CEBPB presence on DNA for the E2F1-DP1–DNA association, according to the contact analysis, appears less significant indicating that the cooperativity is stronger in the WT case.

Protein–DNA interaction energies

We proceed with the analysis of protein–DNA interaction energies. The accurate calculation of the protein–DNA interaction energy is difficult due to the massive negative charge of the DNA backbone. Nevertheless, for a qualitative comparison, we calculate the interaction energies including electrostatic and van der Waals, for specific and nonspecific protein–DNA contacts along all protein-bound trajectories using GROMACS energy analysis tool (Supplementary Figs S15–S20). The interaction energies for the specific and nonspecific contacts follow the trends seen in the dynamic contacts maps (Supplementary Figs S7–S10). The interaction energy distributions for the specific contacts are bimodal (Supplementary Fig. S19), which can be coupled to the flickering behaviour of the long-side chain residues. The nonspecific interaction energy distributions (Supplementary Fig. S20) are normal and broad, which reflect the large variations in the number and contact strength of most of the nonspecific contacts. When we combine the interaction energies for specific and nonspecific contacts (Table 2), it appears that cytosine methylation contributes to a stronger TF-DNA complexation. However, the differences are small with overlapping standard deviations. Furthermore, estimating interaction energies where small changes in contacts, potentially due to different conformational substates exploited by the flickering residues, can lead to overestimated changes in interaction energies. For instance, Arg289 of CEBPB1 can fluctuate between interactions with C30w, A31w and G32c, which explains most of the fluctuations in the interaction energies (Supplementary Figs S15, S16 and S19). In the WT-CEBPB–E2F1-DP1–DNA system, Arg289 shows less sampling of the Arg289-A31w state, which leads to a lower mean value of the specific interaction energy. Taken together, the analysis again indicates that TF-cooperativity greatly affects the E2F1-DP1–DNA binding. We also calculate interaction energies following the MMGBSA/MMPBSA approach (Supplementary Table S2). We exclude calculations of conformational entropies due to the size of the systems. The MMGBSA/MMPBSA energies follow the same trends, again highlighting a cooperative dependency on CEBPB for strong E2F1-DP1–DNA interactions and no deducible impact of DNA methylation.

Table 2. Protein–DNA interaction energies (kcal/mol) including standard deviations, calculated between all DNA–TF-dimer complexes

Note: For the E2F1-DP1–CEBPB–DNA trajectories, the provided interaction energies are for the protein dimer in bold.

Abbreviations: ME, methylated; WT, wild type.

Difference in mean values with respect to a: WT-E2F1-DP1–CEBPB–DNA, b: WT-CEBPB–E2F1-DP1–DNA, c: ME-E2F1-DP1–CEBPB–DNA, d: ME-CEBPB–E2F1-DP1–DNA.

To further address the impact of DNA methylation and TF-cooperativity, we look at configurational entropies calculated with the Schlitter’s formula (Harris et al., Reference Harris, Gavathiotis, Searle, Orozco and Laughton2001) (Table 3) for the entire DNA sequence, the E2F1-DP1-site (b.p. 6–15), and the CEBPB-binding site (b.p. 26–35) for all systems. The configurational entropy of the entire DNA sequence decreases upon the TFs binding and methylation, indicating that the molecule becomes more globally rigid. However, the rigidity is not equally distributed. Looking at the configurational entropies of the specific sites, we observe interesting mechanistic effects. First, for the WT systems, we observe that the CEBPB binding decreases the configurational entropy within the E2F1-DP1-binding site (~0.3 kcal/mol), whereas the binding of E2F1-DP1 increases the configurational entropy within the CEBPB-binding site (~1.0 kcal/mol). This indicates that CEBPB slightly reduces the global conformational flexibility of the E2F1-DP1-binding site, which may facilitate the binding of E2F1-DP1. However, the changes are small, thus, caution should be taken. Second, the configurational entropy of the CEBPB-binding site is about the same (~5 kcal/mol) when CEBPB is bound in the absence or presence of E2F1-DP1, whereas the configurational entropy of the E2F1-DP1-binding site is further reduced when E2F1-DP1 binds in the presence of CEBPB (~4 kcal/mol). This indicates that CEBPB stabilises the binding of E2F1-DP1. Third, for the ME-systems, we observe for unbound ME-DNA relative to unbound WT-DNA that the configurational entropy is reduced within the E2F1-DP1-binding site (~2 kcal/mol) and increased within the CEBPB-binding site (0.7 kcal/mol). These results indicate that DNA methylation may also have an impact on the binding of CEBPB. Fourth, the binding of CEBPB to ME-DNA increases the configural entropy of the E2F1-DP1-binding site (~1 kcal/mol) and the binding of E2F1-DP1 increases the configural entropy of the CEBPB-binding site (0.1 kcal/mol relative to ME unbound DNA and 0.8 kcal/mol relative to WT unbound DNA). This suggests that upon methylation the order of the TF dimers binding to DNA may break.

Table 3. Configurational entropies [TS (kcal/mol), for T = 300 K] for the DNA backbone P atoms

Note: Entropies have been derived for the windows 1,000 ns, last 750 ns, and last 500 ns. Values in parenthesis are the differences in the mean values with respect to naked DNA.

Abbreviations: ME, methylated; WT, wild type.

Changes in DNA helical and groove parameters

Next, we analyse changes in DNA helical and groove parameters for the WT versus ME cases and solo- versus enhanceosome systems (Fig. 5 and Supplementary Figs S21–S27). Previously, we reported that the binding of proteins to their REs results in small conformational adjustments in DNA helical parameters, which facilitates formation of stable specific contacts (Hörberg et al., Reference Hörberg, Moreau, Tamás and Reymer2021). The changes are reflected in reducing the bimodality/multimodality and narrowing of the helical parameters’ distributions, indicating that proteins stabilise a particular DNA substate. The analyses show that all CEBPB-bound DNA exhibit similar changes in helical and groove parameters, with respect to naked DNA, within the CEBPB RE region (b.p. 27–34). This agrees with the protein–DNA contact analyses that show nearly identical CEBPB–DNA contact networks in all systems (Supplementary Fig. S12). In addition, the CEBPB binding introduces changes in shift, slide, and grooves parameters of the four flanking b.p. on each side of the CEBPB-RE (b.p. 23–26 and 35–38), further in the linker region at C_W20 and in the E2F1-DP1 RE at C_W10 (Fig. 5). The CpG20–21 step is not involved in any interactions with the E2F1-DP1 dimer but is a soft dinucleotide step, conveniently located in the middle of the linker. The positive shift and the deeper major groove at CG10, in the presence of CEBPB, stabilise the hydrophobic contacts formed by Tyr170 of DP1, which secures the deep placement of the DP1 recognition helix into the major groove, thus explaining the nature of the observed TF-cooperativity (Fig. 4). In the ME-E2F1-DP1–CEBPB–DNA case, the effects from the CEBPB binding are smaller, due to the changes in helical shift coupled to DNA methylation (Fig. 5).

Fig. 5. Average values for shift, major groove width and depth (in Å) for wild-type (WT) and methylated (ME) DNA alone and in complex with either CEBPB/E2F1-DP1 dimers, or the complete CEBPB–E2F1-DP1 enhanceosome. The WT values are depicted with bold lines, ME – with dashed lines.

In contrast, the changes in DNA helical and groove parameters induced upon the binding of the E2F1-DP1 dimer within its RE (7–14 b.p.), the flanking regions (b.p. 3–6 and 15–18), the linker, and further in the CEBPB RE (b.p. 28, 30, 33–34) are system specific. Similarly, to the CEBPB–DNA system, the changes induced by the E2F1-DP1 dimer happen at C_W20. The changes in helical parameters within the CEBPB RE, induced by the E2F1-DP1 binding, do not significantly alter the CEBPB–DNA specific contacts network, as mentioned above. Methylation of CpG steps leads to a decreased twist and an increased roll (Supplementary Figs S23 and S24), which through DNA backbone BI → BII transitions, induce changes in shift and slide, and consequently in major groove depth and width (Fig. 5 and Supplementary Figs S21, S22, S26 and S27). Focusing on shift, where we observe the opposite sign of shift between the WT and ME cases for the E2F1-DP1-bound systems for the last three b.p. steps of E2F1-DP RE (TTTCCCGCG). The effect is more obvious when we analyse the b.p. displacement in the x-axis (Supplementary Fig. S25). This explains the rearrangement of the DNA contacts formed by the E2F1-DP1 dimer, and the appearance of the new specific contacts by Lys163 and Asn164 of DP1 for the ME-E2F1-DP1–DNA and ME-E2F1-DP1–CEBPB–DNA systems, respectively (Fig. 4). For the linker region, in addition to the change of the shift sign for the ME b.p., we also observe the narrowing of the shift distributions (Supplementary Fig. S21B). This indicates that DNA methylation stabilises different from the WT case conformational substates, which explains the damping of the allosteric communication between the two TF dimers.

Previous studies have shown that cooperativity between two TFs, which do not exhibit direct protein–protein interactions, arises due to DNA-mediated allosteric signals transmitted by one TF to the binding site of the other TF in a time-dependent manner (Balaceanu et al., Reference Balaceanu, Pérez, Dans and Orozco2018). Thus, we next analyse long-distance correlations of DNA groove parameters and helical translational parameters with and without a time lag for the entire trajectories (0.1–1.1 μs) (Fig. 6 and Supplementary Figs S28 and S29). Our results show, in agreement with the previous studies, that the binding of CEBPB, but not E2F1-DP1, results in long-distance correlations of DNA b.ps. parameters (Balaceanu et al., Reference Balaceanu, Pérez, Dans and Orozco2018). The strongest impact of the CEBPB binding on long-distance correlations is seen for the shift parameter and the major groove width (Fig. 6). The signal is translated from C_W32 in the middle of the CEBPB1 RE, over to G_W21G_W22A_W23 and then from T_W19C_W20 in the linker region to the E2F1-DP1 RE b.p. 10–14 (TTTCCCGC). However, the correlation coefficients are small (<0.3), therefore caution should be taken. The addition of a time lag has not improved the long-distance correlations. Furthermore, the long-distance communication between the two binding sites increases when both TFs are bound to DNA. This could indicate that cooperativity and communication first arise when one TF is bound to DNA and the other slides and approaches its binding site. The long-distance correlations somewhat increase (<0.35) when calculated for the trajectories intervals when the flickering residues form specific contacts with DNA. However, we believe that taking the complete trajectories provides a more complete picture of the process. The described long-distance correlations observed within the linker and E2F1-DP1 RE regions upon the CEBPB binding are missing when DNA becomes ME (Supplementary Figs S28 and S29), suggesting that methylation hinders the propagation of the allosteric signal.

Fig. 6. Correlation coefficient maps without time lag for b.p. shift and major groove width along the WT and ME sequences for DNA alone and in complex with either CEBPB/E2F1-DP1 dimers, or with both TFs dimers CEBPB–E2F1-DP1. Correlation coefficients are calculated for the entire (0.1–1.1 μs) trajectories. In all panels, the E2F1-DP1 response element is marked with yellow and the CEBPB response element is marked with magenta.

Ion populations

Previous studies have shown that counterion populations within DNA grooves are sequence-specific: depending both on sequence-specific flexibility of the DNA backbone and electronegative groups in DNA bases, and may impact the TF–DNA association (Hud and Polak, Reference Hud and Polak2001; Dans et al., Reference Dans, Faustino, Battistini, Zakrzewska, Lavery and Orozco2014; Pasi et al., Reference Pasi, Maddocks and Lavery2015). Thus, we also investigate the counterion populations within the major and minor grooves for the naked and protein-bound WT- and ME-DNA systems, using Canion program. First, for the WT systems, we observe an elevated K+ concentration at C10 and C14 positions in the absence of CEBPB (Supplementary Fig. S30A). Second, we observe a slight increase in K+ molarity in the major groove of the E2F1-DP1 RE region for naked and CEBPB-bound ME-DNA, in comparison to the corresponding WT systems. In addition, the analyses show that cytosine methylation affects the radial distribution of K+ ions within the major groove; the ions are shifted by 0.5–1.0 Å closer to DNA helical axis with respect to WT-DNA (Supplementary Fig. S28B). These observations agree with the conclusions derived through dynamic contacts maps showing that the recognition helix of DP1 moves out of the major groove in the WT-E1F2-DP1–DNA case; and in the ME-E1F2-DP1–DNA case, the helix slides downwards reducing the specificity of the DP1–DNA binding.