Structure folding of RNA kissing complexes in salt solutions: predicting 3D structure, stability, and folding pathway.

RNA kissing complexes are essential for genomic RNA dimerization and regulation of gene expression, and their structures and stability are critical to their biological functions. In this work, we used our previously developed coarse-grained model with an implicit structure-based electrostatic potential to predict three-dimensional (3D) structures and stability of RNA kissing complexes in salt solutions. For extensive RNA kissing complexes, our model shows great reliability in predicting 3D structures from their sequences, and our additional predictions indicate that the model can capture the dependence of 3D structures of RNA kissing complexes on monovalent/divalent ion concentrations. Moreover, the comparisons with extensive experimental data show that the model can make reliable predictions on the stability for various RNA kissing complexes over wide ranges of monovalent/divalent ion concentrations. Notably, for RNA kissing complexes, our further analyses show the important contribution of coaxial stacking to the 3D structures and stronger stability than the corresponding kissing-interface duplexes at high salts. Furthermore, our comprehensive analyses for RNA kissing complexes reveal that the thermally folding pathway for a complex sequence is mainly determined by the relative stability of two possible folded states of kissing complex and extended duplex, which can be significantly modulated by its sequence.

INTRODUCTION

RNA complexes are functionally important in gene regulation and expression (Baker et al. 2003; Atkins et al. 2011). An RNA kissing complex is formed through base-pairing between the complementary sequences of two loops (e.g., hairpin loop), and can play critical roles in RNA dimerization and packaging processes in retroviruses (Chen and Varani 2001; Miller and White 2006; Dubois et al. 2018). For example, the dimerization of two RNA copies induced by kissing complex has been proposed to be the essential initiation of HIV-1 viral replication (Laughrea and Jetté 1996; Paillart et al. 1996). Moreover, the formation of RNA kissing complex between untranslated mRNA and antisense RNA executes regulatory functions in both prokaryotes and eukaryotes (Kolb et al. 2001; Gamache et al. 2017). For example, the kissing complex between the ColE1 plasmid-specific transcripts RNA I and antisense regulators RNA II can mediate plasmid replication regulation (Marino et al. 1995; Kolb et al. 2000). In addition, kissing complexes are a kind of fundamental tertiary motifs in RNAs such as ribozymes and riboswitches, which can stabilize RNA spatial tertiary structures and provide ligand-binding/protein-recognition sites for activation (Andersen and Collins 2001; Piganeau et al. 2006; Bouchard and Legault 2014b). Generally, the biological functions of RNA kissing complexes could strongly depend on their three-dimensional (3D) structures and stability (Baker et al. 2003; Atkins et al. 2011). Furthermore, as RNA tertiary structure motifs with high charge density, the formation of RNA kissing complexes induced by inter-strand base-pairing/stacking is very sensitive to solution temperature and ion conditions (Draper et al. 2005; Woodson 2005; Koculi et al. 2007; Chen 2008; Wong and Pollack 2010; Lipfert et al. 2014; Mustoe et al. 2014). Therefore, the knowledge of 3D structures and stability of RNA kissing complexes in ion solutions is highly required to fully understand their biological functions.

Experimental methods such as X-ray crystallography, NMR spectroscopy, and newly developed cryo-electron microscopy have been used to obtain RNA 3D structures (Rose et al. 2017). However, it is still time-consuming and expensive to experimentally derive high-resolution 3D structures of RNAs (Hajdin et al. 2010; Pyle and Schlick 2016; Schlick and Pyle 2017). As a result, there are still limited available structures of RNA kissing complexes deposited in the Protein Data Bank (PDB) (Rose et al. 2017). Simultaneously, as complementary methods, some computational models have been developed in recent years to predict 3D structures of RNAs (Shapiro et al. 2007; Sim et al. 2012; Cragnolini et al. 2015; Somarowthu 2016; Miao and Westhof 2017; Miao et al. 2017; Sun et al. 2017). The fragment-assembly methods such as MC-Fold/MC-Sym pipeline (Parisien and Major 2008), FARNA (Das and Baker 2007), Vfold3D (Cao and Chen 2011a), RNAComposer (Popenda et al. 2012), and 3dRNA (Zhao et al. 2012; Wang et al. 2015, 2017) can efficiently predict 3D structures of RNAs, based on given secondary structures. However, due to the lack of 3D structure templates and inexactly predicted secondary structures, it is still difficult for them to make reliable predictions on the 3D structures of kissing complexes from sequences (Zuker 2003; Dimitrov and Zuker 2004; Rehmsmeier et al. 2004; Andronescu et al. 2005; Bernhart et al. 2008; Huang et al. 2009; Shi et al. 2014b; Yamasaki et al. 2014). Although RNA2D3D can be used to manually construct the 3D structure of a kissing complex (Martinez et al. 2008), its reliability strongly depends on the expert knowledge of users. In parallel, some coarse-grained (CG) models (Hyeon and Thirumalai 2011; Zhang et al. 2012; He et al. 2013; Kim et al. 2014; Bian et al. 2015; Dawson et al. 2016; Li et al. 2016; Boudard et al. 2017; Jain and Schlick 2017; Sieradzan et al. 2017; Uusitalo et al. 2017) such as iFold (Ding et al. 2008), NAST (Jonikas et al. 2009), SimRNA (Boniecki et al. 2015), HiRE-RNA (Cragnolini et al. 2013), RACER (Xia et al. 2013), and oxRNA (Šulc et al. 2014) have been developed to predict RNA 3D structures by involving experimental thermodynamic parameters (Xia et al. 1998) or/and knowledge-based statistical potentials (Tan et al. 2019). However, the structures of kissing complexes have not been involved in these 3D structure prediction models.

