Global importance of RNA secondary structures in protein-coding sequences.

MOTIVATION: The protein-coding sequences of messenger RNAs are the linear template for translation of the gene sequence into protein. Nevertheless, the RNA can also form secondary structures by intramolecular base-pairing.
RESULTS: We show that the nucleotide distribution within codons is biased in all taxa of life on a global scale. Thereby, RNA secondary structures that require base-pairing between the position 1 of a codon with the position 1 of an opposing codon (here named RNA secondary structure class c1) are under-represented. We conclude that this bias may result from the co-evolution of codon sequence and mRNA secondary structure, suggesting that RNA secondary structures are generally important in protein-coding regions of mRNAs. The above result also implies that codon position 2 has a smaller influence on the amino acid choice than codon position 1. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

1 Introduction

RNA secondary structures are involved in the regulation of mRNA function and fate. In many cases, functional RNA secondary structures are located in the 5'- or 3'-untranslated regions (Wilkie ), like the secondary structures involved in regulation of histone mRNA translation (Marzluff ), the iron response elements (Anderson ) or the gamma interferon inhibitor of translation element that is involved in limiting the cellular immune response (Mukhopadhyay ). In many viral RNAs, RNA secondary structures in the untranslated regions regulate the viral life cycle (Liu ). Functional RNA secondary structures in viral RNAs have also been reported in the protein-coding regions, like for the MS2 bacteriophage coat protein mRNA (Ball, 1973), for cis-replication elements (CREs) involved in the replication of plus strand RNA viruses (Fricke ; Liu ) or retroviruses (Konecny ) and for the RNA encapsidation signals of betacoronaviruses (Fosmire ) and Hepatitis B Virus (Junker-Niepmann ). However, in the protein-coding regions of cellular mRNAs such RNA secondary structures are known to exert functions only in specific cases, like the selenocysteine insertion sequence involved in the incorporation of the amino acid selenocysteine into proteins (Donovan and Copeland, 2010).

The relationship between codon usage and RNA secondary structures was discussed frequently (Ball, 1973; Carlini et al., 2001; Fitch, 1974; Gu ; Konecny ; Mao ; Oresic ; Shabalina ). The question arises if mRNA protein-coding regions in general have functional RNA secondary structure elements and if there is a relationship between the evolution of protein sequences and the nucleotide sequences forming RNA secondary structures. Three hypotheses have been considered (Ball, 1973): (i) codon sequence and secondary structure are independent, i.e. only the codon sequence is under selection pressure solely due to the needs of the amino acid sequence in the functional protein, and the most stable RNA secondary structure is determined only according to the requirements of the codon sequence; (ii) natural selection between synonymous codons which encode the same amino acid could permit the optimization of RNA secondary structures (Fitch, 1974) and (iii) a strong selection pressure for functional RNA secondary structures could drive the choice also of non-synonymous codons in the mRNA, leading to the choice of another amino acid according to the requirements of the RNA secondary structure (Konecny ).

2 Materials and methods

2.1 Prerequisites

For a given nucleotide sequence of length n, we labeled each nucleotide position i in a coding region with the corresponding codon index , calculated as . We defined for each base-pair a base-pair index . We further group the base-pair indexes into base pair/secondary structure index classes c1={1:1, 2:3, 3:2}, c2={2:2, 1:3, 3:1} and c3={3:3, 1:2, 2:1}.

We used a pairwise t-test with Benjamini and Hochberg P-value adjustment as significance test. The P-values are encoded as following: ***, ** and *.

2.2 Data

We downloaded all genomes and annotations of the NCBI RefSeq genome database (May 2016, ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq) and (in case of multiple isoforms) randomly selected one isoform of each protein-coding gene sequence. We discarded all sequences with . We analyzed mitochondrial and chloroplast mRNAs separately. Due to runtime reasons, we reduced the dataset. We randomly selected 200.000 CDSs of each group (if available), with a balanced number of CDS per species per group. The mitochondrial and chloroplast mRNAs are downloaded from the NCBI nucleotide database with the following search term: biomol_genomic[PROP] AND refseq[filter] AND mitochondrion[filter] or with the filter chloroplast[filter] (July 2017, www.ncbi.nlm.nih.gov/nuccore). The mitochondrial mRNAs are separated by the taxa plants, protists, fungi and animals.

2.3 RNA secondary structure prediction of CDS

