The Dynamic Conformational Landscapes of the Protein Methyltransferase SETD8

Elucidating conformational heterogeneity of proteins is essential for understanding protein functions and developing exogenous ligands for chemical perturbation. While structural biology methods can provide atomic details of static protein structures, these approaches cannot in general resolve less populated, functionally relevant conformations and uncover conformational kinetics. Here we demonstrate a new paradigm for illuminating dynamic conformational landscapes of target proteins. SETD8 (Pr-SET7/SET8/KMT5A) is a biologically relevant protein lysine methyltransferase for in vivo monomethylation of histone H4 lysine 20 and nonhistone targets. Utilizing covalent chemical inhibitors and depleting native ligands to trap hidden high-energy conformational states, we obtained diverse novel X-ray structures of SETD8. These structures were used to seed massively distributed molecular simulations that generated six milliseconds of trajectory data of SETD8 in the presence or absence of its cofactor. We used an automated machine learning approach to reveal slow conformational motions and thus distinct conformational states of SETD8, and validated the resulting dynamic conformational landscapes with multiple biophysical methods. The resulting models provide unprecedented mechanistic insight into how protein dynamics plays a role in SAM binding and thus catalysis, and how this function can be modulated by diverse cancer-associated mutants. These findings set up the foundation for revealing enzymatic mechanisms and developing inhibitors in the context of conformational landscapes of target proteins.


Introduction
Proteins are not static, but exist as an ensemble of conformations in dynamic equilibrium 1 .
Characterization of conformational heterogeneity can be an essential step towards interpreting function, understanding pathogenicity, and exploiting pharmacological perturbation of target proteins 2-4 . Conventional efforts to map functionally relevant conformations rely on biophysical techniques such as X-ray crystallography 5 , nuclear magnetic resonance (NMR) 6 , and cryoelectron microscopy 7 , which provide static snapshots of highly-populated conformational states.
While complementary techniques such as relaxation-dispersion NMR can resolve a limited number of low-population states, they are incapable of providing detailed structural information 8 .
By contrast, molecular simulations provide atomistic detail---a prerequisite to structure-guided rational ligand design---and insight into relevant conformational transitions 1 . The emergence of Markov state models (MSMs) has shown the power of massively distributed molecular simulations in resolving complex kinetic landscapes of proteins 9,10 . By integrating simulation datasets with MSMs, functionally relevant conformational dynamics as well as atomistic details can be extracted 10 . Recently, MSMs have been used to identify key intermediates for enzyme activation 11,12 and allosteric modulation 13 . However, these approaches are limited by the number of seed structures and timescales accessible by molecular simulations (generally microseconds for one structure) relative to the reality of complicated conformational transitions (up to milliseconds for multiple structures) 14 . To overcome the limitations of individual techniques, we envisioned an integrated approach that combines simulation with experiment to characterize conformational landscapes of enzymes and elucidate their functions with the consideration of dynamic conformations.
Protein lysine methyltransferases (PKMTs) comprise a subfamily of posttranslational modifying enzymes that transfer a methyl group from the cofactor S-adenosyl-L-methionine (SAM) 15 . PKMTs play epigenetic roles in gene transcription, cellular pluripotency, and organ development 16,17 . Their dysregulation has been implicated in neurological disorders and cancers 18,19 . SETD8 (SET8/Pr-SET7/KMT5A) is the sole PKMT annotated for monomethylation of histone H4 lysine 20 (H4K20me) 20,21 and many non-histone targets such as the tumor suppressor p53 and the p53-stabilizing factor Numb 22,23 . Disruption of endogenous SETD8 leads to cell cycle arrest and chromatin decondensation, consistent with essential roles for SETD8 in transcriptional regulation and DNA damage response [24][25][26] . SETD8 has also been implicated in cancer invasiveness and metastasis 27 . High expression of SETD8 is associated with pediatric leukemia and its overall low survival rate 28 . As a result, there is enormous interest in elucidating functional roles of SETD8 in disease and developing pharmacological agents to perturb this target 29-31 .
Given the essential roles of conformational dynamics in enzymatic catalysis 1,32 and our current limited knowledge of conformational landscapes of PKMTs, we envisioned leveraging an integrated experimental-computational approach to characterize dynamic conformational landscapes of SETD8 and its cancer-associated mutants with atomic resolution. To access previously-unseen, less-populated conformational states of SETD8 to seed massively parallel distributed molecular dynamics (MD) simulations, we envisioned trapping these conformations with small-molecule ligands. Here we solved four distinct crystal structures of SETD8 in alternative ligand-binding states with covalent SETD8 inhibitors and native ligands. With the aid of these new structures, we generated an aggregate of six milliseconds of explicit solvent MD simulation data for apo-and SAM-bound SETD8. Using a machine learning approach to select features and hyperparameters for MSMs via extensive cross-validation, we identified 24 kinetically distinct metastable conformational states of apo-SETD8 and determined how the conformational landscape is remodeled upon SAM binding. We then validated these conformational landscapes with stopped-flow kinetics and isothermal titration calorimetry by examining SAM binding, characterizing rationally-designed SETD8 variants with increased catalytic efficiency, and resolving multiple timescales associated with transitions among these conformers. The resulting model furnishes unprecedented key insights on how these dynamic conformations play a role in catalysis and how cancer-associated SETD8 mutants alter this process.

