Goldilocks and RNA: where Mg2+ concentration is just right

Abstract Magnesium, the most abundant divalent cation in cells, catalyzes RNA cleavage but also promotes RNA folding. Because folding can protect RNA from cleavage, we predicted a ‘Goldilocks landscape’, with local maximum in RNA lifetime at Mg2+ concentrations required for folding. Here, we use simulation and experiment to discover an innate and sophisticated mechanism of control of RNA lifetime. By simulation we characterized RNA Goldilocks landscapes and their dependence on cleavage and folding parameters. Experiments with yeast tRNAPhe and the Tetrahymena ribozyme P4–P6 domain show that structured RNAs can inhabit Goldilocks peaks. The Goldilocks peaks are tunable by differences in folded and unfolded cleavage rate constants, Mg2+ binding cooperativity, and Mg2+ affinity. Different folding and cleavage parameters produce Goldilocks landscapes with a variety of features. Goldilocks behavior allows ultrafine control of RNA chemical lifetime, whereas non-folding RNAs do not display Goldilocks peaks of protection. In sum, the effects of Mg2+ on RNA persistence are expected to be pleomorphic, both protecting and degrading RNA. In evolutionary context, Goldilocks behavior may have been a selectable trait of RNA in an early Earth environment containing Mg2+ and other metals.

Here we document and characterize Goldilocks behavior of RNAs, with local maxima of chemical lifetimes bounded by conditions of lower lifetimes ( Figure 1A). A Goldilocks landscape is a continuum of conditions, some of which are just right for RNA persistence. We predicted Goldilocks landscapes for RNA because Mg 2+ directly increases RNA cleavage rates by one mechanism and indirectly decreases cleavage rates by a different mechanism. Mg 2+ , the most abundant divalent cation in cells (10), degrades RNA by catalyzing in-line attack of the 2'-oxygen on the backbone phosphate (7,(11)(12)(13)(14)(15)(16). Mg 2+ can also protect against degradation, by facilitating RNA folding (17)(18)(19)(20)(21)(22). The mechanism of protection involves converting conformationally flexible RNA to more structured RNA (23), which is less likely to adopt the geometry required for cleavage (24). A Goldilocks landscape can arise when a given factor acts as a double-edged sword that differentially degrades and protects. Here, we investigated RNA Goldilocks behavior through simulation and experiment.
In simulations, Goldilocks behavior is observed under a broad variety of parameters that influence RNA folding and cleavage. Goldilocks landscapes are influenced by folding mechanism, Mg 2+ -dependency of folding, and folded and unfolded cleavage rate constants. Goldilocks landscapes are not accessible to RNAs that do not fold or unfold.
In our experiments, Goldilocks behavior is observed for well-established model RNAs, including yeast tRNA Phe (25,26) and the Tetrahymena Group I Ribozyme P4-P6 domain (27,28). An experimental comparison of the Goldilocks landscapes of tRNA and P4-P6 RNA indicates that the Goldilocks phenomena is retained even as the landscape is influenced by sequence and chemical modification. In experiments, Goldilocks peaks were observed when RNA is ∼95% folded; a control RNA that does not fold, rU 20 (polyuridylic acid 20-mer), does not display Goldilocks behavior.
Goldilocks behavior of RNA suggests intrinsic sophistication, allowing ultrafine control of chemical lifetime by a variety of inputs (1,29). RNA chemical lifetimes can be tuned by Mg 2+ -mediated shifts into and out of Goldilocks peaks and by remodeling Goldilocks landscapes via sequence and chemical modification. Goldilocks behavior of RNA is consistent with its selection in a primordial world of stringent and conflicting evolutionary demands.

