Comparative Visualization of the RNA Suboptimal Conformational Ensemble In Vivo.

When a ribonucleic acid (RNA) molecule folds, it often does not adopt a single, well-defined conformation. The folding energy landscape of an RNA is highly dependent on its nucleotide sequence and molecular environment. Cellular molecules sometimes alter the energy landscape, thereby changing the ensemble of likely low-energy conformations. The effects of these energy landscape changes on the conformational ensemble are particularly challenging to visualize for large RNAs. We have created a robust approach for visualizing the conformational ensemble of RNAs that is well suited for in vitro versus in vivo comparisons. Our method creates a stable map of conformational space for a given RNA sequence. We first identify single point mutations in the RNA that maximally sample suboptimal conformational space based on the ensemble's partition function. Then, we cluster these diverse ensembles to identify the most diverse partition functions for Boltzmann stochastic sampling. By using, to our knowledge, a novel nestedness distance metric, we iteratively add mutant suboptimal ensembles to converge on a stable 2D map of conformational space. We then compute the selective 2' hydroxyl acylation by primer extension (SHAPE)-directed ensemble for the RNA folding under different conditions, and we project these ensembles on the map to visualize. To validate our approach, we established a conformational map of the Vibrio vulnificus add adenine riboswitch that reveals five classes of structures. In the presence of adenine, projection of the SHAPE-directed sampling correctly identified the on-conformation; without the ligand, only off-conformations were visualized. We also collected the whole-transcript in vitro and in vivo SHAPE-MaP for human β-actin messenger RNA that revealed similar global folds in both conditions. Nonetheless, a comparison of in vitro and in vivo data revealed that specific regions exhibited significantly different SHAPE-MaP profiles indicative of structural rearrangements, including rearrangement consistent with binding of the zipcode protein in a region distal to the stop codon.

Introduction

Ribonucleic acid (RNA) 3D structures are the result of remarkably complex interaction networks that together create emergent biological functions (1, 2, 3, 4). Although crystal structures reveal these networks with atomic detail, these remain static snapshot models of the conformations existing in the cellular environment (5). RNAs, particularly highly structured RNAs such as ribosomal RNA, exist in multiple conformations, many of which are likely to affect their function(s) (6, 7, 8). Thus, when describing RNA structure, it is more accurate to discuss an ensemble of conformations instead of a single structure (7, 9, 10, 11). However, significant biophysical challenges remain, whether at the secondary or tertiary structural level, including visualization of the ensemble of RNA conformations and identification of essential functional elements within the entire ensemble (9, 12, 13, 14).

The challenge of visualizing an RNA secondary structure ensemble is easily illustrated by the Vibrio vulnificus adenosine deaminase (add) adenine riboswitch (Fig. 1) (15, 16, 17, 18). Typically RNA is represented as a single structure, but, for the riboswitch, at least two structures are required for function: the on-conformation and the off-conformation (Fig. 1 A) (16, 18, 19). These two structures interchange, with the off-conformation favored without the adenine ligand, and the on-conformation stabilized by binding adenine (17, 18, 20). Thus, in solution the RNA exists as an ensemble of conformations that interchange (1, 8, 10, 21, 22, 23). In visualizing such an ensemble, two salient aspects should be highlighted to understand function: 1) the structural similarity and difference between the two conformations and 2) the relative abundance of each conformation in the ensemble.

Figure 1

The conformational states of the V. vulnificus add adenine riboswitch. (A) The accepted structures for the bound and unbound states of the riboswitch are determined by crystallography and NMR (54). The unbound state represses translation, and the bound state activates translation (20, 54). (B) The map of conformational space explores five possible structure clusters for the riboswitch. The representative structure is the cluster medoid structure. The orange cluster represents the translation off-conformation, and the purple cluster represents the translation on-conformation, as confirmed by crystallography and NMR (90). To see this figure in color, go online.

Defining structural similarity requires a representation that captures biologically important structural features of the RNA to facilitate clustering of highly similar conformations. From these clusters, it is then possible to determine the relative abundance of the conformations, which reflects their relative thermodynamic weights in the Boltzmann ensemble. We therefore aim to create a visualization based on a sampling of conformational space like the one illustrated for the add riboswitch (Fig. 1 B), which was stochastically sampled from the Boltzmann ensemble. In Fig. 1 B, we illustrate a map of conformational space, in which each square represents a cluster of similar conformations based on a nested feature vector that we define below. This representation is particularly interesting as it reveals several aspects of the add riboswitch conformational ensemble that are not apparent when considering only two structures (Fig. 1 A). First, this visualization suggests that there are more than two classes of conformations in the add riboswitch conformational ensemble. Second, the on- and off-conformational change is conveniently captured along dimension 1. The methods we describe below provide a robust approach for identifying specific dimensions that capture biologically informative structural differences, such as those in Fig. 1 B.

