Universal cold RNA phase transitions

Significance Life exists even in the extreme cold, yet we know little about how RNA functions at low temperatures. We have recently found unexpected RNA properties at near-to-zero temperatures, specifically a phase transition to a cold RNA phase that implies a hidden, altered RNA biochemistry. We have investigated cold RNA biochemistry using single-RNA force spectroscopy. At low temperatures, we find that sequence-independent contributions of RNA–water interactions outweigh sequence-dependent base pairing, leading to misfolding of fully complementary hairpins. RNA in the cold may have profound implications for understanding the cold adaptation of RNA biochemistry in present-day psychrophilic biota and may have shaped the evolution of a primordial RNA world.

at low temperatures, showing that fully complementary RNA hairpins unexpectedly misfold below a characteristic glass-like transition temperature T G ∼ 20 • C, adopting a diversity of compact folded structures.This phenomenon is observed in both monovalent and divalent salt conditions, indicating that magnesium-RNA binding is not essential for this to happen.Moreover, misfolding is not observed in DNA down to 5 • C.These facts suggest that the folded RNA arrangements are stabilized by sequence-independent 2 ′ -hydroxyl-water interactions that outweigh sequence-dependent base pairing.Cold RNA misfolding implies that the FEL is rugged with several minima that kinetically trap the RNA upon cooling, a characteristic feature of glassy matter (25).RNA folding in rugged energy landscapes is accompanied by a reduction of RNA's configurational entropy.A quantitative descriptor of this reduction is the folding heat capacity change at constant pressure, ∆C p , directly related to the change in the number of degrees of freedom available to the RNA molecule.Despite its importance, ∆C p measurements in nucleic acids remain challenging (26,27,28).We carry out RNA pulling experiments at low temperatures and show that ∆C p abruptly changes at T G ∼ 20 • C, a manifestation that the ubiquitous non-specific ribose-water interactions overtake the specific Watson-Crick base pairing at sufficiently low temperatures.

RNA misfolds at low temperatures
We used a temperature-jump optical trap (Sec. 1, Methods) to unzip fully complementary Watson-Crick RNA hairpins featuring two 20bp stem sequences (H1 and H2) and loops of different sizes (L = 4, 8, 10, 12 nucleotides) and compositions (poly-A or poly-U) (Sec.2, Methods).Pulling experiments were carried out in the temperature range 7 − 42 • C at 4mM MgCl 2 and 1M NaCl in a 100mM Tris-HCl buffer (pH 8.1). Figure 1A shows the temperaturedependence of the force-distance curves (FDCs) for the dodeca-A (12nt) loop hairpin sequence H1L12A at 4mM magnesium.At and above room temperature (T ≥ 25 • C), H1L12A unfolds at ∼ 20 − 25pN (blue force rips in dashed grey ellipse), and the rupture force distribution is unimodal (Fig. 1B, leftmost top panel at 25 • C), indicating a single folded native state (N).
Misfolding can be characterized by the size of the force rips at the unfolding events, which imply a change in the RNA molecular extension, ∆x.The value of ∆x is obtained as the ratio between the force drop ∆f and the slope k s of the FDC measured at the rupture force f r , ∆x = ∆f /k s (inset of left panel in Fig. 1B). Figure 1B shows ∆x versus f r for all rupture force events in H1L12A at four selected temperatures.Two clouds of points are visible below 25 • C, evidencing two distinct folded states, the native (N, blue) and the misfolded (M, red).
A Bayesian network model (Sec.3, Methods) has been implemented to assign a probability of each data point belonging to N or M (color-graded bar in Fig. 1B).At a given force, the released number of nucleotides for N and M (n N , n M ) is directly proportional to ∆x (Sec.S2, Supp.Info).To derive the values of n N and n M , a model of the elastic response of the singlestranded RNA (ssRNA) is required.We have fitted the datasets (∆x, f r ) for N and M to the worm-like chain (WLC) elastic model (Sec.4, Methods) using the Bayesian network model, finding n N = 52(1) (blue dashed line) and n M = 46(1) (red dashed line) for the number of released nucleotides upon unfolding the N and M structures.Notice that n N matches the total number of nucleotides in H1L12A (40 in the stem plus 12 in the loop), while M features 6nt less than n N .These can be interpreted as remaining unpaired in M or that the 5 ′ − 3 ′ end-to-end distance in M has increased by ∼ 3nm, roughly corresponding to 6nt.

