Comparative analysis of codon usage patterns in chloroplast genomes of six Euphorbiaceae species

Genomes and coding sequences

The complete chloroplast genomes of Euphorbia esula, Hevea brasiliensis (Tangphatsornruang et al., 2011), Jatropha curcas (Asif et al., 2010), Manihot esculenta (Daniell et al., 2008), Ricinus communis (Rivarola et al., 2011), and Vernicia fordii (Li et al., 2017) with gene annotations were downloaded from the NCBI GenBank database. The number of raw CDSs of six Euphorbiaceae species was 85, 84, 84, 83, 86 and 85 respectively (Table 1; Table S1). In order to avoid sampling errors, each CDS in the chloroplast genomes of six Euphorbiaceae species should follow certain rules, i.e., the number of bases in each CDS should be the fold of three; the length of sequence encoding gene must be ≥ 300 bp; high-quality sequences with identified bases, i.e., containing only A, T, G and C bases; each CDS contains proper initiation codon (ATG) and termination codons (TAG,TGA and TAA); and sequences without an intermediate stop codon (He et al., 2016; Li et al., 2016; Zhang et al., 2007). We used Perl scripts written by our team to filter the CDSs according to the five rules mentioned above and simplify the names of CDSs replaced with numbers to avoid miscalculation. The GC content of the first, second and third codon positions (GC1, GC2, GC3) and the average GC content of three positions were calculated by Perl script.

Analysis of Relative synonymous codon usage (RSCU) and Relative synonymous codon usage frequency (RFSC)

RSCU value for a particular codon refers to the ratio of its actual usage frequency to expected frequency when it is used without bias. The RSCU was calculated as Eq. (1): (1)${\displaystyle RSCU=\frac{x_{ij}}{\sum _j^{n_i}{x}_{ij}}{n}_i}$ where x_ij represents the frequency of codon j encoding for the i th amino acid, and n_i represents the number of synonymous codons encoding the i th amino acid (Sharp & Li, 1986). If the RSCU value of one codon equals 1 that reflected no codon usage bias and is used equally with other synonymous codons. However, strong positive codon usage bias could be observed for RSCU value >1. In contrast, RSCU value <1 displayed negative codon usage bias that is used less frequently than other codons (Sharp & Li, 1987).

The RFSC value is equal to the ratio of the actual observed number of one codon to the number of all synonymous codons (Sharp & Li, 1986). The RFSC was calculated using Eq. (2): (2)${\displaystyle RFSC=\frac{x_{ij}}{\sum _j^{n_i}{x}_{ij}}}$ where x_ij represents the frequency of codon j encoding for the i th amino acid. The RFSC value of a codon exceeding 60% or 0.5 times higher than the average frequency of the synonymous codons indicates high-frequency codon (Zhou, Tong & Shi, 2007).

Comparative analysis of codon usage frequency

In order to deeply analyze the codon usage patterns of six Euphorbiaceae plant species, codon usage bias data of Arabidopsis thaliana (http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species=3702), Populus trichocarpa (http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species=3694), Escherichia coli (http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species=199310) and Saccharomyces cerevisiae (http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species=4932) downloaded from Codon Usage Database were compared with the codon usage frequencies of six Euphorbiaceae plants. Furthermore, we calculated the ratio of codon usage frequency for six Euphorbiaceae plant species to four model organisms. When the ratio is ≥ 2 or ≤ 0.5, it indicates the difference of codon usage bias between two organisms is greater (Pan et al., 2013).

Analysis of ENc-plot

ENc (effective number of codons) value is used to analyze the degree of deviation of codon usage from the random selection that depicted the degree of imbalanced use of synonymous codons in genes or genomes of the specific species. The range of ENc value is 20-61. The smaller the ENc value, the stronger the codon usage bias and vice versa (Wu et al., 2018). When ENc value is ≤ 35, the codon usage of genes or genomes has very significant bias (Mensah et al., 2019). GC3s value refers to the ratio of G and C content at the third position of one codon to the total number of gene bases excluding Met and Trp. Using GC3s value as abscissa and ENc value as ordinate to make the ENc-plot, the results revealed the influencing factors of codon usage patterns of genes or genomes, and the relationship between gene base composition and codon usage bias (Wright, 1990). The expected values of ENc were calculated according to the Eq. (3): (3)${\displaystyle ENc=2+S+\frac{29}{S^2+{\left(1-S\right)}^2}}$ where S denotes GC3s (Wright, 1990; Zhang et al., 2007). When mutation pressure plays an important role in the formation of codon usage patterns, ENc value lies on or around the expected curve. However, when codon usage is affected by natural selection and other factors, ENc value is far lower than the expected curve (Wright, 1990).