Simulation of RNA lifetime
We mathematically modeled the effect of [Mg 2+ ] on folded and unfolded RNA cleavage rate constants which both contribute to an observed cleavage rate constant k obs . Lifetime is the reciprocal of k obs .
For a two-state model, the fraction folded is f f and the fraction unfolded is f u . We used the Hill equation to describe extent of folding (30), although any model that reasonably describes RNA folding can be used (Eq. 1a and 1b): K D is equivalent to [Mg 2+ ] at the folding midpoint (21) and n is the Hill coefficient, which reflects the cooperativity of folding.
RNA cleavage is a second-order phenomena in which the rate depends on [RNA] and [Mg 2+ ] (7,14). The observed pseudo-first-order rate constant (k obs ) is proportional to the second-order rate constant (k) and [Mg 2+ ] (Eq. 2) (31): 2+ (2) To model RNA lifetime with a two-state folding model, we used two cleavage rate constants: k f for folded RNA and k u for unfolded RNA. Folding offers protection from cleavage, and therefore k f < k u . For an RNA that can occupy two states, k obs is a convolution of cleavage contributions based on fractional occupancies and rate constants for each state (Eq. 3, which is an extension of Eq. 2 with differential cleavage based on folding): RNA lifetime is the reciprocal of the observed cleavage rate constant (k obs ) (Eq. 4) (32): Initial folding parameters, K D = 0.022 mM Mg 2+ and n = 4.1, for native yeast tRNA Phe were obtained from previous work (18). Initial values of relative rate constants were set to k u = 1 and k f = 0.2 t rel −1 [Mg 2+ ] rel −1 , consistent with changes in RNA cleavage rates upon conversion of single strands to duplex (33,34). The software GraphPad Prism 8 was used for simulations.
To model the contribution of folding intermediates to Goldilocks behavior, we modified the two-state model using the approach of Shelton et al. (18) (Supplementary Eq. 1). The transition from unfolded to intermediate is described by the terms K D1 and n 1 , and the transition from intermediate to fully folded is described by K D2 and n 2 . We initialized the three-state equation with K D1 at 1 [Mg 2+ ] rel and k u /k f at 5. We set both n 1 and n 2 to 4.1 and K D2 for the second transition to 2 [Mg 2+ ] rel . k i was varied.
Our simulations and analysis of experimental data assume that RNA folding is fast relative to cleavage. We assume that Mg 2+ binding to RNA is not rate-limiting and that cleavage operates over a fixed folding ensemble. These assumptions are based on experiment and theory. For folding, k = 10 4 to 10 −4 s −1 (35)(36)(37)(38)(39)(40)(41)(42). For cleavage k = 10 −4 to 10 −7 s −1 (7,12,14,34). To either cleave or fold RNA, Mg 2+ must first associate with the RNA. Diffuse association (or diffuse 'binding') of Mg 2+ with RNA has a rate constant of k ≈ 10 10 s −1 , near the rate of diffusion (43), whereas k ≈ 10 5 s −1 for specific binding involving first shell coordination (44,45). Increasing RNA length might attenuate Goldilocks behavior if rate of folding were decreased sufficiently (46) such that folding and cleavage occur on the same timescale.

RNA preparation
Yeast tRNA Phe was purchased from Sigma-Aldrich (R4018). rU 20 was purchased from Integrated DNA Technologies. T7-transcribed, stabilized P4-P6 RNA was produced as in Athavale et al. (47). Background Mg 2+ was ). An RNA that shifts between unfolded and folded states shifts between unfolded and folded lifetimes, to establish a Goldilocks peak. Goldilocks behavior requires conversion from unfolded to folded and a lower cleavage constant of folded vs. unfolded RNA (k f < k u ). (B) A two-state reaction mechanism. U is unfolded RNA and F is folded RNA. removed from the RNAs by dialysis in 180 mM NaCl and 50 mM HEPES buffer pH 7.1 using a 10 kDa MWCO filter.

Circular dichroism
Extent of folding was quantified by CD spectroscopy. A solution of 10 M tRNA, 8 M P4-P6 RNA, or 10 M rU 20 in 180 mM NaCl and 50 mM HEPES buffer pH 7.1 was added to a cuvette with a 0.1 cm path length. Spectra were accumulated on a Jasco J-815 spectropolarimeter with scan rate of 200 nm/min, bandwidth of 3 nm, and data pitch of 0.2 nm from 220 to 350 nm at 65 • C. The RNAs were titrated with small volumes of concentrated known MgCl 2 solutions. CD spectra were blank-subtracted and smoothed with a moving average. For yeast tRNA Phe and P4-P6 RNA, fraction folded was plotted using the theta value at the wavelength that maximized the difference between spectra (260 nm for tRNA and 260.6 nm for P4-P6 RNA). Theta values were baseline-corrected for the effects of dilution, evaporation, and cleavage over the course of CD data acquisition. Yeast tRNA Phe showed little cleavage during CD data acquisition. Some P4-P6 RNA cleavage was observed, consistent with the longer CD acquisition time and greater RNA length (Supp. Figure 1).