RNA flexibility at low-T promotes misfolding
To characterize the ssRNA elasticity, we show the force-extension curves versus the normalized ssRNA extension per base in Fig. 2A for H1L12A at different temperatures.Upon decreasing T , the range of forces and extensions becomes wider due to the higher unfolding and lower refolding forces.Moreover, a shoulder in the force-extension curve is visible below 32 • C (see also Fig. S3, Supp.Info), indicating the formation of non-specific secondary structures.A similar phenomenon has been observed in ssDNA (35).The force-extension curves (triangles and circles in Fig. 2A) at each temperature were fitted to the WLC model, with persistence length l p and inter-phosphate distance d b as fitting parameters (Fig. 2B and Eq.( 1) in Sec. 4, Methods).Only data above the shoulder has been used to fit the WLC (Sec.S1, Supp.Info).
The values l p and d b show a linear T -dependence (red symbols in Fig. 2C) that has been used for a simultaneous fit of the ssRNA elasticity at all temperatures (blue lines in Fig. 2A).Over the studied temperature range, l p (Fig. 2C, left panel) increases with T by a factor of ∼ 2.5, whereas d b (Fig. 2C, right panel) decreases by only ∼ 20%.The increase of l p with T is an electrostatic effect (34) that facilitates the bending of ssRNA at the lowest temperatures, promoting base contacts and misfolding.

Cold RNA misfolding is a universal sequence-independent phenomenon
The ubiquity of cold misfolding is due to the flexibility of the ssRNA rather than structural features such as stem sequence, loop size, and composition.To demonstrate this, we show results for another five hairpin sequences in Fig. 3A with different stem sequences and loop sizes.To assess the effect of loop size, three hairpins have the same stem as H1L12A but tetra-A, octa-A, and deca-A loops (H1L4A, H1L8A, H1L10A respectively).A fourth hairpin features a dodeca-U loop (H1L12U) to avoid base stacking in the dodeca-A loop of H1L12A.The fifth hairpin, H2L12A, has the same loop as H1L12A but features a different stem.Except for H1L4A, all hairpins misfold below T = 25 • C, as shown by the emergence of unfolding events at forces above 30pN (blue rips in the black dashed ellipses in Fig. 3A) compared to the lower forces of the unfolding native events ∼ 20(grey dashed ellipses).Figure 3B shows the Bayesianclustering classification of the different unfolding trajectories at 7 • C and 25 • C, in line with the results for H1L12A shown in Fig. 1B.The hairpin composition impacts misfolding; while H1L8A, H1L10A, and H1L12A show a single M at 7 • C, H1L12U and H2L12A feature two distinct misfolded states at high (M 1 ) and low (M 2 ) forces (black dashed ellipses for H1L12U and H2L12A in Fig. 3A).
The effect of the loop is to modulate the probability of formation of the native stem relative to other stable conformations.Indeed, H1L4A with a tetraloop has the largest stability among the studied RNAs (36), preventing misfolding down to 7 • C (Fig. 3B).Misfolding prevalence increases with loop size due to the higher number of configurations and low entropic cost of bending the loop upon folding.The ssRNA elastic responses in H1L12A, H1L12U, and H2L12A show a systematic decrease of l p upon lowering T (Fig. S5, Supp.Info) and therefore an enhancement of misfolding due to the large flexibility of the ssRNA.Figure 4A shows the fraction of unfolding events at 7 • C for all hairpin sequences for N (blue), M 1 (red), and M 2 (green).Starting from H1L4A, misfolding frequency increases with loop size, with the second misfolded state (M 2 ) being observed for H1L12U and H2L12A within the limits of our analysis.
Compared to the poly-A loop hairpins (Fig. S6, Supp.Info), the unstacked bases of the poly-U loop in H1L12U confer a larger d b and extension to the ssRNA (red dots in Fig. S5, Supp.Info).
Elastic parameters for the family of dodecaloop hairpins are reported in Table S2, Supp.Info.
The fact that hairpins containing poly-A and poly-U dodecaloops misfold at low temperatures demonstrates that stacking effects in the loop are nonessential to misfolding.
To further demonstrate the universality of cold RNA misfolding, we have pulled the mRNA of bacterial virulence protein CssA from N. meningitidis, an RNA thermometer that changes conformation above 37 • C (37). Figure 4B shows several FDCs measured at 7 • C and 4mM MgCl 2 (inset), evidencing that the mRNA misfolds into two structures (M 1 , red; M 2 , green).

