Understanding the mechanistic basis of non-coding RNA through molecular dynamics simulations

Abstract

Noncoding RNA (ncRNA) has a key role in regulating gene expression, mediating fundamental processes and diseases via a variety of yet unknown mechanisms. Here, we review recent applications of conventional and enhanced Molecular Dynamics (MD) simulations methods to address the mechanistic function of large biomolecular systems that are tightly involved in the ncRNA function and that are of key importance in life sciences. This compendium focuses of three biomolecular systems, namely the CRISPR-Cas9 genome editing machinery, group II intron ribozyme and the ribonucleoprotein complex of the spliceosome, which edit and process ncRNA. We show how the application of a novel accelerated MD simulations method has been key in disclosing the conformational transitions underlying RNA binding in the CRISPR-Cas9 complex, suggesting a mechanism for RNA recruitment and clarifying the conformational changes required for attaining genome editing. As well, we discuss the use of mixed quantum-classical MD simulations in deciphering the catalytic mechanism of RNA splicing as operated by group II intron ribozyme, one of the largest ncRNA structures crystallized so far. Finally, we debate the future challenges and opportunities in the field, discussing the recent application of MD simulations for unraveling the functional biophysics of the spliceosome, a multi-mega Dalton complex of proteins and small nuclear RNAs that performs RNA splicing in humans. This showcase of applications highlights the current talent of MD simulations to dissect atomic-level details of complex biomolecular systems instrumental for the design of finely engineered genome editing machines. As well, this review aims at inspiring future investigations of several other ncRNA regulatory systems, such as micro and small interfering RNAs, which achieve their function and specificity using RNA-based recognition and targeting strategies.

Keywords: CRISPR-Cas9; Genome editing; Group II intron; RNA splicing; Spliceosome.

Fig. 1

Conformation of the Cas9 protein as apo form (left, (Jinek et al., 2014)) and bound to the RNA (right, (Jiang et al., 2015)). Cas9 is shown in molecular surface, highlighting the protein domains with different colors. The RNA (yellow) is shown as ribbons. The protein is composed by a recognition lobe (REC), which includes the REC1–3 regions, and by a nuclease lobe (NUC), including the HNH and RuvC catalytic domains and the PAM interacting (PI) region. Representative snapshots along the conformational change from the apo Cas9 (“open” conformational state) to the RNA-bound Cas9 (“closed” conformation) are shown. During the conformational transition, the REC1–3 regions collectively move in opposite directions with respect to each other (as indicated using arrows), with the exposure toward the solvent of a R-rich α-helix (magenta) (Palermo et al., 2017a). In its apo form, Cas9 adopts an auto-inhibited conformation, in which the REC and NUC lobes assume an “open” conformational state. Upon RNA binding, a dramatic conformational transition is observed within the REC lobe, which leads to the formation of a “closed” conformational state of the protein (Fig. 1).

Fig. 2

(A) Potential of mean force (PMF) computed with respect to the E945–D435 distance and to the Root Mean Square Deviation (RMSD) from the X-ray structure of the apo Cas9 (Palermo et al., 2017a). The E945–D435 distance has been employed in Förster Resonance Energy Transfer (FRET) experiments to describe the conformational changes of the protein during RNA binding (Sternberg et al., 2015), allowing the validation of the identified pathway. The PMF reveals three local minima, corresponding to the crystallographic apo (M1) and RNA-bound (M2) states, and to an intermediate state (M3), characterized by the solvent exposure of the arginine residues in the R-rich α-helix. (B) From left to the top right: the electrostatic potential of the Cas9 protein is computed along the conformational transition from the apo protein up to RNA binding (Baker et al., 2001). The positive (blue) and negative (red) potential field strengths are captured at an iso-density of ± 200 kT/e. This shows the formation of a positively charged cavity, upon the exposure of the R-rich α-helix in the pre-RNA bound conformation. Adapted with permission from Palermo et al. (2017a). Copyright 2017 National Academy of Sciences.

Fig. 3

(A) Conformation of the catalytic HNH domain in the X-ray structure of the pre-activated CRISPR-Cas9 complex (5F9R.pdb) (Jiang et al., 2016). The catalytic residue H840 (magenta) locates at ∼19.4 Å from the cleavage site (CS) on the target strand (TS). (B) Close-up view on the HNH domain docked at the cleavage site on the TS in the activated CRISPR-Cas9 complex. (C) Distance of the catalytic residue H840 from the cleavage site (CS) over ∼400 ns of conventional MD (cMD, gray) and over three replicas of Gaussian accelerated MD (GaMD, warm colors) (Palermo et al., 2017a). A black bar indicates the value of the distance in the 5F9R X-ray structure (Jiang et al., 2016).

Fig. 4