In Fig. 1 B, we purposely did not indicate the relative abundance of conformations in each conformational cluster; each square is equal in size. The relative weight of these clusters depends on the underlying thermodynamic parameters of the energy model. Given a nearest-neighbor energy model, it is now computationally efficient to rapidly sample the Boltzmann suboptimal ensemble (24, 25, 26, 27). Furthermore, the nearest neighbor model can be extended to empirically include experimental structure probing data, particularly selective 2′ hydroxyl acylation by primer extension (SHAPE) data (28, 29). Inclusion of SHAPE data is relevant because the RNA structure is readily probed under different experimental conditions. For example, the add riboswitch can be probed with and without the ligand that causes a structural rearrangement (15, 30, 31). As we will show below, the visualization proposed in Fig. 1 B accurately captures this biologically important rearrangement when combined with SHAPE-informed structure probing.

Although visualizing riboswitch ensemble conformations is one important goal of our work, the main motivation for improving the ability to visualize and interpret RNA conformational ensembles stems from our studies of messenger RNA (mRNA) folding in vitro versus in vivo. Quantitative comparison of these two conditions effectively enables us to deconvolute the effect of the cellular environment on mRNA folding. The structural ensembles of these long and flexible RNAs tend to be far more complex than the structural ensembles of riboswitches. As such, we require tools that enable “sorting the forest from the trees” to understand these large and complex molecules. We present here an experimental high-resolution comparison of SHAPE data for the human β-actin mRNA that reveals specific regions in which the RNA folds differently in vitro versus in vivo. We show how these visualizations enable interpretation of the complex rearrangements of the mRNA conformational ensemble that occur in the cell, thereby obtaining meaningful biophysical and biological insight into the specific structure function relationships of the specific messenger. Together, these novel data and methods, to our knowledge, establish a robust approach for interpreting chemical and enzymatic probing data in the context of conformational ensembles.

Materials and Methods

Generating structures for the map of conformational space

Our strategy for establishing a conformational map of an RNA ensemble is illustrated in Fig. 2. Beginning with the RNA sequence (Fig. 2 A), we compute its partition function (probability of basepairing (32, 33, 34, 35)) and the partition functions of all AtoU, UtoA, CtoG, and GtoC single point mutant sequences (Fig. 2 B). These point mutations are experimentally determined to be maximally disruptive of structure (32). The purpose of stochastic sampling of multiple single point mutant sequences is to generate a more diverse ensemble of structures from which to build a visualization space. This strategy converges faster and generates a more diverse ensemble than traditional stochastic sampling of a single sequence, as can be seen in Fig. S7. The sum over the rows in the partition function is the basepairing probability, P, for each nucleotide with every other nucleotide, x (Eq. 1). Our goal is to generate an ensemble of diverse possible conformations and establish a representative 2D map for visualization. Thus, single point mutants with the highest ensemble Shannon entropy (H), as defined by Eq. 1, are selected for further analysis. This definition of Shannon entropy has previously been used to interpret RNA structure (36, 37, 38), and computes the entropy based on the 1D basepairing probability vector. Alternative definitions of Shannon entropy could potentially be used to compute the Shannon entropy from the full 2D basepairing probability matrix, the thermodynamic structural entropy, or algorithm computation (33, 34, 37). In the first pass, we eliminate the lowest 25% Shannon entropy mutants (Fig. 2 C) (36, 37, 38). In a second filter, we perform hierarchical clustering of the basepairing probability P(x) vectors based on their Euclidean distance (39) to identify the most divergent partition functions (Fig. 2 D). We then perform Boltzmann stochastic sampling on the two most divergent partition functions (Fig. 2 E), and create nestedness feature vectors from the sampled structures (Fig. 2 F; Fig. S1), to generate a map of conformational space using metric multidimensional scaling (40) (Fig. 2 H). We iteratively add additional Boltzmann ensemble samples of divergent single point mutant sequences until the map of conformational space converges (Fig. 2 G), as follows:

Figure 2

Building the map of conformational space. The map explores the possible structural space for an RNA sequence and its single point mutants. (A) A single point mutant was created for every position in the RNA. We used only mutations that were expected to lead to the largest changes in structure based on experimental observations from the mutate-and-map experiments (AtoU, UtoA, CtoG, and GtoC) (32, 91). (B) The partition function was generated for the wild-type and single point mutants using established structure prediction methods (22, 49, 50, 92). (C) The RNAs were ranked by Shannon entropy, and the top 75% were retained to filter for individual RNAs with more diverse ensembles (36, 37, 38). (D) We collapsed the partition function for each of the remaining RNAs into their basepairing probabilities, and performed hierarchical clustering on the probabilities (39). This clustering selects the most diverse RNA subsets. (E) We selected the most distant RNA and sampled 1000 structures according to their Boltzmann probability (5). (F) We used data abstraction to identify the number of unique structure clusters. This data abstraction is further described in Fig. S1. (G) We repeated steps (E) and (F) until the number of structure clusters converged. (H) The structure clusters are projected into 2D space using metric MDS. By minimizing the stress function for the Euclidean distance matrix, MDS optimizes the positioning of the structure clusters (40, 41). To see this figure in color, go online.