RNA misfolds into stable and compact structures at low temperatures
The Bayesian analysis of the force rips has permitted us to classify the unfolding and refolding trajectories into two sets, N ⇀ ↽ U and M ⇀ ↽ U (Figs. 1B and 3B).We have applied the fluctuation theorem (38,39) to each set of trajectories of H1L12A to determine the free energies of formation of N and M from the irreversible work (W ) measurements at 7 • C (Sec. 5, Methods and Sec.S4, Supp.Info).In Fig. 4C, we show ∆G 0 estimates for N (blue) and M (red), finding ∆G N 0 = 38(9) kcal/mol and ∆G M 0 = 30 (10) kcal/mol in 4mM MgCl 2 (empty boxes).
Within the experimental uncertainties, ∆G 0 for N is higher by ∼ 5 kcal/mol than for M, reflecting the higher stability of Watson-Crick base pairs in N. Notice that the Mfold prediction for N (∆G N 0 = 47 kcal/mol, black dashed line) overestimates ∆G 0 by 10 kcal/mol.
We have also examined the distance between the folded and the transition state x ‡ in H1L12A to quantify the compactness of the folded structure.We have determined x ‡ from the rupture force variance σ 2 using the Bell-Evans model, through the relation S5, Supp.Info).We find that average rupture forces for N and M decrease linearly with T , whereas σ 2 values are T -independent and considerably larger for M, σ 2 M ∼ 50σ 2 N , giving x ‡ M = 0.7(4)nm and x ‡ N = 4.8(6)nm (Fig. S10, Supp.Info).Therefore, M featuring a shorter x ‡ and a more compact structure than N.

