Recognition of 5-Hydroxymethylcytosine by the Uhrf1 SRA Domain

Uhrf1 binds DNA substrates containing hydroxymethylated CpG sites

Using a newly established DNA binding assay [22], [23] as well as electrophoretic mobility shift assays, we investigated the DNA binding activity of Uhrf1, its SRA domain (SRA^Uhrf1) and the MBD of MeCP2 (MBD^MeCP2) to methylated and hydroxymethylated DNA in direct competition (Figure 1, Supplementary Figure S1; note that all supplementary information can also be found in the Combined Supporting Information File S1). We found that the Uhrf1 constructs bind 5mC and 5hmC containing substrates with similar affinities independent of whether one or both cytosine residues of the palindromic CpG site were modified. Control experiments performed with hemimethylated DNA in competition with either unmethylated substrates or substrates containing no CpG site showed that the observed binding activity to methylated and hydroxymethylated DNA is indeed specific (Supplementary Figure S2). In stark contrast to Uhrf1, we found that MBD^MeCP2 clearly discriminates between methylation and hydroxymethylation, which is in accordance with previous reports [21], [24].

Molecular dynamics simulations of SRA:DNA complexes with 5mC and 5hmC

To investigate the binding mode of the SRA domain to DNA containing 5mC or 5hmC, we performed molecular dynamics simulations for both SRA:DNA complexes. Consistent with the in vitro DNA binding data, modeling of an additional hydroxyl group into the complex structure of the Uhrf1 SRA domain with DNA containing hemimethylated CpG sites revealed no spatial constraints for accommodation of the flipped 5hmC nucleotide within the binding pocket (Figure 2). Based on these initial models of the bound conformation, we performed molecular dynamics simulations for a time interval of 57 ns and monitored the RMSD and RMSF values (Supplementary Figures S3 and S4). In both systems equilibrium was reached after 20 to 30 ns. To assure evaluation of equilibrated systems, we continued the equilibrium simulations for another 27 ns and used only the last 10 ns for subsequent interaction energy analysis [25]. To evaluate the stability of the flipped nucleotides within the binding site, we monitored the occurrence and stability of all hydrogen bonds in the vicinity of the binding site with respect to the progress of the simulations (Figure 3).

Before starting the simulations, all water molecules from the X-ray structure were removed and new water molecules were placed by the setup solvation algorithm of NAMD [26]. Therefore, no water molecules were present in the vicinity of the flipped nucleotides at the beginning of the simulations. Interestingly, in both simulations, water molecules from the water-filled simulation box moved into the nucleotide binding site within the first couple of nanoseconds (Figures 3C and 3D, hydrogen bonds 14 to 18). During the remainder of the simulation time, one water molecule was stabilized within the binding site by formation of distinct hydrogen bonds with protein and DNA. Notably, the position of this water molecule in the 5mC complex corresponds to that of a conserved water molecule in the experimental structure (Supplementary Figure S5), confirming the stability and accuracy of our simulations.

Despite the presence of a conserved water molecule in the binding pockets of both complexes, the corresponding hydrogen bond networks showed interesting differences. In the 5mC complex, this water molecule forms hydrogen bonds with the phosphodiester group of the methylated nucleotide as well as with the SRA residues I454 and G453, thereby bridging the DNA backbone:protein interaction (Figure 3A–C, hydrogen bonds 14–16, Figure 4A). Furthermore, direct hydrogen bonds between the 5mC DNA backbone and the protein are formed involving residues G453, S486, and R489 (hydrogen bonds 1–4).

The hydrogen bond network of the 5hmC complex is more stable compared to the 5mC complex (Figure 3D, compare with 3C). Most prominently, one additional and very stable hydrogen bond is formed between the conserved water molecule and the hydroxyl group of the 5hmC nucleotide (hydrogen bond 17). This interaction seems to specifically stabilize the hydrogen bonding network between the DNA backbone and the binding pocket residues G453, S486, and R489 (hydrogen bonds 1–4). Interestingly, these hydrogen bonds have been previously identified to be important for DNA binding [14] and possibly stabilize the flipped conformation of the nucleotide within the binding site. In addition, the hydrogen bond network within the protein involving residues V466 and G453 as well as residues T484 and D474 is stabilized in the 5hmC complex (hydrogen bonds 11–13).