Since RNA-mediated functions are often coupled to the RNA structural stability (Chen and Varani 2001; Baker et al. 2003; Atkins et al. 2011), some computational models have been used to predict thermodynamics of RNAs. The VfoldCPX/VfoldThermal with involving thermodynamic parameters and special loop entropy can make reliable predictions on thermodynamics of RNAs including kissing complexes at a secondary structure level (Cao and Chen 2011b; Cao et al. 2014). The structure-based (Gö-like) TIS model can well predict the thermodynamics of small RNAs including pseudoknots (Denesyuk and Thirumalai 2013; Hori et al. 2016; Denesyuk et al. 2018), but the model could not predict 3D structures of RNAs solely from the sequences. Although some other models such as iFold (Ding et al. 2008), HiRE-RNA (Cragnolini et al. 2013), oxRNA (Šulc et al. 2014), and NARES-2P (He et al. 2013) can give melting curves of RNAs, there is still a lack of experimental validation for the predictions and these models are still required to be tested for predicting the stability of RNA kissing complexes.

Very importantly, metal ions, especially divalent ions such as Mg2+, can strongly affect the 3D structures and stability of RNAs due to the polyanionic nature of RNAs (Gregorian and Crothers 1995; Weixlbaumer et al. 2004; Lorenz et al. 2006; Vo et al. 2009; Singh et al. 2011; Sehdev et al. 2012; Kim and Shapiro 2013; Li 2013). Thus, it is also necessary for a predictive model to involve the effect of monovalent/divalent ions to predict the 3D structures and thermodynamics of RNAs in physiological environments containing monovalent and divalent ions (Chen 2008; Sun et al. 2017). However, the effect of ions on the structures and stability of RNAs, especially the effect of divalent ions, is seldom involved in the existing 3D structure prediction models. The effect of ions has not been involved in the fragment-assembly methods due to the limitations of templates at various ion conditions, and the existing physics-based models such as iFold (Ding et al. 2008), HiRE-RNA (Cragnolini et al. 2013), NARES-2P (He et al. 2013), and RACER (Xia et al. 2013) also did not involve the effect of ions on 3D structures, although the effect of ions may be incorporated in such kinds of models through certain implicit or explicit modeling. Especially, the involvement of divalent ions in 3D structure predictions of RNAs is a challenge due to the requirement of structure prediction efficiency and the intrinsic treatment complexity of strong correlations of divalent ions (Tan and Chen 2007; Sun et al. 2017; Xi et al. 2018; Hori et al. 2019). To predict 3D structures and stability of RNAs in ion solutions, we have recently developed a CG model with involving an implicit electrostatic potential and sequence-based thermodynamic parameters (Shi et al. 2014a). The model has been validated by making reliable predictions on 3D structures and stability for RNA hairpins, double-stranded RNAs and RNA pseudoknots, and the involvement of the electrostatic potential in the model can visibly improve the prediction accuracy for the 3D structure and stability of the RNAs in extensive monovalent/divalent ion solutions (Shi et al. 2014a, 2015, 2018; Jin et al. 2018). However, the model should be further validated for predicting 3D structures and stability for the very important tertiary motifs of RNA kissing complexes with tertiary contact of inter-strand base-pairing/stacking in monovalent/divalent ion solutions. Moreover, although the coaxial stacking interactions have been found in the folded structures of RNA kissing complexes (Chen and Varani 2001), it is still lacking a deep understanding on its contribution to the 3D structures and stability of RNA kissing complexes. Furthermore, there can be two possible folded structures—kissing complex and extended duplex for certain complex sequences, while the assembly mechanism for RNA kissing complexes of different sequences is still not fully clear. Therefore, for RNA kissing complexes, to predict the 3D structures/stability in monovalent/divalent salt solutions, and to analyze the contribution of coaxial stacking to 3D structures/stability and thermally folding pathways are still challenging and significant for understanding their structure-related biology.

In this work, we used our previously developed CG model with a new modification on coaxial stacking potential to predict the 3D structures and thermal stability of RNA kissing complexes in monovalent/divalent ion solutions and to analyze the contribution of coaxial stacking to the 3D structure/stability and thermally folding pathway of RNA kissing complexes. As compared with the extensive experimental data, the present model can make reliable predictions on 3D structures and thermal stability for extensive RNA kissing complexes, and can capture the effects of monovalent/divalent salts with high accuracy. Our further analyses also revealed the important contribution of coaxial stacking to the 3D structure and stronger stability of kissing complexes than the kissing-interface duplexes at high salts. Furthermore, our comprehensive analyses on the thermally folding pathway showed that, for a complex sequence, the folding pathway is mainly determined by the relative stability of the two possible folded states of kissing complex and extended duplex, which can be significantly modulated by its sequence.

RESULTS AND DISCUSSION

In this section, we first made predictions on 3D structures for extensive RNA kissing complexes in monovalent/divalent ion solutions. Afterward, we used the present model to predict the stability of RNA kissing complexes with extensive sequences, and further investigated the effect of monovalent/divalent ions on the stability for the HIV DIS type RNA kissing complex. Furthermore, we made comprehensive analyses on the contribution of coaxial stacking to the 3D structures and stability for RNA kissing complexes. Finally, the thermally folding pathways of RNA kissing complexes of typical sequences were analyzed thoroughly to reveal the assembly mechanism of RNA kissing complexes.

Predicting 3D structures of RNA kissing complexes

As described in Materials and Methods, for each RNA kissing complex, two random chains were generated from the sequence and separately placed in a simulation box with the size determined by the strand concentration (e.g., 1 mM). Afterward, the MC simulation with a simulated annealing algorithm was performed from a high temperature to the target temperature (e.g., 25°C). As an example, for the HIV-2 TAR-TAR* kissing complex (PDB id: 1KIS), as shown in Figure 1C,D, the total energy of the simulation system decreases gradually along the annealing process, and simultaneously, the total number of the formed base pairs increases. Since the melting temperatures (∼60°C–80°C) of individual RNA hairpins are generally much higher than those of their kissing complexes (∼30°C–50°C at experimental strand concentrations) (Gregorian and Crothers 1995; Huang et al. 2009; Li and Tinoco 2009; Cao and Chen 2011b), the folding process of the HIV-2 TAR-TAR* kissing complex includes two steps: the formations of two hairpins and native kissing complex, respectively (see Fig. 1C,D).

FIGURE 1.

