Improved model to predict the free energy contribution of trinucleotide bulges to RNA duplex stability.

Trinucleotide bulges in RNA commonly occur in nature. Yet, little data exists concerning the thermodynamic parameters of this motif. Algorithms that predict RNA secondary structure from sequence currently attribute a constant free energy value of 3.2 kcal/mol to all trinucleotide bulges, regardless of bulge sequence. To test the accuracy of this model, RNA duplexes that contain frequent naturally occurring trinucleotide bulges were optically melted, and their thermodynamic parameters-enthalpy, entropy, free energy, and melting temperature-were determined. The thermodynamic data were used to derive a new model to predict the free energy contribution of trinucleotide bulges to RNA duplex stability: ΔG°37, trint bulge = ΔG°37, bulge + ΔG°37, AU + ΔG°37, GU. The parameter ΔG°37, bulge is variable depending upon the purine and pyrimidine composition of the bulge, ΔG°37, AU is a 0.49 kcal/mol penalty for an A-U closing pair, and ΔG° 37, GU is a -0.56 kcal/mol bonus for a G-U closing pair. With both closing pair and bulge sequence taken into account, this new model predicts free energy values within 0.30 kcal/mol of the experimental value. The new model can be used by algorithms that predict RNA free energies as well as algorithms that use free energy minimization to predict RNA secondary structure from sequence.

Introduction

RNA has a more expansive functional role than merely serving as an intermediate between DNA and protein. For example, recent studies have investigated the catalytic role of ribozymes,[1] the role of snRNAs in mRNA splicing,[2] and the control of gene expression by riboswitches,[3] siRNAs,[4] and miRNAs.[5] As more functions of RNA are elucidated, greater interest develops in investigating RNA structure. Programs such as RNAstructure,[6]mfold,[7] and the Vienna package[8] predict RNA secondary structure from sequence. These predictive algorithms use a nearest neighbor model based on the thermodynamic parameters of various structural motifs. Current algorithms correctly predict ∼73% of base pairs.[9] Any improvement to the nearest neighbor model could improve secondary structure prediction from sequence. Accurate RNA secondary structure prediction can provide clues about the tertiary structure of a given sequence.

One common naturally occurring secondary structure motif in RNA is a bulge. Bulges consist of one or more consecutive unpaired nucleotides on one side of an RNA duplex. Bulges greater than one nucleotide are assumed to prevent the base pairs adjacent to the bulge from stacking on one another.[9] Bulges serve in a variety of functions including viral replication,[10] ligand binding,[11] feedback regulation,[12] intron splicing,[13] gene expression,[14] as well as the formation of RNA tertiary structures.[15] Thermodynamic data have been collected on all combinations of single nucleotide bulges and adjacent base pairs, and algorithms predicting the free energy contribution of this motif have been developed from these data.[16−18] In contrast to the single nucleotide bulge motif, the thermodynamic parameters of the trinucleotide bulge motif have not been well characterized. This lack of data is unfortunate since trinucleotide bulges are prevalent in nature and have significant biological functions. Trinucleotide bulges have been discovered in all three domains of life. Most notably, two trinucleotide bulges form a bulge-helix-bulge motif that is conserved throughout all archaea and acts as the site of recognition and cleavage by splicing endonucleases.[19] Another notable trinucleotide bulge occurs in the HIV-1 TAR RNA hairpin fragment and is crucial for the interaction between the TAR RNA and Tat protein. Since this interaction is vital to mRNA synthesis, and ultimately, to the replication of the HIV-1 virus, this trinucleotide bulge has been widely studied as a potential drug target.[10]

Despite the natural prevalence of trinucleotide bulges, few studies have characterized the thermodynamic parameters of this motif.[20,21] In fact, the current predictive model is based on free energy values for only six trinucleotide bulges, and it attributes a 3.2 kcal/mol free energy penalty to all trinucleotide bulges, regardless of sequence.[9] The aim of this study was to determine whether a sequence-independent model is the most accurate approach for predicting the free energy contribution of trinucleotide bulges. Reported here are the thermodynamic parameters for 19 trinucleotide bulges, most of which occur frequently in nature. These results can be incorporated into predictive models, allowing experimental rather than predicted values to be used for these naturally occurring bulges. For trinucleotide bulges whose thermodynamic parameters have not yet been characterized, an improved model for predicting their free energy contribution was derived.

Materials and Methods

Compiling and Searching a Database for RNA Trinucleotide Bulges