In-line cleavage of RNA
We approximated native ionic strength at 180 mM NaCl, which required elevated temperature (65 • C) to observe Mg 2+ -dependent folding (48). Solutions of 10 M yeast tRNA Phe , 3.8 M P4-P6 RNA, or 10 M rU 20 in 180 mM NaCl, 50 mM HEPES buffer pH 7.1, and variable MgCl 2 were incubated at 65 • C for 48 hour. Reaction mixtures were separated by electrophoresis on 6% or 7% urea PAGE gels run at 120V for 1 hour and stained with SYBR Green II. Stained gels were digitized with an Azure 6000 or Typhoon FLA 9500 Imaging System. Band intensities were quantified using AzureSpot Analysis software to determine the amount of total intact RNA present.

Sequencing and fragment analysis of P4-P6 RNA cleavage products
RNAs were analyzed by capillary electrophoresis (SeqStudio, Applied Biosystems) following the manufacturer's protocol. Data were analyzed and peaks were aligned in MATLAB. A final concentration of 10 ng/L uncleaved 'fresh' P4-P6 RNA with 0.4 M of FAMlabeled reverse transcription primer that binds to the 3' end (5'-AGCTTGAACTGCATCCATATCAACA-3', Integrated DNA Technologies), 5 mM DTT, 1× first strand buffer (Invitrogen), 0.5 mM each dNTP, and 2 mM of either ddATP, ddCTP, ddGTP or ddTTP was reverse transcribed by SuperScript III (Invitrogen) to create fragments for sequencing. The P4-P6 RNA cleavage reactions containing 5, 10 and 15 mM Mg 2+ were reverse transcribed in parallel with omission of the ddNTPs for identification of cleavage sites. 1 l of each reverse transcription was combined with 1 l GenefloTM 625 size standard ROX ladder (CHIMERx) and 20 l HiDi (Applied Biosystems). Samples were resolved on a SeqStudio instrument using the Fra-gAnalysis run module.

Simulations reveal Goldilocks behavior
We investigated RNA Goldilocks behavior for RNAs that fold in response to Mg 2+ . The simplest model (Figure 1) allows two states (folded and unfolded) and two cleavage rate constants; k u is the cleavage rate constant of an unfolded RNA and k f is the cleavage rate constant of a folded RNA. The observed rate constant shifts from k u when the RNA is fully unfolded, to a weighted average of k u and k f when the RNA is partially folded, to k f when the RNA is fully folded. This model allows a Goldilocks peak of chemical ]. The Hill coefficient n modulates the sharpness of the Goldilocks peak without a substantial change in its position; a larger n gives a sharper peak. k u modifies the slope of the Goldilocks peak on the low [Mg 2+ ] side, and k f modifies the slope on the high [Mg 2+ ] side. The ratio of k u to k f modulates the intensity of the peak. Goldilocks peaks are absent for RNAs that (i) do not fold, (ii) are always folded, (iii) do not change cleavage rate constant upon folding, or (iv) transition very gradually between differential cleavage realms with varying [Mg 2+ ] (Supplementary Figure 2).
Here we define Goldilocks peak intensity as the ratio of the lifetime at the local maximum to the lifetime at the preceding minimum. Positions of maxima and minima were determined from the simulated lifetime derivative across [Mg 2+ ] and solving for [Mg 2+ ] where slopes are zero. Returning each of those [Mg 2+ ] values back into the original lifetime equation solves for the maximum and minimum used for Goldilocks peak intensity. In experiment, the low and high lifetime datapoints were used directly for the ratio.