Projection of the map of conformational space

Our projection is based on the representation proposed in the RNAshapes abstraction that captures whether a stem or stack element exists, ignoring the size of that element (35). Biologically, significant variation is observed in stem length but stack elements are generally more conserved (33, 34, 35). Thus, we expect that basing our projections on this distance metric will capture important structure/function features in the ensembles. Our representation counts the number of inner loops and stacks and then positions that count according to the location of the outermost stack in the nestedness feature vector (Fig. S1). Stems and stacks with fewer than three basepairs can be optionally ignored to simplify particularly complex ensemble visualizations. We determine the nestedness representation for every structure in the map of conformational space and collapse the structures into clusters based on unique nestedness representations (Fig. 2 F). Metric multidimensional scaling (MDS) projects the structure clusters into 2D space by finding the positioning of points in 2D space that best recapitulates the original Euclidean distances calculated from the structure cluster representations (Fig. 2 H) (40, 41). MDS calculates the Euclidean distance matrix for n-dimensional data, d. Initial positions for the data points, x, are set in 2D space, i and j. The initial positions are determined from projection onto the first two eigenvectors from eigen decomposition on the distance matrix. Based on this configuration, MDS evaluates the stress function in Eq. 2 (40, 41). The data points are reconfigured in the direction of steepest descent. This process is repeated to minimize the stress function (40, 41). Minimization of the stress function finds the configuration with the smallest residual sum of squares when compared with the original distance matrix (40, 41). As a result, MDS yields a 2D embedding of the data points (used for visualization) that optimally reflects the pairwise distances between data points as computed within the original n-dimensional data, as follows:

Projection of the wild-type RNA

For the wild-type RNA, we recommend generating 1000 structures using Boltzmann-weighted stochastic sampling (42, 43, 44, 45) (Fig. S2 A). SHAPE data can be included to direct the ensemble prediction (28, 46). Each structure from the wild-type ensemble was converted into our nestedness representation. We then compute the frequency of structures that belong to each structure cluster (Fig. S2 B), and these frequencies are then scaled as relative diameters for each bubble in the resulting plot (Fig. S2 C). Bubbles are colored according to their 2D distances, where groups of similar clusters are closer on the viridis color scale in matplotlib (47). If any sampled structure in this plot does not match existing clusters in the map, the structure is added to the closest cluster based on Euclidean distance. As described here, the wild-type RNA is projected onto the map of conformational space. For comparison between mutants or between the same RNA in different environments, the same map of conformational space is used (as opposed to recomputing a new map for every comparison). This results in a stable space for projecting new ensembles of interest. In the interactive visualization output for EnsembleRNA, we include a measure of diversity for each structure cluster to allow the user to get a sense of the similarity of the structures clustered. This measure compares the frequency of the most common structure (maximum cluster frequency) and the average Jaccard similarity (48) between the binary representations of structures (minimum cluster correlation). Thus, if every structure in a cluster is unique, the value is 1; otherwise diversity is the fraction of nonunique structures.

EnsembleRNA package and webserver

