Non-uniqueness of factors constraint on the codon usage in Bombyx mori

Abstract

Background

The analysis of codon usage is a good way to understand the genetic and evolutionary characteristics of an organism. However, there are only a few reports related with the codon usage of the domesticated silkworm, Bombyx mori (B. mori). Hence, the codon usage of B. mori was analyzed here to reveal the constraint factors and it could be helpful to improve the bioreactor based on B. mori.

Results

A total of 1,097 annotated mRNA sequences from B. mori were analyzed, revealing there is only a weak codon bias. It also shows that the gene expression level is related to the GC content, and the amino acids with higher general average hydropathicity (GRAVY) and aromaticity (Aromo). And the genes on the primary axis are strongly positively correlated with the GC content, and GC3s. Meanwhile, the effective number of codons (ENc) is strongly correlated with codon adaptation index (CAI), gene length, and Aromo values. However, the ENc values are correlated with the second axis, which indicates that the codon usage in B. mori is affected by not only mutation pressure and natural selection, but also nucleotide composition and the gene expression level. It is also associated with Aromo values, and gene length. Additionally, B. mori has a greater relative discrepancy in codon preferences with Drosophila melanogaster (D. melanogaster) or Saccharomyces cerevisiae (S. cerevisiae) than with Arabidopsis thaliana (A. thaliana), Escherichia coli (E. coli), or Caenorhabditis elegans (C. elegans).

Conclusions

The codon usage bias in B. mori is relatively weak, and many influence factors are found here, such as nucleotide composition, mutation pressure, natural selection, and expression level. Additionally, it is also associated with Aromo values, and gene length. Among them, natural selection might play a major role. Moreover, the “optimal codons” of B. mori are all encoded by G and C, which provides useful information for enhancing the gene expression in B. mori through codon optimization.

Background

Codon usage bias refers to differences of the occurrence frequency of synonymous codons in coding DNA. It is considered to be a product of mutation pressure and/or natural selection [1-4], and accounts for accurate and efficient translation, as well as mutation–selection–drift [5]. Codon bias analysis has been introduced into both prokaryotes and eukaryotes, such as Escherichia coli (E. coli), Arabidopsis thaliana (A. thaliana), and human beings [6-9], showing that codon bias has a high correlation to gene length, gene function, hydrophobicity of proteins, and the content of iso-acceptor tRNAs in genomes [9-12]. Hence, the analysis of codon usage can be used to study organism evolution and improve protein expression level [13-15].

The domesticated silkworm, Bombyx mori (B. mori), is a well-studied lepidopteran model system with rich genetic and molecular information of morphology, development, and behavior [13]. So far, the draft sequence for the genome of B. mori has been determined [16], and most studies of B. mori focus on the cloning, expression, and characterization of some genes or application as the bioreactor [17-19]. As we know, the analysis of codon usage is a good way to understand the genetic and evolutionary characteristics of B. mori. It can also help us to study the relationship between expression levels and codon usage bias since highly-expressed genes need abundant ribosomes and matching tRNAs for efficient translation. We have reported the codon usage bias of the mitochondrial genome in B. mori recently [20], however, the codon usage bias in the whole nuclear genome of B. mori is not well investigated in detail. Considering its great potential for expressing foreign proteins as a bioreactor, the codon usage bias of B. mori was examined here for codon optimization of genes.

Results and discussion

B. mori reveals a weak codon bias

As shown in Additional file 1 and Table 1, the GC content for the total 1, 097 genes varies from 29.5% to 69.5%, with a mean value of 46.43%. The GC content of the total genes is distributed mainly between 40% and 50% (Figure 1). The greatest differences of GC content are found in the first and the third codon positions (51.92% and 48.40%, respectively), where most neutral mutations occur [21].

Table 1 Means and standard deviations of GC, GC1, GC2, GC3, GC3s, ENc, CAI, A3s, T3s, C3s, G3s, Gravy, and Aromo of codons from Bombyx mori

Figure 1

The distribution of GC contents in the CDS of Bombyx mori.

The effective number of codons (ENc) in B. mori ranges from 30.06 to 61.00, with an average of 53.12. As shown in Additional file 1, among the 1, 097 genes, only 5 genes reveal a high codon bias (ENc < 35). It indicates that B. mori exhibits a general random codon usage, without strong codon bias. Similarly, the relative synonymous codon usage (RSCU) values of 59 sense codons also support the conclusion that B. mori has a weak codon bias. As shown in Table 2, approximately half of the codons (28/59), denoted in bold lettering, are frequently used, such as GCU and AGA which encode Ala and Arg, respectively.