The RNA glassy transition
The ubiquity of the cold RNA misfolding phenomenon suggests that RNA experiences a glass transition below a characteristic temperature T G where the FEL develops multiple local minima.
Figure 5A illustrates the effect of cooling on the FEL (40,41).Above 25 • C, the FEL has a unique minimum for the native structure N (red-colored landscape).The projection of the FEL along the molecular extension coordinate shows that N is separated from U by a transition state (TS) (top inset, red line).Upon cooling, the FEL becomes rougher with deeper valleys, promoting misfolding (green and blue colored landscapes).The distance from M to TS is shorter than from N to TS, reflecting that M is a compact structure (bottom inset, green and blue lines).
The glassy transition is accompanied by the sudden increase in the heat capacity change (∆C p ) between N and U below T G ∼ 20 • C for H1L12A and H1L4A.∆C p equals the temperature derivative of the folding enthalpy and entropy, ∆C p = ∂∆H/∂T = T ∂∆S/∂T and can be determined from the slopes of ∆H(T ) and ∆S(T ) (Sec.6, Methods and Sec.S8, Supp. Info).We observe two distinct regimes: above T G (hot, H) and below T G (cold, C).While ∆C H p ∼ 1.5 • 10 3 cal mol −1 K −1 is similar for both H1L12A and H1L4A (parallel red lines in Indeed, the higher flexibility of the U-loop in H1L12U enhances bending fluctuations and misfolding compared to the stacked A-loop in H1L12A (Fig. 4A).Cold RNA misfolding has also been reported in NMR studies of the mRNA thermosensor that regulates the translation of the cold-shock protein CspA (42), aligning with the CssA results of Fig. 4B.Cold RNA misfolding should not be specific to force-pulling but also present in temperature-quenching experiments where the initial high-entropy random coil state further facilitates non-native contacts (43).We foresee that cold RNA misfolding might help to identify misfoldon motifs, contributing to developing rules for tertiary structure prediction (44,45).
Most remarkable is the large ∆C C p values for H1L12A and H1L4A below T G ∼ 20 • C (293K), which are roughly 4-5 times the high-T value above T G , implying a large configurational entropy loss and a rougher FEL at low temperatures.The increase in ∆C p below T G (dashed grey band in Fig. 5B) is reminiscent of the glass transition predicted by statistical mod-els of RNA with quenched disorder (46,47).As ∆C p = C U p − C N p , we attribute this change to the sudden reduction in C N p and the configurational entropy loss upon forming N (25).Both hairpins show maximum stability ∆G 0 at T S ∼ 5 • C (278K) where ∆S 0 vanishes (Fig. 5B).
The value of T S is close to the temperature where water density is maximum (4 • C), with low-T extrapolations predicting cold denaturation at T C ∼ −50 • C (220K) for both sequences.This result agrees with neutron scattering measurements of the temperature at which the RNA vibrational motion arrests, ∼ 220K (6, 24).We hypothesize that T S ∼ 5 • C and mark the onset of universal phase transitions determined by the primary role of ribose-water interactions that are weakly modulated by RNA sequence, a result with implications for RNA condensates (48,49) and RNA catalysis (50).The non-specificity of ribose-water interactions should lead to a much richer ensemble of RNA structures and conformational states and more error-prone RNA replication.Cold RNA could be relevant for extremophilic organisms, such as psychrophiles, which thrive in subzero temperatures (51).Finally, misfolding into compact and kinetically stable structures might help preserve RNAs in confined liquid environments such as porous rocks and interstitial brines in the permafrost of the arctic soil and celestial bodies (52,53).This fact might have conferred an evolutionary advantage to RNA viruses for surviving during long periods (54) with implications on ecosystems due to the ongoing climate change (55).The ubiquitous sequence-independent ribose-water interactions at low temperatures frame a new paradigm for RNA self-assembly and catalysis in the cold.It is expected to impact RNA function profoundly, having potentially accelerated the evolution of a primordial RNA world (56,57).In a pulling experiment, the molecule is tethered between two polystyrene beads through specific interactions with the molecular ends (62).One end is labeled with a digoxigenin (DIG) tail and binds with an anti-DIG coated bead (AD) of radius 3µm.The other end is labeled with biotin (BIO) and binds with a streptavidin-coated bead (SA) of radius 2µm.The SA bead is immobilized by air suction at the tip of a glass micropipette, while the AD bead is optically trapped.The unfolding process is carried out by moving the optical trap between two fixed positions: the molecule starts in the folded state, and the trap-pipette distance (λ) is increased until the hairpin switches to the unfolded conformation.Then, the refolding protocol starts, and λ is decreased until the molecule switches back to the folded state.
The unzipping experiments were performed at two different salt conditions: 4mM MgCl 2 (divalent salt) and 1M NaCl (monovalent salt).Both buffers have been prepared by adding the salt (divalent or monovalent) to a background of 100mM Tris-HCl (pH 8.1) and 0.01% NaN 3 .The NaCl buffer also contains 1mM EDTA.The pulling protocols have been carried out at a constant pulling speed, v = 100nm/s.We sampled 5-6 different molecules for each hairpin and at each temperature, collecting at least ∼ 200 unfolding-folding trajectories per molecule.

RNA synthesis
We synthesized six different RNA molecules made of a 20bp fully complementary Watson-Crick stem, ending with loops of different lengths (L = 4, 8, 10, 12 nucleotides) and compositions (poly-A or poly-U).The hairpins are flanked by long hybrid DNA/RNA handles (∼ 500bp).Further details about the sequences are given Fig. S1 and Table S1, Supp.Info.
The RNA hairpins have been synthesized using the steps in Ref. (63).First, the DNA template (Merck, Township, NJ, USA) of the RNA is inserted into plasmid pBR322 (New England Biolabs, NEB, Ipswich, MA, USA) between the HindIII and EcoRI restriction sites and cloned into the E. coli ultra-competent cells XL10-GOLD (Quickchange II XL site-directed mutagenesis kit).Second, the DNA template is amplified by PCR (KOD polymerase, Merck) using T7 promoters.The RNA is obtained by in-vitro RNA transcription (T7 megascript, Merck) of the DNA containing the RNA sequence flanked by an extra 527 and 599 bases at the 3 ′ -end and 5 ′ -end, respectively, for the hybrid DNA-RNA handles.Finally, labeled biotin (5 ′ -end) and digoxigenin (3 ′ -end) DNA handles, complementary to the RNA handles, are hybridized to get the final construct.

Bayesian clustering
We use a mixture hierarchical Bayesian model (probabilistic graph network) to classify unfolding events as either emanating from a native or a misfolded initial folded state.The model is a soft classifier, giving each trace a probability (score) to belong to a given state.The model is described in Sec.S3, Supp.Info.

ssRNA elastic model
The ssRNA elastic response has been modeled according to the worm-like chain (WLC), which reads where l p is the persistence length, d b is the interphosphate distance and n is the number of bases of the ssRNA.More details on the WLC model and the fitting method used to derive its parameters can be found in Sec.S1, Supp.Info.

Free energy determination
Given a molecular state, ∆G 0 (N ) is the hybridization free energy of the N base pairs of the folded structure when no external force is applied (f = 0).∆G 0 (N ) is obtained from the free energy difference, ∆G λ , between a minimum (λ min ) and a maximum (λ max ) optical-trap positions where the molecule is folded and unfolded, respectively.Thus, one can write where ∆G el (λ) is the elastic energy upon stretching the ssRNA between λ min and λ max .The latter term can be computed by integrating the WLC (Eq.( 1)).As unzipping experiments are performed by controlling the optical-trap position (not the force), this requires inverting Eq.(1) (Sec.S1, Supp.Info.).We used the fluctuation theorem (64) (FT) to extract ∆G(λ) from irreversible work (W ) measurements.This is computed by integrating the FDC between λ min and λ max , W = λmax λ min f dλ (inset in Fig. S8).Given the the forward (P F (W )) and reverse (P R (W )) work distributions, the FT reads where the minus sign of P R (−W ) is due to the fact that W < 0 in the reverse process.When the work distributions cross, i.e.P F (W ) = P R (−W ), Eq.( 3) gives W = ∆G(λ).Let us notice that the FT can only be applied to obtain free-energy differences between states sampled under equilibrium conditions.However, pulling experiments at low T are carried out under partial equilibrium conditions, with misfolding being a kinetic state.It is possible to extend the FT to our case by adding to ∆G(λ) from Eq.( 3) the correction term k B log (ϕ i F /ϕ i R ), where ϕ U(R) is the fraction of forward (reverse) trajectories of state i = N, M (38).Given ∆G(λ), the free energy at zero force, ∆G 0 (N ), is computed from Eq.( 2) by subtracting the energy contributions of stretching the ssRNA, the hybrid DNA/RNA handles, and the bead in the trap.The first two terms are obtained by integrating the WLC in Eq.( 1)), while the latter is modeled as a Hookean spring of energy ∆G b (x) = 1/2k b x 2 , where k b is the stiffness of the optical trap.

