Studying Disease-Associated UBE3A Missense Variants Using Enhanced Sampling Molecular Simulations.

Missense variants in UBE3A underlie neurodevelopmental conditions such as Angelman Syndrome and Autism Spectrum Disorder, but the underlying molecular pathological consequences on protein folding and function are poorly understood. Here, we report a novel, maternally inherited, likely pathogenic missense variant in UBE3A (NM_000462.4(UBE3A_v001):(c.1841T>C) (p.(Leu614Pro))) in a child that presented with myoclonic epilepsy from 14 months, subsequent developmental regression from 16 months, and additional features consistent with Angelman Syndrome. To understand the impact of p.(Leu614Pro) on UBE3A, we used adiabatic biased molecular dynamics and metadynamics simulations to investigate conformational differences from wildtype proteins. Our results suggest that Leu614Pro substitution leads to less efficient binding and substrate processing compared to wildtype. Our results support the use of enhanced sampling molecular simulations to investigate the impact of missense UBE3A variants on protein function that underlies neurodevelopment and human disorders.

Introduction

The UBE3A gene encodes a Homologous to the E6AP Carboxyl Terminus (HECT) domain-containing E6 ubiquitin ligase that conjugates ubiquitin to target proteins for degradation.[1] The importance of UBE3A in humans is reflected in the finding that copy number variation, insertion–deletion variants, as well as single nucleotide variations—such as frameshift, nonsense, and missense variants—in UBE3A are associated with neurodevelopmental disorders including Angelman Syndrome (AS) and Autism Spectrum Disorder (ASD). Notably, loss-of-function UBE3A variants are associated with AS, while mutations that lead to an increase in UBE3A expression are associated with ASD.[1] Approximately 5% of AS is attributable to missense UBE3A variants, yet many of these are documented as “variants of uncertain significance”[2] as their causative impacts are unclear.[3] A recent study reported that protein stability, mislocalization, and impaired E3-ubiquitin ligase activity were significant features for disease-associated missense variants,[4] while other studies have defined critical residues for PKA-dependent protein degradation, as well as the role of UBE3A in neuronal differentiation.[2,5−7] The impact of missense variation on the UBE3A protein at a molecular level remains to be better understood.

Structurally, UBE3A contains an Amino-terminal Zinc-binding domain of Ubiquitin Ligase E3A (AZUL) domain and a HECT domain. The HECT domain comprises an N-terminal lobe (N-lobe) and a C-terminal lobe (C-lobe), the latter of which mediates catalytic activity.[8] X-ray crystallographic structures of the complete HECT domain of UBE3A both in native and E2-bound (UBCH7) form[9] and its C-lobe in a domain-swapped form[10] have collectively revealed an unusual “open” form structure in a trimeric arrangement, with N- and C-lobes in extended conformation, in sharp contrast to the more commonly observed monomeric “closed” structures for most HECT domains, where these lobes tightly interact.[11] Beyond UBE3A, a small range of HECT domains have been determined in complex with interacting partners. Of note is the structure of NEDD4-2 bound to an E2 (UBCH5B) and ubiquitin,[12] likely representing an active or preactive state for this enzyme class; this structure reveals a closed conformation distinct from other closed states,[13] whereby the β-sheet within the C-lobe interacts tightly with E2. Notably, these data altogether suggest conformational flexibility within the N-lobe and between the N- and C-lobes;[11] this flexibility is likely relevant to UBE3A function but it has so far not been well studied, despite the potential for such insight to facilitate our understanding of missense variation in the UBE3A function.

Here, we report a child with AS harboring a novel maternally inherited c.1841T>C (p.Leu614Pro) missense variant located within the HECT domain of UBE3A. Notably, pathogenic and likely pathogenic missense variants documented in AS flank this variant, including c.1810G>C (p.Glu604Gln), c.1823A>G (p.Gln608Arg), and c.1864A>G (p.Asn622Asp).[14] We used computational modeling approaches to determine that a Leu614Pro substitution alters the conformational dynamics of UBE3A to disrupt its ubiquitination function.

Results and Discussion

p.Leu614Pro Missense Variant in UBE3A in an Individual with AS