We performed an RNA secondary structure prediction with RNAfold (Lorenz ) (using the default parameters) for all species and each CDS. We used a sliding window approach with a window size of 50 nt and a step size of 5 nt to be fast (for runtime reasons) and to move the edge of the window over all possible frames. For the predicted RNA secondary structure of each window, we counted the frequency of the base pair/RNA secondary structure classes. For the viral CDS, we additionally used window sizes of 25, 50, 100, 200 and 400 nt to exemplary exclude a window specific bias. We are aware of the fact that RNAfold produces false positive secondary structure predictions, but we assume false positives to be equally distributed over all three classes.

2.4 Sequence shuffling

The di-nucleotide shuffling of viral CDS was performed with uShuffle (Jiang ). For the mono-, codon-, and di-codon-nucleotide shuffling, we used a simple in-house script (available on request). The codon and di-codon-nucleotide shuffling will persevere the codon frequency of the underlying mRNA sequences.

3 Results and discussion

We provide a global-scale analysis of the crosstalk between amino acid coding requirements and RNA secondary structure formation in protein-coding sequences. Using RNAfold (Lorenz ), we analyzed the occurrence of three different base pair/secondary structure classes, based on the arrangement of amino acid codons within RNA secondary structures. Each arrangement of valid base-pairs (A:U, U:A, G:C, C:G, U:G, G:U) between the codon positions 1, 2 or 3 of one codon and the codon positions 1, 2 or 3 of an opposing codon in RNA secondary structures can be grouped into the following classes: c1={1:1, 2:3, 3:2}, c2={2:2, 1:3, 3:1} and c3={3:3, 1:2, 2:1} (Fig. 1). Thereby, these sequence stretches do not need to have a length in multiples of three nucleotides to belong to the respective class, as illustrated in an artificial RNA structure (Fig. 1B). Consequently, these three base pair/RNA secondary structure classes cover virtually all possible RNA secondary structures. The frequency of these base pair/secondary structure classes was counted on a global scale in the protein-coding regions of mRNAs from viruses, bacteria, archaea, plants, fungi, invertebrates, non-mammalian vertebrates and mammalia.

Fig. 1.

(A) Three base pair/RNA secondary structure classes that represent all possible base-pairs between positions 1, 2 and 3 of one codon with positions 1, 2 and 3 of an opposing codon in RNA secondary structures. The classes are based on the following groups c1={1:1, 2:3, 3:2}, c2={1:3, 2:2, 3:1} and c3={1:2, 2:1, 3:3}. (B) An artificial example RNA secondary structure with the different classes corresponding to A. Shifts between the different classes are possible, introduced through loops and bulges

As a result, we found base-pairs/RNA secondary structures of class c1 are significantly under-represented globally in the mRNA protein-coding regions in all taxa of life (Fig. 2).

Fig. 2.

Distribution of class c1, c2 and c3 within all selected coding sequences (sCDS, see methods for selection procedure) except mitochondria and chloroplasts for each taxonomic group. #species – number of different species; #CDS – number of coding sequence (CDS); #sCDS – number of selected CDS; #nt/sCDS – number of nucleotides per selected CDS

In class c1, two opposing nucleotides which each are in position 1 of two opposing codons would be required to pair. However, the nucleotide in codon position 1 largely defines the exact nature of an amino acid and by that its function in the protein. Thus, this arrangement would require two specific amino acids in a protein to co-evolve to allow RNA secondary structure formation, which appears highly unlikely. Conversely, the fact of an under-represented class c1 is interpreted as a selection pressure to allow RNA secondary structures in mRNA coding regions. At the same time, pairing of the nucleotides in position 2 of a codon with a nucleotide in position 3 of another codon in class c1 may be highly likely due to the degenerate nature of triplett position 3.

To ensure the under-representation of class c1 is not an artifact of the RNA folding algorithm, we performed several tests. The under-representation of class c1 is not triggered by a single specific base-pair since all base-pairs of class c1 are under-represented (Supplementary Fig. S1). The bias is not affected by the folding range of our method (see Supplementary Material and Supplementary Fig. S2) and occurs not by chance (Supplementary Fig. S4). We can also find the same effect if we use other folding algorithms, such as CONTRAfold (Do ), a secondary structure prediction method based on conditional log-linear models (data not shown). Furthermore, if we consider the complete ensemble of RNA secondary structures, we find the same bias (see Supplementary Material and Supplementary Fig. S3). The effect disappears after artificial shuffling of the underlying RNA sequences with a word size of one or two, but remains with a shuffling of word size three or six (conserving the underlying codons and di-codons, Supplementary Fig. S4), showing the effect being not specific to the nucleotide distribution (like GC content) but specific for the underlying codon distribution. Moreover, if we build artificial codons but keep the nucleotide distribution per codon position constant, the effect remains the same (Supplementary Fig. S5). This demonstrates that the global nucleotide distribution per codon position but not the codon itself is important.