PR2-plot analysis

PR2-plot is used to analyze the composition of four bases at the third position of codon encoding amino acids. It is a graphical analysis based on A3/(A3 + T3) as ordinate and G3/(G3 + C3) as abscissa (Sueoka, 1999). The distribution of points around the center point (A = T, C = G) shows the degree and direction of the base deviation. It was generally believed that the proportion of A/T and C/G is balanced in degenerate codons of genes or genomes upon single mutation pressure (Xiang et al., 2015).

Analysis of Neutrality plot

Neutrality plot (GC12 vs. GC3) was performed to investigate the extent of influence between mutation pressure and natural selection on the patterns of codon usage (Sueoka, 1988). GC12 represents the average value of GC contents at the first and second positions of codon while GC3 is the GC content at the third position. GC3 was calculated excluding the three termination codons (TAA, TAG and TGA) and the three codons for Ile (ATT, ATC and ATA). Meanwhile, two single codons for Met (ATG) and Trp (TGG) were also excluded in all three patterns (Sueoka, 1988). GC12 and GC3 of chloroplast genomes in six Euphorbiaceae species were calculated by Perl script. The slope of the plot regression was zero indicates no effects of directional mutation pressure (complete selective constraints). The slope 1 depicted that the codon usage bias is completely affected by directional mutation pressure representing complete neutrality (Sueoka, 1988; Wen et al., 2017).

Correlation analysis (COA)

The codon usage variations in chloroplast genomes of six Euphorbiaceae plant species were investigated with correspondence analysis based on RSCU using CodonW (Version 1.4.2; Mensah et al., 2019). Correspondence analysis was performed to compare the usage patterns of 59 codons (excluding codons encoding Met, Trp and three termination codons), and the results produce a series of orthogonal axes that can be used to present the codon usage variation in chloroplast genomes of six Euphorbiaceae plant species. The distribution of genes can be drawn according to the synonymous codon usage of the genes in a multidimensional space of 59 axes, followed by the maximum fraction of gene variations, thus the main sources of codon usage variation were analyzed (Xiang et al., 2015). Based on the results of codon usage variation, correlation analysis between axis 1 and codon usage indices including codon adaptation index (CAI), the GC content at the third codon position of synonymous codons (GC3s; Zhang et al., 2007) and the total number of amino acids (L_aa; Wright, 1990) were carried out by SPSS (Version 23). The value is negative means a negative correlation. CAI value is widely used to evaluate the gene expression level and ranges from 0 to 1. The larger the CAI value, the stronger the codon usage bias, otherwise, the weaker the codon usage bias (Sharp & Li, 1986).

Characteristics of codon usage bias

Sources analysis of variation in codon usage

ENc-plot

The distributions of ENc and GC3s of chloroplast genomes of six Euphorbiaceae plant species were similar (Fig. 1). Only few points lie in close proximity to the curve, however, a majority of genes with lower ENc values than expected values lay below the curve (Fig. 1). Analysis about points indicated that codon usage bias of chloroplast genomes was affected slightly by the mutation pressure, but natural selection and other factors play the major role (Wright, 1990). Previous researches suggested the codon usage bias of chloroplast genomes of Populus alba (Zhou, Long & Li, 2008), Poaceae (Zhang et al., 2012), Asteraceae (Nie et al., 2013) were influenced by combined the effects of mutation pressure, natural selection and other factors.

PR2-plot

It was an efficient way to reflect the mutation pressure by analyzing the points representing values for G3/(G3 + C3) and A3/(A3 + T3) distributed around the central spot (A = T, C = G). It was revealed that the AT-bias is 0.488, 0.485, 0.482, 0.485, 0.487 and 0.484 for Euphorbia esula (Fig. 2A), Hevea brasiliensis (Fig. 2B), Jatropha curcas (Fig.  2C), Manihot esculenta (Fig. 2D), Ricinus communis (Fig. 2E) and Vernicia fordii (Fig. 2F), while the GC-bias is 0.499, 0.507, 0.509, 0.499, 0.504 and 0.501, respectively. Thus, T/C bias at the third position of codons of chloroplast genes was observed in Euphorbia esula and Manihot esculenta, however, T/G-bias was observed in other four Euphorbiaceae plant species. As a whole, the usage frequency of A/T and G/C in six chloroplast genomes was unbalanced that was not only affected by the mutation pressure, but also with natural selection and other factors. Similar studies have also been reported for the codon usage of chloroplast genomes of Asteraceae (Nie et al., 2013) which illustrated that purines were used more frequently than pyrimidine in the chloroplast of Asteracceae. The analysis of PR2-plot only reflected the factors that influenced codon usage pattern, hence, further analyses are needed to explore the extent of the influencing factors between mutation pressure and natural selection.

