Insights into channel modulation mechanism of RYR1 mutants using Ca²⁺ imaging and molecular dynamics

The RYR plays a key role in Ca²⁺ release from the SR during excitation–contraction coupling of skeletal muscle (Ríos and Pizarro, 1991; Iino, 1999; Hwang et al., 2012; Rebbeck et al., 2014; Ríos, 2018). Mutations in the RYR1 gene (RYR1) cause severe muscle disorders (Robinson et al., 2006; Snoeck et al., 2015; Witherspoon and Meilleur, 2016; Jungbluth et al., 2018; Pancaroglu and Van Petegem, 2018), such as malignant hyperthermia (MH; OMIM accession no. 145600) and central core disease (CCD; OMIM accession no. 117000). MH is a serious complication resulting from general anesthesia through commonly used inhalational anesthetics. The typical symptoms include a rapid increase in body temperature and induction of a hypermetabolic state with skeletal muscle rigidity. CCD is a rare, nonprogressive myopathy that presents in infancy and is characterized by hypotonia and proximal muscle weakness. An important feature of CCD is its close association with susceptibility to MH. Exposure to caffeine or halothane has been shown to accelerate the CICR response in the SR of striated muscles (Endo, 2009). It is widely believed that MH mutations cause hyperactivation of CICR, resulting in abnormal Ca²⁺ homeostasis in skeletal muscles. On the other hand, CCD mutations are of two different types: (1) gain-of-function CCD mutations found in the N-terminal (NT) and central regions (Robinson et al., 2006; Murayama et al., 2015) and (2) mutations found in the C-terminal region that result in loss of function of the channel (i.e., excitation–contraction uncoupling; Avila and Dirksen, 2001; Dirksen and Avila, 2004).

Genetic analyses of patients with MH and related diseases have revealed >300 mutations in RYR1 (Lanner et al., 2010). There are three main mutation clusters (hotspots), which are found in the NT region (RYR1 amino acid residues 35–614), the central region (2129–2458) in the cytoplasm, and the C-terminal region (4637–4973) near or within channel-forming segments of the RYR1. To understand the mechanism of the pathogenesis of MH, comprehensive study of Ca²⁺ homeostasis associated with structural properties of these mutants is necessary. However, due to lack of such studies, the mechanism remains largely unknown.

RYR1 exists as a huge (∼2 MDa) homotetrameric channel with a large cytoplasmic structure that is involved in the regulation of the channel pore, as well as six transmembrane segments at its C-terminus, forming a cation channel pore. The crystal structure of the NT region (1–559; hereafter referred to as the N-terminal domain [NTD]) of RYR1 was first solved by Van Petegem’s group (Tung et al., 2010). The NTD consists of three domains (A, B, and C) and forms a cytoplasmic vestibule (boxed area in Fig. 1 A), and it contains most of the proposed mutation hotspots underlying NT diseases (Fig. 1 B). The same group further performed structural analysis of a series of NT disease-associated mutant RYR1s and demonstrated that mutations affecting interdomain hydrogen bonds/salt bridges alter related locations of domains (Kimlicka et al., 2013a). However, the travel distance of each domain caused by the mutations was rather small (<2.5 Å), possibly due to crystal packing forces that do not occur under physiological conditions. Thus, evaluating the strength of hydrogen bonds/salt bridges between domains by MD simulation may help to improve understanding of the domain reorientations. Recently, the structures of RYR1 at near-atomic resolution were solved by single-particle cryo-EM by three groups (Efremov et al., 2015; Yan et al., 2015; Zalk et al., 2015; Willegems and Efremov, 2018). The structures revealed a complex architecture, involving a superhelical scaffold in the cytoplasmic region and a channel region of the voltage-gated ion channel superfamily. Subsequently, the RYR1 structure in an open state was solved by many different groups (Bai et al., 2016; des Georges et al., 2016; Wei et al., 2016; Willegems and Efremov, 2018). Although three structures of RYR1 in the open state were revealed, the gating mechanism of RYR is still not fully understood. The NT region located in the cytoplasm is ∼100 Å away from the cation channel pore. Therefore, it remains difficult to understand how a single disease-associated mutation at the NT region influences the opening of the channel pore.

To elucidate the molecular basis of the pathogenesis of MH, we performed functional analysis of RYR1 carrying disease-associated mutations at the NT region using a heterologous expression system in human embryonic kidney (HEK) 293 cells, which we established previously (Murayama et al., 2015). We studied eight MH-associated mutations (C36R, Q156K, R164L, G249R, G342R, R402C, R402H, and Y523C) and two MH/CCD mutations (R164C and Y523S) in the NT hotspot by live-cell Ca²⁺ imaging. Furthermore, multiple sets of MD simulation were performed to investigate the structure of the MH/CCD mutants. Since the quality of the initial models is important to obtain reliable predictions by MD simulation (Law et al., 2005), the atomic structure of the RYR1 NTD (1–559) solved by x-ray crystallography was used as the template molecule for all the simulations. Live-cell Ca²⁺ imaging in RYR1 mutants revealed abnormalities in Ca²⁺ homeostasis in terms of caffeine sensitivity, ER Ca²⁺ content, and resting cytoplasmic Ca²⁺ concentration ([Ca²⁺]_cyt). Statistical analysis of MD simulation revealed that changes in Ca²⁺ homeostasis were correlated with changes in the specific hydrogen bonds/salt bridges. These results will clarify the functional significance of the NT region of RYR1. They might also provide clues to why an increase in the channel pore open probability occurs in mutants. Ultimately, this study may lead to the development of new therapeutic drugs for MH/CCD patients.