Table 2 Codon usage of Bombyx mori genes (363,313 codons)

In addition, most of preferentially used codons end with A/U (A/U-ended: G/C-ended=18:10). This phenomenon was also found in many other AT-rich species, such as Pichia pastoris (P. pastoris), Saccharomyces cerevisiae (S. cerevisiae), Kluyveromyces lactis (K. lactis), and Plasmodium falciparum (P. falciparum) [22,23].

Effects of nucleotide composition in shaping codon bias

Correspondence analysis of the RSCU values was used here, which removes the variation caused by the unequal usage of amino acids (although the degrees of freedom are reduced to 40 [24]), generating a first axis that explains 24.51% of the data inertia. The second axis explains 7.46%, while the next two axes respectively account for 4.02% and 3.39% of the data (Figure 2). Moreover, multivariable correlation analysis was introduced here to study the relationship between relative codon bias and nucleotide composition (Table 3).

Figure 2

The relative and cumulative inertia of the first 20 factors from a correspondence analysis (COA) of the RSCU values.

Table 3 Correlation coefficients between the positions of genes along the first two major axes with index of total genes’ codon usage and synonymous codon usage bias

Although the first axis can’t explain the whole variation, there is an obvious positive correlation between the first axis and G3s, C3s, and GC3s (r=0.343, 0.439, and 0.446, respectively, p < 0.01). However, the correlations between the first axis and A3s or T3s are negative (r=−0.444 and r=−0.367, respectively, p < 0.01). Then all the genes were classified into three categories by their GC content (GC < 45%, 45% ≤ GC < 60%, and GC ≥ 60%). As shown in Figure 3A, the position of each gene was marked along the first two major axes. Interestingly, the genes of GC < 45% are scattered at the left side of the first axis, while most of the genes with GC ≥ 60% are located at the right side of the first axis. The genes whose GC contents range from 45% to 60% are found in the middle of the plot. Additionally, almost all the ribosome genes are located in the range of GC ≥ 60%, implying that the expression level might be related with the GC content in B. mori.

Figure 3

Correspondence analysis of RSCU for the total genes in Bombyx mori. A) Distribution of the total genes in Bombyx mori on the plane corresponding to the coordinates on the first and second principal axes was shown in Panel A. Orange triangles, purple triangles and gray triangles, indicate genes with a GC content higher than or equal to 60%, more than or equal to 45%, but less than 60% and less than 45%, respectively. Additionally, blue squares indicate the coordinates of ribosome on the first and second principal axes. B) Distribution of codons on the same two axes was shown in Panel B. Codons ending with A, U, C and G are shown in red, green, yellow, and blue, respectively.

On the overall consideration of Tables 1 and 3, it seems that the genes containing lower GC3s and GC content values tend to distribute at the left side of the first axis. Thus, we speculated that G/C-ending codons could be clustered at the positive side whereas A/U-ending codons gather at the negative side of first major axis. The corresponding distribution plot of synonymous codons ending with different bases along the two axes was implemented under the above mentioned assumption. The result indicates that the separation of codons on the first axis reflects the difference between the frequencies of A/U and C/G ending codons, while that on the second axis represents the frequency differences between A/G and U/C ending codons (Figure 3B), which is consistent with the above-mentioned hypothesis.

On the other hand, the ENc values show no significant correlation with the first axis (r=0.055, p > 0.05) or GC3s (r=0.037, p > 0.05) values, but a significant positive correlation with the second axis (r=0.184, p < 0.01) (Table 3).

The results above suggest that nucleotide composition has an effect on separating the genes along the first major axis, however, it might be not the main factor in shaping the codon bias.

GC3s plays a minor role in shaping the codon bias of B. mori

ENc-plot is an effective tool to study the codon usage patterns, and it was used here to explore the influence of GC3s on the codon bias of B. mori. As shown in Figure 4, most genes are located below the expected ENc-plot curve while only a small number of genes lay on or above the curve. It indicates that the conditional mutation might be a factor in shaping the codon bias but not the unique one.

Figure 4

The ENc plotted against GC3s. ENc denotes the effective number of codons, and GC3s denotes the GC content on the third synonymous codon position. Black boxes, and blue triangles indicate ribosome genes and total genes, respectively. The red solid line represents the expected curve of positions of genes when the codon usage was only determined by the GC3s composition.