A previously compiled database containing 1349 RNA secondary structures was used to identify common naturally occurring trinucleotide bulges and closing pair sequences. The database contained structures for 484 tRNAs, 309 5S rRNAs, 223 large subunit rRNAs, 123 small subunit rRNAs, 100 group I introns, 91 signal recognition particles, 16 large subunit rRNAs, and 3 group II introns.[22] This database was searched for trinucleotide bulges and their adjacent base pairs. In this work, G-U pairs were considered to be canonical base pairs.

Design of Sequences

RNA duplexes were composed of one strand that was 11 nucleotides long and a second strand that was 8 nucleotides long; the bulge and adjacent pairs were centered in the duplex with three additional Watson–Crick pairs on each side. In order to prevent end fraying during melting, the terminal pairs of each duplex were G-C pairs. The sequence of the trinucleotide bulge and nearest neighbor pairs was selected based on those that appeared most frequently in the secondary structure database described above. Prior to the experiments, all sequences were checked for possible competing unimolecular or alternative bimolecular folding. It was assumed that all of the duplexes studied here formed a trinucleotide bulge, as this was the fold with the lowest free energy fold. Any possible competing folds were at least ∼2 kcal/mol less stable.

RNA Synthesis and Purification

The RNA oligonucleotides were ordered from Integrated DNA Technologies (Coralville, IA). The synthesis and purification procedures were standard and have been described previously.[23]

Optical Melting Experiments and Thermodynamics

Optical melting experiments were performed in 1 M NaCl, 20 mM sodium cacodylate, and 0.5 mM Na2EDTA (pH 7.0). Optical melting curves were obtained, and the thermodynamic parameters for the entire duplex were determined as previously described.[24,25] Each duplex was melted at least nine times, using a different concentration each time, to ensure that the total oligonucleotide concentration range was at least 50-fold. The free energy contribution of the canonical pairs was determined using the nearest neighbor method (a 0.45 kcal/mol penalty was applied per each A-U pair adjacent to the bulge, as it was considered a terminal A-U pair).[26,27] This calculated free energy value was subtracted from the experimental free energy of the entire duplex, obtained from the 1/Tm vs ln CT plots, in order to determine the free energy contribution of the trinucleotide bulge.[26] A reference duplex was not used because multiple different reference duplexes would have been needed. In addition, the nearest neighbor method is quite accurate for canonical pairs[26,27] and was used since the bulge parameters will ultimately be used in conjunction with nearest neighbor parameters as part of predictive models.

Linear Regression and Trinucleotide Bulge Thermodynamic Parameters

This study investigated the thermodynamic parameters of 16 trinucleotide bulge and closing pair combinations which occurred frequently in the secondary structure database. Data for 3 additional bulges[28] were added to these 16, and data for all 19 bulges were used during data analysis. These three additional bulges were studied previously by our laboratory,[28] but these particular bulge sequences were not found in the secondary structure database described above. Two previous studies did report thermodynamic parameters for some trinucleotide bulges;[20,21] however, due to non-two-state melts[20] and different buffer conditions,[21] these data were not included in this analysis. The experimental free energy contribution of the trinucleotide bulge was used as a constant when doing linear regression using the LINEST function in Microsoft Excel. Multiple parameters were tested, and the parameters that yielded the greatest predictive accuracy accounted for the closing base pair and the number of purine and pyrimidine nucleotides in the bulge.

Results

Database Searching

In the database of 1349 RNA secondary structures, 410 trinucleotide bulges were found. The first data set in Table 1 presents the frequency and percent occurrence when bulge and closing pair sequence are specified. When the data is analyzed in this fashion, 122 unique bulge and closing pair combinations were identified. The top 20 most frequent bulge and corresponding closing pair sequences found in the database are shown in Table 1. These bulges and corresponding closing pair sequences account for ∼65% of all trinucleotide bulges identified in the database search. Many of the bulges identified do not occur frequently, as 62% of the unique bulge and closing pair combinations found each account for <0.25% of the total number of trinucleotide bulges. The most frequent occurring bulge and corresponding closing pair sequence is , which accounts for ∼12% of all of the bulges identified. All of the trinucleotide bulges with this bulge and closing pair sequence were discovered in 16S rRNA, suggesting that this sequence is conserved in this ribosomal subunit. However, the thermodynamic parameters of this bulge and closing pair combination were not characterized in this study due to possible competing conformations. Similarly, other bulge and closing pair combinations found frequently in the secondary structure database were also not studied here due to possible competing structures. For example, the most frequent bulge found in the database, 5′GAG3′, (Table 1) was not studied here. Based on the sequence of the bulge and nearest neighbors, it is possible that a different trinucleotide bulge could form, 5′GGA3′, resulting from the 3′ U in the bottom strand forming a G-U pair with the 3′-most G in the top strand. This ambiguity in the bulge sequence is akin to the ambiguity described for single nucleotide bulges.[16−18] Due to the ambiguity of the bulge sequence, this bulge was not studied; however, this study did investigate 7 of the top 10 most frequent bulge and corresponding closing pair sequences found in the database.

