Ensemble Switching Unveils a Kinetic Rheostat Mechanism of the Eukaryotic Thiamine Pyrophosphate Riboswitch.

Thiamine pyrophosphate (TPP) riboswitches regulate thiamine metabolism by inhibiting the translation of enzymes essential to thiamine synthesis pathways upon binding to thiamine pyrophosphate in cells across all domains of life. Recent work on the Arabidopsis thaliana TPP riboswitch suggests a multi-step TPP binding process involving multiple riboswitch configurational ensembles and that Mg2+ dependence underlies the mechanism of TPP recognition and subsequent transition to the expression-inhibiting state of the aptamer domain followed by changes in the expression platform. However, details of the relationship between TPP riboswitch conformational changes and interactions with TPP and Mg2+ ¬¬in the aptamer domain constituting this mechanism are unknown. Therefore, we integrated single-molecule multiparameter fluorescence and force spectroscopy with atomistic molecular dynamics simulations and found that conformational transitions within the aptamer domain's sensor helices associated with TPP and Mg2+ ligand binding occurred between at least five different ensembles on timescales ranging from µs to ms. These dynamics are orders of magnitude faster than the 10 second-timescale folding kinetics associated with expression-state switching in the switch sequence. Together, our results show that a TPP and Mg2+ dependent mechanism determines dynamic configurational state ensemble switching of the aptamer domain's sensor helices that regulates the stability of the switch helix, which ultimately may lead to the expression-inhibiting state of the riboswitch. Additionally, we propose that two pathways exist for ligand recognition and that this mechanism underlies a kinetic rheostat-like behavior of the Arabidopsis thaliana TPP riboswitch. Published by Cold Spring Harbor Laboratory Press for the RNA Society.

INTRODUCTION

Riboswitches, first reported by Henkin (Grundy and Henkin 1998) and later named by Breaker and coworkers (Barrick and Breaker 2007), are mRNA elements that bind metabolites to regulate neighboring transcription, translation, or splice machinery (Mironov et al. 2002; Nahvi et al. 2002; Winkler et al. 2002a,b; McDaniel et al. 2003; Sudarsan et al. 2003; Garst et al. 2011). Typically, riboswitches are composed of two distinct functional domains, the ligand-sensing aptamer domain and the transcription, translation, or splice regulating expression platform (Fig. 1A). An effector molecule, or ligand, binds to the aptamer domain causing structural changes in the expression platform domain to regulate gene expression. In the ligand-free state, aptamer domains can adopt multiple three-dimensional conformations, allowing for target recognition and stabilization of the ligand in the bound state (Winkler and Breaker 2003). Further, riboswitches require high selectivity in target ligand recognition to elicit the appropriate regulatory response.

FIGURE 1.

Model of TPP binding to the aptamer domain of the Arabidopsis thaliana TPP riboswitch. (A) Linear schematic of the TPP riboswitch showing the expression platform (black) and the aptamer domain, highlighting the locations of the donor (green) and acceptor (red) fluorophores in the pyrophosphate (purple) and pyrimidine (orange) sensor helices, respectively, and the P1 switch helix (gray). (B) Primary and secondary structure schematic of the aptamer domain of TPP riboswitch used in all experiments and simulations. Five Watson–Crick base-paired (black) helices (P1–P5) connect through three bulges (J2/3, J2/4, and J4/5) and two loops (L3 and L5). The two sensor helix arms of the aptamer domain (P4/P5, purple and P2/P3, orange) sense the pyrophosphate (pluses) and pyrimidine (ring structure) ends of the TPP ligand (yellow), respectively, and a coordinating Mg2+ ion (blue circle). The P1 switch helix emerges from the sensor helices to form a three-way junction (J2/4), and the expression platform (not shown) extends beyond the 3′ end of the switch helix. We truncated the sequence of L3 to match the crystal structure and added donor (Alexa 488 in position 56, green) and acceptor (Cy-5 in position 25, red) fluorophores (See Materials and Methods for details on riboswitch design; Thore et al. 2006) to L5 and near L3, respectively. (C) Schematic of TPP riboswitch aptamer domain “ON” (top) and “OFF” (bottom) states. Current models of the TPP riboswitch mechanism suggest that the aptamer domain is in the “ON state” (switch helix open, top) when the sensor helix is open and the fluorescent markers (green circle for the donor, red circle for the acceptor, see Materials and Methods) are far apart (low FRET), and it is in the “OFF state” (switch helix closed, bottom) when the sensor helix is closed, and fluorescent markers are close together (high FRET), upon TPP binding between J4/5 and J2/3. Figure partially based on Thore et al. (2006) and Serganov and Nudler (2013).

The thiamine pyrophosphate (TPP) riboswitch, also known as Thi-box riboswitch, regulates the transport and synthesis of thiamine and its phosphorylated derivatives (Manzetti et al. 2014) in response to changes in TPP concentration. TPP is the active form of vitamin B1, an essential coenzyme in many archaea, bacteria, and eukaryotes (Rodionov et al. 2002; Kubodera et al. 2003; Sudarsan et al. 2003; Bocobza et al. 2007). Moreover, the TPP riboswitch has been postulated as an important therapeutic target for novel antibacterial drugs (Blount and Breaker 2006). For example, pyrithiamine pyrophosphate (PTPP, an antimicrobial agent) down-regulates essential thiamine production in Bacillus subtilis and Aspergillus oryzae, demonstrating the potential use of riboswitches as targets for antibacterial or antifungal drug design (Sudarsan et al. 2005).

Similar structures and regulatory responses to TPP have been found for the TPP riboswitch in multiple organisms, including E. coli ThiM riboswitch, A. oryzae ThiA riboswitch, and A. thaliana ThiC riboswitch (Sudarsan et al. 2003; Serganov et al. 2006; Thore et al. 2006). The TPP riboswitch aptamer domain consists of five base-paired helices (Fig. 1B, P1 to P5; Thore et al. 2006) and comprises two sensor helix arms, each composed of two stacked helices connected through a bulge, and a switch helix (P1). The stacking of helices P2 and P3 forms the pyrimidine sensor helix, and the stacking of helices P4 and P5 forms the pyrophosphate sensor helix domain (Ali et al. 2010). The TPP riboswitch's expression platform and aptamer domain comprise a secondary structural switch that regulates the expression machinery. The exact mechanism of switching varies for TPP riboswitches from different organisms. However, they all bind TPP and regulate its production (e.g., modulates ribosome binding to the expression platform in E. coli [Winkler et al. 2002a] or controls mRNA processing through alternative splicing in A. thaliana [Serganov and Nudler 2013]). This expression platform connects to the aptamer domain at the 3′ end of the P1 switch helix, which extends beyond at the junction of the sensor arms (aptamer domain pictured in Fig. 1A,B).

The long-standing model of the TPP riboswitch's mechanism involves two mutually exclusive structural conformations in the aptamer domain (Fig. 1B,C) corresponding to the “ON” and “OFF” states of the riboswitch's expression platform that enable and inhibit the expression of enzymes essential to thiamine synthesis pathways, respectively. The model suggests that the expression platform's OFF state corresponds to the aptamer domain's TPP-bound state, while the ON state corresponds to the ligand-free state (Winkler et al. 2002a). Further, the ON/OFF behavior in E. coli ThiM riboswitch (Schwalbe et al. 2007; Haller et al. 2013) is coupled to the formation of the switch helix, P1, with the formed helix corresponding to the OFF state and unpairing of P1 corresponding to the ON state (Thore et al. 2006). However, implicit RNA dynamics suggest that the aptamer domain may sample multiple ensembles of configurations rather than a small number of discrete conformational states. In this manuscript, the term “configurations” refers to subensembles of rapidly sampled structural states with similar characteristics that comprise a “conformation.” Both binary ON/OFF and ensemble-like tuning behavior have been observed for various riboswitches, in some cases with differing behavior for the same riboswitch in different organisms (Baird et al. 2010). Additionally, models with several discrete states have been proposed to explain riboswitch functions, such as the coupling of ligand and temperature sensing of the adenine-sensing riboswitch (Reining et al. 2013).