In addition, we generally also find a less pronounced but visible under-representation of the pairing of codon position 2 with position 2 of an opposing codon in class c2 (Supplementary Fig. S1). This illustrates also codon position 2 to be important, likely because Uracil in position 2 defines the hydrophobic character of an amino acid. A change of the hydrophobic versus hydrophilic nature of amino acids could result in general misfolding of proteins. However, our results shown here imply that codon position 2 has a smaller influence on the amino acid choice than codon position 1.

Although changes in codon position 2 are generally believed to be more important for protein structure and function than changes in codon position 1 (Berleant ; Pohlmeyer, 2008), our results challenge this general assumption. In this context, it is important to note that in previous studies in which the physico-chemical properties of amino acids like size, hydropathy and charge were considered, charged amino acids were usually classified by simply having a charged side chain or not. However, the polarity of charge and the effects of changing this polarity were not taken into account. As a consequence, changes in amino acid charge from positive to negative or vice versa were usually assigned only low penalties, other than changes in amino acid hydropathy and size (Berleant ; Goodarzi ; Miyata ; Nemzer, 2017).

However, our results suggest that codon position 1 is globally more important than codon position 2. Therefore, we inspected the genetic code to find out which changes in amino acid properties can occur after changes in codon position 1 or position 2. We found that charge reversal (i.e. the change of a positive to a negative charge or vice versa) only occurs upon single nucleotide exchanges in codon position 1 (G to A, G to C or vice versa) but never upon single nucleotide exchanges in codon position 2. We conclude that the selection pressure against charge reversal—which may destroy salt bridges in proteins that are important for the protein’s structure and function, like the His-Asp salt bridge that mediates the Bohr effect in hemoglobin (Russu )—may have been largely underestimated in previous studies.

Globally, structured mRNAs have a more pronounced base-pair bias than less structured mRNAs. Compared to the low structured mRNAs (Fig. 3A), we found a slight increase of cytosin and guanin content at codon positions 1 and 2 in the structured mRNAs but a strong increase at codon position 3 (Fig. 3B). Based on the higher G/C content, the G:C base-pairing was significantly increased. There seems to be a consistent trend over all organisms that the increase in G/C content among the base-pairs of one class occurs to a large extent in those base-pairs which are involved in a pairing with a codon position 3 (Supplementary Fig. S1). The largest differences were found in 1:3, 3:1, 2:3, 3:2 and 3:3 pairings, whereas no big changes were found in 1:1, 2:2, 1:2 and 2:1 pairings (Supplementary Fig. S6). This shows that an increase in G/C content in more structured mRNAs occurred mainly in the wobble position 3, leaving the amino acid sequence largely unchanged. In turn, this suggests that the co-evolution of amino acid sequence and RNA secondary structures may have been mainly driven by the amino acid sequence. Conversely, it appears difficult to estimate at this global scale if the RNA secondary structure influences the amino acid sequence by the choice of non-synonymous codons. Even though the pronounced base-pair bias in all taxonomic groups supports the presence of structured mRNAs in vivo, up to now we are aware of only a few functional RNA structures in mRNA coding regions (Ball, 1973; Fosmire ; Fricke ; Konecny ; Liu ; Nassal ).

Fig. 3.

