Evolutionary selection on synonymous codons in RNA G-quadruplex structural region

Xu Yuming Qi Ting Gu Wanjun Lu Zuhong

(School of Biological Science and Medical Engineering, Southeast University, Nanjing 210096, China)(State Key Laboratory of Bioelectronics, Southeast University, Nanjing 210096, China)

Abstract：To investigate how synonymous codons have been adapted to the formation of ribonucleic acid (RNA) G-quadruplex (rG4) structure, a computational searching algorithm G4Hunter was applied to detect rG4 structures in protein-coding sequences of mRNAs in five eukaryotic species. The native sequences forming rG4s were then compared with randomized sequences to evaluate selection on synonymous codons. Factors that may influence the formation of rG4 were also investigated, and the selection pressures of rG4 in different gene regions were compared to explore its potential roles in gene regulation. The results show universal selective pressure acts on synonymous codons in rG4 regions to facilitate rG4 formation in five eukaryotic organisms. While G-rich codon combinations are preferred in the rG4 structural region, C-rich codon combinations are selectively unfavorable for rG4 formation. Gene’s codon usage bias, nucleotide composition, and evolutionary rate can account for the selective variations on synonymous codons among rG4 structures within a species. Moreover, rG4 structures in the translational initiation region showed significantly higher selective pressures than those in the translational elongation region.

Key words：ribonucleic acid (RNA) structure; G-quadruplex; synonymous codons; evolution; selection

mRNAs of protein-coding genes in eukaryotic organisms carry the genetic information to encode amino acid sequences and contain multiple regulatory and structural signals[1]. The pivot provision enabling mRNA to hold these regulatory functions is the redundancy of the genetic code that allows for many “silent” mutations at synonymous codon sites[2]. Synonymous nucleotide mutations do not change the amino acid sequences of the encoded proteins; however, they can confer dramatic differences to the structure and functions of mRNA[1-3]. Much evidence has shown that synonymous codons are selected for optimized RNA stability[4-5], proper nucleosome positioning[6], efficient mRNA splicing[7-8], correct microRNA targeting[9], efficient translation initiation[10-12]and elongation[13], and proper protein co-translational folding[14-15].

The RNA G-quadruplex (rG4) is a non-canonical super-secondary structure around the G-rich sites in RNA sequences, which consists of the stacking of G-quartets formed by the G-G Hoogsteen hydrogen bonding[16]. Several experimental studies have shown that rG4 is ubiquitous in untranslated regions and protein-coding sequences (CDSs) of eukaryotes[17]. RG4 structures have been shown to perform many diverse and vital functions in a wide range of biological processes, such as pre-mRNA splicing[18], alternative polyadenylation[19], mRNA localization[20-21], microRNA targeting[22], and translational regulation[23-28]. Due to the presence of the 2α-hydroxyl property, the rG4 structure is more stable than DNA G-quadruplex[29]. In a recent study, Mirihana Arachchilage et al[30]have demonstrated that the most stable G4s appeared to be significantly under-represented within the CDS using specific synonymous codon combinations. However, several key problems regarding synonymous codon usage, rG4 formation, and evolution remain unaddressed.

In this study, the evolutionary choices of synonymous codons were evaluated around putative rG4 structures (pG4) at the whole transcriptome scale in multiple eukaryotic species. The searching algorithm G4Hunter was applied to detect rG4 structures in the protein-coding sequences of mRNAs in five eukaryotic species. Then，the native sequences forming rG4s were compared with randomized sequences to evaluate the selection of synonymous codons. Factors that may influence the formation of rG4 were also investigated, and the selection pressures of rG4 in different gene regions were compared to explore its potential roles in gene regulation. These analyses may help describe the evolutionary selections acting on synonymous codons in protein-coding regions.

1 Materials and Methods

1.1 Data

The nucleotide sequences and exonic structures of all protein-coding genes were downloaded in five eukaryotic species, includingH.sapiens(GRCh38.p13),M.musculus(GRCg6a),G.gallus(GRCm38.p6),D.rerio(GRCz11), andD.melanogaster(BDGP6.28), using Ensembl BioMarts (release 97)[31]. Only protein-coding genes with a coding sequence of more than 150 nucleotides were included. Furthermore, miRNA target sites in the protein-coding regions of human and mouse genomes were downloaded from the miRDB database[32-33].