We also estimated the difference between the observed and the expected ENc values using the plot of the frequency distribution of (ENCexp-ENCobs)/ENCexp in total genes (Figure 5). There was a similar single peak for each kind of genes. Peaks located within the 0 ~ 0.1 range of (ENCexp-ENCobs)/ENCexp values suggest that most actual ENc values are smaller than the ENc values from their GC3s. It is consistent with the results depicted in Figure 4, which shows that the difference in codon bias is dependent upon the differences in GC3s, thereby providing further evidence that GC3s works as a conditional mutational bias.

Figure 5

Frequency distribution of (ENCexp-ENCobs)/ENCexp.

Natural selection influences the codon bias as a major role

Although ENc plot can quantify the codon usage bias of synonymous codons, it is not sufficient to easily distinguish the main determinant factor between natural selection and mutational pressure within a species [25]. Therefore, a neutrality plot was implemented here.

The neutrality plot shows that the genes have a wide range of GC3 value distributions, ranging from 19.7% to 93.8% (Figure 6). There is a significant positive correlation between GC12 and GC3 (r=0.394, p < 0.01), suggesting that the effect of directional mutation pressure is present at all codon positions. Moreover, the slope of the regression line of the entire coding sequence is 0.1452. The results reveal that the effect of directional mutation pressure is only 14.52%, while the influence of other factors, for example natural selection, is 85.48% [26]. Accordingly, mutation bias only plays a minor role in shaping the codon bias, whereas natural selection probably dominates the codon bias.

Figure 6

Neutrality plot analysis of the GC12 and that of the third codon position (GC3) for the entire coding DNA sequence of Bombyx mori. GC12 stands for the average value of GC content in the first and second position of the codons (GC1 and GC2). While GC3 refers to the GC content in the third position. The solid line is the linear regression of GC12 against GC3, R²=0.1534, P < 0.001.

Codon usage bias in B. mori has a high correlation to aromaticity and gene length

In order to assess the relationship between the codon usage bias and hydrophobicity or aromaticity or gene length in B. mori, correlation analysis was performed. It could be observed from Table 3 that neither the Gravy values nor the Aromo values have significant correlation with GC3s. However, the Aromo values exhibit strongly positive correlation with the ENc values (r=0.100, p < 0.01), while the GRAVY values do not. The results indicate that the Aromo values are associated with the codon usage bias of B. mori.

The data in Table 3 also reveal that the gene length is positively correlated with the ENc values (r=0.079, p < 0.01), suggesting that gene length has a high correlation to the codon usage bias and might be also one of the factors contributing to the codon usage bias in genes.

Effects of gene expression level

To explore the relationship between codon bias and gene expression level, correlation coefficients were calculated between the codon adaptation index (CAI) values and several other characteristics of the genes, including their position along the first major axis, the nucleotide composition, and the ENc values. Ribosome genes sequences were selected as the reference of highly expressed genes [15].

The results indicate that CAI, which represents gene expression level, shows significant negative correlation with the gene length (r=−0.148, p < 0.01), GC2 (r=−0.081, p < 0.01), A3s (r=−0.444, p < 0.01), T3s (r=−0.061, p < 0.05), Gravy (r=−0.170, p < 0.01), Aromo (r=−0.140, p < 0.01), and ENc (r=−0.210, p < 0.01). However, CAI shows obvious positive correlation with the first axis and the other nucleotide composition indices (i.e. GC, GC1, GC3, GC3s, C3s, and G3s, as shown in Table 3). The results above indicate that both nucleotide composition and gene expression levels are the major factors in shaping the codon usage bias of B. mori.

To statistically measure the relationship between the index of amino acid composition in B. mori and their codon bias, the correlation coefficients between the positions of the genes along the first four major axes with their indices of amino acid usage were analyzed using Spearman’s rank correlation analysis method and shown in Table 4.

Table 4 Correlation coefficients between the positions of genes along the first four major axes with index of total genes’ amino acid usage

The first four axes generated by the correspondence analysis explain 40.31% of the amino-acid variation. And the first axis accounts for 13.90% of the variation in amino-acid usage (Figure  7 ). The genes on these axes are all highly correlated with CAI, GRAVY score and Aromo value. The principle factor is negatively correlated with CAI (r=−0.216, p < 0.01), and is positively correlated with the GRAVY score and Aromo value (r=0.327, p < 0.01; r=0.208, p < 0.01, respectively). The second axis accounts for 10.95%, and is also correlated with the three indexes (r=−0.208, p < 0.01; r=0.545, p < 0.01; r=0.594, p < 0.01, respectively).

Figure 7

The relative and cumulative inertia of the first 20 factors from a correspondence analysis (COA) of the amino acid usage frequencies.

