Comparative analysis of codon usage patterns in chloroplast genomes of five Miscanthus species and related species

Genomes and sequences selection

The complete cpDNAs of M. floridulus (NC_035750.1), M. sacchariflorus (NC_028720.1), M. sinensis (NC_028721.1), M. x giganteus (NC_035753.1), M. transmorrisonensis (NC_035752.1), Sorghum bicolor (NC_008602), Saccharum spontaneum (NC_034802.1) with gene annotation were downloaded from the NCBI GeneBank database. The number of raw sequence coding for amino acids in protein (CDS) of above seven species was 106, 122, 122, 106, 106, 84 and 76 respectively (Table 1). To avoid sampling bias, the CDS sequences were screened from genome-wide data by python scripts (https://github.com/shexuan/codon_analysis) according to the following principles: (1) CDS contains initiation codon (ATG), termination codons (TAA, TAG or TGA) and without intermediate stop codons in the sequences; (2) the number of bases in each CDS must be the fold of three (3) the length of the sequence of CDS should be ≥300 bp (Wright, 1990; Zhang et al. 2007). After filtration, the CDS number, the base composition at the first/second/third site of codons (GC1/GC2/GC3) and average GC, as well as the total amino acids encoded by each CDS were calculated.

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

Relative synonymous codon usage value of a codon (the number of codon occurrences in a gene divided by the number of codon appearances expected under the same codon usage) is the ratio of its actual frequency of utilization to the expected usage frequency without bias. The RSCU was calculated as Eq. (1): (1)${\displaystyle \mathrm{RSCU}=\frac{x_{ij}}{\sum _j^{N_i}{x}_{ij}}{n}_i}$ where x_ij represents the frequency of codon j encoding the i th amino acid, and n_i represents the number of synonymous codon encoding the i th amino acid (Sharp & Li, 1986). If the RSCU value of a codon is equal to 1, the codon is used without bias, whereas a RSCU value greater than 1 reflects a significant codon usage bias (Sharp et al., 1993).

The RFSC value refers to the proportion of the actually observed number of a codon in the number of all synonymous codons. The RFSC were calculated using Eq. (2): (2)${\displaystyle \mathrm{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 high-frequency codon was screened based on the results of RFSC in all codon. The screening principles were as follows: the RFSC > 60% of one codon; or the RFSC of a codon exceeds the average frequency of synonymous codon by 0.5 times (Zhou et al., 2007).

Determination of optimal codons

The effective number of codons (ENC) can be applied to describe the extent of deviation of codon usage from the random selection, which reflects the degree of unbalanced use of synonymous codon in genes. The ENC value range from 20 (each amino acid uses only one synonymous codon) to 61 (Each synonymous codon is equally used), which is inversely proportional to the codon bias (Wright, 1990). The ENC value in each species was calculated by CodonW software and then 10% of the CDS with remarkably high and low expression levels were filtered out according to the ENC value. The RSCU of each codon was obtained from the sequence files of the high and low groups according to the cusp function of emboss (https://www.bioinformatics.nl/emboss-explorer/). Optimal codons were determined by ΔRSCU method. Specifically, the average RSCU values of the two amino acid groups were computed and subtract subsequently (ΔRSCU). The codon will be identified as the optimal codon through comparing the high and low group of the same codon ΔRSCU (>0.08) and RSCU value (high group > 1, low group < 1) (Romero, Zavala & Musto, 2000).

Comparative analysis of codon usage frequency

The ratio of codon usage frequency is one indicator of codon usage bias among species. To further explore the codon usage patterns in the seven species of Miscanthus and their relatives, codon usage bias data of four model species including Escherichia coli (http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species=199310); Saccharomyces cerevisiae (http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species=4932); Populus trichocarpa (http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species=3694) and Arabidopsis thaliana (http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species=3702), which have been used as the most common heterologous expression hosts were downloaded from the Codon Usage Database. Subsequently, the codon usage frequencies of the seven species in this study were compared with the above four model organisms. When the ratio is ≥2 or ≤0.5, it suggests that the codon bias difference between the two organisms is significant, whereas other values outside of this range represent a lack of significance (Pan et al., 2013).

Analysis of ENC-plot

The proportion of G and C content at the third position of a codon to the total number of gene bases are defined as GC3s. ENC-plot is plotted with ENC values as ordinate and GC3 value as abscissa, which can be used to analyze the codon usage characteristics of each gene and to explore the relevance between gene base component and codon preference (Wright, 1990). ENC values are located on or near the expected curve, when mutation pressure plays a key role in the formation of codon usage patterns. Conversely, when the use of a codon is constrained by natural selection, the ENC value will be well below the prospective curve (Wright, 1990).

PR2-plot analysis

In G3/(G3+C3) as the abscisic and A3/(A3+T3) as the ordinate graphic mapping (PR2-plot), which is performed to explore the composition of the four bases at the third nucleotide position of each codon (Sueoka, 1999; Sueoka, 1995). The pattern of splashes around the central spot (A = T, C = G) indicate the extent and orientation of the base offset.

Analysis of Neutrality plot

Neutrality analysis is used to exploring the degree of impact between natural selection and mutation pressure on the mode of codon usage (Sueoka, 1988). GC12 indicates the mean GC content at the first and second nucleotide positions of the codon, while GC3 represents the GC content of the third site. GC content at the third nucleotide position of a codon was counted post eliminating the Codon Met (ATG) and Trp (TGG). GC3 was counted post the eliminating of the three stop codons (TAA, TAG and TGA) and three codons (ATT, ATC and ATA) of Ile (Sueoka, 1988). Both the GC12 and GC3 values of the seven cpDNAs were counted by Python scripts (https://github.com/shexuan/codon_analysis). If the gradient of the curve regression is 0, indicating that there is no impact of mutation pressure. Gradient 1 represents complete neutrality, which describes that codon usage preference is completely influenced by mutation pressure (Sueoka, 1988).

Correspondence analysis of codon usage

The variations of codon usage in the seven analyzed cpDNAs were investigated based on the correspondence analyses (COA) using CodonW (Anue et al., 2019). The usage patterns of 59 codons (exincluding Met, Trp and three termination codons) were compared and all genes can be embeded into a 59-dimension hyperspace, in which each dimension corresponds to the synonymous codon usage of the gene (Xiang et al., 2015). Therefore, the major trends (Axis 1) of these axes in the 59-dimensional hyperspace can be used to determine the maximum fraction of genetic variation, indicating the major sources of codon usage variation. In addition, according to the results of COA, the correlation index between Axis1 and codon usage exponent, including the GC content of codons, GC3s, codon adaptation index (CAI) and the total numbers of amino acids in the encoded polypeptide (L_aa) were computed by python scripts package (https://github.com/shexuan/codon_analysis). CAI value is widely applied to assess gene expression levels, ranging from 0 to 1. Specifically, the larger the CAI value is, the stronger the codon usage preference is, and vice versa (Sharp & Li, 1986).

Characteristics of codon usage bias

RSCU and RFSC

The cp genomes of the seven Miscanthus and related species have 30 common codons (RSCU > 1) with 28 codons ending with the nucleotides A/T (93.3%) (Table S1). Therefore, almost the majority of codons of the seven plants species analyzed (RSCU > 1) are likely to end with A/T. The variation ranges in the RSCU values were close in the seven cpDNAs, i.e., 0.31–1.93 in M. floridulus, M. giganteus and M. transmorrisonensis, 0.32–1.94 in M. sacchariflorus and M. sinensis, 0.32–2.01 in Saccharum spontaneum and 0.33–2.04 in Sorghum bicolor, respectively (Table S1). In addition, the maximum and the minimum RSCU values belonged to TTA and CTG which encode Leu, indicating the vitally positive bias. Furthermore, the pattern of codon usage were summarized in the seven Miscanthus and related species (Fig. 1). Specifically, the high-frequency codons of seven Miscanthus and related species possess strong common base and share a total of 21 high-frequency codons (Table S1). Furthemore, Saccharum spontaneum and Sorghum bicolor were determined to possess two more high-frequency codons, specifically CGT and TAA, than the other five Miscanthus species.

Determination of optimal codons

The ENC values of each CDS were ranked and 10% of genes from both ends were selected to establish high and low expression gene banks respectively. The RSCU values and ΔRSCU values in the two expression libraries were calculated and are listed in Table S2. According to the values of ΔRSCU, the optimal codons in the seven assessed species were determined as follows (Table 2).

Codon usage frequency

The codon usage frequencies of the seven cpDNAs were compared with four model species including Escherichia coli, Saccharomyces cerevisiae, Arabidopsis thaliana and Populus trichocarpa (Table S3). The Results of our analyses indicated that there is little divergence in the codon usage frequencies among the seven assessed plant species with Saccharomyces cerevisiae, Arabidopsis thaliana and Populus trichocarpa, have 9 to11 (accounting for 14.1%–17.2% of total codons), 13 to 14 (20.3%–21.9%), 12 to 13 (18.8%–20.3%) different codons, respectively (Table S3). However, the codon usage frequencies of the seven species with Escherichia coli were comparatively higher (27 different codons). The results indicated that the codon frequency difference between Miscanthus species and Arabidopsis thaliana, Populus trichocarpa and Saccharomyces cerevisiae was the lowest, while was the largest with Escherichia coli. Based on above results, it was optimal to select Saccharomyces cerevisiae, Arabidopsis thaliana and Populus trichocarpa as heterologous expression hosts for Miscanthus and related species. Furthermore, the results indicated that TAG is a different termination codon in usage frequency when comparing the seven assessed plants with the four model species (Table S3).

Source analysis of variation in codon usage

ENC-plot

The ENC and GC3s of the seven analyzed cpDNAs were plotted. It can be seen from Fig. 2 that the ENC values of most genes were lower than expected values and lie below the standard curve. The results of ENC-plot analysis suggested that codon usage preference of the seven cpDNAs is mainly influenced by natural selection and other factors, while mutation pressure was determined to play only a minor role.

PR2-plot

PR2-plot is an efficient method to indicate the influence of mutation pressure by investigating the composition of A, T, C and G at the third nucleotide position of a codon. Our results revealed that the AT-bias is 0.464, 0.463, 0.463, 0.463, 0.463, 0.465 and 0.463 for M. floridulus, M. giganteus, M. sacchariflorus, M. sinensis, M. transmorrisonensis, Saccharum spontaneumand and Sorghum bicolor, while the GC-bias is 0.512, 0.513, 0.516, 0.515, 0.512, 0.515 and 0.518, respectively (Fig. 3). Therefore, T/G-bias was observed in all seven assessed species. When considered together, codon usage bias of A/T and G/C in the seven cp genomes was unbalanced, indicating that the base composition of the seven analyzed cpDNAs is not only influenced by mutation pressure, but also by natural selection.

Neutrality plot

The distribution range of GC12 (the mean GC content at the first and second nucleotide positions of a codon) and GC3 is relatively concentrated, in which the range of GC12 is 0.3272 to 0.5469, and the range of GC3 is 0.179 to 0.512 (Fig. 4). No significant correlation was found for GC1 with GC2 (r1 = 0.157, r2 = 0.168, r3 = 0.128, r4 = 0.127, r5 = 0.161, r6 = 0.155, r7 = 0.140), GC1 with GC3 (r8 = 0.092, r9 = 0.079, r10 = 0.055, r11 = 0.049, r12 = 0.100, r13 = 0.242, r14 = 0.202) and GC2 with GC3 (r14 = 0.054, r14 = 0.053, r15 = −0.014, r16 = −0.020, r17 = 0.063, r18 = −0.032, r19 = −0.032), which suggested mutation pressure only contributes a minor role in the codon usage preference. In addition, the regression coefficient (slope of neutrality plot) was 0.006 to 0.198, indicating that the correlation between GC12 and GC3 is not significant, and the composition of the first two bases may be different from the third base of the codon. These results demonstrated that the codon usage patterns of cp coding sequences in the seven species are mainly affected by natural selection.

Correspondence analysis (COA)

Correspondence analysis is used to explore the variations of codon usage among the analyzed cpDNAs. In the current study, RSCU-based COA was used to compare the usage patterns of 59 codons, which produced a series of orthogonal axes, reflecting the trend of change of codon usage in the seven Miscanhus and related plant species cpDNAs. The first four axes accounted for 36.2%, 38.2%, 38.5%, 38.3%, 36.8%, 40.7% and 40.9% of the overall changes, while the first axis proportion to 14.3%, 15.9%, 15.8%, 15.8%, 14.6%, 15.6% and 9.7% of the total variation in seven species respectively (Table 3). Axis 1, responsible for ∼10% of total variation, was the main source of variation, indicating that the codon usage was influenced by multiple factors. In addition, the relationship between axis 1 and axis 2 was visualized to explore the effects of GC content on codon usage bias (Fig. 5). Genes with different GC content are plotted as different colors, red with GC% < 45% and blue with 45% ≤ GC%< 60% (Fig. 5). In order to determine the factors leading to gene dispersion along axis 1 and axis 2, the correlation index were computed on axis 1 with CAI, CBI, Fop, GC3, GC and L_aa (Table 3). As can be seen from the results in Table 3, axis 1 for M. floridulus, M. sacchariflorus, M. sinensis, M. transmorrisonensis, Saccharum spontaneum and Sorghum bicolor possessed a remarkable correlation with GC3s (p ≤ 0.01), which indicated the base composition stemming from mutation pressure was the main factor impacting codon usage preference.