Previous work on TPP riboswitches has focused on the TPP ligand recognition process. Typically, TPP binds to TPP riboswitches with high affinity (up to 20 nM) (Winkler et al. 2002a; Yamauchi et al. 2005) because the stacking of the phosphate sensor helix over the pyrimidine sensor helix allows recognition of the pyrimidine ring and pyrophosphate group of TPP (Edwards and Ferre-D'Amare 2006; Serganov et al. 2006; Thore et al. 2006). However, the situation is complicated because an Mg2+ ion coordinates TPP binding (Baird and Ferre-D'Amare 2010; Burnouf et al. 2012). A further complication arises because a dynamic equilibrium exists between multiple aptamer domain conformational states, even when bound to TPP (Ali et al. 2010; Baird and Ferre-D'Amare 2010; Baird et al. 2010). Additionally, the sensor helix formation at the junction with the expression platform is necessary to further stabilize TPP binding to the aptamer domain and transition to the OFF state in the E. coli TPP riboswitch (Haller et al. 2013). Moreover, recent experiments suggest a two-step TPP binding process, going from a weakly bound intermediate conformation to the tightly bound final state (Anthony et al. 2012; Savinov et al. 2014).

In the TPP riboswitch, it is thought that the two sensor helices close upon recognition of the pyrimidine and pyrophosphate groups of TPP (Thore et al. 2006). However, the conformation of the aptamer domain upon target recognition before it closes fully remains unclear. It has been suggested that J4/5 kinks to facilitate the Mg2+-coordinated binding of the TPP's phosphates, while the J2/3 bulge quasi-simultaneously kinks to accommodate the pyrimidine base of TPP on the other domain (Thore et al. 2006). However, this binding mechanism remains somewhat speculative because these configurational transitions are thought to be rather subtle and fast (submillisecond); thus, they are challenging to characterize. Nonetheless, there are some reports of slow (second timescales) dynamic transitions, which would require that the transitions have high energetic costs (Haller et al. 2013). Similar riboswitches exhibit slow dynamics, such as the cotranscriptional folding and dynamics of the E. coli thiM riboswitch, which has 68% identity to the Arabidopsis thaliana TPP-riboswitch, that were revealed by single-molecule FRET (smFRET) (Uhm et al. 2018). smFRET has also been used alongside nuclear magnetic resonance (NMR) to reveal complex, ligand-dependent conformational dynamics in the adenine-sensing riboswitch (Warhaut et al. 2017).

To resolve the apparent mismatch between fast and slow mechanisms, as well as to probe the ON and OFF conformations and potential intermediate states, we integrated single-molecule multiparameter fluorescence spectroscopy (Hamilton and Sanabria 2019), including time-correlated single-photon counting (TCSPC) (Lakowicz 2007), calculation of fluorescence parameters for single-molecule events (burst-integrated) (Eggeling et al. 2001), and fluorescence fluctuation methods using a FRET-labeled Arabidopsis thaliana TPP riboswitch aptamer domain (Felekyan et al. 2012, 2013), with all-atom discrete molecular dynamic (DMD) simulations (Ding et al. 2008b) and optical tweezer force spectroscopy measurements. We found that the Arabidopsis thaliana TPP riboswitch's aptamer domain exhibits Mg2+ and TPP-dependent ensemble switching that give rise to two distinct binding pathways in the transition to the known OFF state with both Mg2+ and TPP bound. Further, we only observe the TPP-bound high-FRET conformation at extreme conditions of saturating Mg2+ and TPP, while in moderate and physiological conditions, OFF-like states and other intermediates are transiently populated at submillisecond timescales. Together, our results show that the TPP riboswitch's aptamer domain is a conformational switch capable of tuning function through ligand-dependent dynamics as a kinetic rheostat rather than a binary, all-or-nothing switch.

RESULTS

Submillisecond dynamics of Arabidopsis thaliana TPP riboswitch's aptamer domain

We placed the FRET-labeled aptamer domain of the TPP riboswitch in various buffer conditions, and we observed changes in the shape and position of the two-dimensional MFS histograms (Fig. 2). In the apo buffer (absence of Mg2+ and TPP), we observed a narrow, unimodal distribution centered in the low-FRET regime (Fig. 2A). Even at extreme saturated Mg2+ concentrations (1 M, Fig. 2B), the population was unimodal and centered in the low-FRET regime, with a slight decrease in the mean and increase in the skew of the fluorescence-weighted lifetime population as compared to the apo buffer results according to MFS alone. The population still lies along the static FRET line.

FIGURE 2.

The dynamical nature of the sensor helices depends on the presence of TPP and Mg2+. Multidimensional smFRET populations for the aptamer domain in (A) apo, (B) Mg2+, (C) TPP, and (D) Mg2+ and TPP buffers (see Materials and Methods for buffer compositions; total number of single-molecule events [bursts] are N = 94,999, N = 5,198, N = 8418, and N = 21,836, respectively). The 1D histograms indicate the number of single-molecule bursts corresponding to binned values of the FRET indicators, while the normalized 2D contours relate the two parameters to each other. The FRET indicators shown are the ratio of donor to acceptor fluorescence (F) and the average donor fluorescence lifetime 〈τ〉 for each PIE-selected 1:1 donor-to-acceptor stoichiometry burst (see Materials and Methods). Darker contours correspond to a higher density of bursts, and static FRET lines (black line, Supplemental Information Table S2) describe the expected relationship between the FRET indicators for nonexchanging, static FRET populations according to Förster theory. Large F and large 〈τ〉 correspond to the low FRET regime (long interdye distance), and small F and 〈τ〉 correspond to the high FRET regime (short interdye distance). Note that two static FRET lines are needed in Mg2+ and TPP conditions (D) due to acceptor quenching (see Supplemental Information).

In the TPP buffer, we found a broader distribution of the FRET measurement population that tailed toward lower F values (Fig. 2C), which indicates the presence of more closed, higher-FRET states that we did not observe in either the apo or even at extreme Mg2+ concentrations (Fig. 2A,B). Additionally, in the TPP buffer, we found that the highest density of FRET measurements was shifted slightly to the right of the static FRET line (Fig. 2C), which indicates dynamic averaging of FRET states within the ms measurement timescale (Sisamakis et al. 2010; Hamilton and Sanabria 2019; Sanabria et al. 2020) that was not observed in the apo or Mg2+ buffers (Fig. 2A,B). Together, these results suggest that, in the presence of TPP alone, the sensor helices of riboswitch's aptamer domain exhibit conformational dynamics between open, low-FRET states and more closed, higher-FRET states at timescales faster than the ms observation time.

In the Mg2+ and TPP buffer, we found two distinct FRET populations, which was particularly evident in the one-dimensional projection of the 〈τ(〉 distribution (Fig. 2D). One of these populations was broad (∼3.5 nsec FWHM) and centered around 〈τ〉 = 3.1 nsec (medium-FRET), and the other was narrow (∼0.3 nsec FWHM) and centered below 〈τ〉 < 0.5 nsec (high-FRET) (Fig. 2D). This splitting of the FRET populations (Fig. 2D), as compared to the results in other buffers (Fig. 2A–C), suggests that the exchange processes between open and closed conformations of the sensor helices occur at timescales similar to or longer than the msec timescale (Hamilton and Sanabria 2019; Sanabria et al. 2020) in either saturating Mg2+ or TPP buffer conditions. This result suggests only the presence of both saturating TPP and extreme Mg2+ concentrations slows transitions away from the high-FRET closed state of the sensor helices to a timescale greater than the msec. In any other condition (Fig. 2A–C), there is not a measurable long-lived (>msec) population of closed aptamer domain sensor helices. Additionally, the broadening of the lower-FRET population in TPP and the shift of its peak to medium-FRET upon the addition of Mg2+ suggests that Mg2+ may enable the sensor helices to sample other, intermediate states along the closing pathway.

Together, these results suggest that the Arabidopsis thaliana TPP riboswitch's aptamer domain sensor helices are in an open state (with L3 and L5 loops far apart, Fig. 2A,B) in the absence of TPP, can transiently be found in a more closed state with TPP alone (with L3 and L5 loops closer together, Fig. 2C), and populate a long-lived (>msec timescale) closed state (predicted to be a conformation with the L3 and L5 loops close together like the Mg2+ and TPP bound structural model) only in the presence of saturating TPP and Mg2+ concentrations (Fig. 2D; Thore et al. 2006). Our results also indicate that the aptamer domain's sensor helices respond to its two ligands, Mg2+ and TPP, differently (Fig. 2B,C), which suggests two independent closing pathways. These results do not conclusively determine whether the sensor helices are dynamic, in the apo or Mg2+ alone conditions (Fig. 2A,B), as opposed to the TPP alone condition, in which we observed a shift away from the static FRET line (Fig. 2C). Nonetheless, the observed differences in the shape of the FRET populations (Fig. 2) suggest a dynamic sampling of conformational states of the sensor helices in all measured conditions (Sisamakis et al. 2010; Hamilton and Sanabria 2019). To probe the dynamics of the TPP riboswitch's aptamer domain sensor helices further, we use TCSPC and fFCS, detailed in the following sections.

The stabilizing effect of extreme Mg2+ concentrations on the aptamer domain

To further address the role of Mg2+ ions ability to coordinate conformational changes in the aptamer domain sensor helices, we measured the donor fluorescence lifetime at increasing concentrations of Mg2+ (0–0.5 M, Fig. 3A) in the absence of TPP. We found that the mean donor fluorescence lifetime (〈τ〉) population distributions were nearly invariant over this broad range of Mg2+ concentrations (Fig. 3A). In particular, we found that the expected values and the standard deviations of the mean donor fluorescence lifetime population distributions were not significantly different over this range of Mg2+ concentrations (Fig. 3D, left). This result suggests that the sensor helices stay open in the absence of TPP, even for increasing Mg2+ over a large, nonphysiological range of concentrations.

FIGURE 3.

Increasing concentrations of Mg2+ in the presence of TPP induce a population of closed sensor helix states in the aptamer domain. Average fluorescence lifetime distributions of the donor in the presence of acceptor (〈τ〉) for FRET-labeled TPP riboswitch aptamer domains (A) with increasing concentrations of MgCl2 starting from 0 mM to 1 M, (B) with increasing concentrations of TPP starting from 0.3 mM to 4.8 mM, and (C) with increasing concentrations of MgCl2 starting from 5 mM to 1 M, all at a fixed TPP concentration of 4.8 mM. (D) Mean donor fluorescence lifetimes calculated from the data in A–C. The standard deviations of the distributions are shown as error bars.

We measured the mean donor fluorescence lifetime of the FRET-labeled TPP riboswitch aptamer domain at increasing concentrations of TPP (0–4.8 mM, Fig. 3B), but in the absence of Mg2+, and found that the donor fluorescence lifetime population broadened with increasing TPP concentration. In particular, we found that the expected values of the mean fluorescence lifetimes were relatively constant while the standard deviations of the mean donor fluorescence lifetime population distributions increased concomitantly with TPP concentration (Fig. 3D, middle). This broadening was consistent with our previously noted observation of a broader mean donor fluorescence lifetime in the TPP buffer as compared to the apo buffer (Fig. 2A,C, respectively).

Because we had found that the FRET population deviated from the static FRET line in the TPP buffer (Fig. 2C), which is consistent with dynamic averaging effects and indicates that the riboswitch's sensor helices are unable to be stabilized in the closed state in the absence of Mg2+, we investigated the effect of Mg2+ concentration on the ability of the aptamer domain to be in the high-FRET, closed state by measuring the donor fluorescence lifetime of the FRET-labeled riboswitch aptamer domain at increasing concentrations of Mg2+ in TPP-saturating (4.8 mM) conditions. We found that at a 5 mM concentration of Mg2+, which already is notably high compared to the expected physiological Mg2+ concentration of about 1 mM in Arabidopsis thaliana (Bose et al. 2013), the mean donor fluorescence lifetime distribution of the TPP riboswitch's aptamer domain (Fig. 3C, top) is not significantly different from the case of saturating TPP and no Mg2+ (Fig. 3B, bottom). When the concentration of Mg2+ is increased to 25 mM, we found a broadening of the mean donor fluorescence lifetime distribution due to an increase in the frequency of observations at the higher-FRET tail (Fig. 3C). It is only as the Mg2+ concentration is further increased to extreme limits (to 100, 500, and 1000 mM, Fig. 3C), that we find a significant high-FRET population, indicated by the population peak at lower τ values (Fig. 3C), corresponding lower expected values (Fig. 3D, right), and an increase in the standard deviation of the total population (Fig. 3D, right). Together, these results suggest that Mg2+ stabilizes a closed state of the sensor helices only at nonphysiological concentrations of Mg2+, which we expect is most likely due to the screening effect of Mg2+ against the repulsive negative charges on the RNA backbone and the electronegative polarity of TPP. In total, these results strongly suggest that the TPP riboswitch aptamer domain's sensor helices exist in a highly dynamic state that transiently samples closed-like configurations at physiological concentrations of magnesium.

There are at least five distinct ensembles of the aptamer domain

To identify configurational ensembles of the aptamer domain that persist longer than the duration of the fluorescence lifetime of the FRET-labeled riboswitch aptamer domain (ns timescales), we compared the fluorescence decay histograms of donor-only reference aptamer domains, which we generated using time-correlated single-photon counting (TCSPC, Materials and Methods), with the decay histogram of the FRET-labeled riboswitch aptamer domains in the same buffer conditions. We found that there was a significant population of configurational ensembles with no-FRET (gray bars) compared to the corresponding fraction with donor fluorescence quenched through FRET (colored bars) in each buffer condition (Fig. 4A). The fraction of molecules exhibiting FRET increased from 27% in the apo buffer to 41% in saturating Mg2+ and TPP conditions (Fig. 4A). This result is consistent with the multidimensional smFRET data, in which we observed a population shift toward high-FRET (Fig. 2) and a broadening of the 〈τ〉 distribution (Fig. 3) in Mg2+ and TPP buffer as compared to apo buffer.

FIGURE 4.

Identification of at least five sensor helix conformational states as a function of ligand concentration. (A) Bar plot represents the fraction of the population of molecules exhibiting FRET (color-filled bars) and no-FRET (gray bars) as determined by fluorescence decay from TCSPC FRET experiments. (B) Interdye distance distribution with a two-Gaussian distribution fitting model. a is a semi-open/open ensemble, b is an apparent intermediate ensemble, c is a second intermediate ensemble that is more compact than b, and d is the closed ensemble of the aptamer domain's sensor helices. Each sample is fitted with distances corresponding to two states plus a no-FRET/long distance state as well as the intrinsic variance (width) for each state (Table 1; Supplemental Information Table S5A,B). Widths of the bars represent the uncertainties in the mean distances for each Gaussian-distributed state.

TABLE 1.

Interdye distances associated with the conformational states of the aptamer domain's sensor helices identified in the population of FRET-exhibiting TPP-riboswitches

We calculated the interdye distances of the FRET-exhibiting population of TPP riboswitch aptamer domains in each buffer using the functional model (Supplemental Information, Equation 4). We identified a total of four distinct TPP riboswitch aptamer domain sensor helix ensembles (a, b, c, and d, Fig. 4B), with at least two at each buffer condition undergoing dynamic interconversion between those ensembles. In each of the buffers that contain ligand (TPP, Mg2+, or both) the more populated FRET-exhibiting state was the more closed state [smaller mean interdye distance, 〈R(1)〉, Fig. 4B], while in the apo buffer the more open state (larger mean interdye distance, 〈R(2)〉, Fig. 4B) was more populated. Further, in the ligand-containing buffers, the no-FRET population was smaller compared to the apo buffer. In the apo and Mg2+ buffers, we found that the FRET-exhibiting states 〈RDA〉s were similar to one another (states a and c, Fig. 4B and Table 1), but the relative population fractions switched, with the more closed state c being more populated in Mg2+ and the more open state a being more populated in apo buffer (Fig. 4B). The change in relative state populations suggests a change in the exchange dynamics between these ensembles but does not reveal whether this is due to faster transitions from a to c or slower transitions from c to a, as the MFS histograms appear unimodal for both conditions. Thus, the difference is likely due to subtle changes in the submillisecond regime. In both buffers containing saturating concentrations of TPP (TPP and Mg2+ and TPP buffers), where the fraction of the population exhibiting FRET was 41% (Fig. 4A), we found that 78% of the FRET-exhibiting ensemble had a mean interdye distance of 〈〉 = 25.4 ± 1.0 Å (state d, Fig. 4B). The other 22% of the FRET-exhibiting ensemble had 〈〉 = 43.4 ± 1.2 Å in TPP alone (state b, Fig. 4B) and 〈〉 = 62.1 ± 4.5 Å (state a, Fig. 4B) in the Mg2+ and TPP buffer. We also found that the most-populated FRET-exhibiting ensemble observed in both TPP-containing buffers (〈〉 = 25.4 ± 1.0 Å, state d, Fig. 4B) is more closely related to the most-populated FRET-exhibiting ensemble observed in the Mg2+ and apo buffers (interdye distance of 31.7 ± 1.9 Å, state c, Fig. 4B) than the less-populated ensembles in TPP-containing buffers (states a and b). Finally, we found that the open state in the Mg2+ and TPP buffer was more like state a found in the apo and Mg2+ buffers rather than state b found in the TPP buffer (Fig. 4B). A summary of the interdye distances and their statistical uncertainties are shown in Table 1. We found that ensembles d and c (Fig. 4B) agreed within 5 Å of the accessible volume (AV) modeling of the dyes at the corresponding locations in a structural model with the sensor helix closed state, which is thought to be associated with the TPP riboswitch's OFF state (26.9 Å, Supplemental Information Fig. S7; Kalinin et al. 2012).

The no-FRET (Fig. 4A, gray bars) population has an interdye distance >80 Å. These are likely long interdye distance populations with fluorescent labels that are too far apart to be quantified by FRET due to low acceptor sensitization, rather than mislabeled or damaged constructs, because the synthesis of the FRET-labeled riboswitch aptamer domains ensures 1:1 donor-to-acceptor stoichiometry (Materials and Methods) and we only analyze bursts with S values near 0.5 (Supplemental Information Fig. S3). These data suggest that sensor helices are wide open a significant fraction of the time in all buffers, but more so in those without TPP (Fig. 4A). These no-FRET populations likely include aptamer domains that are not entirely folded and a small number that might contain acceptor fluorophores incapable of FRET.

In total, our results strongly suggest that there are at least five distinct ensembles, but only three are resolved in each buffer condition, with one state for which the interdye distance is too long to be quantified with the chosen FRET pair and assigned labeling locations. A sampling of multiple aptamer domain sensor helix configurational ensembles in each buffer suggests both a dynamic aptamer domain closing pathway and a rheostat-like tuning of the aptamer domain subensemble configurational states that might ultimately control TPP riboswitch ON/OFF switching.

The aptamer domain ensemble switching rates span four orders of magnitude

We identified at least five distinct aptamer domain sensor helix configurational ensembles (Table 1; Fig. 4), but the multidimensional smFRET histograms displayed mostly unimodal distributions with varying widths (Figs. 2, 3). To obtain a holistic solution, we postulate that the unimodal distribution in Figures 2, 3, must reflect dynamic averaging occurring between 25 nsec (TCSPC observation window) and the msec-timescale (the diffusion time of the aptamer domain). To test this hypothesis, we used filtered fluorescence correlation spectroscopy (fFCS), which, as opposed to traditional fluorescence spectroscopy methods, uses characteristic fluorescence decays and time-resolved anisotropy decays to probe state-specific exchange processes (Felekyan et al. 2012, 2013). When only two states undergo kinetic exchange, the cross-correlation function between species corresponding to those states shows a decrease in the amplitude of the correlation at the characteristic timescale of exchange between those states. Additional decay terms appear for additional transitions in more complex kinetic schemes. Neither species-specific auto-correlation function exhibits correlation amplitude decreases at these characteristic anticorrelation times (Felekyan et al. 2012, 2013). Thus, we use fFCS cross-correlation curves to probe the conformational transitions in the aptamer domain's sensor helices that are only reported by FRET changes.

We performed a global fit of the species-specific auto- and cross-correlation data (Felekyan et al. 2013) that simultaneously satisfied all the relaxation terms, all the terms corresponding to photophysical effects, and all the diffusion terms in the correlation functions (Materials and Methods, Fig. 5) using the same photon-stream data that we used in the TCSPC experiments (Fig. 4) for each buffer condition (Supplemental Information Fig. S8). We found relaxation times ranging from 100 nsec to msec, with one observed at every decade in time (e.g., apo buffer, Fig. 5A). We calculated the effective transition rate constants associated with each buffer condition by taking the inverse of the relaxation times (Fig. 5A,B for the apo buffer, Supplemental Information Fig. S8 for the other conditions and Supplemental Information Fig. S9 for filter component probabilities) and plotted them for each measured condition (Fig. 5C). These transition rates indicate switching between at least five (number of relaxation times plus 1) different conformations at timescales ranging from µsec to msec, consistent with the hypothesis that the lifetimes derived by time-resolved analysis correspond to heterogeneous ensembles in all buffer conditions.

FIGURE 5.

The sensor helices exhibit four dynamic timescales. Example filtered FCS species auto-correlation (A, sACF) and cross-correlation (B, sCCF) functions from FRET experiments on the aptamer domain in the apo buffer. Timescales of transitions between low-FRET (LF) and high-FRET (HF) states and vice-versa appear as anticorrelation terms in the sCCF. Four distinct state transitions rates (vertical dashed lines), spanning different decades in time were identified in all conditions, as exemplified here by data obtained in the apo buffer. Shaded regions correspond to timescales typical of chain dynamics, local configurational changes, and global dynamics (Bothe et al. 2011; Mustoe et al. 2014), from darker to lighter, respectively. Functional fits are shown as solid black curves, while colored lines represent the raw correlation data. (C) Four relaxation time populations were observed in each set of conditions. Bar height represents the log of the inverse of each correlation time. fFCS fit parameters are listed in Supplemental Information Table S5C.

We compared the state-switching transition rates in each buffer condition. We found slightly slower effective transition rates at all timescales in the Mg2+ buffer as compared to the apo buffer condition (Fig. 5C), suggesting that the population shift toward shorter lifetimes, which were observed using the donor fluorescence decay data (Fig. 3), requires Mg2+. This Mg2+-dependent mechanism appears to slow transitions from state c to state a in saturating Mg2+ conditions, lowering the overall transition rates. Further, increases in the effective transition rates at all timescales upon the addition of TPP (Fig. 5C) indicate that a saturating concentration of TPP lowers the energy barriers between the sampled subensemble configurations. Additional increases in the fast, µsec-timescale transition rates, corresponding to tR(1) and tR(2), in the Mg2+ and TPP buffer compared to only TPP (Fig. 5C) indicates increased subensemble configuration sampling that may correlate with a transition to the closed aptamer state. Little change in the slowest transition rate, corresponding to tR(4), in the Mg2+ and TPP buffer compared to TPP buffer (Fig. 5C) suggests that either large-scale transitions are not significantly altered or additional rates are not captured because they are longer than the diffusion time of the aptamer domain molecules through the confocal volume (compatible with the long-timescale stabilization of a high-FRET state as observed by MFD, through the additional appearance of a low-donor-lifetime peak in Mg2+ and TPP).

Ensemble switching observed with FRET is predicted by replica-exchange DMD simulations

To provide structural insights into the experimentally observed conformational dynamics, we used replica-exchange discrete molecular dynamics (REXDMD) simulations of the TPP riboswitch's aptamer domain under different conditions, including in the absence and presence of TPP and Mg2+ ions (see Materials and Methods). We observed dynamic, spontaneous conformational exchange in each case (Supplemental Information Figs. S9–S12), highlighting the sampling efficiency of the all-atom REXDMD simulations. We computed the probability distribution function (PDF) of the intersensor helix arm distance between dye-conjugated bases G25 and U56 at room temperature (Fig. 6A). In the absence of the TPP and Mg2+ coligands, we found that the sensor helices stayed in open-like state configurations with the highest probability centered at an inter-arm distance of 94 Å and a broad distribution tailing toward shorter inter-arm distances (Fig. 6A). In the presence of Mg2+ at a 38:1 molecular ratio to RNA (63 mM simulation concentration), the sensor helices exhibited a similar trend as in the apo form, and the PDF displayed a shoulder at the inter-arm distance of 85 Å in addition to the peak centered at 93 Å (Fig. 6A). In the presence of TPP only (1:1 molecular ratio to RNA, or 1.6 mM), we found an intersensor helix arm distance PDF with a broad unimodal distribution centered at a distance of 73 Å (Fig. 6A). In the presence of both Mg2+ and TPP (63 and 1.6 mM, respectively), the PDF exhibits two distinct peaks, with the inter-arm distance centered at 25 Å representing a closed-state ensemble and the population with peaks at 85 and 95 Å resembling the open-state ensemble, similar to the apo and Mg2+ conditions (Fig. 6A). These simulations qualitatively agree with the FRET measurements (Fig. 4), which showed the appearance of low-populated closed-like ensemble (state c) in the Mg2+ buffer, a shift toward a closed-like ensemble (state d) and the appearance of intermediate conformations (state b) in the TPP buffer, and a prominent population peak for a closed state (state d) and a simultaneous reduction in the open-like ensemble (state a) in the Mg2+ and TPP buffer.

FIGURE 6.

Atomistic DMD simulations of the aptamer domain find sensor helix conformational states that correlate to the FRET data. (A) The probability distribution function of the sensor helix inter-arm distance measured between G25 and U56. (B–E) The two-dimensional PMFs as functions of inter-arm (G25 and U56) and P1 switch helix/P2 pyrimidine sensor helix costacked distances, measured between U39 and G73, in each simulated buffer condition. The basins correspond to a sensor arm open state with no P1/P2 costacking (α, α′, α″), sensor arm open state with P1/P2 costacking (β, β′), partially closed sensor arm conformational state with P1/P2 costacking (γ, γ′), and sensor arm closed state with P1/P2 costacking (δ).

The simulation trajectories and RNA conformations suggested that a costacking between the P1 switch helix and P2 (Supplemental Information Fig. S10) in the pyrimidine sensor helix was coupled to the opening and closing dynamics of the sensor arms. Hence, we computed the PMF (see Materials and Methods) as a function of both the inter-arm distance, D (measured between the oxygen atoms of the 2′-hydroxyl groups of the labeled G25 and U56 and representing the aptamer domain's sensor helix conformation), and the costacking distance, D (measured between the C1 atoms of the sugar groups in U39 and G73 and representing the P1/P2 costacking configuration), in each buffer condition (apo, Mg2+, TPP, and Mg2+, and TPP, Fig. 6B).

The aptamer domain's sensor helices were mostly open when simulated in the apo buffer (Supplemental Information Fig. S11 for simulation trajectories), with a major basin centered at D = 88 Å, but they did exhibit both closed-like and open-like conformations (Fig. 6B). The major basin was found at D = 8 Å (labeled β, sensor helices in an open state with P1/P2 costacking, Fig. 6B) and had a shoulder corresponding to the loss of the P1/P2 costacking interaction with D = 17 Å (labeled α, sensor helices in an open state without P1/P2 costacking, Fig. 6B). Additionally, there was a weakly populated basin at D = 54 Å and D = 8 Å (labeled γ, sensor helices in a partially closed-like ensemble with P1/P2 costacking, Fig. 6B). Representative snapshots corresponding to the three intermediate basins are shown in Figure 7.

FIGURE 7.

Representative snapshots of the aptamer domain in the apo (α, β, γ), Mg2+ (α′, β, γ′), TPP (α″,β′, γ), and Mg2+ and TPP (α′, β, γ, δ) conditions. The states correspond to labeled basins in Figure 6B. The RNA is shown in cartoon representation, Mg2+ ions as purple spheres, and nucleotide pairs (G25, U56) and (U39, G73) highlighted in the ball-and-stick representation.

Simulations in the Mg2+ buffer (Supplemental Information Fig. S12 for trajectories) showed three basins corresponding to the aptamer domain's configurational ensembles: The predominant basins α′ and β separated by a weak energy barrier and an isolated γ′ basin. Compared to the α basin in apo condition, α had a smaller P1/P2 costacking distance D = 11 Å in the Mg2+ buffer and, compared to the γ basin in the apo condition, γ′ had a smaller inter-arm distance D = 29 Å in the Mg2+ buffer (Fig. 6B and snapshots in Fig. 7).

Compared to the apo and Mg2+ buffers, simulations of sensor helices in TPP alone (Supplemental Information Fig. S13 for trajectories) resulted in a broader and shorter unimodal inter-arm distance distribution centered around 69 Å and a slightly shorter P1/P2 costacking distance D = 7 Å (β′, Fig. 6B, snapshots in Fig. 7), but a very weakly populated, partially closed sensor helix conformational state strongly resembling the γ ensemble in the apo buffer. We also observed a weakly populated ensemble like the open sensor helix conformational state without P1/P2 costacking (α basin, Fig. 6B, snapshots in Fig. 7).

In the Mg2+ and TPP buffer simulations (Supplemental Information Fig. S14 for trajectories), the aptamer domain was mostly in a closed-like ensemble with sensor helix inter-arm distance D = 26 Å and P1/P2 costacking distance D = 7 Å (δ-basin, Fig. 6B and a corresponding conformation in Fig. 7), but also showed weakly populated ensembles corresponding to α′, β, and γ basins. Hence, both TPP and Mg2+ shifted aptamer domain conformations toward compact states with shorter sensor helix inter-arm and P1/P2 costacking distances, but the extent of this compaction and resulting intermediate ensembles were distinct. These findings are consistent with the FRET observations in Figure 4.

Further, we calculated the 2D PMF with respect to the U39–G73 costacking distance and the distance between the C1′ atoms of C6 and G73 (denoted as C6.C1′ and G73.C1′), the base pair within the switch sequence that are closest to the sensor helices and costacking site between P1/P2 (Supplemental Information Fig. S15A–D). In all conditions, the lowest energy PMF region corresponded to a C6.C1′–G73.C1′ pair distance of 10.9 Å with the two bases coplanarly oriented (Supplemental Information Fig. S15E). However, in apo and TPP conditions, the aptamer domain exhibits a distinct shift in the population toward longer U39–G73 costacking distances and loss of coplanar alignment, or buckling deformation, of the C6 and G73 base pair (Supplemental Information Fig. S15F). Therefore, our simulations suggest that the loss of P1/P2 costaking results in the deformation and subsequent destabilization of the P1 switch helix. The decrease in the buckling behavior of C6–G73 in Mg2+ and TPP condition suggests a stabilizing effect of Mg2+ and TPP for base-pairing in the P1 switch helix.

Sampling conformational ensembles in the unfolding pathway

To confirm that the rapid dynamics observed with smFRET do not correspond to unfolding transitions of RNA helices (P1–P5), we proceeded to measure the unfolding of the aptamer domain of the TPP riboswitch under force (Materials and Methods) in the same buffer conditions as in the smFRET experiments. Using a single-bead optical tweezer assay (Materials and Methods; Supplemental Information Fig. S16), we stretched and unfolded the TPP riboswitch aptamer domain. We considered the aptamer domain to be fully unfolded once the full extension of the aptamer domain was reached, at which point the stretching behavior matches that of a worm-like chain polymer and the complete loss of secondary structure is assured (Fig. 8A). We found that the fraction of molecules that fully unfolded depended on the buffer conditions. In buffers containing TPP (TPP and TPP and Mg2+ buffers, Materials and Methods), the TPP riboswitch's aptamer domain showed significantly more full unfolding events than under apo or Mg2+-only conditions (Fig. 8B). These results are consistent with the smFRET results both in the apo and Mg2+ conditions, in which ∼70% of the population showed no-FRET and ∼30% exhibited FRET (Fig. 4A), and in the TPP and Mg2+ and TPP buffers, in which the population that exhibited FRET was significantly larger than in the apo and Mg2+ buffers (Fig. 4A). These results suggest that the states observed in the FRET-exhibiting populations (states a–d in Fig. 4) do correspond to dynamic aptamer domain conformational states in which the aptamer domain is predominantly folded, with significant secondary structure (as in Fig. 1).

FIGURE 8.

TPP riboswitch aptamer domain folding and unfolding kinetics. (A) A typical force extension curve (FEC) of the aptamer domain construct, in this case, in the presence of 0.5 mM TPP. The arrow represents an unfolding event, as characterized by a sudden drop in the force corresponding with an increase in the length of the tether. The red FEC represents stretching of a construct with a folded aptamer domain, and the black one represents stretching of the construct with an unfolded aptamer domain. (B) Fraction of the stretched TPP riboswitch aptamer domain constructs showing unfolding events without ligand (apo, N = 26), with 1 M Mg2+ (N = 20), with 0.5 mM TPP (N = 27), and with Mg2+ + TPP (N = 21). The occurrence of unfolding events, that is, the likelihood of finding folded TPP riboswitch aptamer domains, is significantly higher in the presence of TPP and Mg2+ and TPP buffers as compared to either in the presence of Mg2+ alone or in the absence of ligand (apo). The error bars represent the relative standard deviation of counting statistics. (C) Example time-series force data showing rapid transitions between the folded (F) and multiple unfolded states (UF) in the presence of 0.5 mM TPP, as determined by Gaussian mixture model analysis. Data was acquired using the optical tweezer at 1 kHz (gray) and smoothed using a 25 msec moving average filter (green). Sudden increases and decreases in force identify the transitions between the folded (F) and various unfolded states (UF1 and UF2). The relative population of states is shown with Gaussian fits (to the right of time-series data) for higher force (dark green, 20.0 pN) and lower force (light green, 18.6 pN). High-resolution time-series trace (inset) obtained at ∼18.6 pN pulling force. Folded states as fast as 25 msec and relatively long-lived unfolded states were detected at this force using change-point analysis. (D) Histogram of aptamer domain folded state dwell time measured in the presence of 4.8 mM TPP and 0.5 M Mg2+ and at 39.8 pN force, as determined by change point analysis. The characteristic dwell time was τoff = 108.4 ± 6.1 msec (single exponential decay fit constant ± std. error of the fit, N = 203). Histograms at additional forces can be found in Supplemental Information Figure S17. (E) The plot of log10(kon) as a function of the force in the presence of 4.8 mM TPP and 0.5 M Mg2+ where kon is the inverse of τoff (dwell time in the folded state, error bars represent the standard error of the exponential fits). Each data point (circle) represents 100–200 folded to unfolded transitions. The data were fitted using a weighted linear regression (solid line), and the dwell time at zero force was calculated to be τoff = 17.2 ± 2.6 sec (y-intercept ± std. error of the fit).

To study the aptamer domain's P1 switch helix unfolding transition, which is not directly addressable with the confocal fluorescence measurements due to the long-lived (>msec timescale) nature of folding states, we held the aptamer domain in a passive optical trap with force applied and probed 5′ end to 3′ end aptamer domain extension transitions. In the example time-traces (Fig. 8C), the aptamer domain exhibited fluctuations between multiple conformational states (Fig. 8C, inset, for detail). We found that the population of more folded aptamer domains states (F and UF1) was higher when subjected to lower applied tensile forces than higher forces, (Fig. 8C, right, for Gaussian fits of the histogram of state data). In this example case, at 20 pN (Fig. 8C, top), we found that the fully folded state (F), which likely corresponds to a hybridized P1 switch helix (Fig. 1B, bottom), was never observed, the first unfolded state (UF1), which likely corresponds to a dehybridized P1 switch helix (Fig. 1B, top), was observed approximately one-third of the time, and the second folded state (UF2), which likely corresponds to both the P1 and P2 helices being dehybridized (Anthony et al. 2012), was observed approximately two-thirds of the time. At 18.6 pN (Fig. 8C, bottom), we found that F was observed (Fig. 8C, inset), but rarely, UF1 was observed approximately two-thirds of the time, and UF2 was observed approximately one-third of the time. Together, these data show that we can probe the state of the TPP riboswitch aptamer domain's switch helix by using force to tilt the aptamer domain's base-pairing energy landscape and observing 5′ end to 3′ end aptamer domain distance transitions with the optical tweezers.

We collected time-trace data, like the example shown (Fig. 8C), in the TPP and Mg2+ buffer condition at multiple applied tensile loads. We used these data to determine the characteristic time of transitions between the F and unfolded states by fitting single exponential functions to conformational state dwell-time histograms acquired at various applied forces. For example, we found that the characteristic time of the unfolding transition from F to UF1 was τoff = 108.4 ± 6.1 msec at a tensile load of 39 pN for the aptamer domain in the Mg2+ and TPP buffer (Fig. 8D). To determine the unfolding rate at zero force (see Materials and Methods for details), we plotted log10(kon) as a function of force and extrapolated to zero force using Bell's equation (Fig. 8E). This calculation allowed for direct comparison to sensor helices and P1/P2 costacking data collected using smFRET and DMD simulations, which were done at zero applied tensile load. We found that the characteristic unfolding time of P1 switch helix at zero force was τoff = 17.2 ± 2.6 sec in the TPP and Mg2+ buffer (Fig. 8E), which is approximately four orders of magnitude slower than the slowest sensor helix dynamics timescales observed in smFRET (Fig. 5).

The measurements of F state residence times (transitions from F to UF1) were limited by the ability to distinguish events from the noise that is inherent in the optical tweezer system (due to the thermal fluctuation of the beads) in the apo and Mg2+-alone buffers, and they were just barely distinguishable from the noise in the TPP alone buffer, at low enough forces to probe the F to UF1 transitions (<70 msec at forces <20 pN). Additionally, the F to UF1 transition corresponded to a change in length of 7.9 ± 0.3 nm. Assuming an internucleotide distance of 4.6 Å (Himbert et al. 2016) for ssRNA, this extension corresponds to ∼17 nt, which is in agreement with the estimated upper limit for extension upon P1 switch helix unfolding, 8.5 nm for 18 nt (7 nt at 3′ end, 6 nt at 5′ end, and 5 nt between the sensor helices). Further, unfolding of the riboswitch under increasing pulling force indicates an extension of 34.1 ± 0.8 nm for complete unfolding (e.g., Fig. 8A), corresponding to approximately 74 unpaired nucleotides. This is only two fewer than the total number (76 nt) of nucleotides that are either base-paired or on the opposite side of the P1 helix from the DNA handles where pulling occurs in the aptamer domain. These data all suggest that F to UF1 transitions correspond to P1 switch helix unfolding events, and that the switch helix is much more stable in the presence of saturating conditions of TPP and Mg2+ than with either TPP or Mg2+ alone.

Together, the optical tweezer data confirm that the aptamer domain folding is stabilized by TPP (Fig. 8B), just as smFRET data suggested (Fig. 4A). Additionally, once the P1 switch helix is hybridized, it is much more stable in the presence of saturating conditions of TPP and Mg2+ (Fig. 8E) than the sensor helix conformations that exchange on timescales faster than milliseconds (Fig. 5C), as suggested by the DMD simulation in similar conditions. In total, our smFRET, DMD simulation, and optical tweezer results are consistent with a model of rapid transitions between ensembles of sensor helix states that underlie and tune much slower transitions between TPP riboswitch ON-like and OFF-like P1 switch helix conformations (UF1 and F, respectively), rather than a model that relies on the direct coupling between long-lived, binary sensor helix conformational and P1 switch helix base-paired states. They suggest that P1 switch helix stabilization is achieved even without significantly long-lived closed sensor helix conformations, but that it still requires both TPP and Mg2+ to be present.

DISCUSSION

Our results suggest that the aptamer domain of the TPP riboswitch is a dynamic structure that is tuned by the presence of its coligands. Specifically, we found that the aptamer domain: (i) exhibits multiple dynamic processes occurring across many orders of magnitude in temporal scale; (ii) is more likely to sample folded and closed states in the presence of TPP, Mg2+, or both; (iii) differentially responds to Mg2+ and to TPP, suggesting multiple dynamic pathways; (iv) requires nonphysiological conditions of Mg2+ in addition to TPP to generate long-lived, closed sensor helix ensembles; and (v) maintains long-lived closed switch helix states, even in the absence of long-lived, closed sensor helix ensembles. Moreover, aptamer domain dynamics occurs even in the presence of saturating amounts of TPP and Mg2+ that most favor the fully folded, switch helix hybridized, OFF-like conformations (Figs. 6–8), likely due in part to Mg2+ enabling the sensor helices to sample intermediate ensembles (state b). Further, these dynamics are consistent with the small but significant decrease in the effective hydrodynamic radius of the aptamer domain in the presence of saturating TPP (Supplemental Information, Fig. S2). A nondynamic, open sensor helix state would have had a larger expected hydrodynamic radius due to the addition of the coordinating TPP, and a nondynamic, closed sensor helix state would have had a larger decrease in hydrodynamic radius associated with a stable, compact, X-ray crystal structure-like conformation.

We further identified five sensor helix configurational ensembles of the Arabidopsis thaliana TPP riboswitch's aptamer domain as a function of Mg2+ and TPP ligand condition (a, b, c, d, and no-FRET, Fig. 4). Each ligand condition gave rise to a unique subset of two sensor helix configurational ensembles (Fig. 4B) in addition to a significant population exhibiting open, no-FRET (Fig. 4A) and unfolded (Fig. 8B) conformations, which became smaller in the presence of ligands (Figs. 4A, 8B). In the absence of ligand, the FRET-exhibiting aptamer domain population was predominantly in a sensor helix open (a) conformation with a smaller fraction of the population in an intermediate (c) ensemble (Fig. 4B). The presence of excess Mg2+ alone increased the relative population of the dynamic intermediate c ensemble but maintained a smaller population of the open a conformation (Fig. 4B). A saturating concentration of TPP alone reduced sampling of the open a ensemble. Instead, it allowed the sensor helices to sample the closed d state, sometimes exchanging with a second intermediate b ensemble (Fig. 4B). In the presence of both ligands, the sensor helices sampled both the a and d conformations, with the kinetics significantly favoring the closed d state (Fig. 4B), thus suggesting that TPP introduces the closed (d) state but that Mg2+ is required to stabilize it, due to the splitting of the 〈τ〉 distribution in saturating Mg2+ and TPP conditions (Fig. 3). Our computational modeling independently corroborates this ensemble switching scheme by identifying structural ensembles through clustering analysis (Fig. 7) that correlate with the experimentally determined scheme (α, β, γ, and δ basins in the free-energy landscape, respectively, Fig. 6B). Additionally, we found stable, folded (F) switch helix states primarily in the presence of saturating Mg2+ and TPP (Fig. 8). These data suggest that the sensor helix closed d state must be significantly populated, but that sensor helices need only to sample the d state dynamically, for the aptamer domain to persist in the OFF-like state for long (second timescale) times.

Together, our results lead us to propose a new model for the conformational landscape of the TPP riboswitch's aptamer domain as it binds to its Mg2+ and TPP ligands (Fig. 9). The model includes multiple conformational states comprised of configurational subensembles that are accessible in each buffer condition with the relative populations and exchange between these populations being sensitive to ligand concentrations. Moreover, our results suggest the aptamer domain can take one of two independent closing pathways, depending on the order of Mg2+ and TPP binding (Fig. 9). The semi-open and open ensembles (very low and no-FRET, respectively; α ensembles) share an exit to two pathways that go from open-like to closed-like configurational ensembles because the α ensembles exchange rapidly and transitions from the open and open-like ensembles are not mutually resolvable, although the low-FRET and no-FRET populations are resolved via TCSPC. Because the MFD histogram for the apo conditions, in which the α ensemble is most populated, exhibits a compact distribution near the static FRET line (Fig. 2A) and the fluorescence decay histogram indicates multiple lifetime populations rather than a single population, we conclude that rapid dynamics (on the order of µs, tR(1)) dominate the exchange between the α ensembles. Conversely, the Mg2+ and TPP conditions gave rise to two distinct peaks in the 〈τ〉 distribution, indicating that dynamic transitions between the high-FRET, OFF-like ensemble (δ) and the intermediate ensembles take place on timescales comparable to or longer than the diffusion time of the aptamer domains through the confocal volume. Thus, we assigned to this limiting reaction step the slowest timescales observed by fFCS (tR(4)). Further, the Mg2+ and TPP conditions exhibited the fastest correlation times among all the conditions (Fig. 4), indicating a degree of cooperativity between Mg2+ and TPP that favors transitions to the OFF-like δ ensemble. This conclusion is further supported by the observation that the high-FRET ensemble only formed a distinct τ peak in extreme, saturating conditions of Mg2+ and TPP, whereas apo, Mg2+ and TPP buffer conditions favored open or intermediate ensembles. Of the two remaining correlation terms from fFCS, we assigned tR(2) to transitions between the α ensemble and the intermediate γ ensemble (Mg2+ binding) and tR(3) to β ensemble (TPP binding) transitions because the MFD histogram for TPP buffer exhibited a larger deviation from the static FRET line than in Mg2+ conditions, indicating slower dynamics. Additionally, TPP buffer conditions led to a larger shift toward high-FRET configurations from low-FRET configurations, as observed by population fractions from donor fluorescence decay fitting (Fig. 4), which was likely associated with the larger-scale conformational changes observed upon TPP binding.

FIGURE 9.

Conformational landscape of TPP riboswitch aptamer domain dynamics and closing pathway. (A) The completely unfolded state accessible by tweezer experiments is shown along with the conformational states observed by smFRET. The long timescale of unfolding relative to kinetics observed by smFRET indicates that a much larger energy barrier must be overcome in these transitions. States observed by smFRET are highlighted in gray. (B) Zoomed representation of states highlighted in A in the context of Mg2+ and TPP buffers. Two pathways exist from the open/semi-open ensembles (α) to the δ state. Conformational states and cartoon representations are shown along the closing pathway of the riboswitch representing the different interactions that are favored either via presence of Mg2+ or TPP (γ, and β ensembles) where the P1/P2 stacking leads to rapid fluctuations in the paring of nucleotides along the P1 switch helix (shaded in gray in the schematic representation). This is followed by stabilization of the P1 switch helix in the δ ensemble. Labels for the times associated with each of these transitions are assigned based on fFCS and MFD. The names for the observed ensembles correspond to those that qualitatively agree between smFRET experiments and DMD simulations. The apparent energy gaps between the conformations and pathways are based on the relative magnitudes of the observed dynamic timescales, with faster transition timescales indicating lower energy barriers between conformations.

Our model of a highly dynamic TPP riboswitch aptamer domain is corroborated by DMD simulations. In all buffer conditions, ensembles α, β, and γ dynamically exchange, with the relative populations sensitive to the presence of Mg2+ and TPP in a manner similarly observed by TCSPC. However, when both ligands were present, the OFF-like δ ensemble was observed and stabilized (Figs. 6, 7).

In addition to the fast ensemble switching of the sensor helices, we found that the dynamics in the P1 switch helix, even under a constant destabilizing force, are significantly slower (τoff = 108.4 ± 6.1 msec under 39.8 pN force, Fig. 8) than the intrinsic conformational dynamics of the sensor helices (less than a few milliseconds, Fig. 5). Under this relatively high force, we observed dynamic switching between a fully folded state (F), which likely corresponds to base-pairing in the P1 sensor helix, and two less-folded states (UF1, UF2), which likely correspond to opening of the P1 sensor helix (Fig. 8C) and P2 (Anthony et al. 2012), respectively. Using Bell's equation, we calculated that the equilibrium (zero-force) transition rate between the P1 sensor helix folded and unfolded states in the presence of both TPP and Mg2+ (Fig. 8E) was four orders of magnitude slower than the slowest transition rate between structural conformations in the sensor helices, also observed in the presence of both TPP and Mg2+ (Fig. 5C). The observed differences between the configuration timescales of the sensor and switch helices indicates that, while TPP and Mg2+ recognition and binding in the aptamer domain certainly stabilize the P1 switch helix (Figs. 6–8), transitions between the structural configurational ensembles cannot correspond one-to-one to switch helix dehybridization and TPP riboswitch aptamer domain ON/OFF transitions (Fig. 9B). Thus, while our results show that the larger-scale ON/OFF behavior of the switch helix takes place on longer timescales, fast ensemble switching of sensor helix conformation and configuration may allow more subtle biasing of the riboswitch toward ON-like or OFF-like states with a rheostat-like tunability given by the switching kinetics.

By utilizing various techniques that probe multiple timescales of dynamics, we determined that the TPP riboswitch aptamer domain exhibits conformational dynamics spanning several orders of magnitude in time at equilibrium conditions and that these dynamics are tunable by changes in Mg2+ and TPP concentrations. Fast kinetics observed by FRET, fFCS, and DMD explain the ensemble switching between sensor helix conformations that alter the likelihood of, but need not directly induce, switch helix folding/unfolding state transitions that occur at slower timescales, as observed by optical tweezers and suggested by DMD simulations. In simulations and optical tweezer experiments, presence of Mg2+ and TPP correlated with increased stability of and transitions to OFF-like, base-paired switch helix conformations. Because changes in Mg2+ and TPP concentrations also regulate fast structural state kinetics to favor subsets of configurational ensembles, which do not correlate one-to-one with the ON/OFF states, our conclusions support a model in which the TPP riboswitch function can be finely tuned with a kinetic rheostat mechanism rather than a binary open/closed aptamer domain conformation model that directly correlates to the ON/OFF regulation of gene expression.

The folding dynamics of aptamer domains are crucial to riboswitch function, and divalent ion coligand binding to the aptamer domain can facilitate this folding (Savinov et al. 2014). Further, the TPP riboswitch aptamer domain binds with TPP in a step-wise fashion (Anthony et al. 2012). Our data indicates that TPP is able to stabilize a population of closed aptamer domain tertiary structures (Figs. 3, 6, and 7) only through coordinating with a divalent Mg2+ ion. Further, all four observed, distinct timescales occur in the submillisecond range, which is well below the timescales associated with slow, high-energy conformational changes or transcriptional folding kinetics (Figs. 5, 6, and 8). Thus, the short-timescale dynamics of Mg2+ and TPP aptamer domain binding regulate the longer-timescale ON/OFF conformational transitions, giving rise to a stepwise closing pathway. Our results suggest a model of TPP riboswitch behavior that spans timescales from nanoseconds to milliseconds. This model couples TPP and Mg2+ binding to both short-timescale configurational changes in the sensor helices and long-timescale switch helix base-pairing in the aptamer domain, bridging the gap between previously suggested fast and slow binding recognition and ON/OFF tuning mechanisms (Thore et al. 2006; Haller et al. 2013). These results also constitute direct evidence that some riboswitches function more as rheostats or dimmers in which ligand binding does not correlate 1:1 with ON/OFF switching, rather than as discrete ON/OFF switches, as suggested by others (Baird et al. 2010; Venkata Subbaiah et al. 2019). We propose that such rheostatic tuning could allow for fast response times for regulation of TPP concentrations in cells.

MATERIALS AND METHODS

TPP riboswitch design and labeling

The sequence of Arabidopsis thaliana TPP riboswitch aptamer domain was truncated to match the crystal structure (Thore et al. 2006) and to compare our results with prior structural information directly. To determine the best locations for fluorophore labeling and to maximize the resolution of expected conformational changes, we ran a coarse-grained simulation (Ding et al. 2008a) whose trajectories were postprocessed for modeling of the fluorescent markers at all possible pair distances. We calculated interlabel distances at each point in the trajectory by sampling accessible volumes generated for coarse-grained dyes linked in silico to the RNA backbone, as previously described (Kalinin et al. 2012). Based on all possible interlabel distances, we selected positions 25 and 56 for maximum resolvability of distance changes between two identified conformations from the coarse-grained simulations (Supplemental Information [SI] Fig. S1). The aptamer domain was synthesized following the sequence 5′-GGG ACC AGG GGU GCU UGU UCA CAG gCU GAG AAA GUC CCU UUG AAC CUG AAC AGG GuA AUG CCU GCG CAG GGA GUG UC-3′ where two selected modified nucleotides [g(25) = 2′-dG-(N2-C6-amino) and u(56) = rU(2′-O-propargyl)] were incorporated for subsequential fluorophore labeling via orthogonal click chemistry (Cy5-NHS to g and Alexa Fluor 488-azide 5/6 isomer to u). As a control during binding experiments, we used an unmodified g/G and u/U sequence. The oligo was purified before labeling via polyacrylamide gel electrophoresis (PAGE) purification to ensure that the synthesized riboswitch product was of the proper size. Labeling was performed using saturating conditions of Cy5 and Alexa 488 to enhance the likelihood of 1:1 donor-to-acceptor stoichiometry followed by desalting to remove excess, unreacted fluorophores. Labeling postsynthesis also ensures that the initial riboswitch aptamer domain fold is unperturbed by the presence of the fluorophores. We prepared the corresponding donor-only reference sample by only replacing the nucleotide u(56) for donor labeling and maintaining the remaining native sequence during synthesis. Purimex GmbH and Biosynthesis Inc. supplied synthesized and purified riboswitch aptamer. Labeling was performed by Purimex GmbH and by us using Click chemistry as follows. Aptamer domain was premixed to 119 µM into a 0.1 M sodium bicarbonate pH 8.2 buffer prepared in RNase-free water. Separately of 100 mM CuSO4, 200 mM THPTA, and 100 mM sodium ascorbate aliquots were prepared in the same buffer. After mixing and waiting ∼30 min, the aptamer domain was incubated with a fivefold molar excess of Alexa 488, 25-fold molar excess CuSO4, the THPTA mixture, and 40-fold excess sodium ascorbate for 2 h. Then, the sample was buffer exchanged back into a pH 8.2 using centrifugal filter units (Amicon Ultra, 10 kDA MWCO) to ensure suitable pH for Cy5 mono-NHS ester labeling. Cy5 was added to 10-fold molar excess and incubated again for 2 h. The labeled aptamer domain was buffer exchanged into Tris buffer (20 mM Tris, 200 mM NaCl at pH 7.0) just before adding MgCl2 and TPP for size exclusion chromatography (SEC) and single-molecule experiments.

To probe whether the fluorescent labels impact the aptamer's ability to fold or bind to TPP, we used SEC. Briefly, we used an NGC Discover 10 chromatography system (#7880009, Bio-Rad Laboratories, Inc.) and equilibrated an ENrich SEC70 10 × 300 size exclusion column (#7801070, Bio-Rad Laboratories, Inc.) with 20 mM Tris, 200 mM NaCl, 200 mM MgCl2 at pH 7.0. We ran 50 µL wild type (unlabeled) and labeled aptamer domain with varying concentrations of TPP. For monitoring the eluate, we used 260 nm (aptamer domain), 480 nm (Alexa 488), and 650 nm (Cy5) wavelengths (Supplemental Information Fig. S2). For generating a standard curve, we ran a 50 µL sample of a Gel Filtration Standard (#1511901, Bio-Rad Laboratories, Inc.) at the beginning and the end of each day (Supplemental Information Table S1). We found that the labeled riboswitch compacts upon binding to TPP as evidenced by the reduction on the hydrodynamic radius (Supplemental Information Figure S2) following the same binding isotherm as the wild-type riboswitch; confirming that the fluorescent labels do not significantly impact the TPP riboswitch aptamer domain's ability to fold or bind to TPP. Moreover, the determined binding affinity of 17.5 ± 7.6 nM was consistent with previously reported values of 20 nM (Winkler et al. 2002a; Yamauchi et al. 2005).

Multiparameter fluorescence spectroscopy (MFS) for single-molecule FRET experiments (smFRET)

A home-built confocal system was used for the MFS-smFRET measurements and calibration, as described in detail elsewhere (Ma et al. 2017; Hellenkamp et al. 2018; Yanez Orozco et al. 2018; Hamilton and Sanabria 2019). To briefly summarize, we use Pulsed Interleaved Excitation (PIE) (Kudryavtsev et al. 2012), or Alternating Laser EXcitation (ALEX) (Kapanidis et al. 2005), with two diode lasers (Model LDH-D-C-485 at 485 nm and laser LDH-D-C-640 at 640 nm; PicoQuant) operating at 40 MHz with 25 nsec interleaved time. The power at the focal point of the 60× Olympus objective was set to 60 µW for the 485 nm laser line and 23 µW for the 640 nm excitation. The emitted photons were collected from the same objective and then spatially filtered by a 70 µm pinhole. The signal was separated into parallel and perpendicular polarizations with respect to the polarization of the excitation source, then split further into two different spectral windows defined as “green” and “red.” We used four photon-counting detectors, two for the green (PMA Hybrid model 40 PicoQuant) and two for the red (PMA Hybrid model 50, PicoQuant) channels. A time-correlated single-photon counting (TCSPC) module (HydraHarp 400, PicoQuant) in time-tagged time-resolved (TTTR) mode was used for data registration.

For smFRET experiments, samples were diluted to pM concentration in one of four buffers (apo buffer: 20 mM Tris and 200 mM NaCl at pH 7.0; Mg: 20 mM Tris, 200 mM NaCl, and 1 M MgCl2 at pH 7.0; TPP buffer: 20 mM Tris, 200 mM NaCl, and 2 mg/mL TPP at pH 7.0; Mg: 20 mM Tris, 200 mM NaCl, 1 M MgCl2, and 2 mg/mL TPP at pH 7.0). The buffers were charcoal-filtered to remove residual impurities. The concentrations of TPP and Mg2+ in these solutions are high to ensure saturating conditions. The particularly high concentration of Mg2+ is well beyond the physiological concentration of ∼1–3 mM. This concentration was used to demonstrate the extreme concentration necessary to stabilize a population of the aptamer domains in the closed conformation (see Results, Fig. 3). At pM riboswitch aptamer domain concentrations, we observed ∼2 molecules per second in the detection volume. NUNC chambers (Lab-Tek, Thermo Scientific) were pretreated with a solution of 0.1% Tween 20 (Thermo Scientific) in water for 30 min and then rinsed with ddH2O to prevent adsorption artifacts. Data collection times were varied from several minutes up to 10 h. In long-time measurements, an oil immersion liquid with an index of refraction matching water (Carl Zeiss AG) was used to prevent drying out of the immersion water. Control experiments consisted of measuring distilled water to characterize the instrument response function (IRF), buffers for background subtraction, and nM concentration green and red standard dyes (Rhodamine 110, Rhodamine 101, and Alexa 647) in water solutions for calibration of green and red channels.

The data analysis used software developed in the Seidel group and is available at http://www.mpc.hhu.de/en/software (Sisamakis et al. 2010; Sindbert et al. 2011; Kalinin et al. 2012; Dolino et al. 2016; Ma et al. 2017). MFS accounts for many challenges associated with smFRET experiments, including photo-bleaching, cis-trans acceptor blinking, and dynamic and static quenching. We used FRET indicators such as the ratio of donor to acceptor fluorescence (F) and the average fluorescence lifetime of the donor (〈τ〉), integrated over single-molecule events, or bursts (Kühnemuth and Seidel 2001; Widengren et al. 2006; Sisamakis et al. 2010). Details on the correction parameters obtained during calibration and the computed FRET lines can be found in Supplemental Information Table S2. Fluorescence lifetime was used, in part, because it is insensitive to instrumental parameters and is thus independent of correction parameters. Bursts were selected for analysis based on the relative fraction of photons observed from the donor vs. the donor plus the acceptor fluorescence signals following their respective direct excitation pulses (S) in order to ensure 1:1 active donor-to-acceptor stoichiometry of analyzed bursts (Supplemental Information Fig. S3; Supplemental Information Methods; Kudryavtsev et al. 2012). Several standards experiments were performed to quantify instrumental correction parameters and ensure the correct selection of donor-acceptor molecules via SPIE (Supplemental Information Fig. S4). Selection against molecules with either inactive donor or acceptor is accomplished with the parameter SPIE, as discussed in the next section and the Supplemental Information Methods. Further, relative dye concentrations of Cy5 and Alexa488 were measured as a first indicator of sample labeling quality (Supplemental Information Fig. S5).

Time-correlated single photon counting (TCSPC) and filtered fluorescence correlation spectroscopy (fFCS)

To identify stable populations that persist longer than the duration of the fluorescence lifetime of the FRET-labeled riboswitch aptamer domain (ns timescales), we compared the donor-only reference aptamer domain with the FRET-labeled aptamer domain in the same buffer conditions using TCSPC. Time-resolved fluorescence decay curves from photon arrival-time histograms (Supplemental Information Fig. S6, henceforth fluorescence decay histograms) ensure enough statistical rigor through large photon counts (on average >109 photons for our FRET-labeled samples). Therefore, we determined the interdye distance distributions from fluorescence decay data by considering the donor-only and FRET-labeled aptamer domain at saturated conditions of Mg2+, TPP, and Mg2+ and TPP, with the corresponding control in the absence of ligand (apo).

To arrive at the correct number of configurational ensembles, we fit the fluorescence decay histograms with an increasing number of fluorescence decay lifetimes corresponding to Gaussian-distributed FRET distances. Gaussian-distributed distances provide a physical representation of the mobility of the fluorescent labels and intrinsic dynamics of the chain. The mean and variance of the distance distribution are both relevant statistical parameters in smFRET experiments (Gopich and Szabo 2005). Therefore, we fixed the variance according to benchmark studies (Sindbert et al. 2011; Kalinin et al. 2012; Ma et al. 2017; Hellenkamp et al. 2018) and optimized for the mean. Next, we evaluated the best results by visual inspection of the weighted residuals, monitoring the distribution of the autocorrelation of the residuals, and by a statistical improvement of the figure of merit, χ2 (Supplemental Information, Equation 7; Box George 1960; Lakowicz 2007). The model of the fluorescence decay includes a no-FRET population that accounts for inactive (15%–55% of Cy5-labeled molecules; Ha and Tinnefeld 2012) acceptor molecules and the extended conformations observed in the multidimensional histograms beyond what could be resolved by FRET. As expected with the addition of distributions, the figure of merit, χ2, and quality of the fit improves. We selected a model of two Gaussian distributions because the change of χ2 improved from an average of 8.52 over all conditions when using a single Gaussian distribution to 2.42 when using two Gaussian distributed FRET distances. Adding a third Gaussian distribution only marginally improved the average χ2 to 2.05 (Supplemental Information Table S3).

To increase the photon counts and improve the statistics for data analysis, the TCSPC and fFCS measurements were performed using the MFS setup described above, but the power at objective was set for 120 µW for the 485 nm laser line and 39 µW for the 640 nm excitation while labeled samples were measured at nM concentration. The durations of these experiments were 10 h each.

TCSPC and fFCS use the same data set but consist of different analyses. TCSPC-generated fluorescence decays were jointly analyzed, along with a reference sample, to identify configurational ensembles that are stable in the nanosecond timescale with high precision. fFCS, a variation of standard FCS (Elson and Magde 1974; Magde et al. 1974) that probes changes in the fluorescence lifetime (Enderlein and Gregor 2005), spectral information, and anisotropy (Felekyan et al. 2012), was used to probe transitions between selected conformers or pseudoconformers characterized by specific fluorescence decay patterns. Thus, it was possible to identify changes in FRET due to conformational changes that occurred at timescales faster than the diffusion time. Details of the functional forms used to process the data, the statistical tests, and error analysis are provided in the Supplemental Information.

For fFCS, we systematically tested models of increasing complexity through the addition of relaxation times (tR) and finally chose a model by iteratively varying the number of terms in each factor in the fFCS function to find the number of terms necessary to adequately fit the data, as determined by visual inspection of the residuals, and to provide a significant relative improvement in the χ2 compared to the use of one fewer term. The global fit model we used here has four distinct relaxation terms describing conformational transitions and two exponential terms to account for dark state kinetics (in the case of the acceptor, this can be due to cis-trans isomerization, and in the case of the donor, this represents a transition to triplet state) and photobleaching (Felekyan et al. 2012).

Discrete molecular dynamics (DMD) simulations

DMD is an event-driven approach to molecular dynamics (MD) simulations that utilizes discontinuous step-functions instead of continuous functions to model the interaction potentials between atoms (Rapaport 2004). Compared to traditional MD simulations that calculate the forces and accelerations of all atoms, DMD computes the ballistic equation of motion, considering conservation laws for energy, momentum, and angular momentum, only for those atoms involved in a collision event, defined as interactions between two atoms within the distances associated with potential steps. This process reduces the total number of calculations required per timestep, allowing simulation of longer timescales and length-scales in complex biomolecular systems (Ding et al. 2008b; Emperador et al. 2010).

The initial structural refinement of the RNA was performed via DMD with a coarse-grained RNA model. The coarse-grained RNA model uses three beads corresponding to the phosphate, sugar, and base groups of each nucleotide (Ding et al. 2008a). The interaction potentials include base-pairing, stacking, and the effective loop energy, which were tabulated from experimental measurements (Mathews et al. 1999). DMD simulations with the coarse-grained RNA model can successfully predict RNA 3D structures and folding kinetics (Ding and Dokholyan 2006; Miao et al. 2017). Thus, initial labeling design was based on these coarse-grained simulations.

We also performed all-atom DMD simulations to uncover the effects of ligands on RNA dynamics. The atomistic RNA model adopts the united-atom approach to model nucleotides, which includes all heavy atoms and polar hydrogen atoms. Similar to all-atom DMD simulations of proteins (Ding et al. 2008b), bonded interactions, including covalent bonds, bond angles, and dihedrals, were assigned according to a statistical analysis of high-resolution RNA structures in the Protein Data Bank (PDB). We used the Medusa force field (Ding and Dokholyan 2006) to describe the nonbonded interactions, including van der Waals (VDW), hydrogen bonding, solvation, and electrostatics. The VDW parameters were taken from the CHARMM19 united-atom force field, and the distance and angular dependent hydrogen bonds were explicitly modeled by a reaction-like algorithm (Ding et al. 2003). The screened electrostatic interactions between charged atoms were computed using the Debye-Hückel approximation, while the solvation term for implicit solvent simulations was modeled by the Lazaridis−Karplus EEF1 (Effective Energy Function 1) model (Lazaridis and Karplus 1999). Details of the all-atom RNA model in DMD simulations were described previously (Williams et al. 2017), and the method has been tested in the community-supported binding prediction of RNA 3D structures (Miao et al. 2020). The ligand parameters were based on MedusaScore, an extension of the Medusa force field that has been parametrized to recapitulate protein–ligand binding affinities (Yin et al. 2008). For interactions of Mg2+ with RNA, we only considered VDW, solvation, and electrostatics without including the effect of metal ion coordination.

The initial structure of the TPP-riboswitch aptamer domain was obtained from the Protein Data Bank (PDB ID: 2CKY [Thore et al. 2006]). Periodic boundary conditions were imposed in all simulations, and the riboswitch aptamer domain was initially positioned in a cubic box with an edge length of 100 Å. Counter ions were added to maintain the overall charge neutrality of the system. Table S4 of Supplemental Information lists the number of Mg2+ and Na+ ions in the simulated systems. Production runs followed 2000 timestep energy minimization runs, where a DMD timestep corresponded to ∼50 fsec. To efficiently sample the conformational space of the free-energy landscape as well as the dynamics of sensor helix closing pathways, we used replica-exchange DMD (REXDMD) simulations in our study. In the replica-exchange scheme, 16 replicas in consecutive temperature ranges within the total range 260–335 K were simulated via the exchange of temperature between replicas at periodic time intervals, thereby overcoming any local energy barriers during the conformational sampling of the energy landscape. The simulation temperature was maintained using the Andersen thermostat (Andersen 1980). Each replica lasted 4 × 106 timesteps, or ∼200 nsec of simulation time totaling ∼3.2 μsec of simulation time.

Since the smFRET experiments focused on the open and close conformational dynamics of sensor helices in the aptamer and the RNA was expected to maintain the native-like base-pairing at the observed fast timescales, we imposed base-pair constraints in our all-atom REXDMD simulations to decouple the three-dimensional conformational dynamics of tertiary structures from the folding and unfolding dynamics of helices that occur on longer timescales (Ding et al. 2012). Similarly, we also assumed that the TPP stayed in its native binding pocket by imposing distance constraints between TPP and the two stacking nucleotides.

Weighted histogram analysis method (WHAM)

The thermodynamics of TPP riboswitch aptamer domain conformational dynamics were computed with WHAM analysis (Kumar et al. 1992; Feig et al. 2004). The density of states, g(E), was calculated self-consistently by combining the energy histograms from all the 16 replicas, using the last 100 nsec of replica-exchange trajectories. The potential of mean force (PMF), that is, the effective free-energy landscape, was computed as a function of a given physical parameter, x: where kB is the Boltzmann constant, T denotes the temperature of interest, and P(x) corresponds to the probability of finding the aptamer domain with specific values of x. P(x) was computed according to: where β = 1/KBT and P(E|x) is the conditional probability of finding conformations with parameter value x while the energy is in the interval (E, E + dE). The x parameter can be comprised of multiple variables. For example, we evaluated the conformational free-energy landscape of the aptamer domain with respect to both the inter-arm distance, D, as measured between the 2′-hydroxyls of labeled G25 and U56, and with respect to the P1/P2 costacking distance, D, as measured between U39 in P2 and G73 in P1.

Clustering analysis

We used the OC hierarchical clustering program (Barton 1993, 2002, 2004) to group similar configurations of the riboswitch aptamer domain. Depending on an input pairwise root-mean-square distance (RMSD) matrix, the clustering algorithm iteratively combined the two closest clusters into one. The “cluster distance” was calculated based on all pairwise distances between elements of two corresponding clusters. We used the mean of all values to compute the distance between two clusters, and the centroid structure of each cluster was chosen such that it had the smallest average distance to other elements in the cluster.

RNA preparation for optical tweezer experiments

Optical tweezer experiments consisted of the TPP riboswitch aptamer domain flanked at both ends by single-stranded DNA extensions (sequence in the Supplemental Information) annealed to 2 kbp DNA handles, each with a single-strand overhang. The DNA handles were amplified from an arbitrarily selected stretch of the pMAL-c5X vector (N8108, New England Biolabs, Inc.) with primer sequences detailed in the Supplemental Information. The 5′ end of the forward primer of handle 1 was modified with an amine group (NH2) for binding to the N-hydroxysuccinimide (NHS) group of N-hydroxysuccinimide functionalized polyethylene glycol (PEG-NHS) on the coverslip. To ensure a single-strand overhang, we used an autosticky reverse primer (Gal et al. 1999) with the tetrahydrofuran abasic site inserted 35 bases from its 5′ end, blocking the synthesis of the complementary strand for a 35-nt 5′ overhang that anneals to the 5′ DNA extension of the RNA construct. The 5′ end of the forward primer for handle 2 was modified with biotin for linking to streptavidin-coated beads. Three phosphorothioate bonds were inserted 37 nt bases from the 5′ end of the forward primer, inhibiting nuclease digestion (Citartan et al. 2011). Handle 2 was then digested with lambda exonuclease (M0262, New England Biolabs, Inc.), which removes nucleotides from linearized double-stranded DNA in the 5′ to 3′ direction, leaving a 3′ overhang on the complementary strand that anneals to the 3′ DNA extension of the RNA construct (Nikiforov et al. 1994).

The efficiency of RNA construct to DNA handle annealing was determined empirically by observing which dilutions of the constructs gave enough tethers to beads during tweezers experiments (Wen et al. 2007). To minimize the likelihood of finding RNA constructs without both handles, we annealed the handles with the RNA construct in 1:0.25:1 (H1:RNA:H2) molar ratio in PBS (pH 8.3) buffer at room temperature for 2 h.

Flow cells (∼10 μL) were made by cutting a channel (2–3 mm wide) into Parafilm sandwiched between a glass slide and 22 × 22 mm cover glass, precleaned with piranha solution (sulfuric acid and hydrogen peroxide in the ratio of 3:1) followed by a series of ethanol and distilled deionized water washing steps (Gell et al. 2010). The piranha solution, being a strong oxidizing agent, cleaned off most of the organic residues and hydroxylates from the surface to facilitate the subsequent covalent bonding with silanes.

We conjugated the RNA construct to coverslip surface by adapting a standard protocol (Schlingman et al. 2011). Briefly, a 2% (w/v) silane-PEG-NHS (Nanocs, Inc.) solution prepared in pure, dry DMSO was flowed into the flow cell and incubated at room temperature for 2 h. The annealed handle-RNA construct with -NH2 and biotin labels at either end was then flowed into the flow cell and incubated in a humidity chamber for 1.5 h to allow covalent bonding between the NHS and the amine modification on DNA handle. 1.1% (w/v) 1.05 µm diameter streptavidin-coated polystyrene beads (Bangs Laboratories, Inc.) were washed twice in one of the buffers detailed above (i.e., apo buffer, TPP buffer, Mg2+ buffer, or Mg2+ and TPP buffer) plus 0.1% Tween 20. The bead solution was then diluted 20 times and flowed into the flow cell. The flow cell was then rinsed with 1% Pluronic F-127 (BASF Group), and the ends were sealed with silicone grease. The flow channel was left to sit at room temperature for 20–24 h.

Optical tweezers and data acquisition

A custom-built, single-beam optical tweezer with a 1064 nm, ytterbium fiber laser (YLR-10-1064-LP, IPG Photonics Corp.) focused into the sample plane by a CFI60 Plan Apochromat Lambda 60× N.A. 1.4 Oil Immersion Objective Lens (Nikon Instruments Inc.) was used to trap the beads and pull on the handles. Measurements of the change in the construct length resulting from the pulling were made by determining the displacement of the trapped bead from the optical trap center using quadrant photodiode (QPD, QP45-Q-HVSD, First Sensor, Inc.) back focal plane detection (Gittes and Schmidt 1998) as well as the translation of the nanostage (Nano-LP100, Mad City Labs, Inc.) holding the cover glass. Force measurements were made by obtaining calibration parameters (Tolic-Norrelykke et al. 2006), detection sensitivity (β in V/µm) and stiffness (κ in pN/nm), from a nonlinear Lorentzian fit of the power spectrum density of the QPD response to a trapped bead subject to thermal fluctuations. The parameters were used to convert the change in QPD voltage signal due to displacement of a trapped bead from the trap center into a corresponding distance (nm) and then that displacement to a force (pN). A trap stiffness of 0.2–0.3 pN/nm was used for all measurements. The data were acquired using an FPGA (PXI-7854R, NI) and custom-written LabView VIs (NI).

The passive force time-series data were recorded at 1000 Hz for 100 sec at each force and filtered offline using a 25 msec moving average filter (MATLAB, MathWorks Corp.). Under the passive mode force setting, the optical tweezer was used with a fixed trap where both the force and extension vary as the molecule undergoes structural transitions.

The MATLAB function findchangepts was used to partition each of the 100 sec time-series data into various force levels. The findchangepts function minimizes the sum of the residual error of a particular region in the data set from its local mean and returns a step at which the mean changes the most. Depending upon the position of the stage, and hence the average force-extension of the RNA construct, the obtained mean forces were fitted to two or three Gaussian functions (using the MATLAB Gaussian mixture model). The transition forces determined from the multiple-Gaussian fit were used to define the folded (denoted by higher force or shorter extension) and unfolded (denoted by lower force or longer extensions) states. The mean dwell time of a folded state was obtained by fitting the dwell times distribution of the state to a single decay exponential function and calculating the characteristic dwell time constant from the fit. The corresponding rate coefficient (unfolding rate, k) was calculated as the inverse of the mean dwell time. To first order approximations, this rate constant depends exponentially upon the applied force as given by two parameters Bell's equation (Walcott 2008): where F is the applied force, k is the rate at zero force, and k is Boltzmann's constant, T is the absolute temperature, and x is the characteristic distance of the stretch. The logarithm of k plotted as a function of force was fitted with a linear equation to determine the k at zero force from the y-intercept.

For force-extension experiments, the stage was repeatedly translated along a horizontal direction at a step size of 5 nm. At each position, after a wait time of 15 msec to a steady-state condition, data were collected for 2 msec at a 500 kHz rate and averaged before moving to the next position. The unfolding rips were identified by sudden drops in force, which correspond to local maxima in the force-extension curve (FEC).

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.