The female proband is the first child born to a nonconsanguineous Caucasian couple. The pregnancy was conceived by assisted reproductive technology and delivery was by ventouse-assisted vaginal delivery at 40 weeks of gestation. Birth weight was 3.286 kg (25th percentile), length was 49 cm (25th–50th percentile), and head circumference was 34 cm (25th percentile). There were early feeding difficulties and associated poor weight gain that resolved once feeding improved.

Developmental progress was satisfactory for the first 12 months, apart from delayed expressive language. Myoclonic seizures developed at 14 months and intensified in frequency over the next 4 months. Development started to regress from 16 months. Antiepileptic treatment with a combination of clobazam and levetiracetam resulted in a significant reduction in seizure activity with concomitant improvement in developmental progress, although balance remained poor. There is no family history of epilepsy, developmental delay, or other health issues. A number of family members on both sides had delayed childhood expressive language development with no long-term developmental sequelae.

Examination did not identify any dysmorphic features and growth generally had been on or just above the third percentile for age. A mild tremor was noted on neurological assessment. EEGs showed a fairly persistent background abnormality with frequent paroxysms of theta and delta activity that appeared sharply contoured and consistent with the pattern seen in AS. No abnormalities were identified on ophthalmology and audiology assessments, brain MRI, metabolic studies, echocardiogram, and chromosome microarray. Multiplex-ligation-dependent probe amplification (MLPA) studies of chromosome 15q11-13 were unremarkable. Targeted epilepsy gene panel testing identified a variant in UBE3A (NM_000462.4(UBE3A):c.[1841T>C];[1841=] p.[(Leu614Pro)];[(Leu614=)]) resulting in a p.Leu614Pro substitution, reported as uncertain significance. Additional family studies showed this variant to be maternally inherited and located on the grandpaternal copy of chromosome 15 (Figure ) and de novo in the proband’s mother. Thus, the variant was reclassified as likely pathogenic.

Figure 1

Family pedigree for the affected female subject (arrowhead).

Leu614Pro Substitution Influences Conformational Dynamics of UBE3A

We first prepared structures of UBE3A in its putative open and closed/active states. The structure of E2-bound UBE3A (PDB: 1C4Z) was passed through the Protein Preparation Wizard in Schrodinger Suite to provide the putative open state, while the closed (putative ubiquitin-binding) state was generated through homology modeling based on aligning segments of this structure to the NEDD4-2:E2:ubiquitin complex (PDB: 3JW0) (Figure ). Notably, Leu614 is contained within a helix and buried within the C-lobe-contacting subdomain of the N-lobe, suggesting that mutation to proline could impact UBE3A conformational dynamics. Next, we employed enhanced sampling techniques for molecular simulation to explore the conformational dynamics of wildtype UBE3A and its Leu614Pro variant. Using the wildtype structure, adiabatic biased molecular dynamics (ABMD) on the intramolecular contact map was initially performed to determine a low-energy path connecting the putative open and closed states of UBE3A. Five ABMD simulations starting from the putative open state and five starting from the closed state were performed. Frame–frame distance matrices were determined for each of the simulations, with 12 evenly spaced frames selected to represent a low-energy path between open and closed states. From this, only one suitable path was identified (Figures and S1).

Figure 2

Preparation of UBE3A closed-state structure. Multiple copies of the UBE3A:E2 complex (PDB: 1C4Z) were aligned to relevant portions of the NEDD4-2:E2:ubiquitin complex (PDB: 3JW0), in accordance with each colored segment. The closed-form complex was constructed by generating a chimeric model based on the aligned segments. The loop structure between the N- and C-lobes (residues 737–741) was generated in a template-free manner during model building to ensure an appropriate positioning of this loop relative to prepositioned domains. Legend: transparent gray—NEDD4-2 and E2 from PDB: 3JW0; solid gray—ubiquitin from PDB 3JW0; pink—C-lobe-contacting subdomain of N-lobe of UBE3A (residues 497–626 and 703–756) from PDB: 1C4Z; blue-violet—E2-contacting subdomain of N-lobe of UBE3A (residues 627–702) and E2 from PDB: 1C4Z; green—C-lobe of UBE3A from PDB: 1C4Z (residues 742–846); and yellow—location of Leu614 within the C-lobe-contacting subdomain of the UBE3A N-lobe.

Figure 3

Frames representing the path between open and closed states of UBE3A generated via ABMD. Interframe root-mean-squared deviations are shown above each arrow. Position along path (s) noted underneath each frame. E2 is colored gray, while UBE3A is colored from N-to-C terminal in a blue-to-red rainbow in each frame.