Table 1

Frequency of Occurrence of Trinucleotide Bulge Sequences in a Secondary Structure Databasea

Not all bulges found in the database are shown due to space limitations.

Trinucleotide bulge when bulge and nearest neighbor sequence is specified.

Frequency of occurrence in the database search.

Percent out of 410 trinucleotide bulges, the total number found in the database search.

Trinucleotide bulge when only the bulge sequence is specified.

Closing pairs of trinucleotide bulges are specified.

The second data set in Table 1 presents the frequency and percent occurrence when only the bulge sequence is specified. This analysis identified 49 unique trinucleotide bulge sequences in the database, out of a possible 64 bulge sequence combinations. Table 1 shows the top 20 most frequent trinucleotide bulge sequences found, which account for ∼86% of all of the trinucleotide bulges found in the database. The trinucleotide bulge sequences not included in Table 1 individually account for ≤4% of the total trinucleotide bulges found in the database search. The most frequent bulge sequence identified was 5′GAG3′ which accounts for 15% of all the trinucleotide bulges found in the database search. As mentioned earlier, the prevalence of the 5′GAG3′ bulge sequence in this database can be attributable to the fact that it appears to be conserved in 16S rRNA. This thermodynamic study investigates a total of 12 different bulge sequences (not considering the nearest neighbors), which represents 48% of the different types of trinucleotide bulges identified in the database search.

The third data set in Table 1 lists the frequency and percent occurrence when only the closing pair combination is specified. The database search yielded 31 of the 36 possible types of closing pair combinations, and 95% of all of the trinucleotide bulges are represented by the closing pair combinations in Table 1. The most frequent closing pair found in the database was , accounting for 25% of all closing pairs found. As was true with the predominant bulge sequences, this closing pair pattern appears to be conserved in 16S rRNA, but it also appears in other rRNA subunits as well. This thermodynamic study investigates a total of nine unique closing pair sequences, which account for 74% of closing pair sequences found in the database.

Thermodynamic Parameters

Table 2 shows the thermodynamic parameters for the formation of duplexes containing trinucleotide bulges. These thermodynamic parameters are derived from an analysis of individual melt curves and an analysis of the 1/TM versus log(CT) plots. The duplexes in the table are listed in order of decreasing frequency in the secondary structure database.

Table 2

Thermodynamic Parameters for the Formation of Duplexes Containing Trinucleotide Bulgesa

Measurements were made in 1.0 M NaCl, 10 mM sodium cacodylate, and 0.5 mM Na2EDTA pH 7.0.

Frequency of occurrence obtained from database described in Materials and Methods.

The trinucleotide bulge is identified by bold letters. The nearest neighbors and bulge are set apart for easy identification. The top strand of each duplex is written 5′ to 3′, and each bottom strand is written 3′ to 5′.

All values are calculated at 10–4 Moligomer concentration.

Melt data from ref (28).

Contribution of the Trinucleotide Bulge to Duplex Thermodynamics

The contribution of the trinucleotide bulge to duplex thermodynamics (Table 3) was calculated as described in Materials and Methods. The free energy contribution of the measured trinucleotide bulges has a large variance, ranging from 2.3 to 5.6 kcal/mol. The most destabilizing trinucleotide bulge reported here is , while the least destabilizing trinucleotide bulge is .

Table 3

Contribution of the Trinucleotide Bulge to Duplex Thermodynamics

Frequency of occurrence obtained from database described in Materials and Methods.

The trinucleotide bulge is identified by bold letters. The top strand of each duplex is written 5′ to 3′, and each bottom strand is written 3′ to 5′.

The experimental free energy contribution of the bulge calculated as described in the text.

The free energy prediction made by the sequence-independent model (ref (9)).

The difference between the free energy predicted by the sequence-independent model (ref (9)) and the experimental free energy.

The free energy prediction made by the sequence dependent model proposed here.

The difference between the free energy predicted by the sequence-dependent model and the experimental free energy.