(A) Conformation of the three regions of the REC lobe in the activated CRISPR-Cas9 complex (Palermo et al., 2018). REC3 and REC2 show a remarkable outward translation, while REC1 moves toward the HNH domain, suggesting that the REC lobe ensures the catalytic competence of Cas9 by “sensing” (REC3), “regulating” (REC2) and “locking” (REC1) the catalytic HNH domain. Three arrows indicate the conformational change from the pre-activated to the activated state. (B) Correlation analysis including the matrix of the generalized correlations (GC, upper triangle), quantifying the coupled motions among residue pairs, and its coarse representation (lower triangle), accumulating the per-domain correlations (Palermo et al., 2017b). High correlations (green) are observed between the REC-HNH and the HNH-RuvC domains. Adapted with permission from Palermo et al. (2018). Copyright 2017 American Chemical Society. https://pubs.acs.org/doi/10.1021/jacs.7b05313. Further permissions related to the material excerpted should be directed to the American Chemical Society.

Fig. 5

(A) Structure of the Oceanobacillus iheyensis (O.i.) group IIC intron before splicing (4FAQ) (Marcia and Pyle, 2012). The intron (green) and the exon (blue) are depicted in ribbons. (B) Close-up view of the active site, displaying two catalytic Mg²⁺ ions (orange spheres), the water nucleophile and the surrounding RNA. MD simulations have been performed using the QM/MM approach, treating the active site at a QM (DFT/BLYP) level of theory, while the rest of the system in explicit solution is treated at the classical MM Amber force field. Representative snapshots along the phosphodiester bond cleavage: the reactant (R), transition state (TS^‡), deprotonation event (DEP) and product (P). (C) Free energy profiles (ΔF, in kcal/mol) for the chemical step, computed via thermodynamic integration, using as a reaction coordinate (RC) the distance between the breaking (O3′–P_RNA) and forming (O_WAT–P_RNA) bonds (Casalino et al., 2016). Adapted with permission from Casalino et al. (2016). Copyright American Chemical Society. https://pubs.acs.org/doi/abs/10.1021/jacs.6b01363. Further permissions related to the material excerpted should be directed to the American Chemical Society.

Fig. 6

Schematic representation of the two-metal aided mechanism for phosphodiester bond cleavage, as observed in the RNase H protein (A) (Rosta et al., 2011; De Vivo et al., 2008; Palermo et al., 2015a,b) and in group II intron (G2I) ribozyme (B) (Casalino et al., 2016). In the RNase H protein, the reaction proceeds through an associative mechanism S_N2 where the nucleophile is bonded to the phosphorus atom before the leaving group is eliminated. In G2I ribozyme, a dissociative mechanism S_N1 is observed, in which the elimination of the leaving group precedes the nucleophile attacks to an unstable meta-phosphate intermediate. The Reactant (R), Transition State (TS^‡) and Product (P) structures are shown for both mechanisms.

Fig. 7

(A) Intron Lariat Spliceosome from the yeast Schizosaccharomyces pombe (S. Pombe, 3JB9.pdb) (Hang et al., 2015; Yan et al., 2015). The protein components are shown in molecular surface, while the snRNAs and the Intron Lariat are depicted as cartoons using different colors. (B) The per-domain correlations are accumulated and plotted as a histogram (reporting the per-column correlation scores calculated for each pair of spliceosomal components). Correlations are based on Pearson coefficients (Pc) and are shown using blue (0 ≤ Pc ≤ 1) and red (1 ≤ Pc ≤ 0). Boxes are used to highlight the correlations of the Spp42 protein with the snRNAs (U2, U5, U6) and with the intron lariat (Casalino et al., 2018). (C) Displacement (highlighted by a black arrow) of the intron lariat/U2 helix, as revealed from principal component analysis of the MD trajectories. The first principal component is computed and plotted on the 3D structure indicating the direction of large-scale motions (shown using arrows). (D) The “polar tweezers” formed by K364, K366, and R388 residues of Cwf19. The electrostatic potential (ESP) is computed and plotted over the molecular surface of the protein components (blue: positive; red: negative) (Baker et al., 2001). Adapted with permission from Casalino et al. (2018). Copyright National Academy of Sciences.

Fig. 8

Schematic illustration of the GaMD metod (Miao et al., 2015). When the system potential $V(\overrightarrow{r})$ is lower than a threshold energy E (E = V_max), the energy surface is modified by adding a harmonic boost potential that follows Gaussian distribution. The harmonic constant k₀ determines the magnitude of the applied boost potential. By increasing k₀ in the range from 0 to 1, higher boost potential is added to the original potential energy surface.

Fig. 9

Schematic illustration of the electrostatic coupling scheme developed by Rothlisberger and coworkers (Laio et al., 2002a,b), to treat the electrostatic interactions between QM and MM atoms in a mixed QM/MM approach.