Well-tempered metadynamics (WTMetaD) in the position along (denoted as s, with s = 1 representing the putative open state and s = 12 representing the putative closed state) and deviation from (denoted as z, with z = 0 being along the path, and z > 0 being a deviation from the path) the ABMD-derived path was carried out for wildtype UBE3A (Figures A–E and S2). WTMetaD of the wildtype UBE3A bound to E2 revealed three major energetic minima (Figure A). The lowest energy of these corresponds to a state intermediate between the open and closed states, although closer to the closed state (Figure B: s ≈ 8, z ≈ 0.125). The second of these is a near-closed state (Figure C: s ≈ 10.5, z ≈ 0.1). The third is a broad minimum representing at least two distinct types of “twisted” closed states (Figure D: s ≈ 10.5, z ≈ 0.275; Figure E: s ≈ 10.5, z ≈ 0.35), wherein the C-lobe rotates about the loop linking it to the N-lobe. States more than halfway to the putative open state (i.e., s < 6) appear strongly disfavored, as evidenced by no substantial energetic minima in this region and large energetic barriers to access such states from the closed states. Next, we performed WTMetaD on the Leu614Pro variant (Figure F). This variant loosely recapitulates the lowest energy state observed in wildtype (Figure G: s ≈ 7.5, z ≈ 0.15), yet all additional observed states are open (Figure H: s ≈ 1.5, z ≈ 0.15) or “twisted” open states (Figure I: s ≈ 1.5, z ≈ 0.25; Figure J: s ≈ 1.5, z ≈ 0.325), with fully closed and near-closed states strongly disfavored. Energetic barriers for moving between states of Leu614Pro appear substantially lower than in wildtype, as evidenced by broader contouring of the free energy surface.

Figure 4

Well-tempered metadynamics simulations of wildtype and Leu614Pro UBE3A-E2 complexes. (A) Free energy surface for the wildtype UBE3A-E2 complex. (B–E) Structural ensembles associated with low-energy states of the wildtype UBE3A-E2 complex. (F) Free energy surface for the UBE3A-E2 complex. (G–J) Structural ensembles associated with low-energy states of the UBE3A-E2 complex. The approximate locations of the structural ensembles in the free energy surfaces are marked with the corresponding panel identifier and the structures have been aligned with respect to E2 and shown with the same camera position in all depictions. Legend for structural ensembles: gray—E2; blue-to-red rainbow—UBE3A N-to-C terminal.

To further explore the utility of this approach, we studied two additional missense variants documented in ClinVar,[14] namely, Ile755Thr, a likely benign variant that maps to the subdomain of UBE3A within which Leu614 is located, which is identified in an individual with intellectual disability, and Met589Lys, a pathogenic variant that maps to a region within UBE3A that facilitates helix–helix interactions in an analogous role to Leu614 (Figure S3). WTMetaD simulation of Ile755Thr largely recapitulates the closed (Figure S3A–B: s ≈ 10.5, z ≈ 0.225) and “twisted” closed (Figure S3C: s ≈ 10.5, z ≈ 0.1; Figure D: s ≈ 7.5, z ≈ 0.15) states that were observed for wildtype. A low-energy open state (Figure S3A: s ≈ 2.5, z ≈ 0.15) is also observed. Large energetic barriers between the low-energy states of Ile755Thr exist, suggesting limited likelihood of interconversion between these states, with the closed states representing the lowest energy states. In contrast, WTMetaD simulation of the pathogenic Met589Lys variant (Figure S3E) revealed favored states that were intermediate between the open and closed states (Figure S3F: s ≈ 6.5, z ≈ 0.2; Figure S3G: s ≈ 6.5, z ≈ 0.325; Figure S3H: s ≈ 7.5, z ≈ 0.15), failing to recapitulate the closed or “twisted” closed states adopted by wildtype or Ile755Thr.

To summarize, wildtype UBE3A and the Ile755Thr variant feature low-energy closed (putative ubiquitin-binding) states, suggesting normal ubiquitin transfer function, while both variants identified in individuals with AS (Leu614Pro and Met589Lys) lost access to the closed (putative ubiquitin-binding) state and thus may fail to efficiently bind and transfer ubiquitin. Altogether, these findings suggest that our approach could be informative to explore the potential differences between UBE3A variants that are associated with clinical features such as AS and intellectual disability.