Results
Crystal structures of SETD8 associated with hidden conformations. To identify hidden high-energy conformational states of SETD8, we envisioned a strategy of trapping the associated conformers with small-molecule ligands. The development of high-affinity SETD8 inhibitors with canonical target-engagement modes is challenging 29 , and led us to exploit covalent inhibitors 31, 33 . These compounds can overcome the high energy penalties associated with hidden high-energy conformers through the irreversible formation of energetically-favored inhibitor-SETD8 adducts. Our prior efforts led to the development of covalent inhibitors containing 2,4diaminoquinazoline arylamide and multi-substituted quinone scaffolds by targeting Cys311 31, 33 .
Upon further optimization of these scaffolds, we identified MS4138 (Inh1) and SGSS05NS (Inh2) 34 , two structurally distinct covalent inhibitors with the desired potency against SETD8 (Figures 1a, S9). X-ray crystal structures of SETD8 were then solved in complex with Inh1 and Inh2, respectively (Figures 1b,c, S10, S11). Notably, despite the overall structural similarity of the pre-SET, SET, and SET-I motifs, the Inh1-and Inh2-SETD8 binary complexes (BC-Inh1 and BC-Inh2) differ from the SETD8-SAH-H4 ternary complex (TC) [35][36][37] by the distinct conformations of their post-SET motifs. The post-SET motif of TC was characterized by its Ushaped topology with a double-kinked loop-helix-helix architecture, which appears to be optimally oriented for binding both SAM and a peptide substrate (Figure 1c,d) [35][36][37] . In comparison, BC-Inh1 and BC-Inh2 rotate their post-SET motifs by 140 ° and 60 °, respectively ( Figure 1d). Moreover, the post-SET motifs of BC-Inh1 and BC-Inh2 adapt more extended configurations with a less structured loop and a singly-kinked helix, respectively (Figure 1c,d).
Whereas multiple factors may influence the overall conformations, the formation of Cys311 adducts likely made the key contribution to the discovery of these hidden post-SET motif conformers. To reveal additional hidden conformers that are structurally distinct from TC, we also solved crystal structures of SETD8 upon depleting native ligands and obtained structures of the SAM-SETD8 binary complex (BC-SAM) and apo-SETD8 (APO) (Figures 1c, S12, S13). Strikingly, BC-SAM and APO differ from TC by their distinct SET-I motifs in the context of the otherwise similar SET-domain (Figure 1d). Furthermore, the post-SET motif of APO structurally resembles an intermediate state between BC-Inh1 and BC-Inh2 but is distinct from those of BC-SAM and TC (Figure 1d). In contrast to the structurally diverse SET-I (I1-3) and post-SET motifs (P1-4) in these structures, their pre-SET motifs show only slightly altered configuration (Figure 1d). The differences between these structures highlight the conformational plasticity of the SET-I and post-SET motifs. Collectively, these observations provide strong structural rationale for the existence of a highly dynamic conformational landscape of SETD8.
Hidden conformations of apo-SETD8 revealed by structural chimeras. The BC-SAM, BC-Inh1, BC-Inh2, APO, and TC structures can be readily classified into three distinct SET-I configurations (I1-3) and four distinct post-SET configurations (P1-4) (Figure 1d). Given the relative independence between the SET-I and post-SET motifs, we expected that additional combinations of discrete motifs can represent yet-unobserved functionally relevant conformations of SETD8. We thus constructed putative "structural chimeras" of apo-SETD8 containing orthogonal I1-3 and P1-4 in a combinatorial (3×4) manner (Figures 2a, S14). Among the twelve structural chimeras as potential seeds for MD simulations, five were crystallographically-determined conformers (BC-Inh1, BC-Inh2, BC-SAM, TC with ligands removed, and APO), four were new structurally-chimeric conformers, and three were excluded because of obvious steric clashes (Figures 2a, S15). Combinatorial construction of structural domain chimeras using crystallographically-derived post-SET and SET-I conformations. Each conformer is boxed and color-coded with black for five X-ray-derived structures, blue for four putative structural chimeras included as seed structures for MD simulations, and grey for three structural chimeras excluded from MD simulations because of obvious steric clashes. b, Schematic workflow to construct dynamic conformational landscapes via MSM. The five X-ray structures and the four structural chimeras were used to seed massively parallel MD simulations on Folding@home (see Method). Markov state models were constructed from these MD simulation results to reveal the conformational landscape. c-e, Kinetically metastable conformations (macrostates) obtained from kinetically coupled microstates via Hidden Markov Model (HMM) analysis. The revealed dynamic conformational landscapes consist of 24 macrostates for apo-SETD8 (left panel) and 10 macrostates for SAM-bound SETD8 (right panel). c, Kinetic and structural separation of macrostates in a 3D scatterplot.
The X, Y axes represent kinetic separation of macrostates with a log-inverse flux kinetic embedding method (see Methods). The Z axis reports RMSDs of each macrostate to APO (left) or BC-SAM (right).
The relative population of each macrostate of apo-or SAM-bound SETD8 ensembles is proportional to the volume of each representative sphere. d, Cartoon depiction of macrostates in a 2D scatterplot. The relative positions of metastable conformations were derived via the log-inverse flux kinetic embedding (see Methods). The diameter of the corresponding circle in the 2D scatterplot is proportional to the diameter of the respective sphere in the 3D scatterplot above. Equilibrium kinetic fluxes larger than 7.14×10 2 s -1 for apo-and 1.39×10 3 s -1 for SAM-bound SETD8 are shown for interconversion kinetics with thickness of the connections proportional to fluxes between two macrostates. e, Chord diagrams and representative conformers of macrostates. The colors are encoded for the free energy of each macrostate relative to the lowest free energy of the macrostates. The equilibrium flux between two macrostates is proportional to thickness of respective arcs.
Dynamic conformational landscape of apo-SETD8 via Markov state modeling from 5ms MD simulation dataset. With seed conformations prepared as above, we envisioned illuminating the conformational landscape with massively distributed long-time MD simulations and resolving its kinetic features with Markov state models (MSMs) (Figures 2b, S14). We conducted approximately 500×1 µs explicit-solvent MD simulations from each seed and accumulated 5 milliseconds of aggregate data in 10 million conformational snapshots for apo-SETD8 (Figures S16, Table S3). To identify functionally relevant conformational states and their transitions, we built MSMs using a pipeline that employs machine learning and extensive hyperparameter optimization to identify slow degrees of freedom and structural and kinetic criteria to cluster conformational snapshots into discrete conformational states (Figures S17-24, Tables S4, S5) 38 . This approach identified 24 kinetically metastable conformations (macrostates) from an optimized, cross-validated set of 100 microstates (Figures 2c, S25-30, Tables S6, S7).
These macrostates are remarkably diverse, spanning up to 10.5 Å Cα RMSD from APO. To visualize the kinetic relationships between functionally important conformations, dimensionality reduction was used to project the landscape into 2D while preserving log inverse fluxes between states (Figure 2d). The relative populations of these macrostates and their interconversion kinetics were calculated on the basis of their transition fluxes, resolving rare conformational states up to 6 kT in free energy (Figure 2d,e).
The dynamic conformational landscape of SAM-bound SETD8. Given the success in constructing the dynamic conformational landscape of apo-SETD8, we applied the same strategy to SAM-bound SETD8. With the two crystal structures of SETD8 in complex with SAM (BC-SAM and TC) as the seed conformations, we conducted ~ 500×1 µs explicit-solvent MD simulations from each structure and accumulated 1 millisecond of aggregate data (2M snapshots) ( Figure S25). The resulting MSM for SAM-bound SETD8 contained 10 kinetically metastable macrostates arising from 67 microstates ( Figure S31, Tables S8, S9). Similar to those of apo-SETD8, the relative macrostate populations of SAM-bound SETD8 and their flux kinetics were computed and embedded into 3D/2D scatter plots and chord diagram (Figure 2c,d,e). The smaller number of metastable states identified for SAM-bound SETD8 is expected given that SAM binding restricts conformational accessibility.  Fluorescence changes of wild-type and K382P SETD8 traced with a rapid-quenching stopped-flow   instrument within 1 s upon SAM binding. e, SAM-binding ITC enthalpogram of wild-type and K382P SETD8. f, Stepwise SAM-binding of SETD8 in the integrative context of biochemical, biophysical, structural, and simulation data. ITC determines the thermodynamic constant of SAM binding by SETD8.
MD simulations and MSM uncover metastable conformations and interconversion rates of apo-and SAM-bound SETD8 (K apo and K SAM ). Stopped-flow experiments revealed that SETD8 binds SAM via biphasic kinetics. Rate constants uncovered by stopped-flow experiments (k1, k-1, k2, k-2) represent macroscopic rates of SAM binding by SETD8 with multiple metastable conformations. The microscopic behavior of individual metastable states and corresponding rates (k1, k-1, k2, k-2) have not been resolved.
Transition probability matrices (red) and microscopic rate constant matrices (blue) are shown as colored grids. A rigorous mathematical model for these processes is shown in Figure S36. g, Kinetic and thermodynamic constants of wild-type SETD8 and its mutants. For k1, k-1, k2, k-2, data are best fitting values ± standard error (s.e.) from KinTek. For Kd-ITC, data are mean ± s.e. of at least 3 replicates.  (Figure 3d). In the context of the conformational landscape of apo-SETD8, we interpreted the major decrease in fluorescence intensity (fast-phase kinetics) as a consequence of the collective changes of Trp390 from the solvent-exposed hydrophilic environment in apo conformations to the hydrophobic environment in SAM-bound conformations (Figure 3a,c). In contrast, the minor decrease in fluorescence intensity (slow-phase kinetics) reflects the slow conformational changes of Trp390 in the SAM-bound SETD8 conformational ensembles (Figure 3d). With unsupervised global fitting to this two-step model, we obtained forward and reverse rate constants for the fast-and slow-phase kinetics, which are in agreement with conventional fitting to double exponential kinetics 40 (Figures 3d,f,g, S32, Table S10). The k-1 value was also confirmed independently by rapid-mixing stopped-flow dilution of SAM-bound SETD8 41 ("ES+E'S", Figure S33, Table S10). Here the k-1/k1 ratio of 309±6 µM corresponds to the average SAM dissociation constant Kd1 of apo-SETD8 conformers, which is consistent with independently determined ITC Kd of 251±16 µM (Figures 3e,f, S34). In contrast, the large k-2/k2 ratio of 30±11 suggests that the second phase corresponds to a slow equilibrium between ES and E'S with minimal contribution of E'S to the overall SAM dissociation constant Kd (Figure 3e).
The conformational ensembles we identified for apo-and SAM-bound SETD8 demonstrate the statistical nature of its SAM-binding process. Therefore, the observed fluorescence changes and herein determined macroscopic kinetic constants represent an ensemble-weighted average of microscopic behaviors of all species that exist in the solution. A rigorous mathematical description of microscopic kinetics of SAM binding was thus obtained under the consideration of interconversion of the metastable conformational states of apo-and SAM-bound SETD8 ( Figure   S36).
We then proposed to confirm our understanding of functionally-relevant conformations and their thermodynamics by identifying SETD8 variants with increased affinity for SAM. We uncovered a collection of characteristic kink motifs around Lys382 in the post-SET motif of SAM-bound SETD8 conformational ensembles (Figure 3c), while this region is less structured in apo-SETD8 conformational ensembles. We hypothesized that a proline mutation (K382P) could better stabilize the conformational ensembles of SAM-bound SETD8 than apo-SETD8 (Figure 3c,h). We also identified a characteristic a-helix in the SET-I motif, which adapts flexible and diverse configurations in apo ensembles but constrained and structurally distorted configurations in SAM-bound ensembles (Figure 3c). We proposed that the replacement of I293 or E292 adjacent to the a-helix with a flexible glycine should relax this distortion to better stabilize SAM-bound ensembles (Figure 3c,h). We therefore characterized the SAM-binding kinetics and affinities of K382P, I293G, and E292G variants of SETD8 with stopped-flow kinetics and ITC (Figures 3c,d,e,f, S32-34). While exhibiting biphasic kinetics similar to that of wild-type SETD8, the stopped-flow mixing experiment revealed the three variants showed a significant two-fold decrease of Kd,SAM (Figure 3d,e). The stopped-flow data further revealed that the two-fold change of Kd,SAM mainly arises from increased SAM-binding rates k1 with relatively unchanged k-1 (Figure 3g). These results are consistent with independently-determined Kd and k-1 from ITC and stopped-flow dilution, respectively (Figures 3e,f, S33, S34, Table S10).
Collectively, these observations confirm the robustness of our conformational landscape model for apo-and SAM-bound SETD8.