Since water dynamics and to some extent also DNA dynamics can depend on the ion concentration parameters used in the molecular dynamics simulation, we performed a second simulation of the 5hmC complex with a higher ion concentration (Supplementary Figure S6). Consistent to the first simulation with 5hmC, we observed the same overall water dynamics and hydrogen bonding patterns including hydrogen bond formation between the hydroxyl group of the 5hmC nucleotide and the conserved water molecule within the SRA structure. Notably, the stable hydrogen bonding between protein residue S486 and the DNA backbone in the first simulation (hydrogen bonds 2a and 2b) seems to be replaced by a stable hydrogen bond of S486 with the water molecule in the second simulation (hydrogen bond 18), indicating two alternative interaction patterns for the S486 residue in the 5hmC complex (Figures 4B and 4C, compare Figure 3D and Supplementary Figure S6B). In conclusion, these data suggest that stable, water bridged hydrogen bond formation of the hydroxyl group of the flipped 5hmC nucleotide with its surrounding occurs in and stabilizes this DNA:SRA complex.

Similar interaction energies for SRA complexes with 5mC and 5hmC containing DNA

To estimate the binding affinity between the Uhrf1 SRA domain and DNA containing either 5mC or 5hmC, we calculated the respective interaction energies using the linear interaction energy (LIE) approach [25]. To exclude energy contributions due to base-flipping when comparing the interaction of the DNA with the protein (bound state) or with the solvent (unbound state), we simulated the DNA in a flipped state in both cases. We determined the difference between the binding energies of the two complexes (ΔΔG = ΔG_5mC−ΔG_5hmC). We included either i) the whole DNA and SRA structure (ΔΔG = −7.94 kcal/mol) or ii) the flipped nucleotide with its five neighboring nucleotides and the binding pocket of the protein, defined as all residues within a distance of 15 Å from the nucleotide in the starting conformation (ΔΔG = −6.65 kcal/mol). These values suggest that the slight difference in binding affinity is predominantly due to interaction of the flipped nucleotide with the proximal protein residues that form the binding site. Considering the estimated uncertainty of about 3–4 kcal/mol in our calculations, these values indicate that both 5mC and 5hmC containing DNA substrates bind with very similar affinity to the SRA domain of Uhrf1.

Expression constructs, cell culture and transfection

Mammalian expression constructs for enhanced green fluorescent protein (GFP), Uhrf1 (GFP-Uhrf1), the SRA domain of Uhrf1 (GFP-SRA^Uhrf1) and the MBD of MeCP2 (MBD^MeCP2-YFP) were described previously [22], [23], [34]. Note that all constructs encode fusion proteins of either GFP or yellow fluorescent protein (YFP). HEK293T cells [35] were cultured in DMEM supplemented with 50 µg/ml gentamicin and 10% fetal calf serum. For expression of GFP/YFP fusion proteins, HEK293T cells were transfected with the corresponding expression constructs using polyethylenimine (Sigma).

DNA substrate preparation

Fluorescently labeled DNA substrates were prepared by mixing two HPLC-purified DNA oligonucleotides (IBA GmbH, Supplementary Tables S1 and S2) in equimolar amounts, denaturation for 30 sec at 92°C and slow cool-down to 25°C allowing hybridization. After purification by 15% non-denaturing PAGE, DNA substrates were resuspended in binding buffer (20 mM TrisHCl pH 7.5, 150 mM NaCl, 1 mM EDTA, 1 mM DTT).

Pull-down DNA binding assay

In vitro DNA binding assays were performed as described previously [22], [23]. In brief, GFP/YFP fusions were purified from HEK293T extracts using the GFP-Trap® (ChromoTek GmbH) and incubated with two differentially labeled DNA substrates at a final concentration of 200 nM DNA/50–100 nM immobilized protein for 45 min at room temperature in binding buffer. After removal of unbound substrate, the amounts of protein and DNA were determined by fluorescence intensity measurements with a Tecan Infinite M1000 plate reader. Binding ratios were calculated dividing the concentration of bound DNA substrate by the concentration of GFP/YFP fusion on the beads, corrected by values from a control experiment using DNA substrates of the same sequence but with different fluorescent labels, and normalized by the total amount of bound DNA.