Goldilocks behavior in complex models
More realistic RNA folding mechanisms involve intermediate states (18). In a model with intermediates, each intermediate I is associated with specific cleavage rate constant k i (Figure 2A). The simulations reveal that the number, intensities, and proximities of Goldilocks peaks depend on the relative magnitudes of the cleavage rate constants and on locations of the folding transitions in [Mg 2+ ]-space. For a three-state model with two transitions that are fully resolved in [Mg 2+ ]-space, RNA can display two distinct Goldilocks peaks ( Figure 2B). When the transitions overlap in Mg 2+ -space, decreasing k i tends to increase the intensity of the Goldilocks peak at the [Mg 2+ ] where the intermediate population is maximum. Increasing k i depresses lifetime at low [Mg 2+ ] (peak position 1, Figure 2B) and shifts the Goldilocks peak to higher [Mg 2+ ] (peak position 2, Figure  2B). An RNA with a folding intermediate that is cleaved more slowly than the fully folded state (k i < k f ) or cleaved more rapidly than the unfolded state (k i > k u ) can display especially intense Goldilocks peaks. If an intermediate state has intermediate protection, the net Goldilocks peak is less intense than in the absence of an intermediate.

Goldilocks intensity
We show that RNA can inhabit a Goldilocks peak of RNA protection flanked by conditions of lability. The level of protection, given by Goldilocks peak intensity, depends on specific RNA properties. We defined Goldilocks peak intensity as the ratio of the peak maximum to minimum, i.e. the ratio of the protected lifetime to the labile lifetime. Using Goldilocks peak intensity, one can compare and rank various RNAs. We observed, in two-state simulations, that Goldilocks peak intensity increases with increased cooperativity of folding (n) or with increased extent of protection af-forded by folding (decrease of k f relative to k u ) ( Figure 3A). We surveyed values of n and k u /k f to create a Goldilocks intensity map ( Figure 3B). A k u /k f = 3 with an n = 4 produces a modest Goldilocks peak. Increasing either n or k u /k f increases Goldilocks peak intensity. Lowering either n or k u /k f disallows a Goldilocks peak unless there is a compensatory increase in the other parameter. We considered area under the curve as an alternative method for quantifying Goldilocks phenomena. By this method, greater area under the curve would correspond to more extensive protection. Intensity and area together describe the shape of a Goldilocks peak. While intensity indicates the protection an RNA achieves near the peak, area can be concentrated locally or widely dispersed in Mg 2+ -space (Supplementary Figure 3). A local description (intensity) appears to be a more useful comparator.

Experimental observation of a Goldilocks landscape of tRNA
To experimentally investigate Goldilocks behavior, we as-  Figure S4). Random coil yeast tRNA Phe folds to the native L-shaped structure upon addition of Mg 2+ (35,(49)(50)(51)(52)(53)(54). Chemical lifetime of yeast tRNA Phe showed a distinct Goldilocks peak near 3 mM Mg 2+ , where the tRNA was ∼95% folded (Figure 4 and Supplementary Figure 5). The tRNA lifetime was longer at 3 mM Mg 2+ than at either 2.0 mM or at 3.5 mM Mg 2+ . In contrast, rU 20 , which does not fold, did not exhibit Goldilocks behavior (Supplementary Figure 6).
We compared the experimental lifetime data with predictions of our models (Figure 4). The observed yeast tRNA Phe Goldilocks landscape is reasonably fit by a two-state model. The fit and observed Goldilocks peaks are centered at the same [Mg 2+ ]. The cleavage rate constant of the unfolded tRNA is predicted to be 2.7 times greater than that of the folded tRNA. However, the observed Goldilocks peak is sharper and more intense than predicted by the two-state model. A three-state model provides a better fit to the data, especially in the center of the Goldilocks peak. In the threestate model the cleavage rate constant for the unfolded RNA is predicted to be 3.2 times that of the intermediate and 2.2 times that of the fully folded tRNA (k u > k f > k i ). The statistical significance of the improved fit of the threestate versus the two-state lifetime and folding models is indicated by residual errors (Supp. Figure 7). Our results are consistent with previous observations of >2 states for folding of yeast tRNA Phe (18,35,49,55). The fundamental conclusion here, which is the prediction and validation of Goldilocks behavior by RNA, is not dependent on the folding model.
Yeast tRNA Phe , with an intense Goldilocks peak, appears to fold via a protected intermediate. Prior work has shown that yeast tRNA Phe is most compact in intermediate ionic strength (56)(57)(58) suggesting that the folding intermediate is more compact that the native state (however, see reference (59)).  ]. The fraction folded of yeast tRNA Phe (gray circles) was determined by CD. Experimental lifetimes were fit to two-state (black dashed) and three-state (black solid) models. Fraction folded was fit with two-state (gray dashed) and three-state (grey solid) Hill equation models. The three-state model better approximates the lifetime data, with a more intense Goldilocks peak than the two-state model. The tRNA Goldilocks peak is coincident with folding. rU 20 lifetimes decrease monotonically with no Goldilocks peak (dotted black). rU 20   tRNA Phe . P4-P6 RNA has a clear Goldilocks peak that is coincident with RNA folding.
The Goldilocks behavior of P4-P6 RNA is approximated by both the two-state or three-state models. Both models recreate the position in Mg 2+ -space and intensity of the single Goldilocks peak. The peak is produced by the intermediate to folded transition in the three-state model wherein k i is 8 times k f . When constrained to two states, k u is 22   Figure 10) and previous observations that P4-P6 RNA has more than two folding states (60,61).