Effects of key simulation parameters on construction of conformational landscapes.
We  Table   S14). These findings argue for the importance of using multiple structures to construct the 14 landscape within achievable computer time. The seed conformations prepared from ligandtrapped SETD8 structures are essential to discovering the complete conformational landscapes of SETD8.
For simulation time, we observed that the fewer seed conformations of apo-SETD8 were employed, the more computing power (the product between the number of simulation trajectories and the time length per trajectory) was required to reach a comparable level of microstate coverage (Figure 4e, Tables S16, S17). When computing power is fixed, comparable microstate coverages of apo-and SAM-bound SETD8 can be obtained by running either multiple short trajectories or few long trajectories (Figures 4f, S38). The current simulation time (5 ms for apo-SETD8 and 1 ms for SAM-bound SETD8) provides 2-10-fold redundant computing power to map the conformational landscapes of SETD8.  The hub-like macrostate A11 consists of two structurally-distinct microstates with comparable populations (Figures 2d, 5a). One microstate structurally resembles the conformation of APO (I3P3), while the other microstate represents a conformer with the I1P23 feature for its SET-I and post-SET motifs (Figure 5b, Table S6). Rapid conformational interconversions within A11 is consistent with its hub-like character, centered between the two lobes of the dumbbell-like network. Interestingly, macrostates kinetically adjacent to A11 have structurally similar SET-I motifs within each lobe but distinct SET-I motifs between the two lobes (I2~3 for the left and I1~2 for the right) (Figures 2d, 5b). Therefore, A11 is a transitiontype state essential for the conformational fluxes of the macrostates between the two lobes, involved in a key step of conformational changes of the SET-I motif between I1~2 and I2~3.
Conformers in the macrostates A22, A24 and A20 are structurally similar to TC and BC-SAM with slightly different but well-defined SAM-binding pockets, suggesting minimal conformational reorganization of A22, A24, and A20 is required to accommodate the cofactor (Figure 5a,b,e). Interestingly, A22 and A24, whose overall structures are similar to each other (TC-like), rarely interconvert in the apo landscape (Figure 2d). In contrast, the basin-like macrostates (A5, A10), A7, A8, (A12, A13, A16) and A15 do not contain a well-defined SAMbinding pocket (Figure 5a,b,e). Here the conformers in macrostate A12 are similar to APO, the conformers in the macrostate A6 are similar to BC-Inh1, and the conformers in the macrostates A10 are similar to BC-Inh2 (Figure 5d). The structural similarity between the simulated conformers and BC-Inh1/2 strongly argue that the two covalent inhibitors successfully trapped key hidden conformers of apo-SETD8.
Similar to that of apo-SETD8, the interconversion network of the macrostates of SAM-bound SETD8 also displays a dumbbell-like shape with S9 as the hub-like state connecting the two lobes of the network (Figures 2d, 5a). The macrostates S1 and S3-S5 are multi-connected states; S6, S8, and S10 are satellite-like states; S2 and S7 are basin-like states (Figure 5a,b).
Notably, the complexity of the overall conformational landscape of SAM-bound SETD8 is dramatically reduced in comparison with those of apo-SETD8 (Figures 2d, 5a). The conformers in S1, S2, and S10 are structurally similar to those of A20, as well as BC-SAM; the conformers in S4, S6, and S8 are structurally similar to those in A22 and A24, as well as TC (Figure 5c,d).
The structural similarities between these apo and SAM-bound macrostates suggest possible pathways for connecting the two conformational landscapes upon SAM binding. representations of TC with cancer-associated SETD8 mutations highlighted. c, Differential residuecontact maps of cancer-associated SETD8 mutants in reference to wild type apo-SETD8 (grey). d, Representative contacts in the differential residue-contact maps of cancer-associated SETD8 mutants. The contacts of SETD8 mutants with >3-fold gain of contact fraction relative to wild-type SETD8 are listed and color-coded according to the increased magnitude of the contact fraction. e, Cartoon representations of neo-conformations revealed by simulations of SETD8 mutants. f, Differential residue-contact maps of the structurally relaxed a-helix at the SET-I motif of SETD8 A296T mutant. Decrease of contact fraction of SETD8 mutants relative to wild-type SETD8 is colored in blue. g, Enzymatic activities of wild-type and mutated SETD8 determined by an in vitro radiometric assay with H4K20 peptide substrate. Here SETD8 mutants are categorized as the following: red, uncovered neo-conformations (Neo-conf.) with > 90% loss of methyltransferase activity; green, populated inactive conformations (Pop. shift) with partially abolished methyltransferase activity; blue, no dramatic change of differential contact maps with comparable methyltransferase activity with wild-type SETD8; brown, unknown relationship between differential contact maps and methyltransferase activities. Data are mean ± standard deviation (s.d.) of 3 replicates.
Characterization of cancer-associated SETD8 mutants. Sequences from tumor samples retrieved from cBioPortal [42][43][44] contain two dozen point mutations in the catalytic domain of SETD8 (Figure 6a,b, Table S11). We expect that some of these mutations perturb SETD8 function. Because of conformational heterogeneity, it has historically been challenging for in silico approaches to annotate how mutations---in particular those structurally remote from functional sites---affect a target protein on the basis of its static structure(s) [45][46][47] . Here, we envisioned addressing this challenge with the aid of the dynamic conformational landscapes of SETD8. To characterize mutations remote from catalytic sites (around 80% of known mutations), 40 independent microsecond-long MD simulations for each of the cancer-associated apo-SETD8 mutants were conducted with seed structures prepared from one ternary complex (TC) conformer. We then constructed a differential residue-contact map for each variant (Figure 6c,d) and extracted snapshots representing most dramatic conformational deviations from the wild type conformational ensembles (Figure 6e). Remarkably, even with modest simulation time, several cancer-associated mutants displayed neo-conformations that were not observed in the 5 ms wild-type dataset and cannot be predicted from static X-ray crystal structures. Strikingly, all of the neo-conformations display distinct reorganizations at the SET-I motif (Figure 6e). For instance, a single point mutation A296T, ~16 Å remote from the active site, yields five distinct neo-conformations (Figure 6e). In addition, relative to wild-type apo-SETD8, this mutant populates several conformations with a structurally relaxed a-helix at the SET-I motif ( Figure   6e). C324del, ~20 Å from the SET-I motif, is associated with three neo-conformations and displays the most dramatic changes in the differential contact map (Figure 6d, panel 13). The remote H340D mutation is associated with one neo-conformation as well as more populated conformations containing spatially compressed active sites (Figure 6d, panel 7; 6e). Using in vitro radiometric assays, the A296T and H340D mutants were characterized by loss of the methyltransferase activity on H4K20 peptide substrate (Figure 6g). The failure to purify recombinant C324del also supports the impact of this deletion on SETD8 function. H388Q, which mutates a histidine involved in substrate binding, is also associated with neoconformations as well as loss of the methyltransferase activity (Figure 6e,g). These observations provide potential molecular rationale for how remote mutations can alter the active sites and the SET-I motif---and hence catalysis---via modulating the conformational landscape. Exceptions are T274I, R279W, R279Q, and A368V, which yielded neo-conformations but showed activity comparable to wild-type SETD8 (Figure 6e,g).
The differential residue-contact maps further revealed that remote mutations can alter conformational landscapes by altering populations of pre-existing conformations (Figure 6c,d).
For instance, E257K, G280S, A301V, T309M, E330Q, D352Y mutations populate conformations containing spatially compressed active sites ( Figure S37); E372D populates conformations containing a constrained post-SET motif; R333C populates conformations with reorganized SET motifs adjacent to the peptide binding pocket. All of these mutations showed partial loss of methyltransferase activity (Figure 6g). Notably, these structural alterations are often remote from the corresponding mutation sites (Figure 6b). In contrast, R244S, T274I, and V356I showed no significant conformational change on the basis of their differential contact maps, consistent with their comparable methyltransferase activity to wild type SETD8 (Figure   6g). Likely due to insufficient simulation time (40×1 µs/mutant), R333L and L334P variants, characterized by partial-to-complete loss of the methyltransferase activity (Figure 6g), showed similar conformational landscapes to that of wild-type apo-SETD8. Exploring conformational landscapes is thus an effective strategy to reveal structural alterations associated with majority of remote-site mutations for functional annotation.