Disease-associated mutation in the NT region corresponding to Q156K was introduced into the rabbit RYR1 complementary DNA (cDNA) by inverse PCR. Inverse PCR was performed using a HindIII-SalI fragment (pBS-RYR1cs1) from cDNA cassettes encoding full-length rabbit skeletal muscle RYR1 (pBS-RYR1; Tong et al., 1997, 1999; Nakano et al., 2014). Other disease-associated mutations in the NT region used in this study (C36R, R164C, R164L, G249R, G342R, R402C, R402H, Y523C, and Y523S) have been previously reported (Murayama et al., 2015). In addition, as mutants for subsequent verification, E40A and D61A were also prepared by the same method as above. Stable and inducible RYR1 mutations in HEK293 cells were generated as previously described (Kakizawa et al., 2012; Nakano et al., 2014; Murayama et al., 2015) using the Flp-In T-REx system (Invitrogen). Clones with suitable RYR1 expression were selected and used in subsequent experiments.

Ca²⁺ measurements were performed using previously described procedures (Kakizawa et al., 2012; Murayama et al., 2015, 2016). Briefly, HEK293 cells expressing WT or mutant RYR1 were cultured on glass-bottomed dishes (MatTek). 24–30 h after the addition of doxycycline at the final concentration of 2 µg/ml, cells were loaded with 4 µM fura-2 acetoxymethyl ester in physiological salt solution containing 150 mM NaCl, 4 mM KCl, 2 mM CaCl₂, 1 mM MgCl₂, 5.6 mM glucose, and 5 mM HEPES at pH 7.4. All experiments were performed at room temperature (23–25°C).

Fluorescence images were acquired at 420 nm using an inverted microscope (IX70; Olympus, Japan) equipped with a 40× objective (numerical aperture, 0.95; UPlanSApo/340; Olympus), together with a cooled charge-coupled device camera (Rolera XR; Teledyne QImaging), at a rate of one frame every 1 or 2 s. Excitation wavelengths were 345 nm and 380 nm. Image analyses were performed using IPLab software (BD Biosciences Bioimaging). Regions of interest corresponding to individual cells were selected, and the mean fluorescence intensity of each region of interest minus the background intensity was calculated for each frame. We used the F345/F380 ratio (the mean fluorescence intensity at an excitation wavelength of 345 nm divided by that at an excitation wavelength of 380 nm) to estimate [Ca²⁺]_cyt, as described elsewhere (Kakizawa et al., 2012; Murayama et al., 2015). The K_d (239 nM) for Ca²⁺ was determined via in vitro calibration of fura-2 fluorescence (Grynkiewicz et al., 1985).

The atomic coordinates of the NTD (1–559) of rabbit RYR1 (PDB accession no. 2XOA; Tung et al., 2010) were used as a template for all calculations. Modeling of missing loops in 2XOA (85–96, 126–129, 185–188, 226–228, 324–327, 362–370, 427–435, 456–464, and 500–505) was performed by using the MODELLER program (version 9.12; Fiser et al., 2000; Martí-Renom et al., 2000; Webb and Sali, 2014). Construction of the initial mutant structures (C36R, Q156K, R164C, R164L, G249R, G342R, R402C, R402H, Y523C, Y523S, E40A, and D61A) and additions of hydrogen bonds/salt bridges and other missing atoms were performed using the LEaP module in the AMBER package (version 14). The AMBER package (version 14) with the parm99 force field (Case et al., 2005) was used for MD simulation.

The initial structure of the NTD was immersed in a rectangular box filled with ∼11,500 TIP3P water molecules (Jorgensen et al., 1983). To ensure a physiological ionic concentration of 0.15 M and zero net charge, K⁺ and Cl⁻ ions were added adequately. The time step was 1 fs, and periodic boundary conditions were applied throughout the simulation. The atomic coordinates of the systems were saved for later analysis every 0.5 ps during MD simulations. A 10-Å cutoff distance was used for nonbonded interactions. The SHAKE algorithm (Ryckaert and Ciccotti, 1977)tella was applied to constrain the bond lengths of hydrogen-containing bonds. After energy minimization, the system was gradually heated in 50–60 steps from 100 K to 310 K under constant-volume conditions. Trajectories 50 ns long were obtained for each mutant and WT under a constant pressure (1 atm) and temperature (310 K) after the equilibration of the system. Twenty such trajectories were obtained for each of the initial structures of the WT and the mutants, applying different sets of initial velocities at the start of the simulation run to sample different sets of the dynamic structure of the molecule.