Comparison of yeast tRNA Phe and P4-P6 RNA
Both RNAs display Goldilocks peaks in experiment and in simulation. Both RNAs are best fit to models with more than two states. For the tRNA the Goldilocks peak is sharp, with a half-height peak width of around 1 mM Mg 2+ . This level of cooperativity is associated with a pronounced Goldilocks peak. tRNA is noted for its high level of structure (62,63). Conversely, P4-P6 RNA has low cooperativity, which is associated with a broad Goldilocks peak, with a peak width at half-height of around 5 mM Mg 2+ . Goldilocks peak intensity for tRNA is greater (2.4) than for P4-P6 RNA (1.5). The fit parameters of both RNAs are shown in Table 1.

Site specificity of cleavage
To characterize Goldilocks phenomena at the level of single nucleotide, we quantified cleavage fragments of P4-P6 RNA using a sequencer. We examined the site-specific extent of cleavage under conditions of the Goldilocks peak (10 mM Mg 2+ ), on the partially folded shoulder of the peak (5 mM), and on the fully folded should of the peak (15 mM) (Supp. Figure S11 and Supplementary Data). The average extent of cleavage is less under the conditions of the Goldilocks peak (0.18) than in regions flanking the peak (0.21 pre-peak and 0.20 post-peak). Subtracting the 10 mM Mg 2+ Goldilocks peak condition as a baseline shows nucleotides that experience more cleavage off the peak than on the peak, with larger cleavage values indicating greater cleavage ( Figure 6A). This information is mapped onto the P4-P6 secondary structure (Figure 6B). The number of detectable cleavage sites and dispersion of cleavage intensities are greatest for the partially folded RNA (at 5 mM Mg 2+ ). These sites overwhelm the few nucleotides that are highly cleaved in the folded state (strong negative cleavage). In the partially folded realm, double-stranded RNA shows more uniform extent of cleavage than unpaired RNA (65). Variability decreases when the RNA fully folds, where the conformations of essentially all nucleotides become fixed. Here, cleavage for 5 mM to 10 mM Mg 2+ is most variable, especially in bulge and loop-forming regions. When increasing from 10 mM to 15 mM Mg 2+ a few cleavage hot spots emerge in unpaired regions, where the RNA is most susceptible to cleavage (24). In P4-P6 RNA, the Mg 2+ -binding metal core does not appear to be folded at 5 mM Mg 2+ , as indicated by extent of cleavage ( Figure 6A, B). The metal core folds and is protected at 10 mM Mg 2+ . This region is a representation of the double-edged sword of Mg 2+ ; by associating with Mg 2+ the RNA is protected from Mg 2+ . As expected (66), by its low reactivity relative to other loop or bulge regions, the GAAA tetraloop appears to be fully folded by 5 mM Mg 2+ .
The sequencing data indicate that even though the low resolution P4-P6 RNA PAGE banding patterns remain reasonably constant, relative extent of cleavage at various sites does in fact change. Although the yeast tRNA Phe banding pattern appeared uniform across [Mg 2+ ] in gels, we assume significant differences in locations of cleavage upon folding. Site-specific analysis of yeast tRNA Phe was not possible with our method because of base modifications.