Derivation of the heat capacity change
To derive ∆C p , we have measured the enthalpy ∆H 0 and entropy ∆S 0 of N at different T 's for H1L12A and H1L4A.∆S 0 is obtained from the extended form of the Clausius-Clapeyron equation in a force (61), while ∆H 0 = ∆G 0 + T ∆S 0 .Both ∆H 0 and ∆S 0 are temperature dependent, with a finite ∆C p (Sec.S8, Supp.Info).This has been obtained by fitting the Tdependent entropies to the thermodynamic relation ∆S 0 (T ) = ∆S m + ∆C p log (T /T m ), where T m is the reference temperature and ∆S m is the entropy at T = T m .
where ∆f = f F − f U is the force difference upon unzipping between the force in the folded branch F (f F ) and in the unfolded branch U (f U ), k F eff is the effective stiffness in the folded branch, i.e. the slope of the FDC before unfolding, and x d is the diameter of the folded structure projected along the pulling axis.We used the value of x r determined from Eq.(S4) to asses whether the unfolding events experimentally observed originate from the native state (hairpin) or misfolded state.Given the number of nucleotides in the folded structure, n, the following relation holds: Therefore, different states characterized by different n give different (f U , x r (f U )) distributions, as shown in Fig. 1B and 3B, of the main text.
The advantage of Eq.( S5) is two-fold.First, by assuming that the WLC parameters l p , d b (see Sec. S1) are known, the equation can be applied to infer the number of monomers n in the folded structure, as is done in the Bayesian hierarchical model presented in Sec.S3.By determining n, we can also distinguish whether the RNA has folded into the native state or a misfolded state.Second, assuming n to be known, the equation can be used to determine the value of the WLC parameters l p , d p with a least squares fitting method.For H1L12A, the native state has n = 52, permitting us to derive the ssRNA elastic parameters.