As in E.coli [27] where the most important trend in the amino-acid usage of B. mori is the usage of hydrophobicity, and the second important trend is the usage of CAI followed by the aromatic amino-acid. Taken all these together, it provides strong evidence for the inference that the effective selection of amino-acid for translational efficiency exists in B. mori.

In summary, the codon usage bias in B. mori is in some way or other, affected by nucleotide composition, mutation pressure, natural selection, and gene expression level. Additionally, it is also associated with Aromo values, and gene length. However, natural selection might play a major role in shaping codon usage variation, manifesting itself though weaker codon usage bias. The selection of amino-acid could also affect the translational efficiency in B. mori.

Translational optimal codons of B. mori

In order to give a reference to enhance the expression level of important proteins with codon optimization, a two-way Chi-squared contingency test was used to compare the codon usage of different genes. Finally, the total putative optimal codons of B. mori are listed in Table 5. For the total genes group, the optimal codons all ended by G or C, and all amino acids—excluding Met and Trp—were identified by different numbers of codons. For example, three codons were identified for Ser, and two codons were identified for Ala, Gly, Leu, Pro, Val, Thr, and Arg. The remaining amino acids were identified by one codon.

Table 5 Optimal codons of genes in Bombyx mori

The optimization of codon usage allows improving the translational efficiency of foreign proteins by replacing the codons which are rarely found in the host organism [28], and it has been introduced into many heterologous systems [29-31]. As we found in this study, the optimal codons of B. mori are all ended by either G or C. This phenomenon is interesting and important to enhance the expression level of foreign proteins in B. mori.

Comparison of codon preferences between B. mori and other model organisms

The ratio of codon frequency in B. mori was compared with five model organisms, including A. thaliana, C. elegans, Drosophila melanogaster (D. melanogaster), S. cerevisiae, and E. coli. The codon with a ratio of greater than 2, or less than 0.5, is defined as the indicative codon, of which usage frequency is markedly distinct from that of B. mori. As shown in Additional file 2, there are six and seven codons revealing distinct usage differences between B. mori and D. melanogaster, S. cerevisiae, respectively. However, there are only one, two or three codons with distinct usage between B. mori and A. thaliana, or E. coli, or C. elegans, respectively. It suggests that the discrepancy in codon preferences between B. mori and D. melanogaster or S. cerevisiae is relatively greater than that comparing with A. thaliana, or E. coli, or C. elegans. This finding implies that B. mori might have some advantages in expressing foreign proteins from certain organisms with fewer preferences in codon usage.

Conclusions

After a series of analyses, the codon usage bias in B. mori is found to be weaker. And it is affected by nucleotide composition, mutation pressure, natural selection, and gene expression level. Additionally, it is also associated with Aromo values, and gene length. However, natural selection might play a major role in shaping the codon usage variation. In addition, it is also found that B. mori has a greater relative discrepancy in codon preferences in comparison with D. melanogaster or S. cerevisiae than with A. thaliana, E. coli, or C. elegans.

In summary, our analysis provides insights into the codon usage pattern in B. mori and is of the benefit to express foreign proteins in B. mori as a bioreactor.

Methods

Sequence collection

