Analysis of synonymous codon usage patterns in sixty-four different bivalve species

Data sources

We considered two bivalve mollusk species with a fully sequenced genome (C. gigas and P. fucata) (Takeuchi et al., 2012; Zhang et al., 2012) and 62 other species whose transcriptome has been sequenced using next generation sequencing technologies and deposited in public sequence databases. When both 454 Life Sciences and Illumina-generated sequencing reads were available for a same species, the latter were chosen due to higher throughput and lower rate of sequencing errors. Namely, Illumina reads were used for species Anadara trapezia, Arctica islandica, Argopecten irradians, Astarte sulcata, Atrina rigida, Azumapecten farreri, Bathymodiolus platifrons, Cardites antiquata, Cerastoderma edule, Corbicula fluminea, Crassostrea angulata, Crassostrea corteziensis, Crassostrea hongkongensis, Crassostrea virginica, Cycladicama cumingii, Cyrenoida floridana, Diplodonta sp. VG-2014, Donacilla cornea, Elliptio complanata, Ennucula tenuis, Eucrassatella cumingii, Galeomma turtoni, Glossus humanus, Hiatella arctica, Lampsilis cardium, Lamychaena hians, Mactra chinensis, Margaritifera margaritifera, Mercenaria campechiensis, Meretrix meretrix, Mizuhopecten yessoensis, Mya arenaria, Myochama anomioides, Mytilus californianus, Mytilus edulis, Mytilus galloprovincialis, Mytilus trossulus, Neotrigonia margaritacea, Ostrea chilensis, Ostrea edulis, Ostrea lurida, Ostreola stentina, Pecten maximus, Perna viridis, Pinctada martensii, Placopecten magellanicus, Polymesoda caroliniana, Pyganodon grandis, Ruditapes decussatus, Ruditapes philippinarum, Sinonovacula constricta, Solemya velum, Sphaerium nucleus, Uniomerus tetralasmus and Villosa lienosa (Qin et al., 2012; Ghiselli et al., 2012; Chen et al., 2013; Meng et al., 2013; Gerdol et al., 2014; Fu et al., 2014; Zhang et al., 2014; Pauletto et al., 2014; De Sousa et al., 2014; Zhao et al., 2014; Cornman et al., 2014; Zapata et al., 2014; Prentis & Pavasovic, 2014; González et al., 2015). The 454 Life Sciences sequences were used for Bathymodiolus azoricus, Geukensia demissa, Laternula elliptica, Mimachlamys nobilis, Pinctada maxima, Saccostrea glomerata, and Tegillarca granosa (Clark et al., 2010; Philipp et al., 2012; Egas et al., 2012; Jones et al., 2013; Fields et al., 2014). Details about the data used for the different species are provided in Table S1.

Sequence data were processed as follows: predicted CDS from the fully sequenced genomes of C. gigas (release 9) and P. fucata were retrieved from http://oysterdb.cn and http://marinegenomics.oist.jp/pinctada_fucata, respectively. De novo transcriptome assemblies were performed for all the other 62 bivalve species with the CLC Genomics Workbench (v.7.5, CLC Bio, Aarhus, Denmark) using the de novo assembly tool with “automatic word size” and “automatic bubble size” parameters selected, and setting the minimum allowed contig length to 300 bp.