Discussion
Here we have demonstrated that tight integration of structural determination---using covalent probes and multiple ligand-binding states to trap hidden conformations (Figure 1)---with massively distributed molecular simulations and the powerful framework of Markov state models (Figure 2b) can provide unprecedented insights into the detailed conformational dynamics of an enzyme. The current work demonstrates the merit of an approach that leverages multiple X-ray structures with distinct diverse conformations for MD simulations and machinelearning-based MSM construction to elucidate complex conformational dynamics, and validates the resulting model experimentally with testable biophysical predictions (Figure 3). Previously, individual components of our integrative strategy have been employed to study the dynamics of transcriptional activators 48 , kinases 11,12 , and allosteric regulation 13 . However, it is the first time that these diverse approaches are consolidated explicitly with the goal of illuminating conformational dynamics of an enzyme in a comprehensive and feasible manner. Assessment of key computational parameters concluded that we have utilized sufficient diverse seed structures and simulation time for microstate discovery and thus robust construction of conformational landscapes (Figure 4). Notably, we relied on a unique computational resource---Folding@home---to collect remarkable six-millisecond simulation data (see Method). Without access to Folding@home, contemporaneous progress on developing adaptive Markov state model construction algorithms---where iterative model building guides the collection of additional simulation data 49,50 ---will still allow research groups to achieve this feat on local GPU clusters or cloud resources in the near future. Furthermore, the concept of adaptive model construction can be extended to identify which new structural or biophysical data would be valuable in reducing uncertainty [51][52][53] and producing refined MSMs. The integrated platform and concept formulated via this work can be readily transformed to explore dynamic conformational landscapes of other proteins.
This work represents the first time that conformational dynamics of a protein methyltransferase has been definitively characterized with atomic details. Strikingly, SETD8 adopts extremely diverse dynamic conformations in apo and SAM-bound states (24 and 10 kinetically metastable macrostates, respectively, Figure 2). Interconversions between metastable conformers cover a broad spatio-temporal scale in particular associated with motions of SETD8's SET-I and post-SET motifs (Figures 1,5). In the apo landscape, the general structural features of the X-ray structures of BC-Inh1, BC-Inh2, APO, BC-SAM and TC (Figure 1) are recapitulated by a subset of macrostates (e.g. A6 for BC-Inh1; A10 for BC-Inh2; A12 for APO; A20 for BC-SAM; A22, A24 for TC, 6 of 24 macrostates, Figure 5). Such observation indicates that these X-tray structures trapped in the different ligand-binding states are not ligand-induced artifacts but indeed relevant snapshots of hidden conformations of apo-SETD8. Similarly, a few macrostates in the SAM-bound landscape also recapitulate major structural features of the two cofactor-bound X-ray structures (e.g. S1, S2, S10 for BC-SAM, S4, S6, S8 for TC, 6 of 10 macrostates, Figure 5). Meanwhile, our results also demonstrate that X-ray crystallography alone is insufficient to capture all metastable conformations of SETD8. Remarkably, there is no correlation of overall structural similarity and interconversion rates between metastable conformers. Though the anticipated findings of fast transitions between structurally similar conformers and slow transitions between structurally distinct conformers (e.g. microstates within individual satellite macrostates A17-A24 of apo SETD8; S6, S8, and S10 of SAM-bound SETD8, Figure 5), we frequently observed fast kinetics of transitions between structurally distinct microstates (e.g. microstates within hub-like macrostates A11 and S8; multi-connected states A1-A4, A9, A14, S1 and S3-S5) and vice versa (e.g. macrostates A22 and A24) (Figures   2,5). It is thus interesting to examine how other factors such as specific residue contacts and cooperative long range motions of certain structural motifs play roles on interconversion kinetics.
Functional annotation of the landscapes revealed that the SET-I motif adopts diverse conformations (Figure 2), and its overall configuration is a key feature that differentiates the lobes of the dumbbell-like conformational landscape of SETD8. The conformational dynamics within the hub-like macrostate A11 mainly involves motions of the SET-I motif. Two gain-offunction I293G and E292G variants of SETD8 were designed for relaxing distorted configurations of the SET-I motif upon SAM binding (Figure 3). These findings argue the functional essentiality of the intrinsically dynamic motions of SET-I motif for SETD8 SAM binding and catalysis. Importance of dynamic conformational modulation of the SET-I motif has also been shown for other SET-domain PKMTs. For instance, the SET domains of MLLs and EZH1/2 alone are catalytically inert but active in the presence of binding partners WDR5-RbBP5-Ash2L-Dpy30 (referred as MLL-WRAD) and EED-Suz12 (referred as PRC2), respectively 15 . Recent structural evidence implicated that the formation of these complexes regulates the conformational dynamics of the SET-I motif, which is essential for catalysis 54,55 .
Interestingly, this region has also been exploited by cancer-associated mutants of PKMTs. For instance, NSD2's E1099 is located in its SET-I motif and its E1099K mutant was characterized as a hot-spot cancer mutation with the gain-of-activity of H3K36 methylation 56 . Additionally, many mutations of PKMTs have been mapped in their SET-I motifs, implicating their potential roles in alternation of function ( Figure S39, Table S12).
In contrast to static X-ray structures, dynamic conformational landscapes greatly facilitated the characterization of cancer-associated SETD8 mutants (Figure 6). A significant portion of cancer-associated, loss-of-function SETD8 mutations, though remote from active sites, were revealed to act allosterically through perturbing the SET-I motif and thus catalysis (Figure 6).
We also discovered significant changes in the connective networks and a dramatic decrease in conformational heterogeneity upon SAM binding (Figure 2). This finding highlights how enzyme-ligand interactions reshape conformation landscapes. The conformational landscapes of SETD8 thus provide an unprecedented platform for virtual screening of ligand candidates as inhibitors via exploring different modes of interaction (SAM-competitive, substrate-competitive, covalent or allosteric). Uncovering hidden conformations can thus be essential for developing potent and selective SETD8 inhibitors. The conformations of individual SETD8 microstates can be further explored to derive their thermodynamic, kinetic, and even transition-state parameters in a catalytic cycle. Similar strategies can be generally applied to native or disease-associated PKMTs for functional annotation.  wrote the manuscript.