To evaluate the stability of the hydrogen bonds/salt bridges between interdomains during the MD simulation, we calculated and analyzed a value of the probability of forming the hydrogen bonds/salt bridges using the “ptraj” program in the AMBER package. The initial 20 ns of data were discarded as preequilibration. The VMD program (Humphrey et al., 1996) was used for analysis and visualizations.

As the initial model of the tetramer of NTDs, the cryo-EM structure of RYR1 (PDB accession no. 5GKZ; residues 1–532; Bai et al., 2016) was used. Modeling of missing loops in 5GKZ (81–99, 116–145, 182–192, 223–230, 235–242, 264–270, 320–335, 360–379, 414–437, 450–473, and 490–514) was performed by using the MODELLER program (version 9.12; Martí-Renom et al., 2000; Webb and Sali, 2014). Construction of the initial mutant structures (R402C and R402H) and additions of hydrogen bonds/salt bridges and other missing atoms were performed using the LEaP module in the AMBER package (version 14). The AMBER package (version 14) with the 14SB force field (Case et al., 2005) was used for MD simulations.

The initial structure of the NTD was immersed in a rectangular box filled with ∼95,000 TIP3P (transferable intermolecular potential with 3 points) water molecules (Jorgensen et al., 1983). To ensure a physiological ionic concentration of 0.15 M and zero net charge, K⁺ and Cl⁻ ions were adequately added. The time step was 2 fs, and periodic boundary conditions were applied throughout the simulation. The atomic coordinates of the systems were saved for later analysis every 10 ps during MD simulations. A 10-Å cutoff distance was used for nonbonded interactions. The WT and mutant systems were refined with energy minimization using the steepest descent method: A 5,000-step energy minimization was performed with harmonic constraints (force constant, 10 kcal/mol/Ų) applied to all protein atoms.

The systems were then heated to 310 K in 500 ps by MD with harmonic constraints (force constant, 1 kcal/mol/Ų) applied to α-carbon (Cα) atoms only, after which the systems were equilibrated for 1 ns by MD in the number of particles, system volume, absolute temperature ensemble with the same constraints as in heating. Finally, the systems were subjected to a 60-ns production MD run in the number of particles, constant pressure, constant temperature ensemble. The SHAKE algorithm was used to constrain the bond lengths of hydrogen-containing bonds, which allows a time step of 2 fs for MD simulations. For each WT and R402C/H mutant system, three 60-ns MD trajectories were generated and combined for hydrogen bond/salt bridge analyses (after discarding the initial 30 ns for system equilibration).