To explore the factors that affect the selection for rG4 structure formation, gene codon usage bias, nucleotide composition, and its evolutionary rate were considered. The effective number of codons (ENC) was used to measure codon usage bias, and ENC values for each gene were calculated as suggested by Wright[34]. A lower ENC value indicates stronger codon bias[34]. For each gene, we also calculated nucleotide compositions, including G and C contents. Moreover, theDnandDsvalues of all human and mouse orthologous genes were also retrieved from Ensembl BioMarts[35]. As a popular indicator of selection acting on protein-coding sequences,Dn/Dsquantifies the mode and strength of selection by comparing synonymous substitution rates (Ds) with nonsynonymous substitution rates (Dn)[36]. Generally,Dn/Dsclose to one indicates neutrality; values greater than one are interpreted as positive selection (selection promoting change), and values less than one usually indicate purifying selection (selection suppressing protein change). Here, the ratio ofDnandDsvalues for each gene are calculated and the value ofDn/Dsis used as the measurement of the evolutionary rate of genes in mice and humans.

1.2 rG4 structure in the protein-coding region

To locate rG4 structures in the protein-coding region, the G4Hunter algorithm[37]was exploited to systematically search for potential rG4 forming sites in the protein-coding sequences in all five species. G4Hunter considers G-richness and G-skewness of a given sequence and presents a quadruplex formation propensity score, G4Hscore, as output. G4HunterApps was run with a window size at 25 nt and a cutoff at 1.2[38]to identify all potential rG4 structures (pG4) in all protein-coding sequences. G4Hunter algorithm was chosen with these two parameters, since a comprehensive evaluation of computational methods for rG4 prediction has suggested that G4Hunter has the best performance in predicting G4 structures[39]. 65 562, 42 153, 38 311, 25 194, and 14 184 rG4 structures were identified in protein-coding sequences forH.sapiens,M.musculus,G.gallus,D.rerio, andD.melanogaster, respectively. To validate the observations of predicted rG4 structures, experimentally detected rG4 sites were also downloaded from the supplemental materials of Guo et al[40]. The experimental rG4 data were achieved by high throughput RT-stop techniques in humanHEK293Tcells and mousemESCcells.

1.3 mRNA randomization

If the choice of synonymous codons influences the formation of rG4 structures in coding sequences, the G4Hscore of mRNA sequences in the naturally occurring pG4 region should be statistically different from that of randomized sequences. Thus, synonymous codons in the coding sequence were randomly shuffled, keeping the same amino acid sequence, GC composition, and codon usage bias. For each CDS sequence, the shuffling process was repeated 1 000 times to obtain a set of randomized artificial sequences. The G4Hscore of each 30 nt window in the native CDS sequence and each permutated sequence were calculated using the G4Hunter algorithm[37]. The difference of G4 forming potential between the native sequence and the randomized sequences was determined by calculating theZ-score of the G4Hscore (ZG4S) for each sliding window using the following formula

(1)

Similarly, the difference between the G (or C) compositions of the native sequence and randomized sequences was evaluated. TheZ-score of the G content (ZG) and that of the C content (ZC) for each sliding window can be calculated as follows:

(2)

(3)

2 Results

2.1 Selection of rG4 structures on synonymous codons

The rG4 forming propensity score, G4Hscore[37], was calculated along the mRNA sequences with a sliding window scheme. For each pG4 structure in the transcriptome, a 30 nt window was moved both upstream and downstream from the start position of the pG4 structure with a 30 nt step, and the G4Hscore of the sequence in each of the 13 windows was calculated. To estimate the background distribution of the formation propensity of rG4 structures, the mRNA sequences were randomized by shuffling synonymous codons 1 000 times, and the G4Hscore in the corresponding sliding windows was calculated. The G4Hscore of the real mRNA sequence in a sliding window was compared with that of 1 000 corresponding sliding windows in the shuffled sequences.ZG4Swas calculated to assess the deviation of rG4 formation in the observed mRNA sequence from a random expectation. A positiveZG4Svalue means synonymous codons are selected to facilitate the formation of rG4 structures, while a negativeZG4Svalue indicates a selective pressure that prevents the formation of rG4 structures.