Competing Interests
The authors declare no competing interests.   Table S2. SETD8 mutants were expressed and purified as described above for wild-type SETD8.    74 . COOT and MOLPROBITY were used for interactive rebuilding and geometry validation, respectively 61,64,75 . Data reduction and refinement statistics are summarized in Table   S1. The model was automatically rebuilt with ARP/wARP, manually rebuilt with COOT, and refined with REFMAC. Data reduction and refinement statistics are summarized in Table S1.  Table S1.

BC-SAM (4IJ8
Preparation of SAM-free SETD8. SAM-free SETD8 was prepared as described previously 81  For the calculation of equilibrium constants, the equations of Kd1 = k-1/k1, Keq = k-2/k2, and Kd = Kd1 × Keq/(1+ Keq) were followed. For conventional fittings, the fluorescence data were fitted into eq. 1, in which F is the fluorescence intensity, A1 and A2 are the amplitude of the signal changes for fast and slow phases, respectively, kobs fast and kobs slow are the observed rate constants for two phases, and t is time. The plot of kobs fast and kobs slow versus SAM concentrations were fitted with eq. 2 and eq. 3, respectively, where [S] is the concentration of SAM, ki and k-i are the association and dissociation rate constants for step i (i = 1 or 2), respectively. For individual rate constants, data are best fitting values ± s.e. from KinTek.
Uncertainties of Kd1, Keq, Kd are shown as s.e. calculated by the propagation of s.e. from individual rate constants and dissociation constants, respectively. KinTek.
Methyltransferase assay of cancer-associated SETD8 mutants. The methyltransferase activities of wild-type and cancer-associated SETD8 mutations were characterized by a previously described filter  Table S3 for details). The structural chimeras were constructed with MDTraj 1.7.2 85 by combining the C-flanking domain (residues 377-393) with the rest of the protein from different crystal structures (see SI for details). Using OpenMM 6.3.1 86 , protonation states appropriate for pH 7 were assigned with openmm.app.modeller, which uses intrinsic pKa values to determine the most likely ionization states of individual residues but ensures all models are created with the same protonation and tautomeric state so they can be analyzed collectively. The protein was then energy-minimized and relaxed To evaluate a large set of hyperparameters to achieve optimal featurization, a reduced dataset subsampled to 5 ns/frame intervals (986,464 frames, 10% of the dataset) was used for computational feasibility. The following trajectory featurization choices were assessed: a) all residue-residue distances (calculated as the closest distance between the heavy atoms of two residues separated in sequence by at least two neighboring residues) that cross a 0.4 nm contact threshold in either direction at least once (yielding 6,567 of 12,720 total residue-residue distances); b) a transformed version of (a) used by MSMBuilder 98 to emphasize short-range distances in the proximity of residue-residue contact via eq. 4, with steepness = 5 nm -1 and center = 0.5 nm; c) backbone (phi, psi) and sidechain (chi1) dihedral angles, with each angle featurized as its sine and cosine (yielding 920 total features). logistic_transformed_distance = 1/(1 + exp (steepness × (distance -center))) (eq. 4) To identify the optimal featurization, we used a 50:50 shuffle-split cross-validation scheme to evaluate various model hyperparameters while avoiding overfitting. In this scheme, subsets of the groups of trajectories initiated from the same conformation (RUNssee Table S3 for further explanation) are randomly split into training and test sets of 2,509 and 2,510 trajectories respectively, using scikit-learn 0.9.1 99 . All further steps until scoring were conducted by fitting the model to the training set only, then transforming the test set according to this model. Scoring was based on the sum of the top 10 squaredeigenvalues of the transition matrix (rank-10 VAMP-2 96 ). Model scores are reported below as means with standard deviations over five shuffle-splits.
To evaluate each featurization choice, the data were projected into a kinetically relevant space using tICA 100 , retaining all tICs, at lag times of either 5 or 50 ns, with either kinetic 101 or commute mapping 102 .
Each of the tICA outputs was discretized using k-means clustering into 50, 100, 500, or 1000 microstate clusters (see Table S4 for the summary of options assessed). Featurization was performed using MDTraj