S3 Bayesian clustering
To model the unzipping experiments of RNA at low temperatures, we used a Bayesian network approach (mixture hierarchical Bayesian model).This has two advantages.First, using latent state variables in the model gives posterior distributions for the state of each data point, allowing a probabilistic soft clustering of each unfolding trace, i.e. the probability of the RNA being misfolded or native is assigned to each point.Second, using appropriate likelihood functions in the model gives a range of useful physical parameters, such as the mode and scale parameters of the rupture force distribution of each state.These parameters are related to the force average and variance.The latter gives us estimates of the distance to the transition state, x ‡ , (see Sec. S5), and the weight of each state, native and misfolded, in the total population.
We recall that Bayesian network models posit that the prior distributions of the parameters to be estimated are known.Similarly, the likelihood function to observe each data point given these prior parameters is known.The estimation of the model parameters is then obtained by computing the posterior distribution of the model, given by the Bayes theorem: which is often done in practice with Monte Carlo methods.
In RNA unzipping experiments, the model data points are the pairs (f , x r ) that characterize the rupture force and released extension of each unfolding event in the forward unzipping process.The model core idea is that the extension (x r ) released in an unfolding event depends both on the initial folded state of the molecule (through the number of released monomers n, see Sec.S2) and the rupture force, f , since force distributions are state-dependent, Sec.S5.In practice, we use Eq.(S5) assuming that it is valid down to the presence of experimental noise, characterized as the difference between the r.h.s and l.h.s of Eq.(S5) and which we posit to be Laplace distributed around 0 with precision t (we comment further on this distribution choice at the end of the section).Therefore, for each data point (f i , x r,i ), we have: where the dependence on the number of monomers released in an unfolding event, n, is introduced through the use of the so-called latent variable z i , which captures the initial state of a trajectory for each data point i = 1, .., N .Here, we use the shorthand notation z i = 1 for the native state and z i = 2 for the misfolded state.
The second core idea of the Bayesian classification consists of explicitly modeling the state dependency of the rupture force distribution.We posit that the parameters underpinning the rupture force distribution depend on the latent variable z i .More specifically, we assume that rupture forces are Gompertz distributed with mode M and scale 1/s, and we have therefore set M ≡ M z i and s ≡ s z i with different values for the native/misfolded states.The overall likelihood of observing an experimental point (f i , x r,i ) is obtained by putting all these elements together: The first term on the r.h.s is based on Eq.(S7) as described above, with the shortened notation The second term is given by the Gompertz likelihood mentioned above.The third term, p (z i | ⃗ w) represents the likelihood of the latent variable z i given a weight vector ⃗ w = (w 1 , w 2 ) whose components give the average occupancy of each state.We use for z i the standard conjugate pair of a Categorical distribution for the likelihood p combined with a Dirichlet prior for ⃗ w.The formal specification of the model can then be finally completed by defining appropriate priors for the parameters we want to infer, namely and ⃗ w.As already mentioned, we use for ⃗ w a Dirichlet prior and parameterize both t and s 1 ,s 2 with gamma priors.Finally, we take normal priors for n 1 ,n 2 , and Laplace priors for M 1 ,M 2 .The overall prior is then given by where the model hyper-parameters are made explicit with z i = 1, 2. Hyper-parameters are given by the Greek variables µ 1 , µ 2 , ν 1 , ν 2 , μ1 , μ2 , τ1 , τ2 , φ1 , φ2 , ω1 , ω2 , α, ϕ and ω.We emphasize that while different valid choices of priors could be made, all the priors chosen here were purposefully parameterized to be very flat in order to minimally constrain the posterior space.
Given the likelihood function and our choice of priors, we use Bayes theorem to compute the posterior distribution of the parameters we want to infer: p({z i } N i=1 ; n 1,2 ; w 1,2 ; M 1,2 ; s 1,2 ; σ) ∝ Likelihood × P rior , where we defined for convenience σ = 1/t, the inverse of the precision t.The model with all its priors, likelihood, and variables is schematically summarized in Fig. S7, Supp.Info.We used the R library RJAGS (68) to set up the Bayesian network.Posterior distributions were obtained by running at least three Monte Carlo Markov Chains (MCMC) using the RJAGS library, with a burn-in of 1000 iterations, followed by 5000 iterations.We ran the usual convergence and diagnostics test for MCMCs (Gelmann, chain intercorrelation coefficient) and visually inspected the MCMC noise term to confirm that our simulations converged.We always took the median of the posterior distribution of interest for point estimates (e.g., n 1 , n 2 ).We give some additional details on other important aspects of the fitting procedure: • The rupture forces are modeled as Gompertz-distributed.This is usually a good approximation in practice and even true in the BE model, Sec.S5.Each rupture force distribution (misfolded/native) is then parametrized by a different mode M z i and scale parameter 1/s z i .Note, however, that JAGS/RJAGS does not offer a Gompertz likelihood function by default.Therefore, we need to input the likelihood manually, using Eq.(S15) and the zero trick.
• When designing the model, we initially modeled the noise term in the l.h.s of Eq.(S7) with a more standard Gaussian likelihood.However, we quickly realized that some experimental points could feature large deviations between x i and f −1 WLC (f i , n z i ), deviations which skew/bias the model when assuming normality and lead to overall poor convergence performance of the Monte Carlo Markov Chain (MCMC) simulation.Hence, we choose to use a more robust Laplace likelihood, which is more accommodating when a few large outliers are present.This considerably improved the model's stability.Moreover, the Deviance Information Criterion (DIC) score of the model with Laplace likelihood was lower than with a Gaussian model, giving further confidence in this choice.