DISCUSSION
By simulation and experiment we validated a Goldilocks model of RNA. Local maxima in lifetime are flanked by conditions of greater lability. RNAs can resist Mg 2+mediated cleavage when Mg 2+ folds the RNA. Increasing [Mg 2+ ] beyond the folding threshold increases Mg 2+mediated cleavage. We use a Goldilocks framework to explain how lifetime landscapes are modulated by specific characteristics of diverse RNAs. We predict that Goldilocks landscapes are modulated by monovalent cation concentrations, type of divalent cation, RNA sequence and modification, protein and ligand association, and temperature. RNA that does not fold or unfold cannot access Goldilocks protection. Self-cleaving ribozymes are exempt from Goldilocks behavior because their folding increases rates of cleavage.

Diversity of Goldilocks landscapes
RNA response to [Mg 2+ ] is modulated by RNA sequence and chemical modification. The number, intensity, profile, and position in [Mg 2+ ]-space of Goldilocks peaks depends on RNA sequence and chemical modification. The position of a Goldilocks peak in [Mg 2+ ]-space is determined primarily by the affinity of the RNA for Mg 2+ . A smaller K D shifts the Goldilocks peak to lower [Mg 2+ ].
Goldilocks behaviors of RNA should extend beyond Mg 2+ to species such as Fe 2+ , which also promote both RNA folding and cleavage (47,68). Even farther, the general principles of Goldilocks behavior can be applied to any agent that has differential opposing effects on biopolymer lifetime. For example, protein is cleaved by hydrolysis (69). Protein folding decreases rates of hydrolysis (1) and is often promoted by high water activity (70). This model predicts water-defined Goldilocks phenomena for proteins.

Goldilocks behavior in vivo
Our simulations anticipate some RNAs in vivo may inhabit Goldilocks peaks. Free Mg 2+ in vivo is near 1 mM (71,72) (Table S1). An RNA hairpin ribozyme used as a model is mostly folded at 1 mM Mg 2+ under molecular crowding, mimicking the cytosol (73). If the minimal [Mg 2+ ] required for folding coincides with the in vivo [Mg 2+ ], RNA may occupy a Goldilocks peak in cells. More specific in vivo conclusions remain unresolved thus far because of differences in in vitro and in vivo conditions and limitations in manipulating in vivo [Mg 2+ ] (74). Goldilocks landscapes remain to be evaluated in vivo and in the context of protein and ligand binding. For mRNAs, with lifetimes in vivo of minutes (75)(76)(77)(78)(79), spontaneous cleavage might be insignificant. However, long-lived RNAs (tRNAs, ∼9 hours to days (80)(81)(82)(83); and rRNAs, ∼5 hours to days (84)(85)(86)(87)(88)) might be subject to spontaneous cleavage, governed by the Goldilocks phenomena. Goldilocks behavior could explain in part why cells invest in careful maintenance of Mg 2+ homeostasis (89). It seems likely that a narrow range of [Mg 2+ ] prolongs specific RNA lifetimes in vivo.

Goldilocks and origins of life
RNA can transit between dangerous spaces and safe spaces. Finely controlled metastability, with access to Goldilocks peaks of protection, is most likely an imprint of evolutionary processes (90) during the emergence of RNA on the ancient Earth. Sophisticated internal control of lifetime is an indication of selection -of backbone structure, base mod-Nucleic Acids Research, 2023, Vol. 51, No. 8 3537 ifications, and sequence, all of which modulate Goldilocks landscapes.

DATA AVAILABILITY
All data are available through NAR online.

SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.