(A) Class distribution and (B) nucleotide distribution of each codon positions of all human mRNAs with less stable (mean MFE>−8 kcal/mol over all folding windows, #7104 mRNAs, TOP) and rather stable (mean MFE  −8 kcal/mol over all folding windows, #15022 mRNAs, BOTTOM) RNA secondary structures. The nucleotide distributions in codon positions 2 and 3 of the less stable mRNAs are similar, since base-pairings of 2:3 and 3:2 are likely. In contrast, mRNAs with more stable RNA secondary structure show an increased GC content at codon position 1 and in particular at position 3. This increases the G:C and C:G base-pairings 1:3, 3:1 and 3:3 (Supplementary Fig. S6) and consequently increases the frequency of classes c2 and c3

Interestingly, in nearly all mitochondrial (except animal) and chloroplast mRNAs the RNA secondary structure class c2 and class c1 is under-represented, whereas in animals only class c2 but not class c1 is under-represented (Fig. 4). In class c2 RNA secondary structures, pairing of two opposing nucleotides in position 2 of each codon would be required. However, variation in this position would likely change the hydrophobic versus hydrophilic nature of the amino acids. Thus, the under-representation of class c2 in plastid genomes may be due to the fact that plastid genomes largely encode those highly hydrophobic membrane proteins which are involved in plastid redox metabolism to allow for short regulatory gene expression circuits (Allen, 2015). This effect appears particularly pronounced in animal mitochondria which have lost all protein-coding genes that do not encode respiratory chain complex subunits (Lang ). Here, only class c2 but not class c1 is under-represented (Fig. 4). This suggests that in these animal mitochondrial sequences exclusively coding for those subunits of the mitochondrial respiratory chain complexes that are located in the membrane (Soto ; Wirth ), hydropathy has become the most important parameter for codon selection.

Fig. 4.

Base pair/RNA secondary structure class distribution in mitochondria and chloroplasts

Some viral RNA genomes contain known functional RNA secondary structures. For a proof of principle, we therefore analyzed the 3'-region of the Hepatitis C Virus RNA genome which contains RNA secondary structure elements known to be involved in genome replication (Niepmann ). The results show that in two regions containing functional signals (J7880 and J8880) class c1 is strongly under-represented, whereas in the coding sequence (CDS) region containing the so-called CRE and encoding the hydrophobic transmembrane region of the protein, class c2 is strongly under-represented, with high values in class c1 (Supplementary Fig. S7). This demonstrates that under-representation of RNA secondary structure classes c1 or c2 can actually correlate with functional RNA secondary structures.

As a control, we also considered if the codon frequencies as found in nature contribute to the under-representation of RNA secondary structure class c1. Therefore, we calculated the probability that two codon positions can base-pair by just using the codons taken from the table of the genetic code, while taking into account the actual frequency of each codon as found in each taxonomic group in nature. Without taking the underlying codon frequency into account, the possible base-pair distribution of classes c1, c2 and c3 would be identical (33.3% each). In contrast, when we weight each possible base pair between the codons with the underlying codon frequency, we found a slight under-representation of class c1 (see Table 1, left part). This demonstrates that even the codon frequencies are sufficient to show an under-representation of RNA secondary structure class c1, underlining the global importance of RNA secondary structures in CDSs. Importantly, in the naturally occurring CDSs with the order of codons as required by the amino acids in the proteins, the under-representation of RNA secondary structure class c1 comes with an additional strong bias (Table 1, right part).

Table 1.

Class distribution over all possible base-pairs between the 64 codons of the genetic code, weighted with the natural codon frequency of both codons (left). Class distribution based on the RNA secondary structure as depicted in Figure 2 (right)

	Codon frequency based			mRNA structure based
	c1	c2	c3	c1	c2	c3
Archaea	31.7	33.8	34.4	25.2	37.8	37.3
Bacteria	32.1	33.7	34.1	27.0	36.3	36.7
Fungi	32.5	33.7	33.7	30.1	35.2	34.1
Invertebrates	32.8	33.5	33.5	30.3	33.5	36.0
Plants	32.8	33.5	33.5	31.2	33.9	34.9
Protozoa	32.7	33.5	33.6	29.6	34.0	36.1
Mammals	32.7	33.6	33.5	29.8	35.1	35.0
Vertebrates	32.8	33.6	33.5	30.3	34.6	35.0
Viruses	32.4	33.8	33.6	30.2	34.8	34.8

It was shown that RNA secondary structure is under positive selection in coding regions of bacteria (Katz and Burge, 2003), but the conclusions are controversial. On the one hand, stable structures in mRNAs can slow down ribosomes (Chen ; Tholstrup ; Wen ), arguing against a high number of very stable RNA structures in coding regions. On the other hand, a recent study of yeast mRNAs found a positive correlation of mRNA secondary structure with protein abundance, without affecting mRNA half-life (Zur and Tuller, 2012). For this finding, Mao et al. proposed a possible explanation that appears counter-intuitive at first sight. By computational methods they found that moderate RNA secondary structures may increase ribosomal density and therewith translation rate (Mao ). Their interpretation is that the first ribosomes that encounter the coding region are slowed down by RNA secondary structures. Then, all following ribosomes catch up closely to create a line of migrating ribosomes which will disable reformation of RNA secondary structures and thus translate the mRNA with high efficiency.

Taken together, synonymous codons in mRNA protein-coding regions are globally selected in all taxa of life in a way that the need to form RNA secondary structures in which a codon position 1 is required to pair with the codon position 1 of another codon is reduced. In turn, this means that mRNA CDSs evolved to allow the formation of RNA secondary structures using synonymous codons, underlining the general importance of such RNA secondary structures for the function of mRNAs in cells. However, this does not exclude that in specific cases conserved functional RNA structures may also drive co-evolution of amino acid sequences using non-synonymous codons.

Funding

This work was supported by the German Research Foundation (DFG) [SFB 1021], AquaDiva [INST 275/367-1] and FungiNet [INST 275/365-1].

Conflict of Interest: none declared.

Click here for additional data file.