All structural figures were prepared with PyMOL (PyMOL Molecular Graphics System; Schrödinger). The Coot program was used for analysis of structures (Emsley and Cowtan, 2004). ProFit (Martin and Porter, http://www.bioinf.org.uk/software/profit/; McLachlan, 1982) was used for fitting of structures. As the cryo-EM structures, we used 5TB0 as the closed state (in the presence of EGTA) and 5T9V as the open state (in the presence of caffeine, Ca²⁺, and ATP) of the RYR1 since both structures are the highest resolution in the structures solved by three groups. To create a hypothetical homotetramer of the NTD, monomer structures before or after the MD simulations are fitted to the cryo-EM structure of RYR1 in the closed state (PDB accession no. 5TB0) using the A domain of the NTD. The HBPLUS program (McDonald and Thornton, 1994) was used to identify hydrogen bonds/salt bridges in the atomic coordinates. DynDom (Hayward and Berendsen, 1998) was used to find domain movements.

All statistical analyses of the data were performed using KaleidaGraph software for Windows (Synergy Software). Differences between two groups were analyzed with Student’s t test. P values <0.05 were considered to represent statistical significance. One-way ANOVA with Tukey’s honestly significant difference post hoc test was used when comparing data from three groups.

Fig. S1 shows the results of Western blot analysis of mutant RYR1 channels in HEK293 cells. Fig. S2 shows Ca²⁺ homeostasis in cells expressing other mutant RYR1s. Figs. S3 and S4 show the results of the MD simulations for other mutants. Fig. S5 shows a comparison of the strength (tightness) of the interdomain interactions in the NTD of the WT and mutant RYR1s (C36R, G249R, R402C, and R402H). Fig. S6 shows the results of the MD simulations using the NTD tetramer of the WT. Fig. S7 shows the results of the MD simulations using the NTD tetramer of the R402 mutants. Fig. S8 shows the x-ray structure of the RYRs and inositol-1,4,5-triphosphate receptor (IP₃R) in the NTD. Table S1 shows the disease-associated mutations used in this study. Table S2 shows the disease-associated mutations located within the domains using a database (https://www.cardiodb.org/paralogue_annotation/).

Ten mutant RYR1 channels with mutations at seven different positions in the NT region were prepared (Fig. 1 B; and Tables S1 and S2). All the mutations have been associated with MH susceptibility, and two of them (R164C and Y523S) have also been associated with CCD phenotypes (i.e., MH/CCD; see Table S1). These mutant RYR1 channels were heterologously expressed in HEK293 cells, which do not express endogenous RYRs (Murayama et al., 2015, 2016). RYR1 expression was induced by the addition of doxycycline (2 µg/ml) and confirmed by Western blotting (Fig. S1 A).

We first characterized the caffeine-induced Ca²⁺ transient as an indicator of CICR activity. As the control, HEK293 cells without transfection of the RYR1 gene were used. No substantial Ca²⁺ increase was observed upon addition of 3 mM of caffeine (Fig. S2 A). This result is consistent with our previous results (Kakizawa et al., 2012). As shown in Fig. 2 A and Fig. S2 B, these mutants showed varied responses to caffeine as compared with the WT. Since small but significant Ca²⁺ transients were detected at 0.3 mM caffeine in WT and the responses plateaued at 10 mM caffeine in all cases (Fig. S2 C), peaks of the Ca²⁺ transient at 0.3 and 10 mM caffeine were compared in Fig. 2, B and C. Five of the mutants studied (Q156K, R164C, R164L, R402C, and R402H) exhibited significant increases in response to 0.3 mM caffeine, indicating marked enhancement of caffeine sensitivity, which can be further appreciated by the leftward shift of dose–response relationships (Fig. 2 E and Fig. S2 C). On the other hand, the peak response at 10 mM caffeine, which represents ER Ca²⁺ content, decreased in all 10 mutants (Fig. 2 C). In particular, the mutant Y523S, which is associated with MH/CCD phenotype, exhibited a Ca²⁺ transient only ∼10% of that in WT, despite the expression level being equivalent to that of WT (Fig. S1). While resting [Ca²⁺]_cyt was ~40 nM in cells expressing WT RYR1, most of the disease-associated mutants caused significantly higher levels of resting [Ca²⁺]_cyt, with values ranging between 50 and 90 nM (Fig. 2 D). Interestingly, there was an inverse correlation between resting [Ca²⁺]_cyt and peak Ca²⁺ transient induced by 10 mM caffeine (R² = 0.68), suggesting that an increase in [Ca²⁺]_cyt was attributable mainly to leaks from the ER, in line with the results of our previous study (Murayama et al., 2015; Fig. 2 F).

To address the structural changes caused by disease-associated mutations, MD simulations using NTD monomers of the WT and 10 mutants were performed. Although the resulting movements may be amplified, owing to the lack of adjacent protomers that would be present in the tetramer, they may indicate some structural changes occurring because of the disease-associated mutation. After modeling the missing loops in 2XOA (the atomic coordinates of the NTD of rabbit RYR1), construction of the initial 10 mutant structures, and addition of protons to the initial structures, multiple MD simulations were performed using AMBER (version 14) with the parm99 force field (Case et al., 2005) under physiological conditions. We observed that the molecule remains stable during the simulation, because the RMSD value of Cα atoms throughout the MD simulation was <3 Å (Fig. 3 A and Fig. S3 A).

We next analyzed the structural changes after the MD simulations. First, we compared the MD structures of R402C/H mutants with the structure of the WT after 50-ns simulation (R402C, Fig. 3, A–C; R402H, Fig. 3 A and Fig. S3, B–D). DynDom was used to analyze the domain movements in the R402C and R402H mutants after simulation. Fig. 3 B and Fig. S3 B display representative results from 20 calculations. Surprisingly, the B and C domains together rotated 13.5 ± 0.7 degrees (n = 20) clockwise with respect to the A domain in R402C (Fig. 3 B and Fig. S3 E) or 14.3 ± 1.2 degrees (n = 20) in the R402H mutant (Fig. S3, B and E). The travel distances of the largest moving parts for the R402C and R402H were ∼8 Å and ∼9.5 Å, respectively. Indeed, these distances were much larger than the distance observed in the crystal structure of the mutants (Kimlicka et al., 2013a). In fact, the B and C domains of the R402C and R402H mutants nearly overlapped with those of the WT after 50-ns simulation (Fig. 3 C and Fig. S3 C). Therefore, R402 must have an important role in keeping a tight connection between the A and BC domains. R402, which is located in the C domain, forms hydrogen bonds/salt bridges with D61 and E40 located in the A domain (Fig. 3 D). However, in the case of R402C/H mutants, NTD loses these connections (Fig. 3 F and Fig. S3 D), and this may initiate the rotation of the BC domains with respect to the A domain.

In contrast to the R402C/H mutants, there is no appreciable change in the MD simulations of C36R located at the AB domain interface and G249R buried in the B domain compared with the WT structure after 50-ns simulation (Fig. S4).

The residues for six other mutations (Q156K, R164C, R164L, G342R, Y523C, and Y523S) are located outside the NTD and are not involved in the interactions between domains but are involved in binding to the adjacent NTD unit. Calculations for these mutants were also performed (Fig. S4). However, there is no appreciable change compared with the WT structure after 50-ns simulation (Fig. S4). Since we used a monomer of the NTD in our MD simulation, it seems reasonable that the effect of the mutations located at the surface of NTD was not detected.

One of the most frequently used approaches for demonstrating formations and movements of interfaces in the molecule during MD simulation is evaluation of the interactions between the side chains in the interdomain. The interactions are generally contributed by hydrogen bonds/salt bridges and hydrophobic interactions. Since the mutations we selected are all hydrophilic residues, we focused on analyzing the interaction of hydrogen bonds/salt bridges. Fig. 4 A shows the hydrogen bonds/salt bridges between domains in the crystal structure of the NTD (2XOA) analyzed by HBPLUS. Eight hydrogen bonds/salt bridges (R242-E481, D61-R402, R45-D447, R221-E397, D18-S207, E40-S406, T286-E481, and E40-R402) were found. Since the MD simulation of R402C/H mutants suggested that the interdomain hydrogen bonds/salt bridges may have an important role in exhibiting their structural changes, we decided to analyze all the interdomain hydrogen bonds/salt bridges during MD simulation. To do so, we defined the probability of forming hydrogen bonds/salt bridges. If the specific hydrogen bond/salt bridge is continuously maintained during the simulation, the value is defined as 100%. Since these hydrogen bonds/salt bridges are repeatedly formed and disrupted during the simulation, it is quite natural that the actual value becomes <100%. Therefore, the magnitude of this value corresponds to the stability of the hydrogen bond/salt bridge. For the analysis, 20 independent runs of MD simulation were performed, and the average values were plotted in Fig. 4, B–F. Interestingly, in addition to the eight hydrogen bonds/salt bridges that were found by HBPLUS, we found another important interaction (D61–R283; Fig. 3 E and Fig. 4 B).

As one of the indexes indicating the tightness between domains, the standard number of hydrogen bonds/salt bridges between domains during the simulation was counted. In practice, the probability of the hydrogen bonds/salt bridges between the pair of residues constituting the interaction is considered as the number of hydrogen bonds/salt bridges between each domain. Then, the total number of hydrogen bonds/salt bridges between each domain was calculated. As shown in Fig. S5 A, the standard numbers of hydrogen bonds/salt bridges between the A and C domains and between the B and C domains are 1.7 and 1.5, respectively. However, the standard number of hydrogen bonds/salt bridges between the A and B domains is 0.7 (Fig. S5 A). Therefore, the interactions between the A and C domains and between the B and C domains are considered to be tighter than between the A and B domains. The behavior and the total number of hydrogen bonds/salt bridges between domains in the C36R and G249R mutants are essentially the same as those in the WT (Fig. S5 B). Thus, it is reasonable to assume that the peak values of the Ca²⁺ transient of these two mutants are mild (Fig. 2). In contrast, the R402C/H mutants showed a loss of two important hydrogen bonds/salt bridges between the A and C domains (D61-R402 and E40-R402). Indeed, the standard number of hydrogen bonds/salt bridges between the A and C domains and between the A and B domains decreased drastically (0.8 and 0.8 in the R402C mutant, respectively, and 0.7 and 0.8 in the R402H mutant, respectively; Fig. S5 B). Meanwhile, the standard numbers of hydrogen bonds/salt bridges between the B and C domains in the mutants were almost the same as those of the WT (1.7 and 1.6 in the R402C and R402H mutant, respectively; Fig. S5 B). Thus, the interaction between the B and C domains is considerably tighter than interactions between the A and C domains and between the A and B domains in mutants.

R402, located in the C domain, seems to play an important role in connecting the A and B domains (Fig. 4 G). R402 forms hydrogen bonds/salt bridges with two residues (E40 and D61) in the A domain (Fig. 3 E). D61, in turn, forms a hydrogen bond/salt bridge with R283 located in the B domain, and E40 also forms a hydrogen bond/salt bridge with S406 located in the C domain. We named this B(R283)–A(D61)–C(R402)–A(E40)–C(S406) network around R402 as the B–A–C–A–C network (Fig. 3, D and E; and Fig. 4 G) because this network is critical for the connection of A, B, and C domains.

We considered why the B and C domains rotate with respect to the A domain in the R402 mutants. Fig. 4, H and I, presents schematic diagrams of the movement of the NTD in the WT and R402C/H mutants. There are four hydrogen bonds/salt bridges between the A and C domains in the WT NTD (Fig. 4 H). Among these, the D61–R402 interaction is stronger than the other three (Fig. 4 B). There are three hydrogen bonds/salt bridges between the B and C domains in the WT NTD (Fig. 4 H). Among these, the R242–E481 interaction is stronger than the other two (Fig. 4 B). Meanwhile, there are only two hydrogen bonds/salt bridges between the A and B domains, and both interactions are relatively weak (Fig. 4 B). Therefore, AC and BC domain interactions are considered to be stronger than the AB domain interaction.

In the R402C/H mutants, however, hydrogen bonds/salt bridges between the A and C domains would become drastically weak because two of the four important interactions (D61–R402 and E40–R402) were lost (Fig. 4, B, E, F, and I). The interaction between the A and B domains was weaker than other interactions in the WT, and this property was shared by the R402C/H mutants (Fig. 4, B, E, F, and I). On the other hand, interactions between the B and C domains were conserved in the R402C/H mutants, as in the WT (Fig. 4, B, E, F, and I), such that the B and C domains tend to move together (Fig. 4 I). This may cause side chain conformational changes in E40 (Fig. 3, E and F), R45, D61 (Fig. 3, E and F), and D18 in the A domain. As a result, it would become possible to bend both AC and AB domain interactions. Thus, the B and C domains would make a rotation-like movement with respect to the A domain (Fig. 4 I).

The links identified as critical for establishing the B–A–C–A–C network were verified experimentally and in MD simulations as follows. We prepared two alanine-substituted mutants in the B–A–C–A–C network (E40A and D61A). Because E40 has tight interactions between R402 and S406, and because D61 has tight interactions between R283 and R402, the two mutants (E40A and D61A) are expected to have the same effect as observed in the examination of R402C/H. Two mutant RYR1 channels were heterologously expressed in HEK293 cells in the same way as described above, and expression was induced by the addition of doxycycline (2 µg/ml) and confirmed by Western blotting (Fig. S1 B). Then, measurements of caffeine-induced Ca²⁺ transients of the mutants were performed (Fig. 5, A–D).

As shown in Fig. 5 A, two mutants showed varied responses to caffeine as compared with the WT. As in the analysis of Fig. 2, peaks of the Ca²⁺ transient at 0.3 mM and 10 mM caffeine were compared in Fig. 5, B and C. Both mutants (E40A and D61A) exhibited significant increase in response to 0.3 mM caffeine, indicating marked enhancement of caffeine sensitivity (Fig. 5 B). On the other hand, the peak response at 10 mM caffeine decreased in both mutants (Fig. 5 C). Both mutants caused significantly higher levels of resting [Ca²⁺]_cyt (58–72 nM) than the resting [Ca²⁺]_cyt (40 nM; Fig. 5 D). Indeed, these series of results were consistent with the results seen with R402C/H.

We next performed MD simulations (50 ns) of the two mutants (Fig. 5, E–K). The molecules remained stable during the simulation, as the RMSD value of Cα atoms throughout the MD simulation was <3 Å (Fig. 5 E). DynDom was used to analyze the domain movements in the two mutants after simulation. Fig. 5, H and I, displays representative results from 20 calculations. When we analyzed the movement of the B and C domains with respect to the A domain, the BC domains in D61A and E40A rotated 13.5 ± 1.1 degrees (n = 20) and 8.0 ± 0.9 degrees (n = 20) clockwise, respectively (Fig. 5, H and I; and Fig. S3 E). Indeed, these rotations of the BC domains were consistent with that observed in the R402C/H mutants. Interestingly, each rotation angle that was observed in the four mutants (D61A, E40A, R402C, and R402H) correlated with the measurement results of resting [Ca²⁺]_cyt in each mutant (Fig. 5 D). To address the molecular mechanism of these differences, the boundaries of A, B, and C domains of the three mutants were analyzed (Fig. 5, F, G, J, and K). In D61A, the interaction among A, B, and C domains via R402 (R402–D61 and R402–E41) no longer existed, due to the lack of D61, although two interactions between A and C (R402–E40 and E40–S406) remained intact. Indeed, this circumstance was almost the same as that observed in the R402C/H mutants. In E40A, as two interactions between A and C domains (R402–E40 and E40–S406) were lost, the interactions between A and C domains became weaker. However, the interactions among A, B, and C domains via R402 (R402–D61 and R402–E41) were retained. Therefore, it was considered that the rotation angle would become rather smaller in the E40A mutant compared with the other three mutants (D61A, R402C, and R402H; Fig. S3 E).

The NT region of the RYR1 located at the cytoplasm is ∼100 Å away from the channel pore. Nevertheless, the region contains a cluster of point mutations that cause MH and CCD. How does a single mutation at the NT region have a long-range effect on the opening of the channel pore? To clarify this question, a combination of functional studies and MD simulations of RYR1 channels carrying disease-associated mutations at the NT region were performed in this study.

We found that the B(R283)–A(D61)–C(R402)–A(E40)–C(S406) network is critical for the connection of A, B, and C domains. Indeed, RYR1 B–A–C–A–C network mutants showed the abnormality in the caffeine-induced Ca²⁺ transient. MD simulations showed that the mutations loosen the connection between A and BC domains. The next question is whether the network is conserved in other subtypes. Comparison of the primary structures at the NTD revealed that most residues involved in the network formation are well conserved in RYR2 and RYR3: B(R298)–A(D61)–C(R417)–A(E40)–C(S421) and B(R287)–A(D62)–C(R409)–A(E41)–C(N413) for human RYR2 and RYR3, respectively. In addition, the B-A-C-A-C network is found in the crystal structures of the NTD of RYR2 (Borko et al., 2014; Fig. S6, A and B). There is an another reported crystal structure of the NTD of RYR2 (Kimlicka et al., 2013b). However this area is substantially different. There is a Cl⁻ bound at the interface and no B-A-C-A-C network. These findings suggest that the NTD in other subtypes may also play an important role in channel opening.

We performed MD simulations using the monomer structure of the NTD and found that B and C domains together rotated clockwise with respect to the A domain in the R402C/H mutants. However, because the RYR channel is a homotetramer, consideration of the movement in the tetramer structure is required. We therefore created a hypothetical homotetramer of NTD by fitting monomer structures before or after the MD simulations to the cryo-EM structure of RYR1 in the closed state (PDB accession no. 5TB0) at the A domain of the NTD monomer using ProFit. Fig. 6, A and B, shows the created hypothetical homotetramer before and after 50 ns of the simulation, respectively. The overall structure was the same before and after the simulation, verifying our MD simulation. Hypothetical homotetramers of the NTD of R402C (Fig. 6 C) and R402H (Fig. 6 D) mutants were also created using the same method as for the WT. Interestingly, the size of the central orifice of the NTD is larger in the mutants (Fig. 6 C, R402C, ∼34 Å; Fig. 6 D, R402H, ∼36 Å) than that of the WT (Fig. 6 B, ∼27 Å).

To interpret the findings, we also compared the cryo-EM structures of NTD in the closed state (PDB accession no. 5TB0; Fig. 6 E, left) and the open state (PDB accession no. 5T9V; Fig. 6 E, right). The diameter of the central orifice in the open state (∼29 Å) was larger than that in the closed state (∼24 Å). Importantly, there was essentially no change in the monomer structure of the NTD in the closed and open states (RMSD of Cα, ∼0.7 Å; Fig. 6 F). Upon widening of the central orifice, NTD translationally moves outward from the molecule, as indicated by the arrow in Fig. 6 G. Structural changes for widening of the central orifice might be related to the channel activity. Thus, it may be possible that the structure of the R402 is closer to the structure in the open state, but this mechanism may differ from that in the channel opening in WT.

Recently, two papers describing MD simulations of the NTD have been reported, and they used the NTD tetramer based on the cryo-EM structure (Zheng and Liu, 2017; Xiong et al., 2018). Therefore, we also performed MD simulations of WT and R402C/H mutants using the NTD tetramer derived from cryo-EM structure (PDB accession no. 5GKZ) for 60 ns (Figs. S7 and S8). Unexpectedly, in both the simulations of WT and mutants, the RMSD value did not completely reach the plateau even after the 60-ns calculation (Figs. S7 A and S8 A), whereas the RMSD value of each protomer reached a plateau (∼2 Å) at 25 ns (Figs. S7 A and S8 A). This finding suggests that the structure of each protomer was not broken, but that the fourfold symmetry in the tetramer would be broken in the simulations. In the case of WT, there is not much difference in the overall structure between before (0 ns) and after 60 ns of the simulation (Fig. S7 B). However, each protomer is translationally and rotationally moved after simulation (Fig. S7 C). In the case of the mutants, in contrast, the structure of each protomer after 60 ns of simulation does not match that at 0 ns of the WT (Fig. S8 C). The BC domain rotated clockwise with respect to the A domain, although the angle in each monomer was widely distributed (Fig. S8 C). This observation is consistent with our results with monomers of the NTD (Fig. 3 B and Fig. S3 B). Hydrogen bonds/salt bridges within each protomer of the mutants were essentially the same as the results using the monomer (Figs. S7 D and S8 D). These results indicate that the MD simulations themselves worked well and that what we observed in the simulations using the monomer occurred with the simulations using tetramer.

We assume that the reason for wide distribution of the angle could be the loss of cooperative movement between the monomer NTDs. In the RYR1 homotetramer, there are many communications between each monomer through transmembrane regions and other associating domains, which enable cooperative movement of NTD. However, the structure that we used in our simulations lacks any of these regions. Thus, each monomer moves independently without cooperativity. The ultimate way to resolve this problem could be calculations using the whole molecule. However, owing to the extremely large size of the molecule (a total molecular weight of 2.2 MDa) and limited resources, performing the calculation correctly is beyond the scope of our current research. An additional problem is incompleteness of the current cryo-EM model: ~15% of amino acid residues were not well identified, and ~19% of the whole model has not been built, due to its poor density.

How do the R402C/H mutants increase the open probability of the channel? A clue to understanding this mechanism is comprehension of the linkage between the movement of the NTD and channel pore in the WT structure. Because nearly the entire structures of RYR1 in the closed and open states solved by cryo-EM are available, we analyzed the movements of each domain when moving from the closed state to the open state. Fig. 7 (left) shows the relationship between the NTD and channel pore movement in the WT as represented by a cartoon. Along with the movement of S6, which constitutes a channel pore, the following four movements are linked, as indicated by the arrow in Fig. 7 (left): (1) Central, which is located above S6; (2) Handle–HD1–P2–HD2; (3) SPRY1–P1–SPRY2–SPRY3; and (4) the NTD. Although the exact sequence of the series of movements is not fully clarified, it is certain that the movement of S6 is related to the movement of the NTD. As shown in Fig. 4 I, the position of the C domain moves slightly upward by the clockwise rotation of the BC domains in the R402H/C mutants. Interestingly, in terms of moving up for the C domain, a similar movement takes place in the WT during the opening of the channel pore. If so, the upward movement of the C domain in the R402C/H mutants may lead to the movement of Handle–HD1–P2–HD2, SPRY1–P1–SPRY2–SPRY3, and Central, the latter of which is located above S6. Therefore, this series of movements may ultimately increase the mobility of S6 and raise the open probability of a channel pore (Fig. 7, right).

In our MD simulations, we could not detect any effect of the mutations located around the subunit interfaces (e.g., Q156K, R164C, R164L, G342R, Y523C, and Y523S; Fig. S4). It has been reported that the mutations located around subunit interfaces may lead to reduced instability in the NTD tetramer, which may in turn serve to be structurally conducive to the open state and, simultaneously, a hindrance to the closed state (Zheng and Liu, 2017). These mutations might also affect the series of movements of the respective domains.

IP₃R is a ubiquitous ion channel and localized in the ER/SR, as is RYR. It is responsible for cytosolic Ca²⁺ signaling and has crucial roles in a variety of physiological functions (Iino, 1990; Bezprozvanny et al., 1991; Berridge et al., 2003). In IP₃R, inositol-1,4,5-triphosphate (IP₃) binds to the NT region, which is located far from the channel pore and regulates the opening of the channel pore. This long-range regulation of the channel pore by the NT region in the IP₃R is also common to RYR. Although IP₃R and RYR have many structural similarities, the total number of amino acids in one monomer of IP₃R is 2,700, which is ∼55% of that in RYR. Thus, the structure of IP₃R is considerably simple compared with that of RYR.

The x-ray crystal structures of NT of the IP₃R in apo and IP₃-bound states have already been determined, and the structural resemblance between the NT of the IP₃R and the NTD of RYR1 has already been discussed (Seo et al., 2012). The x-ray crystal structure is composed of three domains: (1) the β-trefoil domain (βTF1) domain acting as a suppressor and (2) βTF2 and (3) armadillo solenoid folds (ARM) domains involved in the binding of IP₃. Surprisingly, the structures of the NT region (IP₃R) and NTD (RYR) are very similar, and the superposition of both structures is in good agreement (Seo et al., 2012; Fig. S6 C). In IP₃R, binding of IP₃ to the binding site causes movement of both βTF2 and ARM domains (corresponding to the B and C domains in the RYR), and the movement is considered to trigger the opening of the channel pore in the transmembrane region. In addition, a similar interaction between AC domains via R402 found in RYR (E40[A]–R402[C]–D61[A]) also exists in IP₃R between βTF1 and ARM (corresponding to the A and C domains in the RYR). However, in the case of the IP₃R, the charges of the three side chains are totally opposite, and those corresponding to E40(A), R402(C), and D61(A) are K127(βTF1), D444(ARM), and R54(βTF1), respectively (Seo et al., 2012). In IP₃R, the structural studies of the full-length receptor by cryo-EM have already been reported (Fan et al., 2015, 2018; Paknejad and Hite, 2018), and these studies also support the structural change of the NT caused by IP₃ being the key to opening the channel pore in the transmembrane region.

In this study, by combining [Ca²⁺]_cyt measurements and MD simulations, we have shown a structural change in RYR1 with gain-of-function mutations at the domain interface of the NTD. The same line of analysis on the interdomain mutations in the NTD may lead to a deeper understanding of the mechanism enhancing channel opening.

Eduardo Ríos served as editor.

We thank Dr. Tadashi Takemori (Tsukuba University) for valuable discussion and technical advice regarding conduct of the MD calculations, Dr. Shigeru Takemori (The Jikei University) for valuable comments, and Mariko Hanano and Shizuko Syoji (The Jikei University) for their technical assistance.

This work was supported by Japan Society for the Promotion of Science Grants-in-Aid for Scientific Research (grants 19H03198 and 19K07306 [T. Yamazawa], JP16H04748 [H. Ogawa], 19H03404 [T. Murayama], 16K08917 [K. Oguchi and H. Oyamada], 19K07105 [N. Kurebayashi], 15H05648 [K. Kanemaru], and 21229004 and 25221304 [M. Iino]), the Platform Project for Supporting Drug Discovery and Life Science Research (Basis for Supporting Innovative Drug Discovery and Life Science Research [BINDS] grants JP18am0101080 and JP19am0101080 [H. Ogawa and T. Murayama]), a grant from the Nakatani Foundation for Advancement of Measuring Technologies in Biomedical Engineering and The Jikei University School of Medicine (T. Yamazawa), and a National Center of Neurology and Psychiatry intramural research grant (29-4) for neurological and psychiatric disorders (T. Murayama).

The authors declare no competing financial interests.

Author contributions: T. Yamazawa and T. Murayama carried out the cell biological experiments. N. Kurebayashi, H. Oyamada, J. Suzuki, and K. Kanemaru provided experimental tools. T. Yamazawa, T. Murayama, K. Oguchi, T. Sakurai, and M. Iino analyzed the data. T. Yamazawa, M. Yamaguchi, and H. Ogawa performed the MD simulation and analysis. T. Yamazawa and H. Ogawa wrote the manuscript. All authors discussed the results and approved the final version of the manuscript.