Data from ref (28).

The average (absolute value) deviation.

Free Energy Parameters for Trinucleotide Bulges

Currently, a free energy penalty of 3.2 kcal/mol is attributed to all trinucleotide bulges regardless of bulge and closing pair sequence.[9] In order to improve prediction, multiple other models were tested, and a model that resulted in a low average deviation from the experimental bulge values with minimal parameters was derived. This model is

As shown in Table 4, ΔG°37, bulge is dependent upon the number of purines and pyrimidines in the bulge sequence. Thus, ΔG°37, bulge is a 4.19 kcal/mol penalty for three bulged purines, a 5.06 kcal/mol penalty for three bulged pyrimidines, a 3.36 kcal/mol penalty for two bulged purines and one bulged pyrimidine, and a 3.56 kcal/mol penalty for two bulged pyrimidines and one bulged purine. ΔG°37, AU is a 0.49 kcal/mol penalty for each A-U closing pair, and ΔG°37, GU is a −0.56 kcal/mol bonus for each G-U closing pair. When using this model, the average magnitude of difference between the predicted values and the experimental values was only 0.30 kcal/mol (Table 3).

Table 4

Sequence-Dependent Model for Predicting the Free Energy Contribution of Trinucleotide Bulges

ΔG°37, trint bulge parameter	free energy contribution (kcal/mol)
ΔG°37, bulgea
3 purines	4.19 ± 0.30
3 pyrimidines	5.06 ± 0.45
2 purine, 1 pyrimidine	3.36 ± 0.26
1 purine, 2 pyrimidines	3.56 ± 0.37
ΔG°37, AUb	0.49 ± 0.30
ΔG°37, GUc	–0.56 ± 0.26

Free energy contribution attributed to the three bulged nucleotides. One of the four values will be applied depending on the purine/pyrimidine composition of the bulge.

Free energy penalty penalty applied for each A-U closing pair of a trinucleotide bulge.

Free energy bonus applied for each G-U closing pair of a trinucleotide bulge.

Discussion

The database search yielded 124 out of 2304 possible combinations of bulge and closing pair sequence combinations, 49 out of 64 possible bulge sequence possibilities, and 33 out of a possible 36 closing pair combinations. It is likely that many of the possible sequence combinations that were not found in our database do exist in nature and would have been found in a larger secondary structure database. There may be a structural explanation for why certain sequence combinations occur more frequently than others, but such an explanation would require an extensive structural study.

The two most frequent bulge sequences in the database, 5′GAG3′ and 5′GAU3′, are similar and account for 15% and 12%, respectively, of the bulge sequences found in the database. While 5′GAA3′ and 5′GAC3′ were also found in the database search, it is interesting to note that these sequences, with the same first two nucleotides, only accounted for 4% and 2%, respectively, of bulges in the database. Similarly, the two most frequently occurring closing pair sequences for trinucleotide bulges in this database and , are similar and account for 25% and 20%, respectively, of the trinucleotide bulges found in the database. However, the closing pair sequences that result from flipping the C-G pair, and , occur far less frequently. The former sequence accounts for ∼0.5% of all trinucleotide bulges found in the database, and the latter was not found in the database. There are no obvious explanations for these observed idiosyncracies.

Thermodynamic Contributions of a Trinucleotide Bulge to Duplex Thermodynamics

All of the trinucleotide bulges included in this study destabilize the duplex. This destabilization is expected, as it is presumed that every trinucleotide bulge disrupts what would otherwise be stabilizing stacking interactions between neighboring base pairs.[9] Also, the presence of the bulge may put strain on the adjacent base pairs, and the bend or kink at the site of the bulge may also destabilize the helix. Recall that the free energy contribution of the bulge ranged from 2.3 to 5.6 kcal/mol. This large variance suggests that the thermodynamic contribution attributed to the bulge may be dependent upon bulge and closing pair sequence, and thus, a sequence-dependent model would be preferred over a model that attributes a constant value to all trinucleotide bulges regardless of sequence.

At first glance, the thermodynamic data appear to suggest a relationship between frequency of occurrence in the database and stability of the bulge. Considering the 16 duplexes in this study that were found in the database, the average free energy contribution of the bulge for the three most frequent was 2.8 ± 0.5 kcal/mol, while the average free energy contribution of the bulge for the three least frequent was 3.5 ± 0.5 kcal/mol. The three additional sequences included in this study did not occur in the database search, and the average free energy of the bulge for these sequences was 5.2 ± 0.5 kcal/mol. From this analysis, less stable bulges appear to be occurring less frequently in the database.