The sliding window analysis was performed in five eukaryotic species, includingH.sapiens,M.musculus,G.gallus,D.rerio, andD.melanogaster. A significantly positiveZG4Svalue was observed in the window of pG4 structures in all five species (see Fig.1). When sliding windows move to the upstream or downstream of the pG4 structures,ZG4Svalues drop quickly in the flank region of pG4 structures and oscillate around zero for sliding windows away from the pG4 structures (see Fig.1). When in vitro experimentally identified rG4 structures in humanHEK293Tcells and mousemESCcells were analyzed by the same procedure, a similar pattern ofZG4Swas found changing along the sliding windows (see Fig.2).

(a) (b) (c)

(d) (e)

(a) (b)

2.2 Selection for rG4 formation on the GC content of codon combinations

To investigate how synonymous codons are selected for rG4 formation,ZGandZCof a 30 nt window were calculated for each pG4 structure. Fig. 3 shows a significant positive correlation betweenZGandZG4Sfor pG4 structures in all five species. In comparison, a weaker but significant negative correlation betweenZCandZG4Sof pG4 structures was also observed in all five species (see Fig.4).

(a) (b) (c)

(d) (e)

(a) (b) (c)

(d) (e)

2.3 Features of the host gene that associate with rG4 selection

Although the mean value ofZG4Sis significantly larger than zero for all exonic pG4 structures in all five species (see Fig. 1), there are substantial variations among different pG4 structures within a single organism (see Figs. 3 and 4). To explore the factors that affect the selective pressures on synonymous codons for rG4 formation, several putative gene-level factors were considered, including the codon usage bias, evolutionary rate, and nucleotide compositions of the host gene where the pG4 structure is located. Analyses showed thatZG4Svalues of pG4 structures in genes with the highest 5% ENC are significantly higher (p=1.2×10-12in humans andp=2.8×10-11in mice) than those genes with the lowest 5% ENC values (see Fig. 5(a)). For the top 5% genes with the highestDn/Dsratio,ZG4Svalues of pG4 structures are significantly larger (p=4.1×10-3in humans andp=4.3×10-6in mice) than those pG4 structures in the bottom 5% genes with the lowestDn/Dsratio (see Fig. 5(b)). When pG4 structures in genes with the highest 5% G content are compared to those with the lowest 5% G content, it showed that pG4 structures located in genes with the highest 5% G content had significantly smaller (p<2×10-16in humans andp<2×10-16in mice)ZG4Svalues (see Fig. 5(c)). Moreover,ZG4Svalues of pG4 structures in genes with the top 5% C content are also statistically smaller (p=8.5×10-14in humans andp=1.9×10-5in mice) than those in genes with the bottom 5% C content.

(a) (b)

(c)

2.4 Different selective pressure for rG4 structures in different gene regions

Other than the features of the host genes,ZG4Sdifferences of rG4 structures in different gene regions were also evaluated. First, pG4 structures were grouped into two categories: those in the translation initiation region (within 70 nt downstream of the start codon) and those out of the translation initiation region. The results showed that pG4 structures in the translation initiation region had significantly higherZG4Svalues than those in the translation elongation region (see Fig. 6). Next,ZG4Svalues of pG4 structures near exonic splicing sites (within 60 nt of splicing sites) were compared with those of pG4s distant from the splicing sites. AlthoughZG4Svalues for pG4 structures near the splicing sites tend to be smaller than those distant from the splicing sites (data not shown), the differences are not statistically significant for rG4 structures in both downstream and upstream flank regions of splicing sites. Finally,ZG4Svalues of pG4 structures in microRNA (miRNA) target sites were compared withZG4Svalues of pG4s out of miRNA target region. No significant differences were observed betweenZG4Svalues of pG4 structures in miRNA target region and those out of miRNA target sites.

(a) (b)

3 Conclusions

1) Synonymous codons are universally selected for the formation of rG4 structures in the protein-coding sequences of the five eukaryotic organisms under investigation. While G-rich codon combinations are preferred in the rG4 structural region, C-rich codon combinations are selectively unfavorable for rG4 formation.

2) The selective pressures acting on synonymous codons are stronger for pG4 structures in genes with lower G nucleotides, less biased usage of synonymous codons, and a lower evolutionary rate. These differences are consistent for pG4 structures in humans and mice.

3) Synonymous codons in some specific gene regions, such as the translation initiation region, were under different selective pressures for the rG4 formation. However, rG4 structures in the miRNA target region and splicing sites flank region did not show obviously different evolutionary selections on synonymous codons.