Neutrality plot

The neutrality plot reflected the narrow range of GC12 (0.30–0.58) and GC3 (0.16–0.42) value distributions (Fig. 3). The correlation between GC1 and GC2 was very strong (r ₁ = 0.530, r ₂ = 0.493, r ₃ = 0.542, r ₄ = 0.511, r ₅ = 0.559, r ₆ = 0.538, p <0.01). However, no significant correlation was found for GC1 with GC3 (r ₇ = 0.143, r ₈ = 0.092, r ₉ = 0.138, r ₁₀ = 0.106, r ₁₁ = 0.030, r ₁₂ = 0.070) or GC2 with GC3 (r ₁₃ = 0.123, r ₁₄ = 0.194, r ₁₅ = 0.257, r ₁₆ = 0.199, r ₁₇ = 0.129, r ₁₈ = 0.184), which indicated mutation pressure had a minor effect on the codon usage bias. Moreover, the slope of neutrality plot revealed that mutation pressure only accounted for 12.90%–25.58% on the codon usage patterns in six chloroplast genomes while natural selection accounted for 74.42%–87.10%. These results demonstrated that natural selection played a significant role in the codon usage patterns.

Correspondence analysis (COA)

Correspondence analysis is a multivariate statistical method to explore the relationship between the variables in samples (Shields & Sharp, 1987). In the current study, the correspondence analysis based on RSCU was used to reveal the main factors affecting the formation of codon usage patterns in the chloroplast genomes of six Euphorbiaceae plant species. The position of the origin represented the average RSCU value for all genes, with respect to axis 1 and axis 2. The first four axes accounted for 36.36%, 35.97%, 33.64%, 35.38%, 36.03% and 34.17% of the overall variation. The first axis accounted for 11.09%, 11.55%, 10.17%, 11.74%, 11.68% and 9.70% of the total variation in six plants respectively. Therefore, axis 1 was the major source of variation, responsible for ∼10% of total variation. This indicated that the codon usage might be not affected by the single factor. To investigate the effects of GC content on CUB, each gene of chloroplast of six Euphorbiaceae plant species was distributed on the plane with axis 1 as the abscissa and axis 2 as the ordinate axes with different colors (Fig. 4). There was only a gene with GC content within 45%–60% plotted as bottle green in six plants, while all the other genes with GC content were lower than 45%.

To identify the factors resulting in the dispersion of chloroplast genes along axis 1 and axis 2, the correlation coefficients were calculated on axis 1 with CAI, GC3s and L_aa (Table 2). Based on the correlation analysis of axis 1 and codon usage indices (CAI, GC3s, L_aa), it was found that axis 1 for Euphorbia esula, Jatropha curcas and Manihot esculenta had a significant correlation with GC3s (p ≤ 0.01), while Hevea brasiliensis, Ricinus communis and Vernicia fordii had correlation with GC3s (p ≤ 0.05), which indicated GC3s is significant for patterns of codon usage (Table 2). Zhou and colleagues (2008) reported that the axis 1 for the codon usage bias of chloroplast genome of Populus alba significantly correlate with GC3s and gene length, which is in line with the findings of Xu et al. (2011), who reported that axis 1 was significantly correlated with GC3s, gene length and hydrophilicity in the chloroplast genome of the Oncidium Gower Ramsey, suggesting the effect of mutation, gene length and expression level.

Table 2:

Correlation analysis of axis 1 and codon usage index of chloroplast genomes of six Euphorbiaceae plant species.

CAI means codon adaptation index; GC3s indicates the GC content at the third codon position of synonymous codons; L_aa is defined as total number of amino acids.

	Euphorbia esula	Hevea brasiliensis	Jatropha curcas	Manihot esculenta	Ricinus communis	Vernicia fordii
CAI	−0.023	−0.114	−0.035	−0.185	0.152	−0.307*
GC3s	0.421**	−0.320*	−0.581**	−0.352**	0.284*	0.324*
L_aa	0.165	−0.059	−0.020	−0.059	0.068	0.216

DOI: 10.7717/peerj.8251/table-2

Notes:

*P < 0.05.

**P < 0.01.