S4 The H1L12A free-energies
Mechanical work measurements were extracted from unzipping data as described in Sec. 5, Methods.The inset of Fig. S8A illustrates the work measured between two fixed positions (vertical lines) for the unfolding (red) and refolding (blue) FDCs in a given N ⇀ ↽ U cycle.Let P → (W ), P ← (−W ) and ∆G denote the work distributions and free energy difference between N or M and U.In Fig. S8A we show P → (W ) and P ← (−W ) for H1L12A above room temperature, where only N is observed.For the work, W , we have subtracted the energy contributions of stretching the ssRNA, the hybrid DNA/RNA handles, and the bead in the trap (Sec.5, Methods).The free energy at zero force, ∆G 0 , has been obtained by applying statistical approaches such as the Bennett acceptance ratio (BAR) method (39).Additionally, we have also determined ∆G 0 using a diffusive kinetics model for the unfolding reaction, the so-called continuous effective barrier analysis (CEBA) (69) (see Sec. S6).In Fig. S9 (right panel), we show ∆G 0 values obtained with BAR (blue circles) and CEBA (red circles) above 25 • C (Table S3) finding compatible results.The value of ∆G 0 agrees with the Mfold prediction (29) at 37 • C (black triangles).However, a large discrepancy is observed for the enthalpy and entropy values if we assume ∆C p = 0, suggesting a non-zero ∆C p (Sec.S7 and main text).

S5 Bell-Evans model
According to the BE model (70,71), the unfolding and folding kinetic rates between the folded (F ) state and the unfolded (U ) state, can be written as where k 0 is the pre-exponential factor, x ‡ (x U − x ‡ ) are the relative distances between state F (U ) and the transition state.∆G F U is the free energy difference between states F and U at zero force.In a pulling experiment, the force is ramped linearly with time, f = rt, with r the experimental pulling rate.The survival probability in the folded state (F) is zipping experiments (73,74).In CEBA, the effective barrier between the folded (F ) and the unfolded state (U ), B(f ), is derived by imposing the detailed balance condition between the unfolding, k F U (f ), and folding, k F U (f ), kinetic rates (see Eqs.(S11)): where k 0 is the attempt rate, B(f ) is the effective barrier at force f , and is the folding free energy at force f .The latter term is given by where ∆G 0 is the folding free energy between F and U at zero force, and the integral accounts for the free energy change upon stretching the molecule in state U (F ) at force f .We can derive two estimates for B(f ) by computing the logarithms of Eqs.(S16a) and (S16b), which give By imposing the continuity of the two estimations of B(f ) in Eqs.(S18), we can measure the folding free energy at force f , ∆G F U (f ).The free energy of the stretching contribution in ∆G F U (f ) Eq.(S17) can be measured from the unfolded branch.Matching Eqs.S18a and S18b permits us to directly estimate the folding free energy at zero force, ∆G 0 .Details can be found in (69,74).c ",$ ~I(J ",$ , d ",$ ) 8~Dirichlet ⃗ R 9 ",$ ~Laplace( U J ",$ , K ",$ ) # ",$ ~Γ( W F ",$ , X G ",$ )     To be compared, the 100/1 equivalence rule between monovalent and divalent salt concentrations has been applied to the sodium results.The error (in brackets) refers to the last digit.

Fig. 5B )
Fig. 5B), ∆C C p differs: ∆C C p = 8(1)•10 3 cal mol −1 K −1 for H1L12A versus ∆C C p = 5.8(4)•10 3 cal mol −1 K −1 for H1L4A (unparallel blue lines in Fig. 5B) showing the dependence of ∆C C p on loop size at low-T .Despite the different ∆C C p values, ∆S 0 = 0 and stability (∆G 0 ) is maximum at T S = 5(2) • C (Fig. 5B, inset) for both H1L4A and H1L12A (vertical black lines in Fig. 5B, main and inset).Finally, the ∆C C p values predict cold denaturation at the same T C ∼ −50 • C for both sequences.The agreement between the values of T G , T S , and T C suggests that cold RNA phase transitions are sequence-independent, occurring in narrow and well-defined temperature ranges for all RNAs.