A Python package (https://www.python.org/), EnsembleRNA, was created for the visualization of RNA structural ensembles. The package produces bubble charts for the map of conformational space and the wild-type RNA, and allows for comparison between structural ensembles. The package is available at http://ribosnitch.bio.unc.edu/software. Supporting Materials and Methods contains additional information on usage, troubleshooting EnsembleRNA, and tutorials.

In vitro SHAPE treatment

SHAPE-MaP experiments were performed in vitro (37). We obtained a clone of β-actin mRNA (SC319328; OriGene, Rockville, MD) and directly PCR-amplified the coding sequence with a 5′ primer containing the T7 promoter for in vitro transcription (Q5 Site-Directed Mutagenesis Kit and T7 RNA Polymerase from New England BioLabs, Ipswich, MA). To remove DNA after transcription, we treated the reaction with TURBO DNase for 15 min at 37°C (ThermoFisher Scientific, Waltham, MA). Standard bead cleanup was performed between each step (Ampure XP; Beckman Coulter, Brea, CA). The transcribed RNA was folded at 37°C in buffer containing 100 mM Na-HEPES, pH 8.0, 100 mM NaCl, and 10 mM MgCl2. One μg of RNA was treated for 5 min with either 10% dimethyl sulfoxide (DMSO) or DMSO containing the RNA modifying agent 1-methyl-7-nitroisatoic anhydride (1M7) at a final concentration of 10 mM.

In vivo SHAPE treatment

We performed in vivo SHAPE-MaP experiments for β-actin in the 1000 Genome Cell Lines GM07037 and GM12003 (37), obtained from the NIGMS Human Genetic Cell Repository at the Coriell Institute for Medical Research (https://www.coriell.org/). Approximately 50,000,000 cells were collected by centrifugation, resuspended in 1 mL of folding buffer (as in in vitro SHAPE protocol) supplemented with 400 U murine RNase inhibitor, and sonicated three times at 10% power for 10 s (Sonic Dismembrator Model 500; Thermo Fisher Scientific). These samples were incubated at 37°C for 10 min, after which either DMSO (10% final concentration) or 500 mM 1M7 in DMSO (final concentration 30 mM) was added for 5 min with three separate additions. RNA was isolated with Trizol reagent (ThermoFisher Scientific), followed by treatment with TURBO DNase and removal of the majority of ribosomal RNA (RiboMinus Eukaryote System v2; Life Technologies, Carlsbad, CA).

SHAPE data collection and analysis

For all samples, we performed reverse transcription with the specialized reverse transcription conditions for SHAPE-MaP and random nonamer primers (37). The transcription reactions were purified via Ampure XP beads (Beckman Coulter) or G50 columns (GE Healthcare Life Sciences, Little Chalfont, Buckinghamshire, UK), and dsDNA was made by second strand synthesis (NEBNext mRNA Second Strand Synthesis Module; New England BioLabs). To prepare libraries, we used the Nextera or Nextera XT kits (Nextera DNA Sample Preparation Kit, Nextera XT DNA Sample Preparation Kit and Index Kits; Illumina, San Diego, CA). Sequencing for the in vitro samples was performed on the HiSeq 2500 (Illumina) as paired-end, 50-read multiplex runs. Sequencing for the in vivo samples was performed on the HiSeq 2500 as paired-end, 100-read multiplex runs. Analysis was performed with the ShapeMapper pipeline (37) using either β-actin mRNA (NM_001101) to align sequences derived from the in vitro samples or the entire genome (hg38) to align sequences derived from the in vivo alignment. The β-actin data are in the file SNRNASM (see Data S1). SHAPE traces for the wild-type V. vulnificus add riboswitch mutate-and-map experiments were obtained from the publicly available RNA Mapping Database (22, 49, 50). To normalize the SHAPE-MaP data, scaled background reactivities were subtracted from the plus reagent reaction reactivities. A multiplier was used to fit the resulting distribution of values to the distribution of values for the normalized reactivities of a reference mRNA.

β-actin RNA structural modeling

An RNA/protein complex was generated from a starting model of two DNA strands bound to the KH34 protein (37). A custom Python script was used to convert the DNA strands to the appropriate RNA nucleotide sequence. The resulting RNA/protein complex was equilibrated by discrete molecular dynamics (DMD) simulations (51, 52, 53) to accommodate the zipcode binding regions of the RNA strands. The remaining regions of the RNA strands were modeled using coarse-grained DMD simulations (42) in which each nucleotide was represented as three pseudo-atoms corresponding to the phosphate backbone, sugar group, and nucleobase. With the replica exchange approach, we efficiently sampled RNA conformations by utilizing replicas of the same RNA system in parallel at different temperatures. Replicas were allowed to exchange simulation temperatures periodically based on a Monte Carlo algorithm. The replica exchange DMD simulations were run for 50 ns with replica temperatures of 0.200, 0.225, 0.250, 0.270, 0.300, 0.333, 0.367, and 0.400 with units kcal/(mol∗kB). Free energy bonuses were incorporated between nucleotides to model the in vivo basepairing interactions. To select the final RNA model, we used a hierarchical clustering analysis based on the pairwise root mean square deviation of the phosphates and the potential energy as determined by the DMD force field. The coarse-grained RNA model was reconstructed to an all-atom model to combine with the KH34 protein system. We then equilibrated the entire RNA/protein complex using all-atom DMD simulations at a temperature of 0.4 kcal/(mol∗kB) and included static constraints on the protein and harmonic constraints on the zipcode binding regions of the RNA strand.

In vitro model

We incorporated the in vitro secondary structure as constraints in coarse-grained replica exchange DMD simulations, using the same settings as those in the in vivo RNA system. We then performed an root mean square deviation-based clustering analysis to determine the centroid and reconstructed an all-atom model at a temperature of 0.4 kcal/(mol∗kB).

RNA dynamics

The dynamics of the 2′ hydroxyl groups of the in vitro and in vivo RNA strands were calculated using the root mean square fluctuation (RMSF) with the Wordom software package (http://wordom.sourceforge.net/) (43). RMSF calculations were performed on three 100-ns DMD simulations at a temperature of 0.4 kcal/(mol∗kB) for both RNA systems. The in vivo system included static constraints on the protein and harmonic constraints on the zipcode binding protein-interacting regions of the RNA. We calculated the mean using 3-nucleotide windows and SD of the RMSF based on the three DMD simulations for each system.

Results

Generating a robust 2D representation of an RNA structural ensemble

Our first goal in creating a visualization of a structural ensemble was to establish a robust and consistent 2D representation of the conformational space of RNA. Traditionally, principal component analysis is used to identify two Eigenvectors for projection (24, 25). One challenge with this approach is that the first three Eigenvectors often fail to capture enough variance to detect major structural elements (44). If a conformation change is predicted, this limitation of principal component analysis makes it difficult to understand the relative differences in the ensemble. A second challenge is determining which structural features to highlight in the representation to capture important biological aspects of the ensemble. Selecting features to highlight requires picking a specific structural distance representation, which can affect the interpretation as much as which Eigenvectors are used for projection. We propose an approach that provides a stable and robust visualization while also capturing important biological features (e.g., the on- and off-conformation of the add riboswitch in Fig. 1 B).

Our approach is summarized in Figs. 2 and S2. We begin by computing the partition function of the wild-type RNA sequence and all single point mutants. From these partition functions, we select the RNAs that are maximally different, as determined by Shannon entropy and hierarchical clustering on basepairing probability (36, 39). From these partition functions, we sample the Boltzmann suboptimal ensemble and use these structures as the basis to build our visualization (25). This strategy effectively allows us to more comprehensively sample the suboptimal ensemble and the strategy does not depend on the approach used to compute the partition function. The visualization creates a stable space for the comparison of structural ensembles using mutations to explore the possible conformations that an RNA may take (Fig. 1 B). Data abstraction identifies clusters of similar structures that likely have similar function. This cluster representation reduces the map size, thereby creating a more accurate and interpretable visualization of secondary structure. Projecting the structure clusters into two dimensions using metric MDS optimizes their distances (40, 41). This approach enables easy interpretation of the visualization, in which clusters that are farther apart are more different. We can project the RNA ensemble of interest onto this space by varying the size of cluster bubbles based on the number of structures that belong to that cluster (Fig. 1 B). Experimental structure probing data can be included to guide the ensemble prediction (45). This method is further described in the Materials and Methods.

Detecting RNA structure change induced by ligand binding

The add riboswitch is found in the 5′UTR of an mRNA that codes for adenosine deaminase (20, 54). This riboswitch forms two distinct conformations that control translation of the adjacent coding region (20, 54). The adenine-unbound conformation represses translation, and the adenine-bound conformation activates translation. Fig. 1 A shows the accepted secondary structures for the unbound and bound states as determined by crystallography and NMR (54). These secondary structures represent only two of several possible conformations that the riboswitch may adopt in the cell (28, 46). Indeed, the map of conformational space produced by our visualization explores a total of five possible structure clusters including the two accepted conformations (Fig. 1 B). This visualization produces a separation in 2D space between conformations that can bind adenine and activate translation and conformations that cannot bind adenine.

The structural difference induced by ligand binding for the add riboswitch is particularly well suited for the application of SHAPE data. Without experimental data to guide structure prediction algorithms, the accepted bound conformation dominates (Fig. 3 A), and differences in structure that result from changes in environment cannot be discerned. However, including SHAPE data in the ensemble prediction algorithms reveals differences in the add riboswitch structure with and without ligand (Fig. 3, B and C). In each ensemble, the respective structure observed in crystallography and NMR dominates. Thus, our visualization approach combined with SHAPE-directed structural modeling captures key structural features of the ensemble (20, 54).

Figure 3

Visualization of bound and unbound states of the V. vulnificus add adenine riboswitch. (A) Projection of the predicted wild-type ensemble without SHAPE data favors the experimentally determined on-conformation (left). However, alternative conformations are still present (right). (B) When the ensemble generation is guided by SHAPE experiments conducted without ligand, off-conformations are favored in the projection (left). Particularly, the experimentally confirmed off-structure is the most populated conformation. (C) When SHAPE data are collected in the presence of ligand, the experimentally confirmed on-conformation (right) is preferred in the projection (left). Both SHAPE data sets (with and without ligand) are publicly available in the RNA Mapping Database (32, 48, 91). To see this figure in color, go online.

Observing regional structure differences in vitro and in vivo

β-actin is a cytoskeletal protein involved in cell motility and structure (15). The advent of high throughput structure probing methods such as SHAPE-MaP has only recently allowed us to collect information on larger RNAs such as the ∼2-kb β-actin mRNA (37). Structure probing data are collected for RNA in the presence of cellular components, e.g., RNA-binding proteins (in vivo), and for free RNA (in vitro) (55). Thus, it is possible to detect structural differences in long mRNAs caused by differences in environments, such as the presence of ribosomes or RNA-binding proteins in the cell (56, 57). Therefore, we performed SHAPE-MaP structure probing experiments on the β-actin mRNA present in in vitro and in vivo environments (Fig. 4).

Figure 4

Comparison of in vitro and in vivo structure for the human β-actin mRNA. (A) We calculated the Pearson correlation in windows between the SHAPE reactivities collected in vitro and in vivo for the β-actin mRNA. For each step of the trapezoid from bottom to top, the window size increases by five nucleotides from 40 to 140. High correlation (white) corresponds to areas that are similar in structure and low correlation (blue) corresponds to areas that are different in structure. The distances from the median SHAPE value for (B) in vitro and (C) in vivo β-actin were calculated in 50-nucleotide windows. Segments with reactivities above the median are less structured than segments with reactivities below the median. The gray panel highlights a region in which the SHAPE reactivity differs between in vitro and in vivo. (D) This difference is seen in the SHAPE traces for in vitro (top) and in vivo trace (bottom). Structure probing was performed using the high throughput SHAPE-MaP technique. Red nucleotides correspond to high SHAPE reactivity, yellow corresponds to medium reactivity, and black corresponds to low reactivity. The ZBP1-binding region (bright blue) and two zipcode binding protein-interacting sites (purple) are labeled above the windowed correlation and at the bottom of the SHAPE traces. The overlay for the SHAPE traces is in Fig. S5. To see this figure in color, go online.

Because we are specifically interested in differences between the two environments (in vivo and in vitro), we compute the windowed SHAPE correlation coefficient between the two data sets and plot this correlation in Fig. 4 A for a range of window sizes (40–140 nucleotides). Overall, we observe high correlation between the two data sets for a majority of the mRNA’s span, with a mean correlation coefficient of 0.88. This result can be seen clearly in Fig. S3, in which we plot raw data for a highly similar window in the coding region of the gene. We begin our structural analysis by performing SHAPE-directed Boltzmann stochastic sampling of nucleotides 200–400, which we identified as having high in vitro to in vivo correlation (Fig. S3). We expect to observe only small changes in the stochastic sampling because the SHAPE data in this region are highly similar. As expected, the visualization confirmed only small differences, but it identified a remarkably complex ensemble with 24 structural clusters (Fig. S4). This result agrees with the high median SHAPE data (Fig. 4, B and C) observed for this region; high median SHAPE is correlated with higher ensemble entropy, i.e., multiple alternative conformations (37).

The region with the lowest correlation is at the 3′ end of the mRNA. The in vitro-probed mRNA was transcribed in the absence of a polyA polymerase, therefore it was not polyadenylated, which likely explains the differences near the 3′ end because the in vivo mRNA is most likely polyadenylated (and 5′-capped). The region of difference we chose to further characterize structurally occurs 3′ of the stop codon. This region in the mRNA contains functional elements known as the Zipcode Protein Binding Protein Sites (ZPBS1 and ZPBS2). Binding of the zipcode binding protein (ZBP1) mediates mRNA localization and translation, hence the name of the protein (58, 59). We used our ensemble visualization approach to characterize the in vivo conformational rearrangements occurring in ZPBS1 and ZBPS2 within the ZBP1 binding region of the mRNA and to understand these rearrangements in the context of this region’s function. The 54-nucleotide region we model below was previously identified as necessary and sufficient for localization of β-actin mRNA to the cell periphery (58, 59). We therefore decided to specifically focus on the ensemble structure of this region.

Boltzmann stochastic sampling for the ZBP1 binding regions in vivo visualized using our approach revealed a shift in the structural ensemble away from the preferred in vitro conformation toward an alternative conformation (Fig. 5, A and B). Nonetheless, the dominant conformation in vitro (Fig. 5 A) is still significantly populated in vivo (Fig. 5 B). Thus, our visualization suggests a more complex ensemble of conformations in vivo. To further understand the structural context of the shift in ensemble, we visualized the secondary structure medoid for each of the largest structure clusters in vivo and in vitro. Although in both conformations the Zipcode Binding Protein Sites (ZBPS) are unpaired, in vivo the dominant confirmation shows ZBPS1 and ZBPS2 in a contiguous unpaired region, consistent with the larger in vivo SHAPE values. Importantly, the SHAPE reagent is not a footprinting reagent and is only minimally affected by nucleotide accessibility (60, 61). Thus, it is not surprising that we observed higher SHAPE values surrounding the ZBPS. In fact, the ZBP1 is divalent, and it has been shown to simultaneously bind the two ZBPS motifs separated by a linker portion of the RNA, although the precise occupancy of the second site is not known (37, 62). Nonetheless, binding to this region is essential for correct β-actin mRNA localization and translational control (63, 64). To accommodate the ZBP1 protein, the RNA likely has to become more open and flexible, consistent with the higher SHAPE data we observed.

Figure 5

Ensemble visualization for in vitro and in vivo human β-actin mRNA. Generation of structures for the β-actin mRNA ensemble was guided by the in vitro and in vivo SHAPE data. We compared the (A) in vitro and (B) in vivo ensembles for the region where SHAPE reactivities were expected to be different. The ensemble visualization reveals a large shift away from the dominant structure in vitro toward a second structure in vivo. We visualized the second structure for the medoid in each of the largest structure clusters. These nucleotides form different structures in vitro and in vivo. The region that differs includes the zipcode region with the two ZBP1 binding sites (purple). (C) The 3D structure for β-actin in vitro was modeled using molecular dynamics simulations without ZBP1. (D) The 3D structure for β-actin in vivo was modeled with the ZBP1 (in gray). For both 3D models, the ZBP1 binding regions are highlighted in purple. Red nucleotides correspond to high SHAPE reactivity, yellow correspond to medium reactivity, and black corresponds to low reactivity in (A–D) and (E). Comparison of SHAPE reactivity (green) and normalized RMSF (orange) for β-actin in vitro largely follow the same pattern. (F) Comparison of SHAPE reactivity and RMSF for β-actin in vivo also largely follow the same pattern. The SHAPE reactivities and RMSF values are averaged across a 3-nucleotide moving window. The RMSF is calculated from the 3D structural models. ZPB1 binding sites for (E) and (F) are boxed in purple. Fig. S6 includes further comparisons between in vitro and in vivo SHAPE reactivity and RMSF. To see this figure in color, go online.

To further understand the in vivo structural rearrangement, we performed molecular simulations of apo and bound mRNA conformations (Fig. 5, C and D). By using the secondary structure as initial constraints, we aimed to estimate the root mean square fluctuations (RMSFs) of the RNA backbone. We show these data for the apo and bound simulations in Fig. 5, E and F, overlaid with the in vitro and in vivo SHAPE data, respectively. We observed qualitative agreements between the experimental SHAPE data and the simulations, suggesting these molecular models captured overall aspects of the conformational ensemble. One important aspect of these comparisons, especially in the case of the in vivo data, is that the SHAPE data are an ensemble average over all of the β-actin mRNA molecules in the cell. Because ZBP1 binding represses translation, some message molecules are likely not bound by ZBP1, a situation that may explain why a shift to multiple conformations is observed in vivo as opposed to observing only the bound conformation. Nonetheless, these data demonstrate the value of visualizing the structural ensemble to explain structure/function relationships in an mRNA.

Discussion

RNA structure is the key component of cellular function in highly specific instances; the ribosome’s unique catalytic core is a prime example of the role of a specific RNA structure in performing protein synthesis (65, 66, 67, 68, 69). Generally, however, the functions of structures in messenger RNAs are poorly understood, except for a few cases, such as the iron responsive element (22, 70) and the histone stem loop (71, 72, 73), in which single structures are essential for function. Other than ribosomal RNA, no RNA larger than 1 kb, including mRNAs, is known to fold into a unique, well-defined conformation (65, 67, 68, 69, 74, 75). Still, although large RNAs do not adopt single conformations, specific regions do fold into complex 3D structures. One example is riboswitches in bacteria (Fig. 1). Although riboswitches are considered to be structured (i.e., they can be crystallized), riboswitches adopt multiple conformations that lead to different functions (7). Because RNAs such as riboswitches have evolved to form multiple conformations to function, it is essential to consider the suboptimal ensemble when considering structure in messenger RNAs (7).

Our approach to visualizing the suboptimal ensemble is designed to resolve some of the longstanding problems with obtaining a stable projection that allows comparisons of ensembles. A priori, this visualization approach requires sampling the entire suboptimal space to identify good principal components. For any biologically relevant RNA, such sampling rapidly becomes computationally intractable because the number of suboptimal conformations increases exponentially with length (76, 77). Thus, our approach is empirical (Fig. 2) and relies on rapid sampling of suboptimal ensembles for single point mutants of the RNA (78). Combined with multidimensional scaling and a shape-based abstraction (25, 76), our maps have the desired properties of stability and they enable comparison of different ensembles. The stability of our visualization and the complexity of unstructured regions are best illustrated in Figs. S3 and S4.

The main biological motivation for our approach is the need to visualize changes in the ensemble caused by environment. Our results on the V. vulnificus add riboswitch leverage the empirical relationship between SHAPE reactivity and the free energy of folding to recapitulate the apo and bound RNA ensembles (Fig. 3). Importantly, the goal of these visualizations is to facilitate the understanding of a complex process by approximating the specific abundance of each conformation in an ensemble. Moreover, we aim to extract biological insight from the ensemble calculation; for the V. vulnificus add riboswitch, our visualization of the ensemble model recapitulates the understanding of this system in an easily interpreted diagram. The riboswitch, a smaller fragment of a larger bacterial mRNA, is a relatively straightforward example. This is not the case for complex full-length eukaryotic mRNAs that tend to be much more highly regulated and structurally sensitive to their environments (38, 79). Whether prokaryotic or eukaryotic, it is clear that mRNAs are in integral part of cellular regulation (38).

Our analyses of a full-length human mRNA in vivo and in vitro revealed some of the complexities associated with interpreting structures in large RNAs. We observed, in both conditions, regions of high (unstructured) and low (structured) median SHAPE (56), results consistent with locally structured regions. Overall, the high similarity between in vivo and in vitro SHAPE data suggests that the mRNA is not globally affected by its environment, but, instead, specific regions are affected by endogenous molecule binding. Local structure is the case for the ZPB1 binding region in the 3′ UTR of β-actin, which we visualized using our ensemble approach (Fig. 5).

A significant result of this analysis is the median windowed SHAPE, which overall appeared higher in vivo relative to in vitro for the ZPB1-binding region. This result may seem counterintuitive, as the ZBP1 would be expected to protect the RNA from the 1M7 reagent. Although protein binding is detectable by SHAPE comparisons in vitro to in vivo (56, 80), SHAPE chemistry is not a traditional footprinting technique (44, 81, 82, 83). Thus, it is likely that the majority of differences in the SHAPE reactivity in this region are due to a conformational rearrangement due to protein binding, and not the footprint of the protein.

Our model (Fig. 5, A and B) successfully reports a shift in the ensemble, but the model does not suggest a totally dominant alternative in vivo conformation. This restriction is in contrast to the add riboswitch, in which ligand excess shifts the ensemble to almost completely the on-conformation (Fig. 3 C). It is important not to overinterpret the relative ratios of the two dominant conformations proposed for the ZBP1-binding region modeled in Fig. 5 B. However, the model is consistent with our expectation of a mixed population of ZBP1-bound and unbound β-actin mRNA. Also, the fact that the ZBP1 has two binding sites and these sites are not always simultaneously occupied (37, 84) is an additional aspect that our model cannot currently describe. Thus, our visualization accurately represents the likely state of the population of β-actin mRNAs in the cell, but still requires biological knowledge to be fully interpretable.

We performed constrained molecular dynamics simulations of the two proposed structural models of β-actin mRNA to determine if the models agreed qualitatively with the SHAPE data. Because SHAPE chemistry measures backbone flexibility (81, 85), we report RMSFs for both models in Fig. 5, E and F. For the ZBP1-binding region between ZBPS1 and ZBPS1, the agreement between the simulation and SHAPE data is better for the in vitro model compared to the in vivo simulation. For the in vivo model, we constrained both ZBPS1 and ZBPS2 to the binding pockets, which explains the low flexibility of ZBPS1 and ZBPS2. The higher SHAPE data for these two binding sites in vivo are consistent with a significant subset of mRNAs being unbound, which agrees with our ensemble model that suggested a further opening of the structure.

In summary, we have developed a computationally based visualization approach that faithfully represents ensemble mRNA populations and the effects of environment on the ensembles. The β-actin mRNA and the V. vulnificus add riboswitch are two well-characterized systems in which ensemble visualization improves the interpretation of environmentally imposed structural differences. By releasing a software package to create these visualizations easily, we encourage the RNA folding community to simulate more than just minimum free energy structures and to explore the suboptimal ensemble for all mRNAs existing in a cell. It is not clear whether suboptimal alternative conformations are a necessary component of RNA function in the cell or a by-product of the rules that govern RNA folding (28, 86, 87, 88, 89). Regardless, structure ensembles are a thermodynamic reality of RNAs and are accommodated as a feature of their function.

Author Contributions

A.L., C.T.W., L.L., D.G., B.W., and N.V.D. designed the experiments. C.T.W. designed the visualization algorithm. L.L. performed SHAPE-MaP experiments. B.W. performed molecular dynamic simulations. A.L., C.T.W., L.L., and B.W. analyzed data. A.L, C.T.W., L.L., and B.W. wrote the manuscript.