However, a more detailed analysis suggests that no relationship between frequency and stability exists. Recall that the new model derived for predicting the free energy of formation of trinucleotide bulges has a parameter that was dependent upon the sequence of the bulge (ΔG°37, bulge). The corresponding values for this parameter suggest that triple pyrimidine bulges and triple purine bulges are more destabilizing than bulges with a mix of purines and pyrimidines. If there is a relationship between bulge frequency and stability, then it would be expected that the most frequently occurring bulges would favor a mix of purines and pyrimidines. From analyzing the top 10 most frequent occurring trinucleotide bulge sequences from the database search (Table 1), bulges with a mix of purines and pyrimidines do not appear to be significantly favored over what would be statistically expected. This supports the conclusion that there is no relationship between frequency and stability. Conversely, there are 15 bulge sequences which did not appear in the database at all, making these sequences the least frequently occurring. If there was a relationship between bulge frequency and stability, then it would be expected that the least frequently occurring bulge sequences would favor all purine and all pyrimidine sequences. The 15 bulge sequences which did not appear in the database do not significantly favor all purine and all pyrimidine sequences over what would be statistically expected, giving further support to the conclusion that there is no relationship between frequency and stability.

Improving the Model Used to Predict Trinucleotide Bulge Thermodynamics

As shown in Table 4, every trinucleotide bulge receives a ΔG°37, bulge penalty depending upon the purine and pyrimidine composition of the bulge. All four possible ΔG°37, bulge penalties have a far greater magnitude than the G-U closing pair bonus. Thus, all trinucleotide bulges will be predicted to be destabilizing, regardless of bulge and closing pair sequence.

Several of the patterns observed here for trinucleotide bulges are consistent with what was reported previously for single nucleotide bulges. A study by Serra et al. (2007) on single nucleotide bulges concluded that the more stable a duplex is without the bulge, the greater the potential for a single nucleotide bulge to destabilize the duplex.[17] This trend was also evident in this study. The average free energy contribution for trinucleotide bulges located within the five most stable duplexes is 4.45 kcal/mol, while the average free energy contribution for trinucleotide bulges located within the five least stable duplexes is 3.32 kcal/mol. This suggests that the free energy contribution of bulges is somewhat dependent upon non-nearest neighbors.

In eq 1, all of the free energy parameters are penalties, except for the free energy contribution of a G-U closing pair, −0.56 kcal/mol. This bonus is very near −0.6 kcal/mol, the bonus attributed by Serra et al. (2007) to single nucleotide bulges with the closing pair sequence .[17] Our study included a total of 10 trinucleotide bulges with one or more G-U closing pairs. Eight of these 10 bulges matched the closing pair orientation that was attributed the −0.6 kcal/mol bonus in the single nucleotide bulge study.[17] Thus, it is not surprising that the bonus discovered in this study is very near the previously discovered bonus for single nucleotide bulges.

The current model used by RNAstructure to predict the free energy contribution of trinucleotide bulges was derived by averaging the measured free energy contributions of six trinucleotide bulges.[6] Unlike that model, the model derived in this study is sequence-dependent and is based upon the measured free energy contribution of 19 trinucleotide bulges, a 3-fold increase in the sample size. On average, the sequence-independent model[9] predicts a free energy that is 0.67 kcal/mol different than the experimental value. On average, the sequence-dependent model proposed here predicts a free energy that is 0.30 kcal/mol different than the experimental value, a more than a 2-fold improvement over the sequence-independent model. Additionally, the standard deviation of the sequence-independent model is 0.92 kcal/mol, which is more than double the standard deviation of the sequence-dependent model, 0.41 kcal/mol (Table 3). Also, the sequence-dependent model works well on bulges that are predicted rather poorly by the sequence-independent model. For example, the prediction made by the sequence-independent model deviates from the measured value by more than 1 kcal/mol for four bulges, and three of these four bulges have measured free energy values which deviate by more than 2 kcal/mol from the predicted value. In contrast, the sequence-dependent model predicts the free energy of all the bulges studied within 1 kcal/mol of the measured free energy.

With the collection of thermodynamic data for 19 trinucleotide bulges, the model for predicting the free energy of formation of trinucleotide bulges was updated. While the current model used by RNAstructure attributes a constant free energy penalty to all trinucleotide bulges, the new model takes both bulge and closing pair sequence into consideration. This adjustment resulted in a sequence-dependent model with half the average error of the sequence-independent model.