Vaibhav Kumar
Shukla‡
a,
Lucas
Siemons‡
a,
Francesco L.
Gervasio
bcd and
D. Flemming
Hansen
*a
aDepartment of Structural and Molecular Biology, Division of Biosciences, University College London, London, WC1E 6BT, UK. E-mail: d.hansen@ucl.ac.uk
bDepartment of Chemistry, University College London, London, WC1E 6BT, UK
cPharmaceutical Sciences, University of Geneva, Geneva CH-1211, Switzerland
dInstitute of Pharmaceutical Sciences of Western Switzerland, University of Geneva, Geneva CH-1211, Switzerland
First published on 27th May 2021
Human histone deacetylase 8 (HDAC8) is a key hydrolase in gene regulation and an important drug-target. High-resolution structures of HDAC8 in complex with substrates or inhibitors are available, which have provided insights into the bound state of HDAC8 and its function. Here, using long all-atom unbiased molecular dynamics simulations and Markov state modelling, we show a strong correlation between the conformation of aromatic side chains near the active site and opening and closing of the surrounding functional loops of HDAC8. We also investigated two mutants known to allosterically downregulate the enzymatic activity of HDAC8. Based on experimental data, we hypothesise that I19S-HDAC8 is unable to release the product, whereas both product release and substrate binding are impaired in the S39E-HDAC8 mutant. The presented results deliver detailed insights into the functional dynamics of HDAC8 and provide a mechanism for the substantial downregulation caused by allosteric mutations, including a disease causing one.
HDACs are traditionally divided into four classes based on sequence similarities. Class I (HDAC-1, -2, -3, and HDAC8), class II (HDAC-4, -5, -6, -7, -9, and HDAC10), and class III (SIRT1-7) have sequence similarity to yeast Rpd3, Hda1, and Sir2, respectively, whereas class IV (HDAC11) shares sequence similarity with both class I and II proteins.2 In cancers, class I HDACs act as major epigenetic players and are linked to deregulated expression or interactions with transcription factors critical to tumorigenesis.9 Histone deacetylase 8, a class I HDAC, is also implicated in other diseases, including X-linked intellectual disability and parasitic infections.3,10,11 Additionally, genetic mutations leading to impaired HDAC8 activity are reported in patients with Cornelia de Lange syndrome (CdLS).7,8,12–16
HDAC8 is a unique class I HDAC in terms of its structure, activity, and regulation. Unlike the other three members of class I, HDAC8 is constitutively active, whereas HDAC1 and HDAC2 show activity that increases dramatically once present in multi-protein complexes, and HDAC3 is inactive in isolation.2 Moreover, HDAC1, HDAC2, and HDAC3 have serine phosphorylation sites in their flexible C-terminal regulatory tail, which upon phosphorylation activates the enzyme. In contrast, HDAC8 has a regulatory serine phosphorylation site at position 39, in the catalytic domain, and its phosphorylation inactivates the enzyme by stabilising an inactive state.2,17 Recently, in HDAC1, HDAC2, and HDAC3, the region corresponding to the S39 site in HDAC8 was shown to act as a binding platform in holoenzyme HDAC complexes, suggesting a general mechanism for regulation across class I HDACs.17–20
The overall fold of HDAC8, as observed in substrate- and inhibitor-HDAC8 complexes, is comprised of an eight stranded β-sheet, which is surrounded by 11 helices and two helical turns.21–26 The active site is formed by seven loops, L1 to L7, connecting these secondary structure elements. The substrate-binding tunnel is formed by residues D101, H142, H143, G151, F152, H180, F208, M274, and Y306 with the active-site Zn2+ located at the base of the tunnel, Fig. 1a and b. These residues are conserved across the class I HDACs, except for M274, which is a leucine in the other class I family members. For class II members HDAC6 and HDAC10 a lysine and a glutamate is present, respectively, whereas a leucine is present in other class II HDAC family members at this position.21
It has been proposed that, in the unbound form, the active-site Zn2+ adopts a 5-coordinate square pyramidal geometry with D178, H180, and D267 and two water molecules.27 Upon substrate or inhibitor binding, one of the water molecules is displaced by the side-chain carbonyl of the substrate/inhibitor. The Zn2+ ion and H142 activate a water molecule for a nucleophilic attack on the carboxyl of the acetylated lysine side chain and the reaction intermediate has been proposed to be stabilised by Y306, H142, H143, and Zn2+. Subsequently, the hydrolysis is completed by H143 that protonates the lysine side chain.27
Following hydrolysis, the acetate product has been proposed to be released via an internal channel adjacent to the substrate-binding tunnel.21,28 This ‘release channel’ is lined by residues I34, F152, Y306 from the top and Y18, I19, Y20, H42, N136 from the bottom, Fig. 1c and d, with Y306, W141 and F152 located at the intersection between the substrate-binding and acetate release channel. The internal acetate release channel has been hypothesised to divide into to sub-channels with R37 and S138 present at the junction and R37 playing an important role in stabilising the negatively charged acetate for effective release. Despite the availability of many crystal structures of HDAC8 in complex with different substrates and inhibitors,21–26 a high-resolution structure of the unbound form of HDAC8, or any other class I HDACs, is not yet available. Also, many aspects of HDAC8 regulation and function still remain elusive,29–32 in particular the mechanism of product release as well as the role of the dynamic loops, L1–L7, for function and how these are connected to the active site.
In this work we performed 10 μs-long unbiased molecular dynamics simulations of wild-type HDAC8 and two mutants I19S and S39E, which are known to allosterically downregulate activity. We observe that dynamics about the χ1 dihedral angle of aromatic residues near the active site, Y306, H143, and F152, orchestrate the conformational sampling of the functional loops of HDAC8. In both mutants, I19S and S39E, the conformational sampling of the functional loops is significantly affected. In I19S-HDAC8 the dynamics of Y306 and the L6 loop is completely inhibited and a ∼20 Å allosteric path from I19 to Y306 is identified, which passes through several sites that have been identified to be essential for activity. The S39E mutant leads to a substantial stabilisation of all the aromatic side chains near the active site, which effectively traps each of the functional loops in one conformation. Overall, the results show how the active site of HDAC8 is connected to the conformational sampling of functional loops and also how mutations affect the sampled conformations and activity.
The I19S and S39E simulations were prepared similar to the wild-type simulations after introducing the mutation. In order to obtain a simulation for the reverse mutation, I19Sr, the final frame from the I19S simulation was taken and S19 was mutated back to isoleucine. The waters were removed and the model was prepared as described above. Finally, the wild-type HDAC8, S39E-HDAC8 simulations were run for 10 μs, the I19S-HDAC8 simulation for 8 μs, whereas the I19Sr-HDAC8 simulation was run for 2.5 μs. The lengths of the simulations were tailored to the characteristic period of the observed conformational changes and the convergence of the sampling was also monitored through Markov state analysis. An equilibrium period of 2 μs, 1.8 μs, and 2.5 μs was used for the wild-type, S39E-HDAC8, and I19S-HDAC8 simulation, respectively. This equilibrium period was not used for any of the analyses presented and the length of the equilibrium period was determined based on the RMSD to the start structure, Fig. S1.† The loop, between K202 to V217 was not considered in the presented analysis, because convergence of this loop could not be fully justified.
HT2 = H{φ,ψ}2 + Hχ12 |
A Hellinger divergence was only considered significant if it was larger than the average divergence obtained for a block-analysis, {{0 μs, 2.5 μs}, {2.5 μs, 5.0 μs}, {5.0 μs, 7.5 μs}, {7.5 μs, 10.0 μs}}, of the wild-type simulation.
All analyses was carried out using MDAnalysis,46 and Python 3.6 or Python 3.8 with the Numpy and Matplotlib libraries.
The major gauche− state (χ1 ≈ 300°; g−) is characterised by a ∼6.2 Å distance between Y306-Oε and the active Zn2+ ion, whereas the minor gauche+ state (χ1 ≈ 60°; g+) is characterised by a substantial longer Y306-Oε to Zn2+ distance of ∼11.5 Å and the tyrosine side chain pointing towards the L6 loop, Fig. 2b. The flips between the gauche+ and gauche− states of Y306 are strongly correlated with a 5-to-7 Å displacement of the L6 loop, wherein M274 is located, Fig. 2c. Loop L6 and M274 form a part of the substrate-binding tunnel and in the Y306(g+) state the tunnel is substantially wider and solvent exposed, thus showing that the gauche− to gauche+ flip of Y306 acts as a switch to open the substrate-binding tunnel, Fig. 3. Further examination of the Y306(g+) state shows that it is closely associated with a hydrogen bond between Y306-OH and I269-CO, Fig. 2b, which implies that this hydrogen bond may stabilise the Y306(g+) state. Moreover, in the Y306(g+) state a salt-bridge between D233 and R353 is fully formed, whereas in the Y306(g−) state this salt-bridge is only partly formed, Fig. S3.† The D233–R353 salt bridge, in turn, stabilises the L6 loop in the Y306(g+) state via additional charge–charge interactions between R353 and Cys275-CO as well as between R353 and Pro273-CO. In addition to the Y306(g−) and Y306(g+) states, there is also a very low-populated and short-lived state with the Y306 χ1 in a trans state (∼0.5% population), where the Y306 side chain points towards the L1 loop.
Of particular interest is that the substrate-binding tunnel is substantially widened and the active site exposed to the bulk solvent when the L6 loop is distant to the Zn2+, Fig. 3 and S4.† Such a state, with the exposed active site and the widened substrate-binding tunnel, will likely be able to rapidly release products. A movie for visually comparing the structure in the “closed” and “open” states, and with a transition from surface to cut-away, is shown in ESI Movie S1.†
Evidence of conformational flexibility around Y306 was previously observed in a 10 ns simulation of Y306F-HDAC8 (ref. 47) and the gauche+ and trans conformations have been observed in crystal structures of related enzymes, thereby providing evidence for the states of Y306 observed in the unbiased simulations of HDAC8 presented here. Specifically, in the crystal structure of acetylpolyamine amidohydrolase (APAH; PDB 3Q9F), gauche− and trans were observed for the corresponding residue, while in H976Y-HDAC4 the gauche− (PDB: 2VQV) and the gauche+ (PDB: 2VQO) were observed.48,49
The active-site residue H143 is involved in the hydrolysis by forming a hydrogen-bond with the acetylated lysine Nε and thus stabilising the reaction intermediate. As observed for Y306, H143 also exists in a dynamic equilibrium between two conformations that are characterised by different χ1 angles, H143(g−) with χ1 ≈ 300° and H143(t) with χ1 ≈ 205°, Fig. 4.
In all crystal structures of HDAC8, bound to a substrate or inhibitor, only the H143(g−) is present.21–26 When the H143 side chain flips from gauche− and trans, the entire L2 loop region (residues 86–103), as well as the L3 loop located between L2 and the active site (residues 142–154), move approximately 1 Å closer to the active site Zn2+, Fig. S5.† The functional L2 loop region contains D101, which is known to stabilise substrate binding.50,51 Moreover, a correlation is observed between the sampling of the χ1 of Y306 and H143, Fig. 4b. Specifically, when the Y306(g+) state is formed, H143 χ1 is stabilised in the trans conformation.
To gain further insight into how the state of Y306 influence the dynamics of H143, and vice versa, a Markov state model was generated from the 10 μs simulation, Fig. 4c (see Methods). The generated Markov model highlights the stabilisation of H143(t) when Y306 χ1 is in a gauche+ conformation, a stabilisation that appears to be mainly caused by a slower rate from H143(t) to H143(g−) and a slightly faster rate from H143(g−) to H143(t). Also, the dynamics of Y306 is essentially frozen when H143 χ1 is in the trans state. Taken together these observations suggest that the motions of the two functional residues, Y306 and H143, greatly influence each other and the motion of each of these aromatic side chains correlate with the conformations of the functional loops.
The χ1 motion of F152 is substantially faster than those of H143 and Y306, with F152(g−) and F152(t) having an average life-time of 1.4 ns and 1.1 ns, respectively, during the 10 μs simulation, Fig. 5. Thus, there is a slight preference for the F152(g−) state (58%) over the F152(t) state (42%). Moreover, the F152(t) state is stabilised twice, for approximately 1 μs each time, Fig. 5a. The F152(g+) state is only present at 0.3%.
Fig. 5 Side-chain motion of F152 correlates with L1 sampling: (a) the F152 χ1 angle as a function of simulation time (grey) and sampled every nanosecond; the red line is a moving average over a window of 20 ns. The F152 side chain samples the gauche− (χ1 ≈ 300°; g−) and the trans (χ1 ≈ 200°; t) states each with an average life-time of about 1 ns. (b) Two-dimensional histogram showing that the χ1 dihedral angle of F152 correlates strongly with the distance between K33 Cα within the L1 loop and F152 Cα; a distance which has previously been used as a measure of the L1 conformation.17 (c) Two-dimensional histogram showing that the distance between the side chains of K33 in L1 and D87 in L2 as well as ability to form a salt-bridge between the two side chains, correlate with the χ1 of F152. In the F152(t) state, a salt-bridge between L1 and L2 is unlikely (5%), whereas in the F152(g+) state there is a 22% chance of a salt-bridge being formed. |
The F152 residue is found in the gauche− conformation in all human HDAC8 crystal structures solved so far, however, in the HDAC8 homolog from Schistosoma mansoni the corresponding phenylalanine residue, F151, is observed in a trans state,52 see Fig. S6,† thereby providing evidence that the substantial amount of trans conformation observed in the unbiased simulations is realistic for HDAC8 without a substrate or an inhibitor bound.
The correlation between the conformation of the L1 loop and the χ1 angle of F152 is striking, Fig. 5b. When F152 is in the trans state, the L1 loop moves approximately 3 Å away from the substrate-binding tunnel leading to an open state. In line with this movement is that the salt-bridge between K33 in L1 and D87 in L2, previously related to function,50 predominantly forms in the F152(g−) state. These findings agree with the crystal structure of the HDAC8 homolog from Schistosoma mansoni, where the corresponding phenylalanine is in trans, and the L1 loop is ca. 1 Å further away from the active site compared to human HDAC8.
Despite the motion about the F152 χ1 angle being fast, the conformation of other functional sites is correlated with the state of the F152 side chain. One example is R37, Fig. S7,† which has previously been suggested to play a central role for the release of acetate. In all HDAC8-inhibitor or substrate complexes, the guanidinium group of R37 forms hydrogen bonds with the backbone carboxyl oxygen atom of two conserved glycine residues (G303 and G305) in the glycine rich loop adjacent to Y306.21–26 These hydrogen bonds are formed throughout the HDAC8 simulation, however motion around the R37 χ1 angle was observed to be affected by the F152 χ1 sampling, Fig. S7.† Despite the ∼10 Å distance between F152 and R37, both residues are in proximity of the L1 loop, suggesting that the motions of F152 and R37 could be correlated with each other via the L1 loop.
In crystal structures of HDAC8–ligand complexes, both Y306 and F152 are present at the junction of the substrate-binding tunnel and the acetate release channel, and changes in these two residues are important for the opening of the acetate release channel.21,28 Despite this, only a weak correlation is observed between Y306 χ1 and F152 χ1 as well as between Y306 and R37, Fig. S8.† Among all the residues in the vicinity of the substrate-binding tunnel, the active site, and the product release channel, χ1 flips are effectively only observed for the four residues described above, i.e., R37, H143, F152, and Y306. For all the other residues near the active site one χ1 conformation dominates, see Fig. S9.†
In the S39E-HDAC8 simulation, the dynamic sampling of Y306 χ1 and the L6 loop is inhibited and effectively only sample the gauche− state and the proximal conformation, respectively, Fig. 6a. Moreover, the F152 χ1 angle is highly stabilised in the trans state throughout the S39E simulations, which in turn leads to L1 conformations consistently distant to the active site and a very low probability of forming the K33 Nζ–D87 Oδ salt-bridge, Fig. 6b–f. A Markov state model similar to Fig. 4c cannot be derived from the S39E-HDAC8 simulation, because the Y306(g+) state is not populated and the Y306(t) state is only present in a few frames. With regards to the kinetics observed for H143 χ1, the H143(g−) and H143(t) states have average lifetimes of 0.023 μs and 0.048 μs, respectively, during the S39E-HDAC8 simulation compared to 0.020 μs and 0.024 μs in the simulation of wild-type HDAC8.
Fig. 7 Simulation of I19S-HDAC8 and comparison with wild-type. (a) The Y306 χ1 angle as a function of simulation time (grey) and sampled every nanosecond in the I19S-HDAC8 simulation. The cyan line is a moving average over a window of 20 ns. Only the Y306(g−) and the Y306(t) states are sampled, with probabilities of 99.9% and 0.1%, respectively. (b) Two-dimensional histogram showing the correlation between the distance between M274 Cα within the L6 loop and Zn2+, and Y306 χ1. The L6 loop only samples the conformation near the active site. (c) Residue-specific Hellinger distances between the {φ,ψ,χ1} distribution in the I19S and wild-type simulations shown on the structure of HDAC8 (PDB: 2V5W).22 A clear path between the site of mutation, I19, and Y306 is seen. (d) A slightly higher preference for a R37 χ1trans conformation in the I19S simulation. (e) and (f) Very similar probability distributions in the I19S- and wild-type simulation for the H143 and F152 χ1. |
I19 and Y306 are more than 20 Å apart, both in the simulations and in all published structures of HDAC8. A central question is therefore how the I19S mutation leads to an inhibition of the χ1 motions of Y306 as well as the motions of the L6 loop. To gain insight into the differences between the I19S and wild-type simulations, Hellinger distances were calculated between the backbone {φ,ψ} and side-chain χ1 distributions in the I19S and wild-type simulations, for each residue. A block analysis was also performed on the wild-type simulation and a Hellinger distance (I19S to wild-type) was only considered significant if it was significantly larger than what was obtained in the block analysis of the wild-type simulations (see Methods). When the significant Hellinger distances (I19S versus wild-type) are shown on the structure of HDAC8, Fig. 7c, a clear path is observed between residue 19 and Y306. It is important to note here that this path appears naturally and has not in any way been imposed during the calculation of the Hellinger distances. The observed path goes through the central β-sheet of HDAC8, including residues I135 and L299, from where the path runs via the glycine loop, G302–G305 to Y306. The glycine loop, G302–G305, has previously been implicated in function and mutating any of these glycine residues to alanine leads to a small change in the loop conformation in the resulting crystal structures in bound form, but causes a substantial effect on HDAC8 activity, with only about 0.3–3% of enzymatic activity retained.47 Thus, not only do the Hellinger distances provide insight into how the effect of the I19S mutation is transmitted though the structure of HDAC8, but also highlight areas where mutations and regulator binding could alter HDAC8 activity due to a change in the L6 dynamics.
Arginine R37 forms interactions with the glycine loop, G302–G305. Although there is no significant change in the hydrogen-bonding between the guanidinium group of R37 and the glycine loop in the I19S-HDAC8 and wild-type HDAC8 simulations, the distribution of R37 χ1 changes slightly between the two simulations, Fig. 7d, consistent with the allosteric path going through the glycine loop. In contrast, only limited changes are observed for the distribution of the side chain χ1 of H143 and F152 and thus of the L1 and L2 loops, Fig. 7e and f. The motion about the F152 χ1 dihedral angle is slightly faster in the I19S-simulation, where the F152(g−) has a lifetime of 1.0 ns compared to 1.4 ns in the wild-type simulation.
A substantial change in the dynamics of the L1 and L6 loops was observed upon investigation of a CdLS down-regulating mutant, I19S, as well as a phosphorylation mimicking mutation, S39E, of HDAC8. In the S39E-HDAC8 simulation all the aromatic side chains near the active site were highly stabilised in one conformation leading to a consistent open form of the L1 loop and a closed conformation of the L6 loop. Recently, Michaelis–Menten parameters, kcat and KM, were reported for the S39E mutant of HDAC8, where changes in both kcat (∼3 fold decrease) and KM (∼2 fold increase) were observed.53 A change in KM upon mutation typically reflects a change in the substrate-binding affinity, while a change in kcat reflects a change in the rate of hydrolysis and/or a change in the rate of product release. These experimental data therefore suggests that mutating S39 to glutamate affects both substrate-binding as well as hydrolysis and/or product release. In the MD simulations presented above, changes in sites related to both substrate-binding (L1; F152), hydrolysis (H143), and product release (L6; Y306, M274) were observed, in line with the experimental findings.
In the simulation of I19S-HDAC8, the backbone conformation of G302–G305 changed slightly and the L6 loop as well as Y306 were trapped in a closed conformation. In available experimental enzyme kinetics data, mutation of Y306 to phenylalanine resulted in a two-fold decrease in the KM and a 290-fold decrease in the catalytic constant, kcat.54 Similarly, upon mutation of glycine residues adjacent to Y306 to alanine, G304A and G305A, a ∼300 and a ∼38 fold decrease in kcat was observed, respectively, with only a 2 fold and a 1.5 fold increase in KM.47 Therefore, the predominant role of Y306, whose conformational sampling is affected by the adjacent glycine residues, is in stabilising the enzymatic transition state and/or in the release of products, kcat, as opposed to substrate-binding, KM.54 In the crystal structure of the G304A and G305A mutants, in complex with an inhibitor, only a minor change in the position of Y306 is observed relative to wild-type structures. Importantly, in these structures there are no significant change in the distance between Y306-Oε and the carboxyl of the inhibitor: 2.4 Å in wild-type HDAC8 and 2.6 Å in G304A-HDAC8 and G305A-HDAC8. This strongly suggests that these mutations have limited effect on transition-state stabilisation. Hence, the glycine-to-alanine mutations likely change kcat by impeding the release of product, which in turn happens via a change in the conformational sampling of Y306. Based on comparison of these experimental biochemical data with the results of the molecular dynamic simulations, it is reasonable to hypothesise that the downregulation caused by the I19S mutation is due to an inhibition of the L6 dynamics and a slight change in the sampling of R37, which results in the inability of I19S-HDAC8 to release the products.
During the wild-type HDAC8 simulation the Y306(g−) state, which has the substrate-binding tunnel formed, is the most stable state with an average lifetime of 1.3 μs. Within the Y306(g−) state, R37, H143 and F152 are in rapid dynamic equilibrium between their gauche− and trans states and we therefore propose that the Y306(g−) state is a substrate-binding-competent state. In all crystal structures of HDAC8 bound to substrate or inhibitor, R37, H143, F152 and Y306 are all found in the gauche− state. A comparison of the HDAC8 frames with the side chain χ1 of R37, H143, F152 and Y306 in the gauche− state with the crystal structures shows some differences for residues forming the active site, the substrate-binding tunnel and the release tunnel. This includes the orientation of W141 and Y306, Fig. S11.†
However, it is very likely that these side chains move slightly in the presence of substrate, to stabilise the substrate for catalysis, including forming a hydrogen-bond between Y306 and the substrate carboxyl oxygen. The Y306(g+) state can be considered as a product-release state, since the substrate-binding tunnel is widely exposed (see ESI Movie S2†). This agrees with the experimental enzyme kinetics data available for Y306F-HDAC8, where only the catalytic rate, kcat, is substantially affected whereas binding of substrate is only slightly affected. A concomitant movement of R37 and F152 χ1 to trans in the Y306(g+) state facilitates the release of acetate. The structure of the binding tunnel in the three states: apo-state, bound-state, and release-state are shown in Fig. 8.
Fig. 8 Change in structure of binding tunnel of HDAC8 during catalysis. The cut-away surface views of representative structures of HDAC8 in the apo-state: {R37(g−), F152(g−), Y306(g−)} (top), release-state: {R37(t), F152(t), Y306(g+)} (bottom left), and bound-state: crystal structure 2V5W22 (bottom right). |
Footnotes |
† Electronic supplementary information (ESI) available: Supporting Fig. S1–S11, supporting Movies S1 and S2. See DOI: 10.1039/d1sc01929e |
‡ These authors contributed equally to this work. |
This journal is © The Royal Society of Chemistry 2021 |