(A) The representation of all-atomistic (left) and CG (right) for the HIV-2 TAR-TAR* RNA kissing complex (PBD id: 1KIS), and the schematic representation for one nucleotide in the present CG model (middle). (B) A general secondary structure representation of RNA kissing complex including three helices (H1, H2, and H3) and four loops (L1, L2, L3, and L4), and the helix H2 is formed by base-pairing between two individual RNA hairpins. The corresponding secondary structure elements in A and B are in the same colors: H1 (red), H2 (green), H3 (blue), L1 and L3 (orange), and L2 and L4 (cyan). (C,D) The time-evolution of the energy (top panel in C), the number of base pairs (bottom panel in C), and the typical 3D structures (D) during the Monte Carlo simulated annealing process. All the 3D structures are shown with PyMOL (http://www.pymol.org).

To better capture the geometry of helical parts, additional MC simulation of further structure refinement was performed based on the conformation at the end of the annealing process at the target temperature (e.g., 25°C), in which the parameters of bonded potentials were used for the base-pairing regions to more sufficiently depict the helical geometry of the stems (see also Parahelical in Supplemental Table S1; Shi et al. 2014a, 2015, 2018; Jin et al. 2018). The predicted structure ensemble from the structure refinement can be evaluated by their root-mean-square deviation (RMSD) values calculated over all the CG beads to the corresponding atoms in the native structures in PDB (Rose et al. 2017), since there is still no reliable scoring function to identify the nearest-native structure from the predicted CG structure ensemble for the all-atom structure conversion (Li et al. 2018; Tan et al. 2019). Since the CG beads were simplified from the structurally fundamental atom groups (phosphate, sugar, and base groups) and are located at the coordinates of key atoms of the groups (P, C4′ and N1/N9 atoms), the predicted CG structures of RNAs keep the important structure features of the RNAs, and consequently, the RMSDs from predicted CG structures can serve as an examination quantity for prediction reliability. For predicted structures, the mean RMSD (the averaged value over the structure ensemble in the refinement process) and minimum RMSD (corresponding to the structure closest to the native one in the refinement process) were used to evaluate the reliability of the predictions (see more details in Shi et al. 2014a, 2015, 2018; Jin et al. 2018). As shown in Figure 2, for 1KIS, the mean and minimum RMSDs are 4.8 and 2.5 Å, respectively, suggesting a reliable precision of the present model.

FIGURE 2.

(A) The predicted 3D structures (ball-stick) with the mean RMSDs and the minimum RMSDs for three sample RNA kissing complexes (PDB ids: 1KIS [left], 2B8R [middle] and 2MI0 [right]) and their corresponding predicted secondary structures (bottom). The mean (minimum) RMSDs for the three kissing complexes are 4.8 Å (2.5 Å), 5.8 Å (2.4 Å), and 4.9 Å (3.0 Å), respectively. (B) The mean (minimum) RMSDs for 18 RNA kissing complexes (from 32-nt to 78-nt) predicted by the present model with/without involving coaxial stacking potential (solid/dashed lines). The 3D structures are shown with PyMOL (http://www.pymol.org).

3D structure predictions for RNA kissing complexes

According to the above process, we used the present model to predict the 3D structures of 18 RNA kissing complexes, each of which has been determined by experimental methods (e.g., NMR and X-ray crystallography); see Supplemental Table S3. For the RNAs measured by NMR, we made predictions at their respective experimental ion conditions; while for the RNAs from X-ray crystallography, we made predictions at 1 M NaCl (i.e., RNAs are nearly full-neutralized) due to the absence of experimental ion conditions.

As shown in Figure 2, for the 18 RNAs with different sequences as well as different lengths of loops, the overall mean (minimum) RMSD between the predicted structures and the experimental structures is ∼5.4 Å (∼2.7 Å), which suggests that the present model can make reliable predictions for 3D structures of RNA kissing complexes with three coaxial stacked helical stems. For visual convenience, the predicted secondary structures for the 18 RNAs were also shown in Supplemental Figure S2. Generally, in RNA kissing complexes, the junction loops L1 and L3 across the deep major groove of the helix, which can cause a visible bending of the helix toward its major groove (Chen and Varani 2001). We further calculated the stem–stem (H1–H3) angles for the 18 RNA kissing complexes using the DSSR software (Lu et al. 2015). As shown in Figure 3, the calculated stem–stem angles are in good agreement with the experimental data (Chang and Tinoco 1997; Lee and Crothers 1998; Mujeeb et al. 1998; Kim and Tinoco 2000; Ennifar et al. 2001, 2006; Baba et al. 2005; Ennifar and Dumas 2006; Kieken et al. 2006; Lebars et al. 2007, 2008; Van Melckebeke et al. 2008; Bouchard and Legault 2014a; Huang et al. 2017), and the mean deviation of the stem–stem angle is ∼6.5° between our predictions and experimental data, suggesting that the present model is also effective in capturing the helix orientations of RNA kissing complexes. For example, the HIV-2 TAR-TAR* kissing complex (PDB id: 1KIS) with the junction loops of L1 = L3 = 0-nt has a sharp bending (stem–stem angle ∼120°) to the major groove of H2, while the HIV-1Lai genomic DIS kissing complex (PDB id: 2B8R) with longer junction loops of L1 = L3 = 3-nt has a slight bending angle (stem–stem angle ∼160°; see Fig. 3A). The sharp bending for the HIV-2 TAR-TAR* complex is attributed to the strong mechanical pulling effect from the very short loops of L1 = L3 = 0-nt, and the slight bending for the HIV-1Lai DIS complex results from the coaxial stacking between H2 and H1/H3 stems and the slight pulling effect of long loops of L1 = L3 = 3-nt.

FIGURE 3.

(A) The predicted 3D structures (ball-stick) of three RNA kissing complexes with different stem–stem (H1–H3) angles (∼120°, ∼140°, and ∼160°). The three helical stems of H1, H2, and H3 are represented in red, green, and blue, respectively. Junction loops L1 and L3 are represented in orange. The long cyan sticks represent the helical axes of the helical stems of H1 and H3, which are calculated by DSSR software (Lu et al. 2015). (B) The stem–stem angles predicted by the present model with/without involving the coaxial stacking potential for the 18 RNA kissing complexes in comparisons with the experimental data. The 3D structures are shown with PyMOL (http://www.pymol.org).

Effect of coaxial stacking on 3D structures

Since coaxial stacking is important for RNA kissing complexes to keep the quasi-straight helix with proper bending (Marino et al. 1995; Chang and Tinoco 1997; Lee and Crothers 1998; Ennifar et al. 2001, 2006; Baba et al. 2005; Ennifar and Dumas 2006; Lebars et al. 2007, 2008; Van Melckebeke et al. 2008; Li and Tinoco 2009; Mundigala et al. 2014), we further used our model to examine the effect of coaxial stacking on the 3D structures of the 18 RNA kissing complexes through removing the coaxial stacking potential in our predictions.

As shown in Figure 2B and Supplemental Table S3, the overall mean/minimum RMSDs (6.9 Å/3.4 Å) of the structures predicted by the model without the coaxial stacking potential are apparently larger than those (5.4 Å/2.7 Å) from the model with the coaxial stacking potential. Furthermore, our analyses show that there is significant deflection in the stem–stem (H1–H3) angles between the experimental structures and the structures predicted by the model without involving the coaxial stacking potential (see Fig. 3B). Such apparent deviations due to removing coaxial stacking potential suggest that the coaxial stacking makes a significant contribution to the global 3D structures and helix orientation at the junction region of RNA kissing complexes, and the inclusion of the coaxial stacking potential could be important for a model to effectively reproduce the 3D structures and helix orientation of RNA kissing complexes.

Monovalent/divalent salt effects on 3D structures

Since ion solutions may strongly influence 3D structure of RNAs (Gregorian and Crothers 1995; Jossinet et al. 1999; Serra et al. 2002; Weixlbaumer et al. 2004; Lorenz et al. 2006; Vo et al. 2009; Leipply and Draper 2011; Singh et al. 2011; Sehdev et al. 2012; Kim and Shapiro 2013; Li 2013; Wu et al. 2015; Zhang et al. 2017; Chen et al. 2019), we have introduced a structure-based implicit electrostatic potential in the model (Jin et al. 2018). Here, to examine monovalent/divalent salt effects on 3D structures, for the 10 RNA kissing complexes determined by NMR, we made two separate structure predictions using the present model at their corresponding experimental ion conditions and at 1 M NaCl, that is, RNAs are nearly full-neutralized (Shi et al. 2015; Jin et al. 2018). As shown in Table 1, the involvement of the experimental ion conditions causes visible improvements for all the sequences in the 3D structure predictions, and such improvement appears even stronger for longer sequences, for example, for the sequence of 2D1B (78-nt), the mean RMSD of the predicted structures decreases by more than 3 Å. It is reasonable, because RNA would be slightly less compact at lower salt due to the polyanionic nature and such effect would be more pronounced for longer RNAs. For the 10 RNA kissing complexes, the overall improvement is ∼1.7 Å in mean RMSD and ∼1.0 Å in minimum RMSD.

TABLE 1.

Ten RNA kissing complexes for structure predictions at experimental salt conditions

Furthermore, we examined the effects of Na+ and Mg2+ on the 3D structures of RNA kissing complexes of two typical RNA kissing complexes, HIV-2 TAR-TAR* (PDB id: 1KIS) and HIV-1Lai DIS (PDB id: 2F4X), in wide ranges of [Na+] and [Mg2+]. Due to the available experimental measurements at various salts, we examined the dependence of stem–stem (H1–H3) angles of the complexes on monovalent/divalent salts. As shown in Figure 4A,B, the predicted stem–stem angles of the two RNA kissing complexes at different ion conditions agree well with the available experimental data (Chang and Tinoco 1997; Mujeeb et al. 1998; Baba et al. 2005; Kieken et al. 2006; Lebars et al. 2008; Van Melckebeke et al. 2008). It is interesting that the dependences of stem–stem angles on ion concentration are opposite for the two complexes: (i) The HIV DIS type RNA kissing complex (L1 = L3 = 2-nt, L2 = L4 = 1-nt) bends to a stem–stem angle of ∼140° toward a major groove at very low Na+/Mg2+concentrations and becomes quasi-straight with stem–stem angle ∼160° at high Na+/Mg2+ concentrations; (ii) The HIV-2 TAR-TAR* complex (L1 = L2 = L3 = L4 = 0-nt) exhibits bent (stem–stem angle ∼140°) at very low Na+/Mg2+ concentrations and becomes visibly more bent (stem–stem angle ∼100°) at high Na+/Mg2+ concentrations. Our analyses revealed that such interesting opposite salt dependences for the two typical complexes are attributed to the competition between coaxial stacking and mechanical pulling of the loop at the junction regions which are oppositely modulated by ions. Generally, coaxial stacking would make a quasi-continuous helix straight while a loop at junctions would cause it bent (Shi et al. 2015). For the folded complexes, ions can enhance the coaxial stacking interaction due to the neutralization for negative charges on RNAs (Chen et al. 2019), which can overcome the (weak) mechanical pulling by the (long) loops of L1 = L3 = 2-nt and promote the structures of HIV-1Lai DIS to be more straight (see Fig. 4C,D). However, for the HIV-2 TAR-TAR* complex with very short loops of L1 = L3 = 0-nt where the loop pulling effect is strong, the increase of salt concentration could reduce the (inter-helix) electrostatic repulsion induced by the loop-constrained inter-helix bending, which would indirectly strengthen the mechanical pulling effect of the very short loops. Such a salt-enhanced effect can be overwhelming the promotion on the coaxial stacking interaction due to higher salt, causing a stronger bending at high Na+/Mg2+ concentrations (see Fig. 4C,D; Supplemental Table S3).

FIGURE 4.

The predicted stem–stem (H1–H3) angles and coaxial stacking ratio of two typical kissing complexes as functions of [Na+] (A,C) and [Mg2+] (B,D): HIV-2 TAR-TAR* (PDB id: 1KIS) and HIV-1 DIS (PDB id: 2F4X) complexes. The coaxial stacking was determined through the method of RNAView (Yang et al. 2003), and the coaxial stacking ratio was calculated as the fraction of the conformations in which the coaxial stacking between H1 and H2 stems and H2 and H3 stems are both formed. The experimental data are taken from Chang and Tinoco (1997), Mujeeb et al. (1998), Baba et al. (2005), Lebars et al. (2008), and Van Melckebeke et al. (2008). The 3D structures in panels C (left, 2F4X at 10 mM Na+; right, 2F4X at 1000 mM Na+) and D (left, 1KIS at 0.1 mM Mg2+; right, 1KIS at 10 mM Mg2+) are shown with PyMOL (http://www.pymol.org).

Predicting thermal stability of RNA kissing complexes in salt solutions

Beyond 3D structure predictions, we used the present model to predict the thermal stability of extensive kissing complexes and to analyze the Na+/Mg2+ dependence of thermal stability of a typical kissing complex. Furthermore, the contribution of coaxial stacking to the stability of kissing complexes was explicitly analyzed through the present model.

Thermal stability of RNA kissing complexes

Based on the predicted structure ensemble at each temperature for a kissing complex, the fractions of unfolded state at different temperatures can be obtained, and afterward, the melting curve as well as the melting temperature Tm of the RNA kissing complex at an experimental strand concentration (e.g., 0.01 mM) can be calculated through Equations 2 and 3 in Materials and Methods. As an example shown in Figure 5, the predicted melting temperature of HIV-1 DIS type RNA kissing complex (sequence of kissing-interface duplex H2: CUAAAC/GAUUUG) at a strand concentration of 0.01 mM is ∼31.2°C, which agrees well with the experimental value (∼32.0°C) (Weixlbaumer et al. 2004).

FIGURE 5.

(A) The time-evolution of the number of base pairs (vertical red lines) and the average fractions of unfolded state (transverse green lines) for HIV-1 DIS type RNA kissing complex (sequence 4 in Table 3) at different temperatures. (B) The fractions of the unfolded state as functions of temperature for the HIV-1 DIS type RNA kissing complex. (Symbols) The predicted data at different temperatures; blue triangle: Φ(T) at high strand concentration of 1 mM; red circle: corrected value is θh(T); green square: θ(T) at experimental strand concentration (0.01 mM) derived through Equation 2 (Jin et al. 2018). (Lines) The fitted melting curve to the predicted data through Equation 3. More details can be found in the Supplemental Material. (C) The first derivative of θ(T) with the temperature at different strand concentrations.

TABLE 3.

Melting temperature predictions for RNA kissing complexes and corresponding kissing-interface duplexes

For a more extensive examination on thermal stability, we made predictions on the thermal stability for 15 RNA kissing complexes with different sequences or at different strand/ion conditions. As shown in Table 2, for all complexes, the predicted melting temperatures are very close to the experimental measurements (Gregorian and Crothers 1995; Muriaux et al. 1996; Weixlbaumer et al. 2004; Lorenz et al. 2006; Lebars et al. 2007, 2008; Van Melckebeke et al. 2008; Sehdev et al. 2012) with an overall mean deviation of ∼2.4°C. This suggests that the present model with the sequence-dependent base-pairing/stacking potential and structure-based electrostatic potential can make successful predictions on the stability for RNA kissing complexes of extensive sequences in ion solutions.

TABLE 2.

Melting temperature prediction for RNA kissing complexes

Effects of monovalent and divalent ions on stability

Since ions can play important roles in stabilizing RNA tertiary structures, we examined the monovalent/divalent ion effects on the thermal stability for a typical RNA kissing complex, HIV-1 DIS type RNA kissing complex (sequence of kissing-interface duplex H2: CCGACC/GGCUGG). As shown in Figure 6A,B, the predicted melting temperatures as functions of [Na+] and [Mg2+] for the RNA complex are in good accordance with the experimental data (Weixlbaumer et al. 2004), suggesting that the present model can almost give the quantitative predictions on the melting temperatures of RNA kissing complexes over broad ranges of Na+/Mg2+ concentrations. It should be noted that with the combination with the TBI model (Tan and Chen 2006), the present model can capture the competition between monovalent and divalent ions on the RNA kissing complex stability. For example, at very low [Mg2+], Na+ dominates the stability of the HIV-1 DIS type RNA kissing complex, while the increase of [Mg2+] enhances the stability significantly until the strong electrostatic neutralization at very high [Mg2+]; see Figure 6B. This is because the bindings of Na+ and Mg2+ are generally anti-cooperative and Mg2+-binding is more efficient in stabilizing RNA kissing complexes with compact tertiary structures (Draper et al. 2005; Woodson 2005; Koculi et al. 2007; Chen 2008; Wong and Pollack 2010; Lipfert et al. 2014; Mustoe et al. 2014).

FIGURE 6.

The predicted melting temperatures Tm (solid lines) and corresponding experimental data (symbols) (Weixlbaumer et al. 2004) as functions of [Na+] (A) and [Mg2+] (B) for the HIV-1 DIS type kissing complex (sequence of helical stem H2: CCGACC/GGCUGG) and its corresponding kissing-interface duplex. The dashed lines in panels A and B are the predictions from the present model without involving the coaxial stacking potential. The 3D structures in panels A and B demonstrate the kissing-interface duplex (left in A and B) and kissing complex (right in A and B), which are shown with PyMol (http://www.pymol.org).

It is interesting to make comparisons between kissing complex and the corresponding kissing-interface duplex H2 on the salt-dependence of stability. As shown in Figure 6A,B, the kissing complex has higher stability than the duplex H2 at high salts, while it has apparently lower stability than the duplex H2 at low salts. This causes that the Na+/Mg2+-dependence of the kissing complex stability is apparently stronger than that of the duplex. The predicted salt (Na+ and Mg2+)-dependence of Tm calculated as ΔTm/Δlog[Na+] in the range of 50–500 mM and ΔTm/Δlog[Mg2+] in the range of 0.01–1 mM are ∼24°C and ∼15°C for the kissing complex, and are ∼10°C and ∼5°C for the corresponding kissing-interface duplex, respectively. The stronger salt (Na+ and Mg2+)-dependence of Tm and the apparently lower Tm at low salts for the kissing complex than those for the corresponding kissing-interface duplex are because the folding of the kissing complex involves a much higher charge buildup than the duplex H2 (Tan and Chen 2011; Xi et al. 2018).

Contribution of coaxial stacking to stability

For a deep understanding on the relative stability of kissing complex (to the corresponding kissing-interface duplex), we would analyze the contribution of coaxial stacking in the following. For four complexes, we made additional calculations through the present model without involving the coaxial stacking potential. As shown in Table 3, the melting temperatures of the kissing complexes decrease by a mean value of ∼2.4°C due to removing the coaxial stacking, suggesting the important role of coaxial stacking in stabilizing RNA kissing complexes. Furthermore, the melting temperatures from the calculations with removing the coaxial stacking potential are still slightly higher than those of the corresponding kissing-interface duplexes by a mean value of ∼1.3°C. This indicates that coaxial stacking makes the major contribution to the higher stability of kissing complex than the kissing-interface duplex (at high salts). However, the slightly higher melting temperatures without coaxial stacking than the corresponding kissing-interface duplexes may be attributed to that the formation of the kissing-interface duplex H2 from two folded hairpins loses less entropy than the isolated duplex from two unfolded single strands (Tan and Chen 2011). Therefore, the higher stability of kissing complex than the corresponding kissing-interface duplex is mainly attributed to the coaxial stacking between H1–H2 and between H2–H3 stems.

Since salt can modulate the contribution of coaxial stacking to the 3D structures of RNA kissing complexes, it is also necessary to examine the effect of coaxial stacking on the Na+/Mg2+-dependence of the kissing complex stability. We performed additional predictions by our model with removing the coaxial stacking potential for the HIV-1 DIS type RNA kissing complexes at various Na+/Mg2+ concentrations. As shown in Figure 6, at low Na+/Mg2+concentrations, the predicted melting temperatures from the model without coaxial stacking potential are very close to those with involving the coaxial stacking potential. It is because the coaxial stacking can be disrupted by the (inter-helix) electrostatic repulsion due to low Na+/Mg2+ concentrations as discussed above. As Na+/Mg2+ concentrations increase, the calculated melting temperatures for the kissing complex by the present model with/without coaxial stacking both increase rapidly, but the Na+/Mg2+-dependence for the melting temperature without coaxial stacking is not as strong as that with the coaxial stacking. This is reasonable. Higher ion concentrations can enhance the charge neutralization to phosphates and consequently can promote the coaxial stacking in the kissing complex, which causes that the kissing complex has visibly higher stability than that without involving coaxial stacking at high salts.

Thermally folding pathways of RNA kissing complexes

Since the metastable structures of RNAs could be important to their functions (Huthoff and Berkhout 2002; Baker et al. 2003; Atkins et al. 2011; Plumridge et al. 2018; Wang et al. 2018), the thermally folding pathways of RNA kissing complexes would be crucial for understanding their functions (Muriaux et al. 1996; Bernacchi et al. 2007; Mujeeb et al. 2007; Tam et al. 2007; Lee et al. 2009; Cao and Chen 2011b; Livshits et al. 2011). To explore thermally folding pathways of RNA kissing complexes, we made further analyses for three typical sequences of RNA kissing complexes, based on the extensive examination of the folding structures in the annealing simulations for the 18 RNA kissing complex sequences used for the structure prediction. According to the base-pairing states of predicted structures at different temperatures, there are mainly four possible states for typical RNA kissing complexes: dimer of kissing complex (K, at least two continuous kissing base pairs between two folded hairpins), dimer of extended duplex (D), monomer of two hairpins (H, at least one hairpin formed), and monomer of two unfolded single strands (S) (Jin et al. 2018). Although RNA helices can assemble to nanoparticles through blunt-end coaxial stacking (Dibrov et al. 2011), the state of two hairpins with coaxial stacking between blunt ends almost never appears in our simulations. The major reason may be that the stability of coaxial stacking of hairpin blunt ends (e.g., ∼ −2 kcal/mol for two G-C ends) is much weaker than that of a kissing interface with base-pairing/stacking and coaxial stacking (e.g., ∼ −15 kcal/mol for a kissing complex with three kissing base pairs) (Walter and Turner 1994; Walter et al. 1994; Xia et al. 1998).

For three typical RNA kissing complexes (HIV-2 TAR-TAR*, HIV-1Lai DIS and HIV-1Mal DIS), the fractions of four states at different temperatures were calculated as functions of temperatures, and were shown in Figure 7 where the sequences and predicted structures of the four complexes can also be found. For the HIV-2 TAR-TAR* kissing complex, as temperature is decreased from a high temperature (e.g., ∼120°C), the fraction of S state gradually decreases, and simultaneously the fraction of H state increases gradually and reaches its maximum value at ∼50°C, indicating that the two hairpins would fold first. With the further decrease of temperature, the fraction of K state increases accompanied with the decrease of the fraction of H state, suggesting the formation of the kissing complex at low temperatures (e.g., ∼20°C). Thus, the thermally folding pathway of the HIV-2 TAR-TAR* kissing complex follows the order of S → H → K (see the illustration in Fig. 7D), which is in accordance with the previous experiments (Ennifar et al. 2001; Lebars et al. 2007, 2008; Van Melckebeke et al. 2008). In contract, for the HIV-1 DIS kissing complexes (e.g., Lai and Mal), the fractions of K and D states both increase accompanied with the decrease of the fraction of H state when temperature is decreased from ∼70°C (see Fig. 7B,C). This suggests that there are two major folding pathways (S → H → K and S → H → D) and one very minor pathway (S → D) for the two HIV-1 DIS kissing complexes, which has been suggested previously by several experiments (Muriaux et al. 1995b; Bernacchi et al. 2007; Mujeeb et al. 2007; Tam et al. 2007; Lee et al. 2009; Livshits et al. 2011). Such coexistence of both kissing complex and extended duplex has been proposed to be necessary for the dimerization process of the HIV-1 RNA genome (Laughrea and Jetté 1996; Paillart et al. 1996). In addition, at low temperature, the fraction of D state is comparable to that of K state for HIV-1 DISLai complex, while for HIV-1 DISMal complex, the fraction of D state is visibly higher than that of K state. In addition, our predictions showed that the thermally folding pathways of the 18 RNA kissing complex sequences used for the structure prediction can be classified into the three typical types described above.

FIGURE 7.

(A) The state fractions of kissing complex K, hairpins H, and single-stranded chains S as functions of temperature for the thermal folding of HIV-2 TAR-TAR* (PDB id: 1KIS). (B,C) The state fractions of kissing complex K, extended duplex D, hairpins H, and single-stranded chains S as functions of temperature for the thermal folding of HIV-1Lai DIS (B, PDB id: 2B8R) and HIV-1Mal DIS (C, PDB id: 2B8S) complexes. In panels B and C, the dashed lines in cyan and yellow represent the fractions of dimmers (K + D) and monomers (H + S), respectively. (D,E) The schematic diagrams show the structural transitions of HIV-2 TAR-TAR* (D) and HIV-1Lai kissing complex (E) along the thermally folding pathways. The ball-stick graphs represent the 3D structure and the corresponding secondary structure of each state. The left down arrow denotes the direction of decreasing temperatures.

Why do RNA kissing complexes with different sequences have very different folding pathways? To examine what dominates the folding pathway of the complexes, we further calculated the free energies of the states for the three typical complexes at different temperatures and 1 M [Na+] using the Vfold2D algorithm (see Supplemental Fig. S4; Cao and Chen 2011b). For the HIV-2 TAR-TAR* kissing complex, the two possible folded states (K and D states) are sufficiently different in stability, and the K state appears significantly more stable than the D state especially at low temperatures with the relative-free energy ΔΔG = ΔG − ΔG between two states of ∼ −15 kcal/mol at 25°C. This should be the reason that the D state never appears in the folding process of the complex and there is only one folding pathway (S → H → K) for the HIV-2 TAR-TAR* complex. However, for the two HIV-1 DIS kissing complexes, the two (K and D) states are comparable in stability with the relative free energy ΔΔG of 0.4–1.6 kcal/mol, which would lead to the two major folding pathways in parallel. Due to the difference between the sequences in the kissing-interface for the two complexes (HIV-1Lai DIS: CGCGCG and HIV-1Mal DIS: GUGCAC), the relative free energy ΔΔG = ΔG − ΔG = ∼0.4 kcal/mol for the HIV-1Lai DIS complex is much smaller than that (ΔΔG = ∼1.6 kcal/mol) for the HIV-1Mal DIS complex, which would result in the slightly higher fraction of K state for the HIV-1Lai DIS complex and the apparently higher fraction of K state (than D state) for the HIV-1Mal DIS complex at low temperatures. In addition, since D state can be formed directly from S state and has slightly higher stability than K state, there is a minor pathway of S → D. Therefore, the above analyses for the three complexes show that the thermally folding pathway of an RNA complex sequence is mainly dependent on the relative stability of the folded states (K state and D state), which can be significantly modulated by the sequence (Cho et al. 2009; Chen et al. 2018; Shi et al. 2018).

Conclusions

The knowledge of 3D structures and thermal stability of RNA kissing complexes in ion solutions is crucial for understanding their biological functions. In this work, we used our CG model with the implicit structure-based electrostatic potential to predict 3D structures, thermal stability, and folding pathways for RNA kissing complexes in monovalent/divalent ion solutions. The following are the major conclusions:

For extensive RNA kissing complexes in monovalent/divalent ion solutions, we successfully predicted the native-like 3D structures with overall mean (minimum) RMSD of 5.4 Å (2.7 Å) from their sequences, and found that the coaxial stacking interaction and salt condition can play important roles in maintaining their 3D structures.

For various RNA kissing complexes at different strand/ion conditions, we made reliable predictions on their thermal stability with an overall mean deviation of ∼2.4°C in melting temperature from the experimental data, and our analyses showed that ions, especially divalent ions, can promote coaxial stacking and consequently stabilize RNA kissing complexes more strongly at higher salts. Such coaxial stacking makes the major contribution to the higher stability of kissing complex than the corresponding kissing-interface duplex at high salts.

Our comprehensive analyses on thermally folding pathways of typical RNA kissing complexes revealed that the folding pathways as well as final folded structures (kissing complex or extended duplex) are mainly determined by the relative stability between different states, which could be significantly modulated by their sequences.

Although the present model can predict 3D structure/thermal stability for RNA kissing complexes with relatively high accuracy and the role of coaxial stacking/thermally folding pathways was analyzed thoroughly, there are still several limitations that can be improved in the future. First, this work only involved hairpin–hairpin kissing-loop complexes in the predictions, and other types of kissing-loop complexes such as hairpin-bulge/internal kissing-loops were not covered. Fortunately, the present model as well as the analyses can be conveniently extended to examine other types of kissing-loop complex structures beyond hairpin–hairpin kissing-loop complexes. Second, the present model only includes the interactions between canonical Watson–Crick and wobble base pairs and consequently ignores the noncanonical base pairs (e.g., A–G and A–A base pair) which could appear frequently in the complicated RNA tertiary structures (Lemieux and Major 2002). The noncanonical base-pairing may make an additional contribution to the 3D structure and stability of some RNA kissing complexes. At present, our CG model has been successfully used to make structure prediction and stability analyses for small RNAs including hairpins, duplexes, and pseudoknots and kissing complexes in ion solutions. However, it is still difficult for the model to reliably predict structures and stability of large RNAs with more complicated structures (e.g., riboswitches and ribozymes) in ion solutions, since these RNAs generally involve a lot of complex segments/interactions such as multibranched loops, noncanonical base pairs, base triple interactions, base-backbone hydrogen-bond interactions and possible specific ion binding. These complicated structure segments/interactions are required to be taken into account in the further development of the model. Finally, the predicted RNA 3D structures from the model are still at CG level, and it is still very necessary to rebuild all-atom structures based on the predicted CG ones for more strict examination and possible practical use.

Nevertheless, the model can well predict the 3D structures and stabilities of extensive RNA kissing complexes and can well capture the effects of ions on their 3D structure and stability. Furthermore, our comprehensive analyses suggested the important contribution of coaxial stacking to the 3D structure and stability, and revealed the assembly mechanism for the folding of RNA kissing complexes. Therefore, this work can also be very helpful for understanding the roles of ions and coaxial stacking in stabilizing the 3D structure and the structure assembly mechanism for the folding of RNA kissing complexes.

MATERIALS AND METHODS

Coarse-grained RNA structure representation

In our model, to reduce computational complexity and retain the major structure features of RNAs, each nucleotide was simplified into three CG spherical beads: phosphate (P) bead, sugar ring (C) bead, and base (N) bead with van der Waals radii of 1.9, 1.7, and 2.2 Å, respectively (see Fig. 1A; Supplemental Fig. S1; Shi et al. 2014a, 2015, 2018; Jin et al. 2018). The three (P, C, and N) beads were placed at the positions of P, C4′, and N1 for pyrimidine (or N9 for purine) atoms, respectively. Each P bead has a charge of −e on its center (Shi et al. 2014a, 2015, 2018; Jin et al. 2018).

Coarse-grained force field

The total potential energy (Utotal) of an RNA CG conformation is composed of two contributions of bonded potential (Ubond) and nonbonded potential (Unonbond) (Shi et al. 2014a, 2015, 2018; Jin et al. 2018): The bond potential in Equation 1 accounts for the energy associated with virtual bonds in the CG structure, and consists of bond-length energy Ub, bond-angle energy Ua, and dihedral energy Ud. The nonbonded potential in Equation 1 describes excluded volume interaction Uexe, base-pairing interaction Ubp, base stacking interaction Ubs, coaxial stacking interaction Ucs, and effective electrostatic interaction Uel between negatively charged P beads. The coaxial stacking interaction between two discontinued RNA double-helical stems has been found of crucial importance in stabilizing the structures of RNA kissing complexes (Marino et al. 1995; Chang and Tinoco 1997; Lee and Crothers 1998; Ennifar et al. 2001, 2006; Baba et al. 2005; Ennifar and Dumas 2006; Lebars et al. 2007, 2008; Van Melckebeke et al. 2008; Li and Tinoco 2009; Mundigala et al. 2014). As illustrated in Figure 1B, a typical RNA kissing complex contains three helical stems (H1, H2, and H3) and four loops (L1, L2, L3, and L4), in which H2 is formed by the interstrand base-pairing/stacking between two hairpin loops. Generally, H1 and H2 as well as H2 and H3 are coaxially stacked to keep the whole kissing complex as a quasi-continuous helix (Chen and Varani 2001). In the present model, a slight modification was made to more accurately describe the stacking interaction between the two ends of the helical stems at the kissing interface, since our previous examination has shown that the Ucs in our previous model underestimates the coaxial stacking interaction range (Jin et al. 2018). Such slight modification on Ucs has been described in detail in the Supplemental Material. In addition, as described previously (Jin et al. 2018), an implicit structure-based electrostatic potential was involved in capturing the contribution of monovalent/divalent salt to RNA 3D structure and stability, based on the combination of counter ion condensation theory (Manning 1978) and tightly bound ion (TBI) model (Tan and Chen 2006, 2010, 2011; Xi et al. 2018). The involvement of the electrostatic potential would enable our model to predict the effect of monovalent/divalent ions and the effect of the competition between monovalent and divalent ions on the 3D structures and stability of RNAs (Shi et al. 2015; Jin et al. 2018). The detailed description of the potential components in Equation 1 can be found in the Supplemental Material.

Simulation procedure

The Monte Carlo (MC) simulated annealing algorithm was used to predict 3D structures of RNA kissing complexes in ion solutions based on the above described CG force field. Following the previously developed simulation method (Jin et al. 2018), two initial unfolded CG chains randomly generated based on sequences were separately placed in a cubic box, the size of which was determined by the strand concentration. The simulation system was gradually cooled from a high temperature (e.g., 120°C) to the target low temperature (e.g., 25°C). At each temperature, the RNA conformations were renovated by intra-strand pivot moves and inter-strand translation/rotation through the Metropolis algorithm until the system reached enough equilibrium (Jin et al. 2018).

For some sequences of RNA kissing complexes, there may be two possible folded states with comparable folding free energies (Huthoff and Berkhout 2002; Cao and Chen 2011b), for example, the HIV-1 DIS sequence can fold to kissing complex and extended duplex which are coexistent in vivo and vitro (Muriaux et al. 1996; Paillart et al. 1996; Bernacchi et al. 2007; Mujeeb et al. 2007; Tam et al. 2007; Lee et al. 2009; Livshits et al. 2011). In this work, we generally implemented multitrajectory simulations to constitute the conformation ensemble at each temperature. By analyzing structure ensemble in the order of interchain/intrachain base pairs, the 3D structures at different temperatures and the thermodynamic properties for a complex sequence can be obtained (Jin et al. 2018).

Calculations for melting temperatures

The stability of RNA complexes generally depends on the strand concentration due to the contribution of translation entropy of RNA strands (Privalov and Crane-Robinson 2018). In order to improve the simulation efficiency for the systems at each temperature, we generally performed the MC simulations for RNA complexes at a relatively high strand concentration (e.g., 1 mM) instead of low experimental strand concentrations Cs (e.g., 0.01 mM) (Jin et al. 2018). Based on the equilibrium conformations at each temperature T, the fraction θh(T) of unfolded state at can be obtained (Jin et al. 2018). After a correction for finite-box-size effect at the high bulk strand concentration (Ouldridge et al. 2010), the fraction θ(T) of unfolded state at an experimental strand concentration Cs can be calculated by Jin et al. (2018): where . Afterward, the melting temperature Tm for a complex can be obtained by fitting the fraction θ(T) to the two-state model (Shi et al. 2014a), where dT is an adjustable parameter. More details about the calculation of melting temperature were given in the Supplemental Material.

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.