Accession numbers for a total of 1,213 reference sequences (RefSeq) of B. mori were obtained from Silkworm Genome Database (ftp://silkdb.org/pub/current/otherdata/Refseq/silkref.seq) (Downloaded on 1-Sep-2014). Coding DNA sequences (CDS) were downloaded from GenBank (http://www.ncbi.nlm.nih.gov). In these sequences, we only chose the CDSs without unidentified bases. To improve the quality of sequences and minimize sampling errors, genes without correct initiation and termination codons or with internal termination codons were ruled out. Additionally, only genes greater than 300 nucleotides in length were used for further analysis. As we only collected the CDSs from nuclear genome, 13 mitochondrial genes were excluded from the analysis. CDSs with gaps were also excluded. Finally, 1,097 CDSs were left for analysis, and each corresponds to a unique gene in B. mori.

Indices of codon usage and synonymous codon usage bias

GC3s is a useful parameter for evaluating the degree of base composition bias, and represents the frequency of either a guanine or cytosine at the third codon position of synonymous codons, excluding Met, Trp, and stop codons.

GRAVY (General Average Hydropathicity) values are calculated as a sum of the hydropathy values of all the amino acids in the gene product divided by the number of residues in the sequence [32]. The more negative the GRAVY value, the more hydrophilic the protein, while the more positive the GRAVY value, the more hydrophobic the protein.

Aromo values denote the frequency of aromatic amino acids (Phe, Tyr, Trp) in the hypothetical translated gene product. The index and GRAVY value have been used to quantify the major COA trends in the amino acid composition of E. coli genes [27].

RSCU (relative synonymous codon usage) is the ratio of the observed frequency of codons relative to the expected frequency of the codon under a uniform synonymous codon usage. The RSCU value would be greater than 1.0 when the observed frequency is larger than the expected frequency [33].

ENc (Effective Number of Codons) values, varying from 20 to 61, are used to measure the magnitude of codon bias for individual genes, though it is worth noting that ENc values are affected by base composition [34]. A value of 20 indicates a gene with extreme bias using only one codon per amino acid, while a value of 61 indicates the absence of bias. In general, a gene is thought to possess strong codon bias if its ENc value is lower than 36 [35].

CAI (Codon Adaptation Index) values are often used to measure the extent of bias toward codons which are known to be preferred in highly expressed genes. With values ranging from 0 to 1.0, the higher the value, the stronger the codon usage bias and the higher the expression level. The set of sequences used to calculate CAI values in this study were the genes coding for ribosomal proteins in B. mori [35], so that it can provide an indication of gene expression level under the assumption that translational selection can optimize gene sequences according to their expression levels. These noted values and parameters were all utilized in this study.

All the indices of total genes and ribosomal protein genes are shown in Additional file 1, respectively.

ENc-plot

The ENc-plot is a general strategy to investigate patterns of synonymous codon usage, where the expected ENc values are plotted against GC3s values. Expected ENc values were calculated according to Equation 1. In genes where codon choice is constrained only by a G + C mutation bias, predicted ENc values will lie on or around the GC3s curve, whereas if other factors such as selection effects are present, the values will deviate considerably below the expected GC3s curve [35].

$$ \mathrm{E}\mathrm{N}\mathrm{c}=2+\mathrm{S}+\left(29/\left({\mathrm{S}}^2+{\left(1\hbox{-} \mathrm{S}\right)}^2\right)\right) $$

(1)

S is the frequency of G + C (i.e. GC3s)

Neutrality plot

A neutrality plot (GC12-GC3) [26] was used to estimate and characterize the relationships amongst the three positions in B. mori codons. A plot regression with a slope of 0 indicates no effects of directional mutation pressure (complete selective constraints), while a slope of 1 is indicative of complete neutrality.

Determination of optimal codons

Based on axis 1 ordination, the top and bottom 5% of genes were regarded as the high and low datasets, respectively. Codon usage was compared using a two-way Chi-squared contingency test to identify optimal codons. The test dataset with the lower ENc values were putatively assigned as highly expressed, and those codons which occur significantly more often (p < 0.01) were defined as optimal codons [24].

Correspondence analysis of RSCU

Correspondence analysis (COA) is a widely used method in the multivariate statistical analysis of codon usage patterns. Since there are a total of 59 synonymous codons (including 61 sense codons, minus the unique Met and Trp codons), the degrees of freedom was reduced to 40 in removing variations caused by the unequal usage of amino-acids while generating a correspondence analysis of RSCU [24].

Software used

Mobyle server (http://mobyle.pasteur.fr), including CodonW (Ver.1.4.4) (http://mobyle.pasteur.fr/cgi-bin/portal.py#forms::CodonW), CHIPS (http://mobyle.pasteur.fr/cgi-bin/portal.py#forms::chips), and CUSP (http://mobyle.pasteur.fr/cgi-bin/portal.py#forms::cusp), were used to calculate useful indices of codon usage bias, such as GC, GC3s (G + C content at the third position of codons), silent base compositions (i.e. A3s, T3s, C3s, and G3s, which indicate the frequency that codons have an A, U, C, or G, respectively, at their synonymous third position), GRAVY values (general average hydropathicity), Armoro values (aromaticity), RSCU (relative synonymous codon usage), and ENc (effective number of codons). Similarly, a COA (correspondence analysis) was also performed.

CAI (codon adaptation index) and gene length were calculated using the CAI calculate server (http://genomes.urv.es/CAIcal). GC1, GC2, and GC3 values were also calculated to determine the GC content at the first, second, and third codon positions, respectively.

Together these indices allow for an assessment of the level to which selection has been effective in shaping codon usage [33]. Codon preferences of other organisms were downloaded from the Codon Usage Database (http://www.kazusa.or.jp/codon) for comparison.

Statistical analysis

Correlations between codon usage variations amongst indices of codon usage were carried out using the multi-analysis software SPSS Version 22.0 (SPSS Inc. software, Chicago, Illinois, USA) and GraphPad Prism 5 (GraphPad Software, San Diego, California, USA).