Enhanced Sampling Approaches to Study the Functional Impact of UBE3A Missense Variation

UBE3A is challenging to study at the biochemical level because while dosage is relevant to its role within cells, steady-state levels of this protein are influenced by its autocatalytic activity.[2,4] While computational approaches, such as PolyPhen2,[15] SIFT,[16] and PROVEAN,[17] can estimate whether missense variants are likely damaging to function, their impact on protein–protein interaction and enzymatic catalysis is difficult to establish through these approaches. Advances in both computer hardware and software[18,19] have underpinned our application of enhanced sampling approaches to study missense variation on UBE3A function; yet, important considerations apply. First, as Protein Data Bank entries increase[20] along with improved quality of single chain protein structure predictions,[21] accurate computational structures of proteins alone or in complex with partners are required, which, ideally, are cross-validated through laboratory studies.[22] Second, a hypothesis-driven approach is required to study the effects of variants, as this dictates the approach employed for investigation. For example, alchemical perturbation-based approaches[23] are invaluable to study missense variants that influence protein–protein interactions and protein stability.[24] However, where large conformational differences are likely to be relevant to protein function, for instance, in association with a receptor activation event,[25] such approaches may be insufficient to characterize missense variants.

Enhanced sampling approaches employed here provide the opportunity to thoroughly explore protein conformation and hence identify likely low-energy states.[26] The use of these approaches is not without its own challenges, which primarily involve defining relevant state variables to facilitate transition between given molecular endpoints and a suitable sampling therein. Further, the particular state variables used and their definition will generally vary depending on the protein being explored and the nature of the transition, implying no one particular approach to applying these types of methods. Nonetheless, when these approaches can be applied, they provide insight into molecular function, as demonstrated here for UBE3A missense variants, including a novel UBE3A p.(Leu614Pro) AS-associated variant described in this study. More broadly, these approaches could be useful as part of an approach to systematically document the functional impact of the spectrum of UBE3A variants in health and disease.

Methods

Medical History and Clinical Data

Informed consent was obtained from the patient/primary caregiver and other family members tested. Chromosome microarray analysis was performed using the Illumina HumanCytoSNP-12 Beadchip with 68 Kb mean resolution and Karyostudio software version 1.4 (genome assembly GRCH37/hg19). MS-MLPA analysis of the chromosome 15q11-13 region was undertaken using the probe mix ME028-C1 (MRC Holland). Massively parallel sequencing and bioinformatically targeted analysis of 292 genes on the Illumina TruSight One Expanded kit identified the UBE3A p.(Leu614Pro) variant. SNP microarray family haplotype analysis was conducted using array information from relevant family members using the above Illumina platform.

Preparation of UBE3A-E2 Complex Structures

The Protein Data Bank structure of E2-bound UBE3A (PDB: 1C4Z)[9] was obtained as the open-form wildtype UBE3A:E2 complex. UBE3A molecules not bound to the E2 were removed, and missing side-chains and loops were added using the Protein Preparation Wizard within Schrodinger Suite 2018-3. For the closed complex, a chimeric homology modeling strategy was pursued, aligning specific portions of PDB: 1C4Z to the NEDD4-2:E2:ubiquitin complex (PDB: 3JW0)[12] and retaining ubiquitin during the build (Figure ; geometry summary statistics from MolProbity[27] in Table S1). The initially built closed-form complex was then refined using Prime Minimization on all atoms, following which the ubiquitin molecule was removed. To generate open- and closed-form variant UBE3A:E2 complexes, variants were introduced in each of the generated wildtype complexes using UCSF Chimera.

Open- and closed-form UBE3A:E2 complexes were parameterized for simulations in AmberTools 14 using the AMBER ff14SB force field.[28] Topologies were converted to GROMACS format using acpype(29) and subsequent steps performed using GROMACS 2018.3.[30] Structures were solvated with TIP3P water[31] in a triclinic box with a minimum distance of 25 Å from complex to box edge. The box was charge-neutralized by the addition of sodium ions, following which steepest descent minimization was performed. Energy-minimized systems were subject to brief (100 ps) equilibrations in the NVT and NPT ensembles, during which harmonic restraints on complex heavy atoms (1000 kJ/(mol nm2)) were employed. Berendsen thermostat with velocity rescaling[32] was used to achieve a target temperature of 310 K and a target pressure of 1 atm.[33] A 2 fs time step was employed in all simulations. Long-range electrostatic interactions were calculated using the Particle Mesh Ewald (PME) method.[34] All covalent bonds were constrained using the Linear Constraints Solver (LINCS).[35] Frames were written at intervals specific to a given simulation, as described below.