Electrophoretic mobility shift assay

For competitive electrophoretic mobility shift assays, equimolar amounts of two differentially labeled DNA substrates (250 nM each) were incubated with increasing amounts of GFP/YFP fusion protein (Supplementary Figure S1), subjected to 6% non-denaturing PAGE and analyzed with a Typhoon scanner (GE Healthcare), which allowed separate detection of DNA substrates and protein by ATTO labels and GFP tag, respectively, using the following laser/filter settings: 532 nm/580 nm (ATTO550), 633 nm/none (ATTO700), 488 nm/520 nm (GFP/YFP).

Molecular dynamics simulations

Molecular dynamics simulations were performed based on the X-ray structure of the Uhrf1 SRA domain with the PDB identifier 3FDE [14], using the program NAMD 2.7b1 [26] and the CHARMM22/27 force field [36], [37]. Binding free energies were estimated using the Linear Interaction Energy (LIE) model [25].

After energy minimization of 50,000 steps, one hydrogen atom of the methyl group of the protein-bound 5-methylcytosine (5mC) residue was substituted by a hydroxyl group using the tool psfgen. CHARMM22 force field parameters were available for 5mC (patch: PRES 5MC2), but not for 5-hydroxymethylcytosine (5hmC). Therefore, a new 5hmC residue was created based on the 5mC parameters and topology. For this purpose, one hydrogen atom of the 5mC methyl group was exchanged by a hydroxyl group. The charges of the hydroxyl group were subsequently set to charges of the hydroxyl group of a serine residue according to the CHARMM27; the charges of the CH₂ group were adjusted accordingly (Supplementary Table S3). After solvation, the 5mC and 5hmC structures were further energy minimized for 50,000 steps. For each structure, two simulations were performed, in which the charges were either neutralized or a salt concentration of 0.5 M was used.

Each simulation was performed using periodic boundary conditions and particle-mesh-ewald summation [38] for long range non-bonded interactions. The non-bonded cutoff was set to 14 Å with a switching/shifting distance of 12 Å. A stepsize of 1 fs was chosen. The systems were heated from 0 to 200 K for 160 ps under constant volume. Harmonic restraints (1000 kcal mol⁻¹ nm⁻²) were applied to all atoms of the complex. The heat up was continued without harmonic restraints from 200 to 300 K for 80 ps under constant pressure conditions, using a Nose-Hoover barostat [39], [40] with a target pressure of 1.01325 bar, an oscillation time scale of 100 fs, and a damping time scale of 50 fs. The temperature was maintained by Langevin dynamics using a damping coefficient of 5/ps. The temperature bath was not coupled to hydrogen atoms. After the heat up procedure, the simulations were continued for 57 ns. During the simulations, all bond lengths were constrained to ideal values using the Shake algorithm [41], [42].

For analysis of the simulation results, all hydrogen bonds formed by the flipped nucleotides and the binding site were identified and monitored throughout the simulations and the occurrence of water molecules in and around the binding site was monitored every 5 ps. In order to estimate the difference in the binding free energy of the two nucleotides, we performed three further simulations in which the protein and the two DNA molecules were simulated separately using the conditions described above. To keep the DNA in the flipped state, we additionally applied harmonic restraints to the whole DNA backbone (atom names: C4′, P, O1P, O2P, O5′, C5′, C3′, O3′). The solvated single protein was simulated for 34 ns and the separated DNA molecules were simulated for 20 ns.

To estimate the binding affinity of the two DNA molecules to the protein, we estimated the binding free energy according to the Linear Interaction Energy (LIE) model [25]:

(1)

(2)

In this approach the binding free energy is approximated by the difference between the interaction energies ΔV^el and ΔV^vdw of the ligand in the protein-ligand complex (bound state) and in solution (unbound state). The <> denotes the average values obtained from the simulation trajectories. According to the linear response approximation the weights α and β were set to 1 and 0.5, respectively. We calculated the DNA-(protein+solvent) (bound state) and the DNA-solvent (free state) interaction energies from the trajectories of the DNA/SRA and the DNA/solvent simulations, using the average energy over the last 10 ns.