Final featurization and microstate number selection.
To determine the optimal number of microstates, we again used variational scoring [94][95][96][97] combined with cross-validation 38 to evaluate model quality. The full dataset (4.931 ms, 0.5 ns/frame, 9,862,657 frames) was separately featurized with the top-scoring feature sets: 6,567 distances (featurization a above) and 920 dihedral angles (featurization c above). As for the featurization selection, we used the 50:50 shuffle-split cross-validation scheme, using the same 5 data splits. All further steps until scoring were conducted by fitting the model to the training set only, subsampled to 5 ns/frame intervals for computational feasibility, then transforming the full training set and the test set according to this model. Data were projected into the tICA space using a lag time of 5 ns.
The tICs were scaled by commute mapping, with subsequent clustering operations using a sufficient number of tICs necessary to explain 95% of the total kinetic content. The tICA outputs were discretized using k-means clustering into 100, 500, 1000, 2000, 3000, 4000, or 5000 microstate clusters (see Table   S5 for the summary of options assessed MSMs were constructed with PyEMMA 2.5.1 103 from the discrete trajectories of the training set using a lag time of 50 ns and subsequently scored on the test set, using the rank-10 VAMP-2 96 score. The highest scoring model ( Figure S21) had dihedral features (featurization c above) and 100 microstates (VAMP-2 = 9.25 (SD = 0.32)). tICA and k-means clustering were refitted to the full dataset subsampled to 5 ns/frame intervals for computational feasibility. Keeping the number of tICs necessary to explain 95% of the total kinetic content resulted in 466 tICs used for k-means clustering. The full dataset was then transformed to give the final discretized trajectories at 0.5 ns/frame intervals. Checking the convergence of the implied timescales validated the choice of the MSM lag time (Figure S22). The Chapman-Kolmogorov test 104 was then conducted on the MSM to validate the self-consistency of the model ( Figure   S23). To aid structural interpretation, the 10 frames closest to each of the 100 microstate cluster centers were extracted from the dataset. To quantify the structural diversity of each macrostate, a sample of 100 frames was drawn from each macrostate with probabilities for each frame given by the observation probability of the frame's microstate from the given macrostate divided by the total number of frames in the frame's microstate. The C-alpha RMSD (after superpose of the SET motifs only) of each frame versus all other 99 frames was then calculated, and the minimum average RMSD over all 100 reference frames was reported ( were assigned using Antechamber 112 , with missing parameters assigned using Antechamber's ParmChk2.

Coarse
The SAM parameter files were then converted from the Amber format to the OpenMM XML format using the conversion script distributed with the openmm-forcefields package 113 . The systems were solvated in rectilinear water boxes with 1 nm padding and neutralized with a minimal amount of NaCl.
This resulted in the final systems containing 34556 atoms (system prepared from 4IJ8) and 35588 atoms (system prepared from 2BQZ). These were energy-minimized and relaxed for 1 fs in the NVT (T = 10 K) ensemble using the OpenMM Langevin integrator with a 0.01 fs timestep, 91 ps -1 collision rate.
Nonbonded interactions were treated with the reaction field method only during minimization (due to its increased stability over PME when steric clashes needed to be resolved following introduction of mutations) 114 at a cutoff of 0.9 nm. The systems were then equilibrated for 10 ns in the NpT (p = 1 atm, T = 300 K) ensemble using the OpenMM Langevin integrator, the PME nonbonded method, a Monte Carlo molecular-scaling barostat with anh updated interval of 25 steps, and packaged with OpenMM 6.3.1 86 as seeds for production simulation on Folding@home 88 . All other force field parameters and simulation settings were as previously described for apo-SETD8.
Folding@home simulations. In total, 1,000 independent MD simulations were generated on Folding@home: 500 each for the two seed structures prepared above. Simulations employed the same settings as for NpT production of apo-SETD8. 99.7% of the generated trajectories (997 trajectories) successfully reached 1 µs each (see Figure S25 for length distribution), resulting in 1.003 ms of aggregate simulation time and 2,005,945 frames. This amount of simulation time corresponds to ~ 46 GPU years on an NVIDIA GeForce GTX 980 processor. Prior to data analysis, the first 25 ns of each trajectory were discarded to allow the trajectories to relax away from the initial equilibrated configurations. One trajectory was shorter than the length being discarded and was removed from the dataset.  Figure S27). Finally, a Hidden Markov model (HMM) was constructed for a 50-ns lag time using 10 macrostates (the minimal number of macrostates to achieve increasing structural separation between distinct SET-I and post-SET configurations was chosen in the same way as for apo-SETD8). As for apo-SETD8, log-inverse fluxes between macrostates were used to construct a two-dimensional representation, and the 10 frames closest to the microstate cluster centers were assigned to the macrostate to which they had the highest fractional membership to aid structural interpretation. Prior to visualization, frames were re-imaged with MDTraj 1.8 85 to ensure SETD8 was centered and the SAM ligand was in the same unit cell. Further, as for apo-SETD8, macrostate C-alpha RMSDs to the homology models derived from all 5 crystal structures were calculated by weighted mean over microstate average RMSDs, and structural diversity was quantified by the reference frame with the minimum average RMSD of each macrostate.

Cancer-associated mutant apo-SETD8 simulations
Preparation of molecular dynamics (MD) simulations. All-atom models of the same 162-residue SETD8 fragment with each of 24 cancer-associated single mutations (including 1 deletion giving a 161residue fragment) identified from the cBioPortal for Cancer Genomics [42][43][44] were prepared in an analogical way to apo-SETD8 using Ensembler 1.0.5. The mutants prepared are summarized in Table S11 (#1-21, 23-25). As we aimed to gain a direct interpretation of the influence of these mutations on the enzymatic activity of SETD8, only a single chain of the TC structure was used as the template. To choose the particular chain out of the 18 TC chains used in the apo-SETD8 simulations (Table S3), the homology models of all the TC chains generated by Ensembler were projected into the apo-SETD8 tICA space using PyEMMA 2.4 103 . The distances between the points in the tICA space were then calculated with SciPy 1.0 and chain "A" of the 1ZKK structure, which had the smallest average distance to all others, was selected.
The modeling procedure was the same as for apo-SETD8, except the appropriately mutated sequences were passed to Ensembler 83 . Briefly, homology models were created by UCSF MODELLER 9.16 84 , protonation states appropriate for pH 7 were assigned with OpenMM 6.3.1 86 , the models were then energy-minimized and relaxed with 100 ps of implicit solvent dynamics. The proteins were then solvated in rectilinear boxes with 1 nm padding and neutralized with a minimal amount of NaCl. This resulted in the final systems containing between 35,185 and 35,208 atoms. These were equilibrated for 5 ns and packaged as seeds for production simulation on Folding@home. All force field parameters and simulation settings were as previously described for wild type apo-SETD8.
Folding@home simulations. In total, 960 simulations were generated on Folding@home: 40 for each of the mutants. Simulations employed the same settings as for NpT production of wild type apo-SETD8. 99.7% of the generated trajectories (957 trajectories) successfully reached 1 µs each (see Figure S29 for length distribution), resulting in the final aggregate dataset of 0.966 ms and 1,931,849 frames. This amount of simulation time corresponds to ~ 44 GPU years on an NVIDIA GeForce GTX 980 processor.
This trajectory dataset without solvent is available via the Open Science Framework at https://osf.io/2h6p4. The code used for the generation and analysis of the molecular dynamics data is available via a Github repository at https://github.com/choderalab/SETD8-materials.

Contact map analysis.
Prior to data analysis, the first 750 ns of each trajectory were discarded to allow for successful metastable transitions out of the wild type kinetic basin. For the remaining frames of each mutant, all residue-residue distances (calculated as the closest distance between the heavy atoms of two residues) for which the two residues are separated in sequence by at least two other residues (yielding 12,720 residue-residue distances) were calculated with MDTraj 1.8 85 . These were converted into binary contact maps by replacing all distances below the 0.4 nm contact threshold with 1's and all other distances with 0's, and casting into a square-form matrix. These were then averaged over all frames of each mutant, yielding one contact map for each mutant. In the same way, the wild type contact map was calculated using the 60 wild type apo-SETD8 trajectories started from chain "A" of the 1ZKK structure (Table S3). The wild type contact map was then subtracted from the mutant contact maps to generate one absolute contact change map for each mutant. In the one case of amino acid deletion, all contact changes corresponding to that residue position were set to zero. Relative contact changes were also calculated by dividing the absolute contact change value by the wild type contact value and taking the modulus of the result. The result was set to zero where the wild type value was zero. Contacts with absolute changes of more than 0.2 and relative changes of more than 3 were selected for further structural annotations.

Extraction of hypothetical new conformations.
For each mutant, of the contacts that showed an absolute fractional change of more than 0.2 or more negative than -0.2, up to 5 contacts with the most positive changes ("positive contacts") and up to 5 contacts with the most negative changes ("negative contacts") were noted. For each trajectory (after discarding the first 750 ns) of a given mutant, every 10th (every 20th if more than 7,500 such frames were present) of the frames that had all "positive contacts" present or none of the "negative contacts" present was extracted. These frames were manually inspected using PyMOL 1.8.4 115 against the representative conformations of the wild type apo-SETD8 microstates (the 10 frames closest to each of the microstate cluster centers) and frames that displayed similar SET-I and post-SET motif configurations to any of the wild type conformations were discarded. For the remaining frames, the C-alpha RMSDs to all frames of the wild type apo-SETD8 dataset subsampled to 5 ns/frame were calculated using MDTraj 1.8 85 and the wild type frames with the lowest RMSD to each mutant frame were extracted. The mutant frames were manually inspected using PyMOL 1.8.4 115 against the extracted wild type frames and further mutant frames similar to the wild type frames were discarded.
The remaining frames for each mutant were clustered based on the manual inspection and their C-alpha RMSDs to all frames of the given mutant (without discarding the first 750 ns of each trajectory) were calculated using MDTraj 1.8 85