Structured RNA molecules form complex 3D architectures stabilized by multiple interactions involving the nucleotide base, sugar and phosphate moieties. A significant percentage of the bases in structured RNA molecules in the Protein Data Bank (PDB) hydrogen-bond with phosphates of other nucleotides. By extracting and superimposing base-phosphate (BPh) interactions from a reduced-redundancy subset of 3D structures from the PDB, we identified recurrent phosphate-binding sites on the RNA bases. Quantum chemical calculations were carried out on model systems representing each BPh interaction. The calculations show that the centers of each cluster obtained from the structure superpositions correspond to energy minima on the potential energy hypersurface. The calculations also show that the most stable phosphate-binding sites occur on the Watson-Crick edge of guanine and the Hoogsteen edge of cytosine. We modified the 'Find RNA 3D' (FR3D) software suite to automatically find and classify BPh interactions. Comparison of the 3D structures of the 16S and 23S rRNAs of Escherichia coli and Thermus thermophilus revealed that most BPh interactions are phylogenetically conserved and they occur primarily in hairpin, internal or junction loops or as part of tertiary interactions. Bases that form BPh interactions, which are conserved in the rRNA 3D structures are also conserved in homologous rRNA sequence alignments.
Structured RNA molecules form complex 3D architectures stabilized by multiple interactions involving the nucleotide base, sugar and phosphate moieties. A significant percentage of the bases in structured RNA molecules in the Protein Data Bank (PDB) hydrogen-bond with phosphates of other nucleotides. By extracting and superimposing base-phosphate (BPh) interactions from a reduced-redundancy subset of 3D structures from the PDB, we identified recurrent phosphate-binding sites on the RNA bases. Quantum chemical calculations were carried out on model systems representing each BPhinteraction. The calculations show that the centers of each cluster obtained from the structure superpositions correspond to energy minima on the potential energy hypersurface. The calculations also show that the most stable phosphate-binding sites occur on the Watson-Crick edge of guanine and the Hoogsteen edge of cytosine. We modified the 'Find RNA 3D' (FR3D) software suite to automatically find and classify BPhinteractions. Comparison of the 3D structures of the 16S and 23S rRNAs of Escherichia coli and Thermus thermophilus revealed that most BPhinteractions are phylogenetically conserved and they occur primarily in hairpin, internal or junction loops or as part of tertiary interactions. Bases that form BPhinteractions, which are conserved in the rRNA 3D structures are also conserved in homologous rRNA sequence alignments.
Structured RNA molecules form compactly folded architectures, superficially resembling folded proteins. This tight packing is remarkable considering that each RNA nucleotide bears a full negative charge. RNA molecules must therefore overcome significant electrostatic self-repulsion to fold compactly. Part of the repulsion is eliminated by mobile as well as structurally bound monovalent and multivalent cations, especially Mg2+ (1). Basic proteins also play important roles in facilitating RNA folding and stabilizing the compact biologically active structures of large RNAs (2). As the negative charge of each nucleotide is concentrated on the anionicoxygen atoms of phosphate groups, they can form very strong hydrogen bonds (H-bonds) with appropriate donors. The RNA bases (A, C, G and U) each present multiple H-bond donors that can interact with phosphate groups and thereby help reduce intra-molecular RNA self-repulsion and stabilize compactly folded, functional structures. In this contribution, we examine the stabilizing interactions between nucleotide bases and the phosphate backbone moieties, which we refer to as ‘base-phosphate’ (‘BPh’) interactions.RNA 3D motifs, ordered arrays of non-Watson–Crick (non-WC) base pairs formed by the nucleotides of hairpin, internal and junction ‘loops' of the secondary structure, mediate most of the tertiary RNA interactions that lead to compact folding (3). Previous work has shown that many recurrent 3D motifs contain specificBPhinteractions, in addition to base pairing (BP) and base stacking (BSt) interactions (3,4). The bases forming BPhinteractions cited in the literature tend to be highly conserved. For instance, many hairpin loops, including the anti-codon and T-loops of tRNAs (5), contain the ‘U-turn' submotif, a sharp bend in the backbone stabilized by a H-bond between the U(N3) imino position and the phosphate of the N+3 base of the loop (6,7). GNRA hairpin loops contain similar motifs, with the imino N1 and/or amino N2 positions of the conserved G forming the BPhinteraction (8). A conserved BPhinteraction also involving G is observed in recurrent sarcin/ricin (S/R) internal loop motifs (9). Conserved GU wobble base pairs were observed to bind anionicoxygen phosphate atoms in the minor groove to facilitate tight packing of helical elements (10). These examples suggest that BPhinteractions may be widespread in RNA structures and may play significant roles in RNA folding.In this article, we empirically identify and categorize BPhinteractions from 3D structures and propose a unified classification and nomenclature for them. We apply ab initio quantum chemical methods to obtain optimized geometries for each empirically identified interaction and to evaluate their intrinsic stabilities. We then use the optimized geometries to define geometriccriteria to identify and classify BPhinteractions in 3D structures. The classification criteria were implemented in the ‘Find RNA 3D'—FR3D (11) software suite (http://rna.bgsu.edu/FR3D). We assess the significance of these interactions in stabilizing RNA 3D structures by identifying the contexts, occurrences and base conservation of these interactions, using 3D structural alignments of the 16S and 23S ribosomal RNA (rRNA) structures of Escherichia coli and Thermus thermophilus (12). After comparing the rRNA 3D structures to identify conserved BPhinteractions, we examine rRNA sequence alignments to determine the base frequencies and substitution patterns at sites where we observe conserved BPhinteractions in the 3D structures. Finally, we compare the contributions of BPhinteractions, WC and non-WC base pairs and BSt interactions to nucleotide conservation, using linear regression.Base-phosphateinteractions may also affect the properties of nucleicacids by more subtle effects, such as the polarization of the nucleobase rings. Such effects are more difficult to capture and their consideration remains for future investigations. Transient BPhinteractions may play a role in RNA catalysis, e.g. by stabilizing transition states or facilitating protonation of leaving groups (13–15).
MATERIALS AND METHODS
Datasets
RNA-containing atomic-resolution X-ray crystal structures of RNA with resolution better than 4.0 Å were downloaded from the Protein Data Bank (PDB, http://www.rcsb.org/pdb/) (16). FR3D was used to add hydrogen atoms to experimental 3D structures by fitting ab initio calculated base geometries containing hydrogens to experimental base coordinates (11). Sequence alignments for 16S and 23S rRNAs were downloaded from the European Ribosomal RNA Database (17), in January 2007 from http://bioinformatics.psb.ugent.be/webtools/rRNA/.
Software
We used MATLAB version 7.5.0.338 (R2007b) for program development and Microsoft Excel for tables. PDB files were analyzed and classified using the ‘FR3D' program (11) available at http://rna.bgsu.edu/FR3D/. To eliminate redundant sequences from sequence alignments, we used the SeqQR program (18) obtained from http://www.scs.uiuc.edu/∼schulten/software/.
Hardware
Data were analyzed using a MacBook (Mac OS X) with an Intel Core Duo running at 2 GHz and with 2 GB of RAM along with a Dell Optiplex GX280 with two Intel Pentium 4 processors running at 3.4 GHz and with 1 GB of RAM.
Selection of non-redundant sequences for BPh analysis
We identified the most complete 16S or 23S sequences for each species in the rRNA sequence alignments obtained from the European Ribosomal RNA Database as previously described (12) and employed the SeqQR program to filter redundant sequences from the alignments, using the following program parameters: (i) sequence identity cutoff of 95%; (ii) gap scale of 0.5 and (iii) norm value of 2 (18). The final sequence alignments comprise 717 16S sequences and 136 23S sequences (12).
Selecting a reduced-redundancy set of PDB files for analysis
The RNA-containing 3D structures deposited with the PDB contain multiple versions of some RNA structures (e.g. 1ffk, 1jj2 and 1s72 are all 3D structures of the Haloarcula marismortui 23S rRNA). We selected a reduced-redundancy set of PDB files for analysis as previously described (12) and these files are listed in ‘Supplemental Materials S4' of Stombaugh et al. (12). Statistical analyses of this dataset in the current paper were restricted to the 159 files having reported resolution 2.5 Å or better.
Initial geometries used in quantum chemical calculations
Representative interactions from 3D structures were identified through the initial search described in the first section under ‘Results', by calculating the H-bond distance between phosphateoxygens and base H-bond donors. A qualified estimate of the H-bonding distance was obtained from the sum of the van der Waals radii of the H-bond donor group and the acceptor atom (19). This means that the distance of H-bond donor and acceptor atoms was selected to be below 3.2 Å for N-H … O and 3.5 Å for C-H … O H-bonds. After identifying the representative interactions, we constructed model geometries of the representative BPhinteractions for quantum mechanical (QM) calculations. Each model consisted of a nucleobase and a phosphate group, with the oxygen atoms representing the O3′ and O5′ atoms of the sugar–phosphate backbone terminated by hydrogen atoms. Thus, the total electrical charge for the studied models was –1. Since the purpose of the QM calculations was to make a general classification of the interaction patterns, we intentionally constructed the simplest possible models. We subsequently performed structural relaxation on these geometries so that the final geometries were not affected by small details of the starting structures.The QM analysis was performed in two steps. First, we optimized the geometries, and then we used the optimized geometries to derive interaction energies to characterize the strength of each BPhinteraction.
Geometry optimizations
Models were optimized using the gradient optimization technique. For each model, we optimized all geometrical parameters (full geometry optimization). The calculations were carried out at the Density Functional Theory (DFT) level of theory using the Gaussian03 program package (20). The density functional wasconstructed using Becke's three-parameter exchange (21) and Lee-Yang-Parr's (22,23) correlation functional (B3LYP). The 6–31G** basis set was used for geometry optimizations, augmented by diffuse d-polarization functions on phosphorus (exponent: 0.0348) as well as on the anionicoxygens of the phosphate group (exponent: 0.0845). This is a standard approach used to describe highly polarizable atoms and anions with electron clouds reaching far from the atomiccenters (24). Our previous studies have shown that the B3LYP/6–31G**-optimized structures compare quite well with reference RIMP2/cc-pVTZ data for H-bonded systems and are entirely sufficient for subsequent high-quality interaction energy calculations (25).Because all systems in the present study have a formal charge equal to –1, all computations were carried out in a dielectriccontinuum using the COSMO solvation model (26–28), setting the dielectricconstant at ϵ = 78.4 (the value for pure water). This method reasonably represents the balance of various stabilizing forces acting in charged intermolecular complexes (29). QM interaction energies are usually calculated in the gas phase (in vacuo). However, for electrically charged systems, the interactions in vacuo are dominated by long-range molecular ion-molecular dipolecontributions. For biomolecules in their normal environments, these long-range interactions are compensated by solvation. Therefore, to make the calculations more relevant to RNA molecules in solution, we included solvation effects in the calculations so as to quench the full expression of the anionic electrostaticinteractions (30). In addition, use of continuum solvent prevents formation of irrelevant H-bonds in the course of geometry optimization and helps us to determine more reliably the relative stability order of the individual types of BPhinteractions.
Interaction energies
Interaction energies were calculated for the optimized geometries. The BPhinteraction energy (ΔEBPh) is used to characterize the strength of the direct intermolecular forces acting between the phosphate and the nucleobase. ΔEBPh is defined in Equation (1).
Here EBPh stands for the electronic energy of the whole system, and EB and EPh are the electronic energies of the isolated subsystems, i.e. nucleobase (B) and phosphate (Ph). The interaction energy is usually derived assuming in vacuo environment, and in this case, (Equation 1) corresponds to the change in enthalpy, without the zero point vibrational energy contribution, of a hypothetical in vacuo dimerization process at a temperature of 0 K. This describes the electronic structure part of the interaction energy, which modern QM methods can derive with very high accuracy (31).In general, the interaction energy (Equation 2) has two components, the Hartree–Fock (HF) term, ΔEHF, and the electron correlation term, ΔEcorr (32). They are calculated consecutively.
The ΔEHF term accounts mainly for the electrostatic (Coulombic) effects. This term also includes a large portion of the polarization (induction) and charge-transfer effects. Besides these attractive contributions, the ΔEHF term also contains the short-range exchange repulsion reflecting the mutual repulsion of the electronicclouds once they start to penetrate each other. The ΔEcorr term includes dispersion attraction and corrections to other terms arising from the correlated motion of the electrons.In this study, we derived the interaction energies in the frame of the COSMO dielectriccontinuum approach that adds the mean-field effect of solvent screening of the electrostatic forces (26–28). The calculations were done at the RIMP2 level of theory using the aug-cc-pVDZ basis set of atomic orbitals. The RIMP2/aug-cc-pVDZ method achieves a close to quantitative accuracy for H-bonded and ionic molecular clusters and is entirely sufficient for the purpose of this study. The RIMP2 method for calculating interaction energies has been validated in previous studies (25,33). These calculations were done as follows: first we computed interaction energies in the gas-phase using the optimized geometries obtained with the COSMO solvation model (see above). We carried out these calculations with and without correction for the ‘basis set superposition error'—BSSE (34). BSSE is a mathematical artifact that arises from using a finite basis set of atomic orbitals. All calculations should be corrected for BSSE, as explained in Sponer et al. (35). The difference of the BSSE-corrected and uncorrected interaction energies gives the BSSE correction (δBSSE). In the next step we computed the BSSE-uncorrected interaction energy at the RIMP2/aug-cc-pVDZ level using the COSMO dielectriccontinuum model (ΔECOSMO). Then, the final BSSE-corrected interaction energy, computed with the COSMO dielectriccontinuum model (ΔECOSMOBSSE), wascalculated using Equation (3).
This represents the desired interaction energy, i.e. ΔEBPh from (1) with the added solvent screening correction. For simplicity, we rename and refer to this quantity as EINT. Thus,
All RIMP2calculations were performed with the Turbomole code (36–38). For the sake of completeness, gas-phase interaction energies derived without the solvent screening are given in the supporting information (Supplementary Data S1).
Modifications to the FR3D program suite
New modules were written for FR3D to find and classify BPhinteractions in RNA 3D structures. These modules have been incorporated in the new release of FR3D, http://rna.bgsu.edu/FR3D/. With this new release, users can launch searches for RNA motifs using as search criteria BPh or near-BPh (nBPh) interactions asconstraints between specified pairs of nucleotides.
Construction of rRNA 3D structural alignments
The FR3D program suite, modified as described above was used to identify and classify BPhinteractions from selected 3D files of the 16S and 23S rRNAs of E. coli—PDBs: 2avy and 2aw4 (39)—and T. thermophilus—PDBs: 1j5e (40) and 2j01 (41). The BPh lists generated by FR3D for homologous rRNA structures were imported into Excel and aligned horizontally to identify conserved BPh positions. The complete BPh 3D alignments for 16S and 23S rRNA are provided in Supplementary Data S2 as MS Excel files to allow the reader to access and manipulate the data as desired. BPhinteractions are listed in the alignments in the order of the number of the nucleotide that H-bonds using its base. Column A in the Excel file provides a running index to restore the original order of the alignment after manipulation. Column B is used to indicate BPhinteractions that are conserved between the E. coli and T. thermophilus structures (an ‘X' indicates conservation). The nucleotide number and base identity for the H-bond donor in each BPhinteraction are listed in columns E and F for the E. coli structure and in N and O for the T. thermophilus structure. BPhacceptor nucleotide data are listed in columns G and H (E. coli) and P and Q (T. thermophilus). Occurrence frequencies from homologous rRNA sequence alignments were computed for each H-bond donor base (columns V-AA) and phosphateacceptor (columns AC-AH) in the BPh alignment using the E. coli structure as the reference. Columns AJ and AK provide comments describing those BPhinteractions that could not be aligned between the E. coli and T. thermophilus 3D structures.
RESULTS
Identification of BPh interactions in 3D structures
We searched the reduced-redundancy set of RNA 3D structures to locate all phosphate groups within H-bonding distance of some nucleobase (12). To visualize the potential interactions, we superposed the interacting bases for each potential BPhinteraction and then plotted the positions of the phosphorous and the phosphateoxygen atoms that are closest to the corresponding H-bond donor atom of the base. Projections of these data onto the plane of the donor base are shown in Figure 1, color-coded according to the distance between the H-bond donor and the nearest phosphateoxygen atom.
Figure 1.
Potential BPh interactions extracted from the reduced-redundancy dataset of RNA 3D structures having resolution 2.5 Å or better. Bases from each instance were superposed and relative locations of phosphorus (black dots) and nearest phosphate oxygen atoms (colored dots) are plotted. Dots representing oxygen atoms are colored toindicate distance to the nearest base H-bond donor. These range from 2.5 Å or less (red) to 4.5 Å (dark blue). Only those phosphate oxygens within 4.5 Å of a base H-bond donor and with bond angle >110° are included.
Potential BPhinteractions extracted from the reduced-redundancy dataset of RNA 3D structures having resolution 2.5 Å or better. Bases from each instance were superposed and relative locations of phosphorus (black dots) and nearest phosphateoxygen atoms (colored dots) are plotted. Dots representing oxygen atoms are colored toindicate distance to the nearest base H-bond donor. These range from 2.5 Å or less (red) to 4.5 Å (dark blue). Only those phosphateoxygens within 4.5 Å of a base H-bond donor and with bond angle >110° are included.The data in Figure 1 show that all four RNA bases can form more than one BPhinteraction. Most interactions involve H-bonding to a single phosphateoxygen, but the C and G bases have interactions in which two oxygen atoms of the same phosphate simultaneously H-bond to distinct base H-bond donors. For C, these simultaneous interactions involve N4 and C5 on the Hoogsteen edge, and for G, they involve N1 and N2 on the WC edge. These H-bond donors can also make non-simultaneous interactions, bringing the total number of distinct potential BPhinteractions to 17. Rather than simply labeling them 1–17, we categorize the interactions as indicated in Figure 1, using labels 0BPh through 9BPh, so that the label indicates, in a consistent way for all four bases, the location of the base H-bond donor group relative to the glycosidic bond and the base edges. The 4BPh and 8BPh labels are used for the simultaneous interactions made by G and C, respectively. There are a few instances in which a single phosphateoxygen H-bonds simultaneously to G(N1) and G(N2). These instances involving bifurcated H-bonding are also classified as 4BPh.Idealized BPhinteractions and the corresponding category labels are shown in Figure 2. While each base presents distinct combinations of H-bond donors and acceptors, certain bases have electropositive functional groups at equivalent positions and therefore can form equivalent BPhinteractions, as shown in Figure 2. BPhinteractions which are formed by more than one base include 0BPh, 5BPh, 6BPh, 7BPh and 9BPh. All bases can form 0BPhinteractions using the purine H8 or pyrimidine H6 proton as H-bond donor. Both A (upper-left panel) and C (upper-right panel) can form 6BPh and 7BPhinteractions using equivalent exocyclic amino groups at the A(N6) and C(N4) positions. Similarly, G (lower-left panel) and U (lower-right panel) can form equivalent 5BPhinteractions using the G(N1) or U(N3) imino groups on their WC edges. The 3BPh, 4BPh and 5BPhinteractions all occur on the WC edge of G and are possibly interchangeable during conformational changes or thermal fluctuations of RNA molecules. Likewise, the 7BPh, 8BPh and 9BPhinteractions all occur on the Hoogsteen edge of C and may also be interchangeable. The other interactions are specific to individual bases: 1BPh, 3BPh and 4BPh are specific to G, 2BPh is specific to A and 8BPh is specific to C. We refer to BPhinteractions that occur on the same edge of the same base as neighboring interactions. When comparing homologous 3D structures, we expect to observe the same or possibly neighboring BPhinteractions at homologous locations where the base forming the BPhinteraction is conserved. When a base substitution occurs, we expect that for the BPhinteraction to be maintained without a large conformational change or disruption of the structure, an equivalent BPhinteraction must form. This suggests BPhinteractions that are crucial for RNA folding, stability and function constrain the sequence variation of homologous RNA molecules. This idea will be examined in detail later in the paper.
Figure 2.
Proposed nomenclature for BPh interactions and superpositions of idealized BPh interactions observed in RNA 3D crystal structures for each base. H-bonds are indicated with dashed lines. BPh categories are numbered 0–9, starting at the H6 (pyrimidine) or H8 (purine) base positions. BPh interactions that involve equivalent functional groups on different bases are grouped together, i.e. 0BPh (A, C, G, U), 5BPh (G, U), 6BPh (A, C), 7BPh (A, C) and 9BPh (C, U).
Proposed nomenclature for BPhinteractions and superpositions of idealized BPhinteractions observed in RNA 3D crystal structures for each base. H-bonds are indicated with dashed lines. BPhcategories are numbered 0–9, starting at the H6 (pyrimidine) or H8 (purine) base positions. BPhinteractions that involve equivalent functional groups on different bases are grouped together, i.e. 0BPh (A, C, G, U), 5BPh (G, U), 6BPh (A, C), 7BPh (A, C) and 9BPh (C, U).
Quantum chemical calculation of optimized geometries and interaction energies for BPh interactions
Quantum calculations were carried out for each of the 17 distinct potential BPhinteractions identified above. As described in the ‘Materials and Methods’ section, model geometries were constructed for each interaction formed by each base using the empirical structures as a guide. The geometries were optimized quantum mechanically and then interaction energies were calculated. PDB files with the optimized geometries can be viewed and downloaded at http://rna.bgsu.edu/FR3D/BasePhosphates. We present the calculated interaction energies and optimal distances in Table 1 and summarize the results base by base in this section.
Table 1.
Energies and H-bond distances for BPh interactions (0BPh-9BPh) by nucleotide (A, C, G, U)
A
C
G
U
EINT (kcal/mol)
H-Bond Distance (Å)
EINT (kcal/mol)
H-Bond Distance (Å)
EINT (kcal/mol)
H-Bond Distance (Å)
EINT (kcal/mol)
H-Bond Distance (Å)
1BPh
−4.3
2.85 (N2)
2BPh
−0.1
3.41 (C2)
3BPh
−5.4
2.80 (N2)
4BPh
−10.1
2.85 (N2)(1)
2.89 (N1)(2)
5BPh
−4.0
2.92 (N1)
−4.2
2.81 (N3)
6BPh
−3.1
2.90 (N6)
−3.5
2.89 (N4)
7BPh
−2.8
2.77 (N6)
−4.8
2.85 (N4)
8BPh
−5.6
2.86 (N4)(1)
3.53 (C5)(2)
9BPh
−0.6
3.35 (C5)
−0.6
3.34 (C5)
0BPh
−1.1
3.22 (C8)
−1.0
3.21 (C6)
−1.1
3.25 (C8)
−1.1
3.23 (C6)
Each nucleotide is represented by two columns:
(1)The calculated interaction energy (kcal/mol);
2Distance (in Å) from H-bond donor to acceptor. The H-bond donor site is given in parentheses.
Energies and H-bond distances for BPhinteractions (0BPh-9BPh) by nucleotide (A, C, G, U)Each nucleotide is represented by two columns:(1)The calculated interaction energy (kcal/mol);2Distance (in Å) from H-bond donor to acceptor. The H-bond donor site is given in parentheses.
Adenine
Adeninecan H-bond to phosphates with the exocyclic amino group (N6) as well as the polarized C2-H2 and C8-H8 groups (Table 1, columns 1–2, and Supplementary Data S3). H-bonding with the less polarized C-H groups is markedly weaker (interaction energies are –0.1 and –1.1 kcal/mol for C2-H2 and C8-H8, respectively) than with the N6-amino group (computed interaction energies are –3.1 kcal/mol for 1H6 and –2.8 kcal/mol for 2H6).
Cytosine
Cytosine forms five different binding patterns with phosphate, utilizing its N4 exocyclic amino group as well as the C5-H5 and C6-H6 groups (Table 1, columns 3–4, and Supplementary Data S4). The strongest binding mode for cytosine occurs on the Hoogsteen edge when the H-bond from the N4 amino group to the phosphate is supplemented by an H-bond from C5-H5. We have found three different types of stable binding that differ as to which phosphateoxygens are involved: (i) N-H binds to an anionicoxygen, while C-H binds to O3′ or O5′; (ii) N-H and C-H both bind to anionicoxygens and (iii) N-H and C-H both bind to O3′/O5′. While the first two variants possess almost identical stabilities (computed interaction energies are –5.4 and −5.6 kcal/mol, respectively), the Bphcontact is noticeably weaker in the third case (interaction energy is −2.5 kcal/mol). We have also attempted to optimize a fourth possible variant, i.e. when C-H binds to an anionicoxygen and N-H to O3′ or O5′. We could not locate the corresponding minimum energy geometry on the potential energy hypersurface, since after a couple of steps the optimization was driven into the first geometry, with the amino group binding to the anionicoxygen. These results indicate that the most stable binding occurs when cytosine N6 H-bonds to an anionicoxygen and C5 H-bonds either to an anionic or esteroxygen.Each of the other four binding patterns is stabilized by single H-bonds between the nucleobase and one of the anionicoxygens of the phosphate.
Guanine
Guanine establishes H-bonds with the phosphate group in five different ways (Table 1, columns 5–6, and Supplementary Data S5), three of which use the exocyclicN2 amino group. In addition, the endocyclic imino group (N1) and the C8-H8 bond are capable of donating H-bonds to phosphateoxygens.Guanine is unique in forming the 4BPhinteraction, which is overall the most stable BPhinteraction by far. The 4BPhinteraction comprises two strong H-bonds between nitrogendonor groups of the base (N2 and N1) and two distinct oxygenacceptors of the phosphate group. The acceptor groups can be the anionicoxygens of the phosphate group or the O3′ or O5′ atoms of the riboses. Depending on the type of acceptor oxygens, the computed interaction energy pattern ranges from –5.3 to –10.1 kcal/mol. In the weakest form of the 4BPhinteraction (with an interaction energy of –5.3 kcal/mol) N2 and N1 of guanineinteract with O3′ and O5′ of the same phosphate group. The interaction is noticeable stronger (–7.8 and –8.3 kcal/mol) when N1 and N2interact with one anionicphosphateoxygen and one O3′ or O5′ atom. Finally, the highest stability is obtained when N1 and N2 both H-bond to anionicoxygens (−10.1 kcal/mol).In addition, we have found a peculiar alternative for 4BPh, a bifurcated binding mode in which N2 and N1 bind to the same anionicoxygen of the phosphate group (Supplementary Data S5). The computed stability of this binding mode is relatively high (–8.8 kcal/mol). Note that this is the only bifurcated binding pattern found for BPhinteractions. For example, we tried to optimize the bifurcated analog of 8BPh for cytosine but this failed since the optimized geometry was identical to 7BPh of cytosine, showing that the bifurcated binding is not stable if the H-bonding ability of the adjacent H-bond donors differ markedly (the 8BPh binding mode of cytosinecombines C–H and N–H H-bonds.). This does not rule out the occurrence of such an interaction in RNA molecules, however, the geometry would have to be constrained by some other interactions, in contrast to the BPhinteractions that are intrinsically stable.Due to the significant stability of the 4BPh binding pattern, we could not optimize the geometry of the G 5BPh binding pattern which involves exclusively N1 of guanine. All optimization attempts resulted in the significantly more stable 4BPhcombined binding mode. To avoid H-bond formation by the N2 amino group during geometry optimization, we replaced the N2-H2 atom with a methyl group. This electronically neutral substitution blocks H-bond formation by the N2 amino group and allowed for optimizing the geometry and calculating the energy of the 5BPh model representing G(N1)-binding.
Uracil
Uracil may bind to phosphates using its N3 endocyclic H-bond donor group on the WC edge or its C5-H5 and C6-H6 groups on the Hoogsteen edge (Table 1, columns 7–8, and Supplementary Data S6). Of these three phosphate-binding modes, the H-bond formed with the N-donor is markedly stronger (with interaction energy of –4.2 kcal/mol) than those having C-donors (interaction energies are –0.6 kcal/mol and –1.1 kcal/mol for C5 and C6, respectively). Since the N3 position of uracil is on the WC edge, while C5 and C6 are on the Hoogsteen edge, no combined binding modes involving an N-H donor are possible for uracil.
Significance of computed BPh interaction energies
The interaction energies reflect the intrinsic binding energy between the base and phosphate in a given geometry and report the difference in electronic energy between the bound and unbound systems. The relationship between the geometry and the energy complements the purely geometrical information obtained experimentally. The interaction energy describes the most interesting component of the binding that comes from the electronic structure, which tells us how the electronic shells of base and phosphatecommunicate with and adapt to each other (42). However, these calculations do not include non-electroniccontributions to binding and thus cannot be directly compared with the free energies of binding. The non-electronic terms, including entropy terms, are highly variable and depend on the specificcontexts of the interactions, while the electroniccomponent reports on the intrinsic strength of the interaction. Since the BPh systems are not electrically neutral, we report the interaction energies derived with continuum solvent screening, to attenuate the long-range ionic electrostatic effects. In summary, the interaction energies provide the basic stability ranking for the individual BPhinteraction types reflecting the strength of the bond formed between the H-bond donor and acceptor. The actual contributions to free energy may be further modulated by the specificcontext of each interaction.We refer the reader to Supplementary Data S7, where we present the results of calculations mapping the potential energy hypersurface at the optimized minima for the 6BPhinteraction of cytosine, to represent nitrogen H-bond donors, and the 2BPhinteraction of adenine, to represent carbon H-bond donors.
Relative stabilities of BPh interactions
The QM calculations show that all of the experimentally identified BPhinteractions, with the exception of 5BPh, are stable in that they occupy minima in the potential energy hypersurface. Moreover, the calculations show that 5BPh is not separated from the adjacent and more stable 4BPhinteraction by a significant energy barrier. One therefore expects that 5BPh will tend to convert to 4BPh unless the overall local arrangement favors the 5BPh. One must bear in mind that the energy difference between the 4BPh and the 5BPh or 3BPhinteractions, results from the fact that 4BPh is associated with two direct N-H … O H-bonds between the base and phosphate while 3BPh and 5BPhcontain only one. Although the COSMO solvation model includes overall solvent screening effects, this treatment is not equivalent to the explicit consideration of the eventual second direct H-bond to the phosphate group provided by another H-bond donor, possibly water, in the case of the 3BPh or 5BPhinteractions. In other words, we have some uncertainty in the quantitative comparison of the strength of 4BPh versus 5BPhinteractions. The saturation of the phosphates by direct H-bonds in these computations is unequal and the energy difference between 4BPh and 5BPh probably is somewhat exaggerated. While our stability order, 4BPh>>3BPh>5BPh is correct, the scale of energy differences is likely reduced when all possible interactions, including directly bound water molecules, are considered.
Comparison of BPh to base pair interaction energies
Comparison of the relative strength of the BPhinteractions to base pairs is not straightforward, asBPhs are electrically charged systems which respond differently to solvent screening than do base pairs, most of which are electrically neutral. Nevertheless, the interaction energy of the strongest BPhinteractions (∼–10 kcal/mol) exceeds in absolute value the interaction energy of the cis WC/WC (cWW) GC base pair calculated with the identical protocol (∼–7 kcal/mol). Therefore, it is reasonable to conclude that the most stable BPhinteractions involve very strong H-bonds and provide stabilization fully comparable with strongly H-bonded base pairs. We note that assessment of the overall thermodynamicscontribution of a given interaction in a complex nucleo-protein structure such as the ribosome is a more delicate issue than is often assumed, as the reference state may be defined in different ways. The bound state is usually compared to the state with entirely unbound base and phosphate in aqueous solution. However, in the folded ribosome, it may be more relevant to compare the stability of the formed BPh with respect to the alternative substates where the interaction is locally disrupted. The free energy balance will depend not only on the direct strength of a given interaction, but also on the cooperative formation of other interactions in its vicinity. The present calculations obviously ignore differences in solvation and cannot capture ion displacement effects (surely an important issue in BPhinteractions) and other context-dependent effects. The stabilizing effect of a given type of BPhinteraction in a given specificcontext will depend on the details of its neighborhood, including, for example, whether direct competitive binding of a positive ion to the phosphate is sterically allowed. Another open issue is competition between direct and water-mediated BPhinteractions. Nevertheless, the intrinsicinteraction energies, ascaptured by QM calculations, reflect the direct (electronic structure) communication between the interacting systems, which is the single most important term contributing to the final energy balance. Their knowledge is fundamental for basic understanding of the interactions. Although such knowledge is not sufficient for unambiguous quantitative free energy predictions in individual cases, it can at least prevent misleading interpretations that sometimes do occur in the structural biology literature when judgments are based purely on observed geometries.Thus, the main purpose of the present QM calculations is to refine and verify the basicclassification of the BPhinteractions and to qualitatively assess the strengths of distinct BPhinteraction patterns relative to each other and to base pair interactions. Further analysis of the stability of such interactions will require in-depth studies specifically dealing with distinct BPh patterns, while considering broader aspect of their structural contexts. The present data nevertheless suggest that BPhinteractions have the potential to provide substantial stabilization that is energetically comparable to well-bonded base pairs.
Base-phosphate classification in RNA 3D crystal structures
The quantum calculations indicate the optimal bond lengths for the H-bonds in BPhinteractions. In RNA 3D structures, the observed bond lengths and angles of non-covalent interactions are rarely optimal; they vary for a variety of reasons including thermal agitation of interacting partners, limited resolution of the experimental data and model-building errors. Here we discuss the variability in the bond length and angle parameters and how we set classification cutoffs for BPhinteractions.Figure 3 displays the data behind Figure 1 in another way, by showing each BPhinteraction as a point on a graph with the horizontal coordinate indicating the distance between the H-bond donor and the nearest phosphateoxygen and the vertical coordinate indicating the H-bond angle, defined as the angle between the donor-H vector and the H-O vector. The data are color-coded to indicate which phosphateoxygen atom is H-bonding. Four disjoint subsets of BPhinteractions are plotted: (i) upper-left, only intra-nucleotide (self) 0BPhinteractions; (ii) upper-right, all other BPhinteractions involving C-donors; (iii) lower-left, imino N-donors and (iv) lower-right, amino N-donors. These data show that the H-bond lengths for BPhinteractions involving carbon atoms are longer than for nitrogen atoms, consistent with the QM calculations (Table 1 and Supplementary Data S7). Furthermore, the data show that certain oxygen atoms more often participate in some interactions than others.
Figure 3.
Hydrogen bond donor to oxygen distance and angle parameters for instances of BPh interactions extracted from 3D structures in the reduced-redundancy list having resolution 2.5 Å or better. Upper left, intra-nucleotide (self) interactions involving C6 and C8 H-bond donors. Upper right, interactions with base carbon atom H-bond donors, excluding self-interactions between a base and its own phosphate group. Lower-left, interactions with imino Nitrogen, lower-right, interactions with exocyclic amino Nitrogen. Red dots correspond to H-bonds to anionic O1P or O2P acceptor atoms, green dots, H-bonds to O5′ acceptor atoms and dark blue dots, H-bonds to O3′ acceptor atoms. Black vertical and horizontal lines indicate implemented thresholds for distances and bond angles for classification as BPh (upper-left lines) and near-BPh interactions (lower right lines). Self-interactions with C6/C8 are too numerous to display all instances, so 4000 were chosen at random for display.
Hydrogen bond donor to oxygen distance and angle parameters for instances of BPhinteractions extracted from 3D structures in the reduced-redundancy list having resolution 2.5 Å or better. Upper left, intra-nucleotide (self) interactions involving C6 and C8 H-bond donors. Upper right, interactions with base carbon atom H-bond donors, excluding self-interactions between a base and its own phosphate group. Lower-left, interactions with imino Nitrogen, lower-right, interactions with exocyclicamino Nitrogen. Red dots correspond to H-bonds to anionic O1P or O2P acceptor atoms, green dots, H-bonds to O5′ acceptor atoms and dark blue dots, H-bonds to O3′ acceptor atoms. Black vertical and horizontal lines indicate implemented thresholds for distances and bond angles for classification asBPh (upper-left lines) and near-BPhinteractions (lower right lines). Self-interactions with C6/C8 are too numerous to display all instances, so 4000 were chosen at random for display.Clear clusters for discriminating ‘interacting' from ‘not interacting' instances are only evident in the scatterplot for intra-nucleotide interactions (upper-left panel of Figure 3). These interactions are by far the most numerous BPhinteractions in RNA structures and are probably most consistently modeled. We use these data as a guide to set cutoffs for BPhinteractions involving carbon H-bond donors. To be classified as a BPhinteraction, we require that the distance from a carbon H-bond donor to oxygenacceptor atom be <4 Å and the bond angle >130°. When these conditions are not met, but the H-bond length is <4.5 Å and the bond angle is >110°, the interaction is labeled ‘nBPh'. We use the same cutoffs for non-self interactions with carbon (upper-right panel of Figure 3). For interactions with nitrogendonors, we use the data in the lower-left panel of Figure 3 to set the H-bond length cutoff to 3.5 Å and the bond angle cutoff 130°, and extend this to the case illustrated in the lower-right panel of Figure 3. Cutoffs for nBPhinteractions with N-donors were set at 4 Å and 110°, as illustrated. Using ‘near' interaction categories softens the sharp divisions produced by a dichotomous classification scheme and has been found to be useful in classifying BP interactions (11). The data in Figure 3 also show that the phosphateoxygen making each type of BPhinteraction is usually the one predicted from the quantum calculations to form the more stable H-bond (see section ‘Quantum chemical calculation of optimized geometries and interaction energies for BPhinteractions’). Thus for N-H donors, the acceptor is almost exclusively an anionicoxygen, while for C-H donors a large proportion of O3′ and O5′ acceptors is observed.
Annotation of BPh interactions in secondary structures
In the previous work (43), we proposed annotating BPhinteractions in secondary structures with the letter ‘P' inside a circle, to indicate the phosphate, connected by a line to a circle, square or triangle symbol, to indicate the base edge of the H-bond donor, consistent with the Leontis–Westhof annotations for base pairs (44), where circles represent WC edges, squares Hoogsteen edges and triangles Sugar edges. Now we propose to add a number within the symbol representing the base edge to designate the BPh type (0–9). Figure 4 shows the annotated secondary structures of three RNA 3D motifs that contain BPhinteractions (two hairpin loops and one internal loop). Figure 4a shows a T-Loop motif (6,7), where U55 forms a 5BPh with A58 (the circle represents the WC edge of U55) and C61 forms a 7BPh with C60 (the square represents the Hoogsteen edge of C61); panel (b) represents a GNRA motif (8), where G159 forms 3BPh with A162 and panel (c) shows a sarcin motif (45), where A889 forms a 6BPh with A908, G890 forms a 4BPh with A907 and G906 forms a 5BPh with U893.
Figure 4.
2D annotations for (a) T-Loop from yeast Phe-tRNA (b) GNRA from T. thermophilus 16S rRNA and (c) Sarcin/ricin motif from T. thermophilus 16S rRNA.
2D annotations for (a) T-Loop from yeast Phe-tRNA (b) GNRA from T. thermophilus 16S rRNA and (c) Sarcin/ricin motif from T. thermophilus 16S rRNA.
Occurrence of BPh interactions in structured RNA molecules
In this section we report the number of BPhinteractions in selected reference structures, the distribution over the different types of interaction, and the locations of BPhinteractions in secondary structures.A large number of intra-nucleotide (self) interactions were identified, almost all of which are 0BPhinteractions involving purine (C8) or pyrimidine (C6) positions as H-bond donors. These interactions stabilize the default ‘anti' configuration of nucleotides, but are of little interest for understanding RNA 3D motifs and tertiary interactions and so were excluded from analysis. However, we included the small number of inter-nucleotide 0BPhinteractions.
Occurrence of BPh interactions by type and interaction energy
We used FR3D to classify the BPhinteractions in the reduced-redundancy subset of the PDB as described in the ‘Materials and Methods’ section. We counted the number of instances of each category of interaction, in structures having resolution 2.5 Å or better and present these data for comparison with the calculated interaction energies, in Figure 5. Generally, the more stable (lower energy) interactions occur more frequently for each base type. In particular, the most frequent interactions by base are the 6BPh of adenine (upper-left panel), the 8BPh of cytosine (upper-right panel), the 4BPh of guanine (lower-left panel) and the 5BPh of uracil (lower-right panel). These are also the lowest energy interactions calculated for each base.
Figure 5.
Comparison of calculated BPh interaction energies (red) and BPh occurrence frequencies (blue) from a reduced-redundancy set of crystal structures with resolution better than 2.5 Å (blue). Each panel represents one of the four nucleotides: Adenine (upper-left), Cytosine (upper-right), Guanine (lower-left) and Uracil (lower-right).
Comparison of calculated BPhinteraction energies (red) and BPh occurrence frequencies (blue) from a reduced-redundancy set of crystal structures with resolution better than 2.5 Å (blue). Each panel represents one of the four nucleotides: Adenine (upper-left), Cytosine (upper-right), Guanine (lower-left) and Uracil (lower-right).When lower resolution structures are included in this analysis (resolution 4.0 Å or better), we find slightly more 3BPh than 4BPhinteractions for G. This suggests that neighboring BPhinteractions are less well resolved in lower resolution structures.
Locations of BPh interactions in relation to the secondary structure
We studied the locations of BPhinteractions in relation to the secondary structure of E. coli 16S and 23S (PDB files 2avy and 2aw4, respectively). We extracted 504 classified non-self BPhinteractions involving 480 distinct bases and the phosphate groups of 473 distinct nucleotides. (Bases and phosphate groups can simultaneously make more than one BPhinteraction.) There are a total of 4488 nucleotides in these two structures and they make 1303 canonical or wobble cWW base pairs (AU, UA, CG, GC, GU or UG) and 690 non-WC base pairs. Thus, the number of BPhinteractions is on the same order of magnitude as the number of non-WC base pairs.Bases making BPhinteractions are fairly evenly distributed along the length of the nucleotide chain, although they tend to occur in small clusters, in which several successive nucleotides form BPhinteractions. Such clusters of BPhinteractions usually occur in multi-helix junctions, as is apparent for 16S rRNA in Figure 6. Secondary structure diagrams of E. coli 16S and 23S rRNA (46), annotated with the BPhinteractions that are conserved between the E. coli and T. thermophilus rRNA structures, are shown are in Figure 6 and Supplementary Data S8–S9, respectively. In these diagrams, each base that acts as an H-bond donor in a conserved BPhinteraction is marked by a semi-transparent red symbol, according to the base edge involved in the BPhinteraction (circle for WC edge, square for Hoogsteen edge and triangle for Sugar edge). A red diamond is used for the Adenine 2BPhinteraction, which straddles the WC and Sugar edges. A green dot placed between consecutive nucleotide letters indicates phosphates that act asBPhacceptors. 1BPhinteractions conserved at the base pair level are annotated with red triangles placed between bases forming the WC base pair.
Figure 6.
BPh interactions conserved between E. coli and T. thermophilus rRNA 3D structures mapped on the 2D structure of E. coli 16S rRNA (46). Red symbols were used to denote the edge used by each base donor (circle for Watson–Crick edge, square for Hoogsteen edge, triangle for Sugar edge and diamond for the Adenine 2BPh which straddles the WC and Sugar edges). The 1BPh interactions that are conserved at the base pair level, are marked by red triangles placed between the bases forming the WC base pair. Green circles denote the locations of phosphate acceptors.
BPhinteractions conserved between E. coli and T. thermophilus rRNA 3D structures mapped on the 2D structure of E. coli 16S rRNA (46). Red symbols were used to denote the edge used by each base donor (circle for Watson–Crick edge, square for Hoogsteen edge, triangle for Sugar edge and diamond for the Adenine 2BPh which straddles the WC and Sugar edges). The 1BPhinteractions that are conserved at the base pair level, are marked by red triangles placed between the bases forming the WC base pair. Green circles denote the locations of phosphateacceptors.To study the contexts in which BPhinteractions occur, we annotated each nucleotide in the E. coli 16S and 23S rRNA structures with the kind of secondary structure element to which it belongs (Helix, Hairpin Loop, Internal Loop or Junction Loop). Table 2 shows the number of each such element that contains one or more BPhinteractions in which both the base and the phosphate forming the interaction belong to the same element. We find that most hairpin, internal and junction loops in the rRNAscontain at least one BPhinteraction, but few helices do so. Indeed, the BPhinteraction is quite rare in helices and only occurs between adjacent, stacked WC base pairs at the ends of helices.
Table 2.
Fraction (%) of secondary structure elements in E. coli 16S and 23S rRNA (2avy and 2aw4) that contain one or more BPh interactions internal to that element
Secondary structure element
Number with BPh Total number
Number of secondary structure elements containing at least one internal BPh interaction
Helices
Hairpin loops
Internal loops
Multi-helix junction loops
Fraction (%) of secondary structure elements in E. coli 16S and 23S rRNA (2avy and 2aw4) that contain one or more BPhinteractions internal to that elementBPhinteractions occur between structural elements (helices, hairpin, internal and junction loops) as well as within them. We define as ‘short-range' all BPhinteractions that are internal to a single element (Table 2) or involve nucleotides in adjacent elements, such as a hairpin loop and its helical stem. We define as ‘long-range' all BPhinteractions between distinct structural elements that are not adjacent in the secondary structure. We present an analysis of short-range and long-range BPhinteractions in the E. coli 16S and 23S rRNA structures according to the locations of base H-bond donors and phosphateacceptors in Table 3. The diagonal entries in Table 3 show the numbers of internal BPhinteractions according to the type of structural element. These numbers are larger than the corresponding numbers in Table 2 because some elements of each type have more than one internal BPhinteraction. This is especially true for multi-helix junctions, where we find that 35 junctions out of total 44 contain 113 BPhinteractions or more than an average of 3 BPh per junction.
Table 3.
Number of BPh interactions occurring between secondary structure elements in E. coli 16S and 23 S rRNA (2avy and 2aw4), according to the range of the interaction
Phosphate location
Helix
Hairpin
Internal
Junction
Total
Short-range base–phosphate interactions
Helix
8
19
25
14
67
Hairpin
3
83
0
0
86
Internal
3
0
69
0
72
Junction
13
0
1
113
127
Total
27
102
95
127
353
Long-range base–phosphate interactions
Base location
Helix
24
5
5
6
42
Hairpin
10
7
9
11
37
Internal
15
2
5
5
27
Junction
15
4
7
4
31
Total
67
19
27
27
143
The data are tabulated according to the location of the base H-bond donor (rows) and phosphate acceptor (columns).
Number of BPhinteractions occurring between secondary structure elements in E. coli 16S and 23 S rRNA (2avy and 2aw4), according to the range of the interactionThe data are tabulated according to the location of the base H-bond donor (rows) and phosphateacceptor (columns).Adjacent short-range BPhinteractions occur mostly between bases located at the ends of helices and adjacent phosphates in hairpin, internal or junction loops. Long-range BPhinteractions occur in every combination, but in largest numbers between phosphates located in helices and bases located in hairpin, internal or junction loops, or in helices and between bases located in hairpin loops and phosphates in internal loops.Ribosomal RNAs are rich in BPhinteractions, having about one BPh for about every eight nucleotides. The average over the entire reduced-redundancy dataset is only 1BPh for every 24 nucleotides. However, this list contains many small RNAs that consist only of short duplexes or single strands that do not contain any junctions and few if any internal or hairpin loops, which are the primary locations at which BPhinteractions are found. Generally speaking, the BPh frequency increases as the number of non-cWW base pairs per nucleotide increases, while the BPh frequency decreases as the number of cWW base pairs per nucleotide increases. Thus 3D structures that consist primarily of helical regions have few BPh; those with many non-cWW base pairs tend to have relatively many BPhinteractions.
Sequence conservation of BPh interactions
The significant energies of BPhinteractions and their widespread occurrence in hairpin, internal and junction loop motifs and tertiary interactions suggest they play an important role in stabilizing RNA 3D structure. We can therefore expect that the locations of BPhinteractions are likely to be conserved in homologous RNA molecules. Furthermore, most BPhinteractions are base specificas to the identity of the donor base. Thus, on the one hand, we may also expect that the H-bond donor base for most BPhinteractions is likely to be conserved or to be substituted only by bases which can make geometrically and energetically equivalent interactions. On the other hand, we do not expect to observe much effect of BPhinteractions on the base identity of phosphateacceptors, since all four nucleotides have the same phosphate group.In this section, we report four tests of the conservation of BPhinteractions in homologous RNA molecules. In ‘Conservation of BPhinteractions between homologous RNA 3D structures’ section, we consider the conservation of these interactions between the 3D structures of E. coli and T. thermophilus 16S and 23S rRNA. In ‘Conservation of H-bond donor bases in BPhinteractions among homologous rRNA sequences’ section, we examine the conservation of the H-bond donor between the 3D structures and multiple sequence alignments of homologous molecules. In ‘Predicting base conservation in sequences using pairwise interactions in 3D structures’ section we compare the contributions of BPhinteractions and other interactions (WC and non-WC base pairs, BSt and protein interactions) to base conservation at homologous sites. Finally, in ‘Conservation of BPhinteractions at the base pair level’ section we discuss conservation of certain BPhinteractions at the base pair level.
Conservation of BPh interactions between homologous RNA 3D structures
Using FR3D, we extracted, classified and aligned all BPhinteractions for the 16S and 23S rRNA structures of E. coli (2avy and 2aw4) and T. thermophilus (1j5e and 2j01), as described in the ‘Materials and Methods’ section. This alignment is provided as an Excel file in the Supplementary Data S2. We found a total of 592 BPhinteractions in the E. coli. rRNA structures (195 in E. coli 16S and 397 in 23S rRNA) and 601 total BPhinteractions in the T. thermophilus structures (203 in T. thermophilus 16S and 398 in 23S rRNA). Since some bases form more than one BPhinteraction, we found a total of 557 distinct bases in the E. coli rRNA and a total of 562 distinct bases in T. thermophilus rRNA acting as H-bond donors in BPhinteractions, or about 13% of bases. These data are summarized in Table 4.
Table 4.
Frequencies of non-self BPh interactions in E. coli and T. thermophilus 16S and 23S rRNAs
16S rRNA 3D structures
23 S rRNA 3D structures
Escherichia coli (2avy) (1j5e)
Thermus thermophilus (lj5e)
Escherichia coli (2aw4)
Thermus thermophilus (2j01)
Total nucleotides
1530
1513
2841
2772
Number of distinct bases involved in BPh interactions
185 (12.1%)
191 (12.6%)
372 (13.1%)
371 (13.4%)
Total number of BPh interactions
195
203
397
398
Number of BPh interactions at corresponding E. coli and T. thermophilus positions
174
338
Conserved bases at corresponding BPh positions
About 13% of all bases in the bacterial rRNA structure from BPh interactions and ∼86% of these interactions are common to the E. coli and T. thermophilus rRNA structures. For the corresponding BPh interactions, the base is ∼95% conserved between E. coli and T. thermophilus.
Frequencies of non-self BPhinteractions in E. coli and T. thermophilus 16S and 23S rRNAsAbout 13% of all bases in the bacterial rRNA structure from BPhinteractions and ∼86% of these interactions are common to the E. coli and T. thermophilus rRNA structures. For the corresponding BPhinteractions, the base is ∼95% conserved between E. coli and T. thermophilus.Comparing the lists of BPhinteractions in the E. coli and T. thermophilus rRNA 3D structures, we found 512 ‘corresponding' BPhinteractions (174 in 16S rRNA and 338 in 23S rRNA) for which the nucleotide positions of both the base H-bond donor and the phosphateacceptor can be aligned in the structures. BPhinteractions were counted as aligned if FR3D detected a BPhinteraction in both structures at aligned positions or if FR3D detected a nBPhinteraction in one structure at aligned positions annotated with a full BPhinteraction in the other structure. Thus, ∼86% ((2 × 512)/1193) of the BPhinteractions are common to the E. coli and T. thermophilus rRNA 3D structures. Furthermore, in 484 out of the 512 aligned BPhinteractions (∼95%), the H-bond donor base is the same in the E. coli and T. thermophilus structures (164 in 16S and 320 in 23S rRNA, Table 4).Thus, <15% of the BPhinteractions found in the E. coli 16S or 23S rRNA structure are not found at the corresponding position of the homologous T. thermophilus structure and vice versa. We have examined each of these instances manually and find that they belong to one of five categories: (i) In the first category, data are missing in one of the two structures due to crystallographic disorder (∼2%). Thus, for these instances one cannot tell at present whether the interaction is conserved. (ii) In the second category, there is a major structural difference between the E. coli and T. thermophilus molecules at the position in question, such as a motif swap or a difference in length of a peripheral element. For such instances, about 3% of the total, no comparison can be made and these are best excluded from analysis. (iii) The third category consists of instances where the structures are modeled differently at the nucleotide level (∼5%). Most such instances involve a difference in modeling of the glycosidic bond conformations of the corresponding nucleotides acting as H-bond donors. For example, A65 is modeled ‘syn' in E. coli 16S but ‘anti' in T. thermophilus 16S and is only able to form the 2BPhinteraction in the E. coli structure. Differences in modeling the glycosidicconformation, which may indicate possible errors in one or both structures, are not infrequent in large RNA structures solved at moderate resolution (12) and can be resolved with better crystal data. Thus, some of these instances may also prove to be conserved BPhinteractions. (iv) In the fourth category, the BPhinteraction is conserved at the base pair level rather than the nucleotide level (∼1% of instances). This case will be discussed in more detail in ‘Conservation of BPhinteractions at the base pair level’ section. (v) The fifth category comprises actual base substitutions that preclude forming equivalent BPhinteraction in both structures. These cases (about 3% of instances) are genuine instances where a BPhinteraction does not appear to be conserved. Even in these cases, however, a conformational change between the two structures can compensate for the base change. An example is position 121 in 16S which is U in E. coli 16S and C in T. thermophilus 16S. Both bases are bulged and both form a BPhinteraction with helix 7. U121 forms a 5BPh with the C234 phosphate and C121 forms a 7BPh with the G236 phosphate.Focusing on the ∼86% of BPhinteractions which can be aligned between the E. coli and T. thermophilus 16S and 23S rRNAs, we estimated the conservation rate of the H-bond donor base between structures. These data are shown in Table 5, which organizes the BPhinteractions according to the base edges to which the H-bond donors belong. The cells along the diagonal correspond to BPhinteractions which are identical as to interaction type and base identity in the E. coli and T. thermophilus rRNA structures. These constitute ∼75% of the BPhinteractions. Examination of the off-diagonal cells shows that those with the highest counts involve neighboring BPhinteractions in which the base is the same in the two structures. These cells are indicated by dotted red lines and correspond to guanine 3BPh, 4BPh and 5BPh neighboring interactions (56 instances) or cytosine 7BPh, 8BPh and 9BPh neighboring interactions (27 instances). Combining these instances of neighboring interactions with the instances of identical interactions on the diagonal, we find 92% of BPh instances have the same base in the E. coli and T. thermophilus rRNA structures making the same or neighboring BPhinteraction. Of the remaining ∼8% of BPh instances from the rRNA 3D structures, 11 instances (∼2%) are cases of equivalent BPhinteractions formed by bases which are different in the E. coli and T. thermophilus structures. These are indicated by solid red outlines in Table 5 and include cases of 0BPh, 5BPh, 6BPh, 7BPh and 9BPh. Another set of instances involves neighboring BPhinteractions with base changes, for which there are 16 instances (∼3%). This leaves only 12 instances (∼2%) that involve changes of base edge (pink background). Thus, the data in Table 5 clearly show that when changes occur, the interacting edge of the base is almost always conserved (yellow cells).
Table 5.
Corresponding BPh interactions observed in the 3D structures of E. coli and T. thermophilus 16S and 23S rRNAs
Diagonal entries (dark green) correspond to identical BPh interactions (same base donor and BPh category). Yellow shaded cells correspond to differences in base or BPh category that preserve the geometry of the interaction. Pink cells indicate differences that do not preserve the BPh geometry.
Corresponding BPhinteractions observed in the 3D structures of E. coli and T. thermophilus 16S and 23S rRNAsDiagonal entries (dark green) correspond to identical BPhinteractions (same base donor and BPhcategory). Yellow shaded cells correspond to differences in base or BPhcategory that preserve the geometry of the interaction. Pink cells indicate differences that do not preserve the BPh geometry.
Conservation of H-bond donor bases in BPh interactions among homologous rRNA sequences
In the previous section we obtained direct evidence for conservation of BPhinteractions by comparing the 3D structures of the homologous but phylogenetically distant E. coli and T. thermophilus rRNA molecules. The BPhinteractions common to these structures are likely to be present in other 16S and 23S rRNAs, because they are observed in both structures and are found within the core regions of the rRNA secondary structures (12). In this section, we examine rRNA multiple sequence alignments to obtain larger numbers of instances of BPhinteractions for more robust statistical analysis of BPh sequence variations. We use the multiple sequence alignments of 717 16S and 136 23S sequences that we previously reported in a study of base pair substitutions (12). For each aligned BPhinteraction identified in the 3D alignment, we examined the column corresponding to the H-bond donor in the multiple sequence alignment and tallied the number of instances of A, C, G and U observed in that column. By restricting our attention to locations of BPhinteractions that are conserved between the E. coli and T. thermophilus rRNA 3D structures we gain assurance that the interaction is also found in other homologous rRNA molecules for which we only have sequences.Figure 7 presents on separate plots the data obtained from the alignments, for each BPh type of interaction observed in the E. coli 3D structures. We first consider the 9BPh subplot. This BPhinteraction is specific to C and U. The blue dots represent instances in which U is the H-bond donor in a 9BPhinteraction in the E. coli 3D structures. Those blue dots clustered near the corner of the graph labeled ‘U' correspond to instances that are also predominantly U in the sequence alignments. Likewise, the yellow dots clustered near the Ccorner represent instances that are predominantly C in the sequence alignments, while the blue and yellow dots along the edge that connects the U and Ccorners of the square indicate instances that vary exclusively between U and C in alignments. Together these account for most instances and provide support for the idea that mutations that interchange U and C are structure-neutral as either U or Ccan support the 9BPhinteraction. There is one instance with about 60% U and 40% G (blue dot on bottom edge of the square) and only a handful of yellow dots with significant fractions of A and G as well asC and U. These may represent regions where the alignments are less certain or sites where the 9BPhinteraction is not conserved in all structures.
Figure 7.
Base variations in BPh interactions observed in E. coli rRNA 3D structures and corresponding columns of bacterial rRNA sequence alignments. The subplot titles show the number of BPh interactions of each type that are aligned between the E. coli and T. thermophilus 3D structures. Each aligned instance in a 3D structure is represented by a dot, colored to indicate which base was present in the E. coli structure: red for A, yellow for C, green for G and blue for U. The location of the dot indicates the percentage of A, C, G and U found in the corresponding column in the multiple sequence alignment described in the ‘Materials and Methods’ section. The precise location of the dot is a weighted average of the four corner locations. Note that the 0BPh category excludes intra-nucleotide self-interactions.
Base variations in BPhinteractions observed in E. coli rRNA 3D structures and corresponding columns of bacterial rRNA sequence alignments. The subplot titles show the number of BPhinteractions of each type that are aligned between the E. coli and T. thermophilus 3D structures. Each aligned instance in a 3D structure is represented by a dot, colored to indicate which base was present in the E. coli structure: red for A, yellow for C, green for G and blue for U. The location of the dot indicates the percentage of A, C, G and U found in the corresponding column in the multiple sequence alignment described in the ‘Materials and Methods’ section. The precise location of the dot is a weighted average of the four corner locations. Note that the 0BPhcategory excludes intra-nucleotide self-interactions.Some BPhinteractions are highly conserved in the sequence alignments. For example, the 3BPh and 4BPhinteractions, which are specific to G (all green dots), are also largely G in the sequence alignments (all dots map close to the G corner). Likewise, the 8BPhinteraction, which is specific to C (all yellow dots) is also largely C in the sequences (all dots close to the Ccorner), with the exception of one instance that is mainly A. The 1BPhinteraction is also specific to G and most instances in sequences are also G, but some sites have significant percentages of C, U or A. This interaction, which occurs on the minor groove, sometimes involves an interesting base pair covariation, which may explain the dots on the diagonal between the G and Ccorners of the 1BPh subplot and will be discussed in more detail in ‘Conservation of BPhinteractions at the base pair level’ section.The 2BPhinteraction is specific to A, but is geometrically similar to the 1BPh or 3BPh formed by G. The 3D structure alignment reported in ‘Conservation of BPhinteractions between homologous RNA 3D structures’ section provides examples of A2BPh to G1BPh substitutions, Table 5. This interpretation may explain the red dots along the A–G edge of the 2BPh subplot, which are instances for which we see significant fractions of G in the sequence alignments.Category 5BPh is only made by G or U and we see some G ↔ U variation, as expected.Category 6BPh is only made by A or C. The 5BPh and 6BPhinteractions, which can be made by A, C, G or U, all involve the WC edge of the base, so it is no surprise to see wider sequence variation in these subplots. Table 5 shows that such substitutions are observed in 3D structures. Even so, most locations show near-total conservation; in each subplot, there are only ∼12 instances out of total 139 5BPh or 6BPhinteractions that are not close to a corner.Finally, the 0BPhinteraction is not base-specific and shows the greatest sequence variation of all.
Predicting base conservation in sequences using pairwise interactions in 3D structures
The results presented in the previous two sections indicate that when a BPhinteraction is observed in one RNA 3D structure, it is very likely to be present at the corresponding location in homologous RNA molecules. Moreover, the bases acting as H-bond donors in BPhinteractions tend to be very conserved in homologous molecules. We can ask whether this conservation should be attributed to the BPhinteraction itself, or to other interactions that the base makes. In other words, are BPhinteractions primarily incidental, or do they have a selective effect on RNA sequence? It is difficult to answer this question conclusively for individual instances of BPhinteractions, but we gain some understanding if we examine more globally the degree of base conservation in RNA molecules in relation to other factors that are known to constrain sequence variation.We mapped the 3D structures of E. coli 16S and 23S rRNA (PDB files: 2avy and 2aw4) to the sequence alignments (described above), to determine for each nucleotide in the 3D structure, which bases appear in the corresponding column of the multiple sequence alignment. We define the conservation of each base in the structure as the percentage of bases (A, C, G or U) in the corresponding column of the alignment that are identical to it. For this calculation, we ignore gaps and symbols indicating incomplete data in the sequence alignment. Thus, if the 3D structure has an A and in the corresponding column there are 200 As, 100 Cs, 50 Gs, 50 Us, 85 gaps and 35 Ns, the conservation percentage is calculated 200/(100 + 200 + 50 + 50) = 50%.We investigated possible factors that might influence nucleotide conservation as defined in the previous paragraph. Presumably the interactions the nucleotide makes in the folded RNA structure have a major effect. We therefore used FR3D to tabulate all the annotated interactions that each nucleotide in the structure makes, including the number of cWW base pairs in which the nucleotide is involved (either 0 or 1), the number of non-cWW base pairs (0, 1, 2 or 3), the number of distinct BPhinteractions in which the base of the nucleotide is the H-bond donor (0, 1 or 2) and the number of interactions in which the phosphate of the nucleotide is the acceptor (0, 1 or 2). We also calculated whether a nucleotide is within 3.7 Å of a protein residue in the 3D structure (1 if so, 0 if not). Protein contacts were found using Swiss PDB Viewer (47), by selecting all nucleotides within a neighboring radius of 3.7 Å of a protein moiety. Finally, to account for the general level of conservation among neighboring nucleotides, for each nucleotide we calculated the average conservation percentage for the nucleotides which base stack on it (‘stacking partners'), or, if there are no stacking partners, the average conservation percentage of the nucleotides immediately before and after it in the chain (‘adjacent nucleotides'). Treating these numbers as the predictors of base conservation, we performed a simple linear regression on these parameters, resulting in the following model:Conservation percentage of base in sequences= 26.3 (constant term)+ 58.1 × conservation percentage of stacking partners or adjacent nucleotides, divided by 100+ 8.0 × number of non-cWW pairs it forms+ 7.4 × number of BPhinteractions in which it is the H-bond donor+ 2.7 × near a protein (1) or not (0)+ 1.8 × number of BPhinteractions in which it is the phosphateacceptor+ 0.3 × number of cWW pairs it formsAll coefficients except the last two are non-zero at the 0.05 significance level. The conservation percentage of stacking partners makes the strongest contribution. It accounts for the general level of conservation of the RNA region in which the nucleotide occurs and for possible base specificity of stacking interactions. With that accounted for, nucleotides which participate in one or more non-cWW base pairs have elevated conservation percentages. Non-cWW base pairs generally exhibit less variation between structures (12). Bases that are H-bond donors in one or more BPhinteractions also have elevated conservation percentages, with roughly the same strength of effect as the number of non-cWW base pairs formed. Proximity of the nucleotide to a protein has a moderate effect. The last two factors are the weakest, and may be statistical fluctuations rather than real effects. We do not expect nucleotides which act asphosphateacceptors to be more or less conserved due to this interaction, since it is not base specific. The very small contribution to conservation due to a nucleotide making a cWW pair is consistent with our understanding that the AU, UA, CG, GCcWW pairs are isosteric and can freely substitute for one another, in the absence of further constraints. Any base specificity of stacking for cWW pairs is accounted for by the conservation percentage of stacking partners.This is a simple but fairly complete model for the dependence of base conservation in homologous RNA sequences on the type and number of pairwise interactions a nucleotide makes. The model indicates that bases that act as H-bond donors in BPhinteractions are subject to sequence constraints comparable to those of bases forming non-WC base pairs.In Supplementary Data S10, we use boxplots to show the statistical distribution of conservation percentages for subsets of nucleotides having all possible combinations of three binary attributes: (i) whether the nucleotide makes a cWW pair, (ii) whether the nucleotide makes at least one non-cWW pair and (iii) whether the nucleotide is the H-bond donor in at least one BPhinteraction. The boxplots provide additional evidence that a base, i.e. the H-bond donor, in a BPhinteraction is more conserved than one that is not.
Conservation of BPh interactions at the base pair level
In previous sections, we examined the conservation of individual nucleotides making BPhinteractions by comparing the corresponding positions in two 3D structures or the 3D structures and homologous sequence alignments. We found some cases where Gs forming 1BPhinteractions in one rRNA structure did not appear to be conserved in the other structure. Closer examination revealed that for these cases the 1BPhinteraction is in fact conserved on the base pair level rather than the nucleotide level. Figure 8 shows an example: G2692 in E. coli 23S makes a 1BPh with G2848 and a cWW base pair with C2717. The corresponding base pair in T. thermophilus hasGC replaced with CG, but the G again makes a 1BPhinteraction. This is possible because when superposing a GCcWW base pair on a CGcWW base pair, the G(N2) amino groups of the Gs coincide in space so that the G(N2) group is positioned to participate in an equivalent 1BPhinteraction with a phosphate located in the same place in the middle of the minor (shallow) groove.
Figure 8.
Conservation of 1BPh interaction at the level of base pairs. The G2692/C2717 cWW base pair of E. coli 23S rRNA (a) corresponds to the C2692/G2717 base pair in T. thermophilus 23S (b). The G in each structure forms a conserved 1BPh interaction with the phosphate of nucleotide 2848, as shown (c).
Conservation of 1BPhinteraction at the level of base pairs. The G2692/C2717 cWW base pair of E. coli 23S rRNA (a) corresponds to the C2692/G2717 base pair in T. thermophilus 23S (b). The G in each structure forms a conserved 1BPhinteraction with the phosphate of nucleotide 2848, as shown (c).We have identified two other instances of GC ↔ CG substitutions in the E. coli and T. thermophilus rRNA 3D structures that preserve 1BPhinteractions. In Table 6, we summarize these cases and four others in which a GCcWW base pair in one structure is substituted by a UG or CA cWW base pair in the other structure and the purine in each pair makes a similar BPhinteraction (G makes 1BPh, A makes 2BPh).Table 6 also presents base pair substitution data from the rRNA sequence alignments for the seven identified instances of base pair level BPhconservation. In six of the seven cases, the base combinations observed in the 3D structures are the combinations most often seen in the sequence alignments as well. These data offer confirmation that the base pair substitutions observed in the 3D structures occur with some frequency in other homologous molecules. As noted previously (12), base pair substitutions occur most often between base combinations that make isosteric base pairs. Thus, GC base pairs are most often substituted by CG, AU and UA (isosteric, red boxes), but also by GU, UG, AC and CA (near isosteric, yellow boxes). This explains the very low counts in the boxes shaded blue and gray.
Table 6.
Substitution tables for identified cWW base pairs in 16S and 23S bacterial rRNAs involved in 1 BPh conservation at the base pair level
For each substitution table the sequence for E. coli and T. thermophilus are given in the header, while the corresponding nucleotide numbers are given in the upper-left corner of each table. The cells are colored according to whether the base pair is isosteric (red), near isosteric (yellow), or non-isosteric to the GC cWW base pair as discussed in (12).
Substitution tables for identified cWW base pairs in 16S and 23S bacterial rRNAs involved in 1 BPhconservation at the base pair levelFor each substitution table the sequence for E. coli and T. thermophilus are given in the header, while the corresponding nucleotide numbers are given in the upper-left corner of each table. The cells are colored according to whether the base pair is isosteric (red), near isosteric (yellow), or non-isosteric to the GCcWW base pair as discussed in (12).We have identified 20 other aligned positions in the E. coli and T. thermophilus rRNA 3D structures where Gs form 1BPhinteractions as well ascWW GC pairs. In these cases both structures have the G in the same position. We compiled base pair substitutions for each of these 20 instances from the sequence alignments and find that thirteen of these show 90% or greater conservation of the GC base pair per se, while the other seven show patterns of variation similar to those shown in Table 6. Including the three cases in Table 6, of which one shows 90% or greater conservation of the GC base pair, there are 23 cases of GCcWW pairs making 1BPh in both structures and nine of them (∼40%) exhibit significant levels of base pair substitutions in the sequence alignment. Moreover, most of the variability involves GC, CG, GU or UG covariations which, as we have seen, preserve the 1BPhinteraction in the minor groove. This phenomenon can help to explain some of the variability in the 1BPh subplot of Figure 7. When a GC in the 3D structure in which the G makes a 1BPh is substituted by a CG or UG base pair in the sequence alignment, the 1BPh subplot of Figure 7 will show non-conservation of the G, but this masks the fact that the BPhinteraction may still be conserved at the base pair level. Figure 9 shows the data in the 1BPh subplot of Figure 7 separated according to whether the G making the 1BPhinteraction also forms a cWW base pair. Those that do make a cWW base pair show higher tendency toward substitution, consistent with the phenomenon described above.
Figure 9.
Base substitution data for G1BPh interactions, as described in the caption for Figure 7. The left subplot shows cases in which the G does not make a cWW base pair, while the right panel shows cases in which the G does make a cWW base pair.
Base substitution data for G1BPhinteractions, as described in the caption for Figure 7. The left subplot shows cases in which the G does not make a cWW base pair, while the right panel shows cases in which the G does make a cWW base pair.
DISCUSSION
The analysis of the 3D structures of E. coli and T. thermophilus 16S and 23S rRNAs indicates that a significant fraction of bases in structured RNA molecules (about 13%) form conserved inter-nucleotide BPhinteractions. Moreover, they are widespread in hairpin, internal and multi-helix junction loops. We find that ∼87% of hairpin loops, ∼80% of junction loops and ∼71% of internal loops of E. coli 16S and 23S rRNA contain one or more BPhinteractions involving nucleotides of the same loop. As many of these motifs are recurrent, it is likely that BPhinteractions are also found widely in other structured RNAs. A significant number of the conserved BPhinteractions are long-range (143 out of 496), which suggests that BPhinteractions also play significant roles stabilizing the tertiary structures of structured RNA molecules.Our analysis revealed 17 unique phosphate-binding sites on the standard RNA bases (A, C, G and U). Based on the reported quantum chemical calculations and consideration of the locations of these binding sites along the base perimeters (edges), we decided to classify them in 10 families. This groups together, under the same designation, phosphate-binding sites that occur at equivalent sites on different bases. Thus, 5BPh designates the imino nitrogen H-bond donors on the WC edges of G and U with a phosphateoxygen. Other BPhinteractions that are shared by more than one base are 6BPh, on the WC edges of A and C, 7BPh, on the Hoogsteen edges of A and C, 9BPh, on the Hoogsteen edges of U and C and the non-specific 0BPh, also on the Hoogsteen edge. The other BPhinteractions are base-specific. BPhinteractions that occur on the same base edge, such as the 3BPh, 4BPh and 5BPhinteractions of G, are called ‘neighboring interactions' and are very conserved.We developed the proposed BPh nomenclature to facilitate integration of structural and sequence data. Comparing the 3D structures of E. coli and T. thermophilus 16S and 23S rRNA we found that at least 86% of the BPhinteractions detected in the homologous E. coli or T. thermophilus structures are common to both structures. AsE. coli and T. thermophilus are distantly related bacteria, it is likely that most of these interactions also occur at corresponding sites in homologous bacterial rRNAs. This degree of conservation over billions of years of evolution provides further evidence that BPhinteractions play significant roles in the ribosome. Comparing the annotations of conserved BPhinteractions in the E. coli and T. thermophilus structures we find in most cases (about 75%) the base and the class of interaction are the same in the E. coli and T. thermophilus structures. Closer examination showed that most of the differences could be accounted for in two ways: in many cases, the base is conserved, but the BPh annotations of the E. coli or T. thermophilus structures indicated neighboring interactions, that is BPhinteractions involving the same edge, such as G 4BPh in one structure and G 5BPh in the other. In other cases, the bases were not the same in the two structures, but the interactions were the same, as in U 9BPh in one structure and a C 9BPh in the other, or the structures had neighboring interactions, e.g. A 6BPh in one structure and G 5BPh in the other. Therefore, we conclude that, taking into account conservative base substitutions and neighboring interactions which allow BPhinteractions to form without drasticchanges in geometry, close to 98% of BPhinteractions that occur at homologous sites are conserved. This number includes 1BPhinteractions in the minor groove of RNA helices that are conserved at the level of the base pair, a novel kind of sequence covariation described here for the first time. A linear regression model developed to compare the influences of BPh to other interactions that are known to constrain sequences shows that bases which form BPhinteractions are strongly constrained, in a manner comparable to bases forming non-WC base pairs.Auffinger et al. showed by analyzing X-ray crystal data and MD simulations that mobile anions can intrude into the first hydration shell of nucleic-acid bases (48). They observed a number of anion-binding sites and inferred or proposed several more. They identified about 13 anion-binding sites along the base edges. They designated binding sites according to the H-bond donor base (ADE, CYT, GUA or CYT) and interacting base edge (WC, S, H and CH, for the WC, Sugar and Hoogsteen/CH edges). The following correspondences can be made between the two classifications: Auffinger et al. GUA_S anion-binding site corresponds to our 1BPhclass; the inferred GUA_WC_S site corresponds to 3BPh; GUA_WCcorresponds to 4BPh or 5BPh; ADE_WC_H and the inferred CYT_WC_CH correspond to 6BPh; CYT_CH corresponds to 7BPh or 8BPh. In addition, Auffinger et al. propose bifurcated (B) anion binding to both amino hydrogen atoms simultaneously (ADE_B, CYT_B, GUA_B) and binding to A or C protonated on the WC edge (ADE+_WC and CYT+_WC).
Role of BPh interactions in stabilizing RNA structures
The quantum computations show that experimentally identified BPhinteractions are intrinsically stable. We did not find it necessary to apply constraints to maintain the observed interaction patterns during geometry optimization, except in the case of G5BPhinteractions, which minimize to 4BPh when not constrained. The QM studies clearly indicate that the observed geometries are energetically stable per se and do not depend on neighboring stabilizing interactions. By contrast, in related QM studies of certain non-WC base pairs we had to apply constraints to maintain the experimentally observed geometries (49). Moreover, BPhinteractions are directional and can stabilize specific RNA geometries. Further, the BPhinteraction energies are large enough to play significant roles in RNA folding and stability, corroborating our findings based on comparative structure and sequence analysis that these interactions are very conserved in rRNA evolution. It is likely that this is also the case for other structured RNAs. We found that for each base the relative abundances of the various BPh-binding patterns roughly correlate with the computed interaction strengths, providing further evidence of their relevance and importance in RNA structure. These observations also indicate that the intrinsic energetics of molecular interactions can contribute to modulating or supplementing the primary geometrical (isostericity) signatures of molecular interactions. Thus, relative energies can also affect the evolution of sequences of otherwise isosteric structures (see ‘Comparison of BPh to base pair interaction energies’ section for a discussion of the limitations of the energy analysis). While the classification of molecular interactions in RNAs was until now based solely on structural information, the present study suggests that complementary energy calculations can help to refine and better understand the classification. In fact, the present energy calculations helped to classify the BPh-binding patterns and clarified the nature of several binding modes.There is evidence from several sources that BPhinteractions effectively ‘solvate' the negatively charged phosphate groups of the RNA backbone. The first line of evidence is the size of the interaction energies themselves. Second, the ‘internal' solvation they provide can modulate or limit the formation of deep pockets of negative electrostatic potentials (ESP). For example, the bacterial 5S Loop E motif and the S/R Loop motif (which is also present in Loop E of eukaryal and archaeal 5S rRNAs) comprise some of the same non-WC base pairs (50). However, the S/R loop contains at least three conserved BPhinteractions, all in the major groove, that are not observed in the bacterial Loop E motif. Consistent with these observations, only bacterial Loop E exhibits a deep negative ESP pocket in its major groove, which is known, moreover, to be a strong binder of cations (51–54). The conserved BPhinteractions in the S/R loop are formed by the conserved purine (usually G) of the AR tHS base pair (R = A or G), the conserved G of the cHS UG base pair and the conserved A of the AM tHH pair (M = A or C). The first two BPhinteractions occur on the WC edge of the base and as we have seen, the strongest BPhinteractions on the WC edge are formed by Gs. G is observed exclusively as the second base in the cHS base pair and is by far most frequent as the second base in the tHS pair. Substitution of A for G in the tHScan be accommodated by a slight adjustment in the backbone to form an A 6BPh, but substitution of U or C does not allow a BPhinteraction to form in this context. In other motifs, such as T-loops, substitution of U for a G, to form a 5BPhinteraction is possible, but requires a small adjustment in the backbone, as shown by superposition of T-loops containing U or G at the site forming this conserved BPhinteraction.In a recent survey of base pair frequencies by geometric base pair family (12), we found that 69% of tHS base pairs are AG, 71% of tHH base pairs are AA and 47% of cHS base pairs are UG. Searching the non-redundant 3D database using FR3D, we find that the nucleotides in these non-WC base pairs also frequently have a BPhinteraction between them (45). The fact that these base combinations can simultaneously form strong BPhinteractions may in part account for their high frequencies.
How these data can be used
This work shows that the presence of BPhinteractions in RNA 3D structures strongly constrains the sequences of homologous RNA molecules. The base acting as H-bond donor in conserved BPhinteractions is very likely to be identical at corresponding positions in homologous sequences, when the sequences are properly aligned. We now provide annotations of BPhinteractions for all RNA-containing PDB/NDB structures, in addition to BP and BSt interactions. Annotations can be accessed on the FR3D website (http://rna.bgsu.edu/FR3D/AnalyzedStructures/) annotations index page, by clicking on the link labeled ‘FR3D' in the second column of the record corresponding to each PDB entry.Knowledge that a particular BPhinteraction occurs at a particular position in a 3D structure also indicates possible base substitutions that may occur in homologous sequences, although with generally small probabilities. If a 1BPhinteraction is observed in the 3D structure and the G H-bond donor also forms a GC or GU WC base pair, the interaction is likely to be conserved at the base pair level, as described in ‘Conservation of BPhinteractions at the base pair level’ section. In homologous sequences, the base pair is most likely to be conserved asGC or GU or to be substituted by CG or UG. If the G forming the 1BPh is not forming a WC base pair, the G is most likely to be conserved as G in homologous sequences, or, with much lower probability, to be substituted by A, forming a neighboring 2BPhinteraction. Likewise the presence of a 2BPhinteraction in a 3D structure indicates the A acting as H-bond donor is most likely to be conserved as A in homologous sequences, or with lower probability, to be substituted by G to form either a neighboring 1BPh or 3BPhinteraction, as the 2BPh straddles the Sugar and WC edges.Gs forming BPhinteractions on their WC edges can, with small probabilities, be substituted by As or Us, depending on the context. As and Cs forming 6BPh on the WC edge show more variability and can be substituted by each other or by G. Cs forming 7BPh or 8BPh are quite conserved but can be substituted by A with small probabilities. Of all BPhinteractions the pyrimidine 9BPhinteraction showed the most regular variations in sequences, almost exclusively between C and U. The inter-nucleotide 0BPhinteraction is dominated by purines, which show some tendency to substitute.Comparison of the E. coli and T. thermophilus 3D structures, which at best are based on moderate-resolution data, shows that the corresponding BPhinteractions on the WC edges of Gs or the Hoogsteen edges of Cs are modeled differently in a number of cases (e.g. G4BPh in one structure and G 3BPh or 5BPh in the other). In higher resolution structures, we find higher fractions of interactions modeled as the more stable G4BPh or C8BPhinteractions. Thus, limited resolution is probably one factor contributing to differences between the structures. Another factor may be the inherent dynamics of BPhinteractions. Base pairs are known to open and close on the milli-second time-scale (55,56). AsBPhinteractions have intrinsic stabilities comparable to base pairs, it is possible that they fluctuate between neighboring interactions (e.g. C 7BPh or 9BPh and 8BPh, or G 3BPh or 5BPh and 4BPh) in response to structural dynamics of the RNA. Thus, the uncertainty in the X-ray structures may in part reflect intrinsic dynamics of the structures themselves at some sites.Given the intrinsic stabilities of BPhinteractions, comparable to base pairs, it is not surprising that the presence of BPhinteractions protects the base from reactive chemical probes. In a landmark study, published long before the 3D structures were known, Noller and coworkers (57) probed E. coli 16S rRNA with chemical probes that react with the WC edges of RNA bases. They detected reactive sites by primer extension assays. Their data, in combination with the FR3D annotations of the 3D structures, show clearly that WC edges of the bases acting as H-bond donor are protected when forming 3BPh, 4BPh, 5BPh or 6BPhinteractions, once the 16S achieves its native fold. The reported protection is complete for most bases annotated with these BPhinteractions in the 3D structures and at least partial for the rest. For Gs, bases annotated as making G 3BPh or 5 BPhinteractions were more likely than G 4BPh to be partly accessible to probes, consistent with the partial exposure of the WC edge in these cases.
CONCLUSIONS
By combining empirical and quantum chemical methods, we have identified optimal phosphate-binding geometries for each RNA base and implemented software in FR3D to find and classify these interactions in 3D structures (http://rna.bgsu.edu/FR3D/BasePhosphates). We have shown that BPhinteractions are highly conserved in homologous rRNA molecules. We have shown that BPhinteractions strongly constrain the sequences of homologous structured RNAs. We suggest that conclusions we have drawn based on studies of rRNA may also apply to other homologous structured RNAs. These data should be useful in the interpretation of chemical probing experiments and should provide better understand of RNA sequence variability, to improve algorithms for aligning homologous sequences and predicting RNA 3D structure and motifs.
SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.
FUNDING
Grants to NBL from the National Institutes of Health (2 R15GM055898-04), from the National Science Foundation (Research Coordination Network Grant No. 0443508), and from Bowling Green State University (Research Capacity Expansion Program, funded by the Ohio Board of Regents Research Incentive Fund); Academy of Sciences of the Czech Republic (grants no. AV0Z50040507 and AV0Z50040702); the Grant Agency of the Academy of Sciences of the Czech Republic (grants No. IAA400550701, IAA400040802 and 1QS500040581); Grant Agency of the Czech Republic (grant GA 203/09/1476) and Ministry of Education of the Czech Republic (grant LC06030). Funding for open access charge: National Science Foundation (Research Coordination Network Grant No. 0443508).Conflict of interest statement. None declared.
Authors: Vojtech Mlýnský; Pavel Banás; Daniel Hollas; Kamila Réblová; Nils G Walter; Jirí Sponer; Michal Otyepka Journal: J Phys Chem B Date: 2010-05-20 Impact factor: 2.991