In all transcriptomes, ORFs (Open Reading Frames) longer than 100 codons were predicted with TransDecoder (http://transdecoder.sourceforge.net). We selected the predicted CDS of C. gigas, and of one representative species for the Imparidentia (R. decussatus), Protobranchia (S. velum) and Palaeoheterodonta (P. grandis) lineages to identify a subset of evolutionarily conserved protein-coding genes with a 1:1 orthology ratio across Bivalvia. This was achieved by performing reciprocal tBLASTx searches (the e-value threshold was set a 1 × 10⁻¹⁰ and only hits displaying sequence identity >50% were considered). This procedure resulted in a selection of 2,846 conserved protein-coding genes, whose orthologous sequences were retrieved in the remaining 60 species. Due to the heterogeneous tissue and developmental stage origin, the different sequencing platforms and depth applied, several of these evolutionarily conserved sequences could not be identified or were fragmented in some transcriptomes. In order to ensure a minimum quality criteria, all the selected species had to display at least 25% of the sequences included in the dataset of evolutionarily conserved genes, with an average length >500 nucleotides. A number of additional transcriptomes derived from publicly available data did not meet such criteria and were therefore not included in our analyses (Table S2).

Codon frequencies and codon usage statistics

The sets of evolutionarily conserved genes retrieved for each species were individually processed with the cusp tool of the EMBOSS package (Rice, Longden & Bleasby, 2000) obtaining codon frequencies and GC composition for each codon position. RSCU values for each individual codon were calculated for each species as described by Sharp and colleagues (1986). The effective number of codons (ENC) for each species was calculated according to Wright (1990) using the EMBOSS chips tool, summing codons over al sequences (Rice, Longden & Bleasby, 2000). The sENC-X values were determined for every amino acid for each species and scaled to a range of values between 0 and 1 according to Moriyama & Powell (1997). EMBOSS chips was also used to calculate ENC for individual genes whenever necessary. We identified a reference set of 50 highly expressed genes for the calculation of Codon Adaptation Index (CAI) based on the average expression in C. gigas digestive gland (SRA:SRX093412), gills (SRA:SRX093414) and hemocytes (SRA:SRX093417) RNA-seq libraries and their inclusion in the above mentioned set of 2,846 genes conserved across bivalves. Gene expression was calculated as TPM (Transcripts Per Million) (Wagner, Kin & Lynch, 2012), with the RNA-seq mapping tool included in the CLC Genomics Workbench 8.5 (Aarhus, Denmark), setting length and similarity fraction parameters to 0.75 and 0.98 and insertion/deletion/mismatch penalties to 3. Orthologous genes were used for CAI calculation in other species. CAI values were computed with CAI calculator 2 (Wu, Culley & Zhang, 2005). The gene expression levels of M. galloprovincialis transcripts were calculated using the digestive gland (SRA:SRX126945-8), gills (SRA:SRX389466) and hemocytes (SRA:SRX389338) RNA-seq libraries (Gerdol et al., 2014; Moreira et al., 2015).

Scatter plots were generated between ENC and the average GC content calculated at the third codon position (GC₃) for each species, between ENC and sENCx and between ENC and CAI; Paerson correlation coefficients and linear regression analyses were computed with R 3.1.0 (http://www.r-project.org).

Hierarchical clustering analysis

The Relative Synonymous Codon Usage (RSCU) values calculated for the 59 informative codons in each species were used to build a tabular file. STOP, ATG (encoding Met) and TGG (encoding Trp) codons were excluded from this analysis. This file was used as an input for Cluster 3 (De Hoon et al., 2004), thus generating a species distance matrix. Hierarchical clustering was performed by using the Euclidean distance as a similarity metric and complete linkage as a clustering method.

CUB varies across bivalve species

The codon frequencies and RSCU values calculated for the 64 bivalve species analyzed are displayed in Tables S3 and S4. As shown in the comparative overview of Fig. 1, RSCU values are usually very similar in closely related species, such as in the case of Mytilida, Ostreoidea and Unionida (A, B and C), but marked differences can be observed in an higher-order comparison (D).

Different species clearly show a different tendency to the preferential usage of specific codons, as exemplified by the average Effective Number of Codons (ENC) (Wright, 1990) values in Table 1. Overall, the observed ENC values range between 40.65 (in the Chinese surf clam M. chinensis) and 56.80 (G. turtoni) across the analyzed species, while the theoretical value is comprised between 21 (if only a single codon is used for each amino acid) and 61 (if all codons are used with equal frequency). Most bivalve species display a weak CUB, using on average over 50 out of the 61 available codons, and only a limited number of bivalve species display an ENC value comparable to that of other invertebrates (the ENC range is 45–48 in Drosophila and nematodes) (Powell & Moriyama, 1997; Mitreva et al., 2006; Vicario, Moriyama & Powell, 2007).

Species clustering based on CUB does not reflect the evolutionary history of Bivalvia

We computed RSCU values for the 59 informative codons of each species to perform a hierarchical clustering of species with Cluster 3, thereby investigating the role of CUB in the evolution of bivalve genomes. The resulting dendrogram is shown in Fig. 2. Although phylogeny and CUB-based clustering are in agreement, in several cases, up to the order level, the six major lineages of Anomalodesmata, Archiheterodonta, Imparidentia, Palaeoheterodonta, Protobranchia and Pteriomorphia expected from molecular phylogeny (Bieler et al., 2014) are hardly distinguishable. This observation is consistent with data previously reported for nematodes by Mitreva and colleagues (2006), who identified a connection between codon distribution and phylogeny only for closely related species, up to the genus level. In bivalves such a relationship seems to extend a bit further, in some cases up to one of the six major evolutionary lineages, as for instance in Palaeoheterodontha, which include the freshwater mussels of the order Unionida and the saltwater clams of the order Trigoniida, or Archiheterodonta, which only comprise four relatively small extant families. While in some cases (e.g., Mytilida and Ostreoidea) all the species maintain a similar usage of synonymous codons (Fig. 1), in others (e.g., Venerida) remarkable differences among species are clearly visible.

In essence, the clustering based on codon usage divides the bivalve species into two largely divergent groups:

(I) The first group is very heterogeneous, comprising 43 species with ENC > 52 (with the exception of L. hians). Two subgroups are detectable: group I-a comprises all the Pectinida and Ostreoidea species (Pteriomorphia), two Anomalodesmata (M. anomiodes and L. elliptica), two Protobranchia (E. tenuis and S. velum) and the Imparidentia G. turtoni, M. arenaria and S. constricta. The species pertaining to this subgroup show a weak CUB, with ENC 54-58 and GC₃ very close to 50% (averaging ∼49%).

The subgroup I-b comprises Unionida and Trigoniida (Palaeoheterodonta), Pinctada spp. (Pteriomorphia, Pterioidea), eight unrelated Imparidentia and the three Archiheterodonta species. Compared to subgroup I-a, the observed CUB is slightly higher, including species with ENC ranging from ∼52 to ∼54, but GC₃ is remarkably lower than group Ia, averaging ∼43%.

(II) The second major group comprises the remaining 20 species, including the Pteriomorphia groups Mytilida, Arcida and Pinnoidea (represented by the lone species A. rigida), together with various Imparidentia species classified as Cardioidea, Cyamiidae, Mactroidea, Sphaeriidae and Venerida, implying that this well-defined group of bivalves with high CUB (ENC ∼40–52) and low GC₃ (30–40%, with the outlier M. chinensis reaching ∼23%) comprises phylogenetically distant species.

Codon usage is biased towards A/T ending codons in Bivalvia

Here we report the most commonly used codon(s) for each amino acid, designing as “preferred codons” all the codons used more frequently than expected (RSCU > 1), and we investigate the correlation between their frequency and overall codon usage bias (negative ENC) across the bivalve species analyzed (Vicario, Moriyama & Powell, 2007).

Figure 3A summarizes the number of bivalve species where a given codon is preferred, evidencing that, despite the difference in overall ENC values, several preferred and avoided codons are shared by most, if not by all, bivalves: this is the case, for example, of ACA (Thr), AGA and AGT (Arg), TCA (Ser), TGT (Cys), CCA (Pro), GCT (Ala) and GGA (Gly). On the other hand, the RSCU values of ACG (Thr), AGC (Ser), TCG and CTA (Leu), TGC (Cys), CCG (Pro), CGC and CGG (Arg), GCG (Val) and GGG (Gly) are always lower than 1, indicating that these codons are avoided in all species. In general, A/T-ending codons appear to be preferred over those ending in G/T, but some notable exceptions exist, including the two C-starting codons encoding the six-fold degenerate amino acids Ser and Arg. While G-ending codons are not uncommon, C-ending codons are almost invariably avoided.

However, when the correlation between codon frequencies and overall CUB of bivalve species is taken into account, the important weight of A/T ending codons on bivalve codon bias becomes evident, (Fig. 3B). Indeed, the high CUB of all the species pertaining to clustering group II (Fig. 2) appears to be mostly resulting from an increased use of A/T-ending codons over those ending in G/C, also explaining the significant negative correlation between ENC and GC₃ across species (see ‘Effects of mutational bias and selection on CUB in Bivalvia’ below). This correlation is significant for most codons with, once again, the exception of the C-starting codons encoding the six-fold degenerate amino acids Leu and Arg. On the contrary, the frequency of G/C ending codons is significantly and negatively correlated with CUB in all cases, with the exception of TTG (Leu).

We explored in detail the contribution of different amino acids to CUB in bivalves by calculating s-ENCx values for each amino acid and correlating this parameter to the overall CUB of each species. This measure is a variation of ENC (Wright, 1990) which is scaled in a range from 0 to 1 for each amino acid independently from the level of redundancy, and which can be used to estimate the relative intensity of CUB across the 18 degenerate amino acids (Moriyama & Powell, 1997). As reported in Table S5 and Fig. 4, the sENC-x values of all amino acids negatively correlate with ENC with significant p-values, including those with relatively low average s-ENCx values. The only exception is represented by Gln, which is likely related to the fact that it is the only two-fold degenerate amino acid to display a strong preference for a G/C-ending (CAG) over an A/T-ending codon (CAA) (see Fig. 3A). Overall, Arg is certainly the amino acid which accounts for the greatest CUB in bivalves, as highlighted by the high average s-ENCx value (0.32), followed by Pro, Thr, Cys, Ala, Gly, Ser and Leu, all characterized by values >0.1.

Effects of mutational bias and selection on CUB in Bivalvia

Although the knowledge of the factors underlying the preferential use of certain codons has remarkably increased over the past few decades, the real contribution of the several potentially contrasting forces involved in shaping CUB still remains a matter of debate in the scientific community.

Mutational bias, intended as the non-randomness of mutational patterns, certainly has a major role in determining codon usage in prokariotes, as well as in many eukaryotes (Muto & Osawa, 1987; Sémon, Lobry & Duret, 2006). A critical parameter tightly linked to mutational bias is genomic GC content, which in turn strongly influences the coding GC content (Knight, Freeland & Landweber, 2001) and GC₃ (Mitreva et al., 2006). Actually, all amino acids, with the exception of the non-degenerate methionine and tryptophan, tolerate A/T to C/G changes in the third codon position and thus, in organisms where mutational bias is dominant over other factors, the occurrence of G/C-ending codons should follow the total genomic GC content. However, even in such cases some deviations are observed, in particular for the 6-fold degenerated Arg and Leu, which also tolerate synonymous mutations at the first codon position (Palidwor, Perkins & Xia, 2010).

In bivalves, average GC₃ appears to be tightly related to ENC (Fig. 5), as these two variables could be strongly correlated by linear regression with R² = 0.92 (p-value =2.86 × 10⁻³⁶). Overall, while GC₃ can assume values as low as ∼23 in M. chinensis, GC₃ and AT₃ contents become even in species displaying weak CUB (and high ENC). Based on these observations and the well-documented correlation between genomic GC content and coding GC₃, one could hypothesize that the bivalves pertaining to the clustering group II possess a similar genomic GC content, lower than all the bivalve species included in the heterogeneous clustering group I, thus explaining their common placement in a well-distinct cluster regardless of the reported phylogenetic distance (Fig. 2). On the contrary, the bivalve transcriptomes of the clustering group Ia show a GC₃ not far from 50%, accompanied by ENC values higher than 54 (Fig. 5), which would suggest that the absence of a strong CUB is linked to a lack of mutational bias.

However, the genomic data currently available for bivalves do not support this view, since the differences in genomic GC content among species are minimal. The fully sequenced genomes of P. fucata and C. gigas possess a relatively low G/C composition (33.69% and 32.33%, respectively) (Takeuchi et al., 2012; Zhang et al., 2014) and partial data from A. farreri and M. galloprovincialis indicate a similar GC content for these two species (35.75% and 31.65%, respectively) (Zhao et al., 2012; Nguyen, Hayes & Ingram, 2014). Therefore, while the A/T-rich nature of most bivalve genomes is consistent with a role of mutational bias towards the preferential choice of A/T-ending codons in bivalves (‘Codon usage is biased towards A/T ending codons in Bivalvia’), other factors might be taken into account to explain the differences in CUB among these marine organisms.

Among these, selection for translational speed and accuracy are certainly among the most relevant, as they can potentially overcome and mask the effects of mutational bias on CUB. This is the case, for example, of Drosophila spp. and of Strongylocentrotus purpuratus, which preferentially use G/C ending codons despite having A/T rich genomes (Vicario, Moriyama & Powell, 2007; Kober & Pogson, 2013). The selection for translational speed is evident in many unicellular and multicellular organisms, where the preferential use of optimal codons by highly expressed genes, with the aim to maximize the rate of elongation during protein synthesis (Marais & Duret, 2001), has been clearly demonstrated (Gouy & Gautier, 1982; Powell & Moriyama, 1997; Duret, 2000). Nevertheless, the correlation between gene expression and CUB is not universally applicable, as the selection for translational speed appears to be weak in some species (Hiraoka et al., 2009).

Natural selection also potentially shapes codon usage to improve translational accuracy, reducing the risk of missense errors during the translational process (Akashi, 1994; Eyre-Walker, 1996). Even though this model of selection predicts a higher CUB in genes encoding longer proteins, a positive correlation between CUB and protein length has only been observed in bacteria, while on the opposite an unexpected negative correlation has been described in a number of eukaryotes, including S. cerevisiae, D. melanogaster and C. elegans (Moriyama & Powell, 1998; Duret & Mouchiroud, 1999). Even though differences in selective constraints between prokaryotes and eukaryotes have been evoked to explain this contradiction, the relationship between CUB and protein length still remains obscure (Marais & Duret, 2001).

The investigation of these two selective forces in bivalves is complicated by the incomplete/fragmented nature of de novo assembled transcriptomes and by the limited gene expression resources available. We correlated CUB with gene expression (in three tissues: digestive gland, gills and hemocytes) and with ORF (protein) length in C. gigas, the most appropriate bivalve species for this purpose due to the completeness of its genome annotation and availability of gene expression data. We observed a significant positive correlation between CUB (negative ENC) and gene expression in the three tissues analyzed, as well as between CUB and ORF (protein) length (Table 2), which would suggest that both selection for translational speed and accuracy are actively shaping codon usage in oyster. However, we also observed a highly significant, negative correlation between GC₃ and gene expression and between GC₃ and protein length, which seem to contradict the mutational bias given by the A/T-rich nature of the oyster genome. This observation matches the results obtained in a previous work conducted with limited expression data based on Sanger EST sequencing, which suggested that translational selection acts as a contrasting force to mutational bias in oyster, effectively counteracting its action in highly expressed genes (Sauvage et al., 2007). Our data further indicate that, besides the selection for translational speed, also the selection for translational accuracy provides a contribution to the selection of G/C-ending codons in oyster.

The contrasting action of mutational bias and selection becomes particularly evident while taking into consideration the correlation between CAI, a directional measure of CUB which is based on a reference set of highly expressed genes (Sharp & Li, 1987), and ENC, which on the other hand is a non-directional measure which does not permit to appreciate the contribution of opposite forces (in this case mutational bias towards A/T-ending codons and selection towards G/C-ending codons). Indeed, in oyster and in all the other species clusterized in group I (see Fig. 2), the scatter in the correlation plot between CAI and ENC appears to be quite relevant (Fig. 6; Paerson correlation coefficients are −0.44 for C. gigas and −0.38 for P. magellanicus). This indicates that the mutational bias towards A/T is counterbalanced by natural selection in favor of G/C-ending codons in a relevant number of genes. However, the significant correlation between ENC and GC₃ at the whole protein-coding transcriptome level (Paerson correlation = 0.16, p-value = 5.42 × 10⁻¹⁴⁸) indicates that the weight of A/T mutational bias is still dominating over that of G/C selection in most oyster genes.

Overall, it is likely that in all species pertaining to clustering group I, whose coding GC₃ sensibly deviates from the frequency expected by genomic GC content, the mutational bias towards A/T-ending codons is countered by the opposite forces of selection for translational speed and accuracy, leading to moderate/low CUB.

On the other hand, the correlation between CAI and ENC in the bivalve species pertaining to clustering group II is highly significant (Fig. 6, Paerson correlation coefficients are −0.87 for M. chinensis and −0.71 for M. californianus). To better interpret this result, we extended the analysis performed in oyster to a representative species of this group, M. galloprovincialis, limiting our calculations to full-length protein-coding transcripts (Table 2). Overall, like in oyster, the significant positive correlation between gene expression and GC₃ over three different tissues indicates that selection for translational speed acts towards G/C-ending codons. However, unlike oyster, this appears to have no effect on CUB on the whole transcriptome scale, as the correlation between gene expression and ENC was not significant. In addition, we observed that, in contrast with oyster, the selection for translational accuracy is biased towards A/T-ending codons (i.e., GC₃ and protein length are negatively correlated). Overall, the M. galloprovincialis average GC₃ is 30.80%, a value which is very close to the GC genomic content, further indicating that mutational bias and selection for translational accuracy, which both act in favor of A/T-ending codons, are largely dominating over the opposite selection for translational speed in mussel.