Identification of a Low-Energy Path between UBE3A States

Adiabatic biased molecular dynamics (ABMD)[36] as implemented within PLUMED 2.5[37] was used to identify a low-energy path connecting open and closed states of UBE3A in the wildtype UBE3A:E2 complex. We adapted an approach for generating a path used previously for investigating G protein-coupled receptor activation,[38] performing five ABMD simulations starting from open conformation and five simulations starting from closed, each for 10 ns in total (recording atomic coordinates, velocities, and energies every 1 ps), each employing a gradually increasing array of 10 force constants over the course of 10 ns (1 × 10–7, 2 × 10–7, 5 × 10–7, 1 × 10–6, 2 × 10–6, 5 × 10–6, 1 × 10–5, 5 × 10–5, 1 × 10–4, 1 × 10–3 kJ/mol nm2; conducted for 1 ns each). The collective variable biased during the simulations was the distance in the atom–atom contact map from the desired target state; this was biased toward zero (i.e., exact match to target state). The contact map considered pairs of atoms that moved together/apart by greater than 2.0 nm between open and closed states, and that occurred within 0.5 nm of each other in either state. The distance of each atom–atom pair considered in the contact map was transformed using the following switching function:

For each of the 10 ABMD simulations, a pairwise distance matrix was generated, describing the distance in coordinates between each frame and every other frame in the simulation. Optimal paths 12 frames in length were determined using a Monte Carlo-derived algorithm (Figure S4). The algorithm was carried out five times on the pairwise distance matrix for each ABMD simulation. The path ultimately used for metadynamics simulations was selected from the collected identified paths, with distances no greater than 15% between frames.

Metadynamics Simulations of UBE3A

Well-tempered metadynamics (WTMetaD)[39] in PLUMED 2.5 was employed to widely explore the conformational ensembles of wildtype, Leu614Pro, Ile755Thr (likely benign), and Met589Lys (pathogenic) UBE3A bound to E2. The collective variables selected for the WTMetaD bias were the position along the ABMD-generated path (s) and deviation from this path (z); a comprehensive description of path collective variables is provided elsewhere.[40] A short unbiased simulation (20 ns) of the wildtype UBE3A-E2 complex was performed initially to establish relevant values for the σ parameters for WTMetaD in s and z (0.15 and 0.0275, respectively) (Figure S5). Simulations starting from the open state of each complex were performed at 310 K for 500 ns, with Gaussian hills 1 kJ/mol in size added every 1 ps and employing a bias factor of 15. Calculation of the metadynamics reweighting factor[41] was enabled (recorded every 1 ps) and the bias was stored on a grid for computational efficiency (updated every 10 ps using grid bin widths of 0.1 in s and 0.005 in z). Atomic coordinates, velocities, and energies were saved every 10 ps. To limit the exploration of deviations from the defined path, an upper wall in z at 0.5 was employed, with a force constant of 100 kJ/mol and a rescaling factor of 0.001.

Free Energy Calculations and Identification of Low-Energy States of UBE3A

From the WTMetaD simulations, free energy surfaces in terms of s and z were calculated for UBE3A and variant complexes with E2, incorporating relevant reweighting procedures[41] and bin widths used during WTMetaD. To avoid bias for the initial state, the first 10 ns of accumulated data was discarded when calculating the free energy surfaces. Free energy was reported relative to the lowest energy value in the determined free energy surface. Convergence of simulations was assessed by monitoring collective variable space over time to verify complete sampling, monitoring Gaussian hill height over time to verify the identification of no additional energetic minima and monitoring the difference between free energy surfaces over time (at 1 ns intervals). Ensembles of structures for each local minimum in the free energy surface within 5 kJ/mol of the global minimum were extracted from the WTMetaD simulations and then clustered using the GROMACS gmx cluster utility, employing the GROMOS algorithm[42] and a 0.2 nm threshold.