Figure 1 :Figure 2 :
Figure 1: cold RNA misfolding.(A) Unfolding (blue) and refolding (red) FDCs from H1L12A unzipping experiments (top-left) at temperatures 7 − 42 • C and 4mM MgCl 2 .The grey-dashed ellipse indicates native (N) unfolding events.Unexpected unfolding events from a misfolded (M) structure appear below 25 • C (black-dashed ellipse) that become more frequent upon lowering T from 17 • C to 7 • C. (B) Classification of N (blue dots) and M (red dots) rupture events at T ≤ 25 • C and WLC fits for each state (dashed lines).The top panels show rupture force distributions at each T .The inset of the leftmost panel shows the parameters of rupture force events (see text).

Figure 3 :Figure 4 :Figure 5 :
Figure 3: Universality of cold RNA misfolding.(A) Unfolding (blue) and refolding (red) FDCs of hairpins H1L4A, H1L8A, H1L10A, H1L12U, and H2L12A at 25 • C and 7 • C. Greydashed ellipses indicate native (N) unfolding events.Except for H1L4A, all RNAs show unfolding events from misfolded (M) structures at 7 • C (black-dashed ellipses).Hairpins H1L12U and H2L12A (featuring a dodeca-U loop and a different stem sequence) show a second misfolded structure at low forces (zoomed insets).Hairpin sequences are shown in each panel.(B) Bayesian classification of the unfolding events for the hairpins in panel (A) at T = 7 • C. The dashed lines are the fits to the WLC for the different states.The top panels show the rupture force distributions.

Figure S4 :B
Figure S4: Multi-T fit of the H1L12A ssRNA elastic response.We simultaneously fit the relation l p = l p (T ) = a 1 T + b 1 , l B = l B (T ) = a 2 T + b 2 on data points at all temperatures.This gives for l p the values a 1 = 0.135 ± 0.006 Å/C, b 1 = 3.5 ± 0.1 Å and for l B the values a 2 = −0.023± 0.002 Å/C, b 2 = 7.37 ± 0.05 Å.The fit is performed over the filled symbols only.The three-dimensional T -force-extension surface is represented in light grey.The black lines plot force-extension cross-sections at a given temperature (red, T = 7 • C; blue, T = 10 • C; green, T = 17 • C; orange, T = 25 • C; yellow, T = 32 • C; brown, T = 36 • C; pink, T = 42 • C).

Figure S6 :
Figure S6: T -dependent ssRNA response for H1L4A, H1L8A, H1L10A, and H1L12A.Results are shown at T = 25 • C (panel A) and 7 • C (panelB).If we plot the force versus the total extension, data for different hairpins do not collapse (insets of A and B).In contrast, upon normalizing the extension per base, the force-extension curves of all hairpins collapse into a master curve (main A and B).Extension is normalized per base by dividing the measured extension by each hairpin's total number of bases.

Figure S7 :Figure S9 :Figure S10 :
Figure S7: The Bayesian classification algorithm.(Left) Specification of the prior (in green) and likelihood (in blue) functions used.The hyper-parameters used are indicated by Greek letters.Γ stands for the gamma distribution, N for the normal distribution.(Right) Probabilistic graph view of the Bayesian network used.Misfolded and native states are represented by the superscript 1, 2 and are encoded in the latent variables z i = 1, 2.We highlight in red the part of the model that depends on z i : the rupture force distribution through its mode M z i and scale 1/s z i parameters and the number of monomers released within an unfolding event n z i .

Table S2 :
Temperature dependence of the persistence length (l p ) [nm] and the interphosphate distance (d b ) [nm] of the different dodecaloop hairpins.The error (in brackets) refers to the last digit.

Table S3 :
∆G 0 [kcal/mol] values of N for H1L12A in the high-T regime measured with BAR and CEBA methods compared to predictions by Mfold.The error (in brackets) refers to the last digit.T [ • C] T [K] State Mg 2+

Table S4 :
∆G 0 [kcal/mol] values of M for H1L12A at T = 7 • C derived with the BAR method.