The CDS genomic sequences of T. urartu accession G1812, wild tetraploid wheat (T. turgidum) accession “Zavitan,” Ae. tauschii accession AL8/78, and common wheat cultivar “Chinese Spring” were downloaded from the Ensembl Plants.¹ The CDS sequences were filtered using Python and the following criteria: (1) the gene is a complete CDS, (2) the length of the CDS sequence is more than 300bp, and (3) the CDS begins with ATG and ends with a termination codon (Liu et al., 2020). The leaf transcriptome data in common wheat and its ancestors were collected from the NCBI² under accessions SRR3274670 (T. urartu), SRR11292316 (T. turgidum), SRR11292312 (Ae. tauschii), and SRR9647023 (T. aestivum). In addition, the sequences of genes for grain size (GS; Li and Yang, 2017), the CCT gene family (Zheng et al., 2017), and the prolamin genes (Ling et al., 2018) identified from previous studies were used. The collected fasta sequences of genes from different families were aligned to the local protein database of other species using BALSTX, which was downloaded from URGI database³ with the expected value (E-value) of the lowest. Meantime, the conserved domains of the gene family were identified using HMMER 3.0.

The codon usage of the T. aestivum genome and its ancestor species were analyzed using the software program CodonW 1.4.2⁴ (Peden, 2000). The parameters for the codon usage calculation were listed as following: Relative Synonymous Codon Usage (RSCU), Effective Number of Codon (ENC), GC content (GC), codon first and second and the third position for GC content (GC1, GC2, and GC3), the average GC content of the first and second codons (GC12) and the content of each base in the third position of the synonymous codon (A3s, T3s, G3s, and C3s).

(1) RSCU:

Among the synonymous codons, the codon with the larger RSCU value has a higher usage probability (Wright, 1990). RSCU values were estimated using the following formula:

(2) ENC:

A value below 35 indicates high codon usage bias and values above 50 indicates low bias (Sharp and Li, 1987). ENC were calculated by the following formula:

Effective Number of Codon-plot and neutrality plot analysis were used to analyze factors affecting the codon preference of T. aestivum and its ancestor species. In ENC-plot analysis, when ENC value is within the allowable range of the standard curve, the codon preference is caused by base mutation, otherwise the codon preference is caused by natural selection (Cai et al., 2006).

Neutral plot analysis was used to investigate the effect of mutational and selection pressure on the codon usage pattern. The neutrality plot was mapped with the scatterplot to analyze the correlation between GC12 and GC3. When the regression coefficient was close to 0, the correlation was low, indicating that the GC composition for the three positions of the codon differed, the GC content is highly conserved, and the codon preference for the species is mainly affected by natural selection (Su et al., 2009; Nasrullah et al., 2015). Otherwise, the codon bias of species is more affected by mutations.

PR2-plot was used to analyze whether the third base of a codon is biased according to the bias rule (Sueoka, 2001). The bias rule states that if A=T and G=C between the two complementary strands of DNA then there is no mutation and selection bias. If otherwise, then the codon is biased.

Correspondence analysis (COA) was generated to analyze codon usage changes in different genes. Correlation analysis is used to determine the main factors affecting codon bias (Jia et al., 2015). In COA, all the genes were distributed in multi-dimensional vector spaces using the RSCU COA function of codonW. The first axis represented the most different codon usage changes with sequential decreases represented in the second, third, and fourth axes.

To determine optimal codons of T. aestivum genome and its ancestor species, genes with the 5% highest and lowest ENC values were screened to establish a low preference group and a high preference group. Among genes from two groups, the RSUC difference was greater than 0.08, and the codon with the RSUC value greater than 1 was determined to be the optimal codon (Liu et al., 2020).

To determine the factors affecting the codon bias of the T. aestivum genome and its ancestor species, we analyzed the number of tRNA anti-codons and the length of genes. Firstly, CDS sequences of common wheat and its ancestor species were divided into four length categories: <1,000bp, 1,000–2,000bp, 2,000–3,000bp, and >3,000bp. Secondly, codon bias analysis for genes of different length was done using CodonW 1.4.2 and the numbers of genes and values of GC3s were calculated. Anti-codons of common wheat and its ancestor species were identified using tRNAscan-SE 2.0 (Chan and Lowe, 2019).

The rawdata of common wheat and its ancestors were used for quality control using fastp (Chen et al., 2018). To analyze the relationship between codon preference and function of expressed genes, the TPM of CDS sequences from common wheat and its ancestors were calculated using Salmon v1.3.0 (Patro et al., 2017). To determine differences in the function of genes with different codon bias, gene ontogeny (GO) enrichment were conducted separately using genes with different codon bias. For gene annotation, coding sequence regions of each gene in wheat and its ancestor species reference genomes were compared against the default database using eggNOG-mapper (Huerta-Cepas et al., 2017, 2019). The CDS sequences of the TPM>1 of all species were divided into high GC bias genes and high AT bias genes according to GC content of each species. The CDS sequences with different codon bias were calculated via GO enrichment analysis using the clusterProfiler R package (Yu et al., 2012).

The codon usage characteristics of the T. aestivum genome and those of its ancestor species, T. urartu, T. turgidum, and Ae. tauschii, were investigated using CodonW1.4.2. Triticum aestivum showed a stronger codon usage bias (ENC 48.6) compared to its ancestor species T. turgidum (ENC 51.6) and Ae. tauschii (ENC 51.2; Supplementary Table 1). ENC values of the coding genes in the A and B genomes were compared. Among the 57,979 coding genes from the A genome of T. turgidum, 3,874 sequences (6.68%) showed high codon usage bias with ENC values less than 35. In genes from the A genome of T. aestivum, there were 5,606 sequences (13.66%) with ENC values less than 35. Sequences with ENC values less than 35 accounted for 6.10% of the T. turgidum B genome and 12.53% of the T. aestivum B genome. Thus, both the A and B genome from T. aestivum showed more codon usage bias than the A and B genome from T. turgidum.

GC3s frequencies in T. aestivum (63.7%) and T. urartu (62.0%) were higher than those in T. turgidum (52.6%) and Ae. tauschii (54.6%; Figure 1; Supplementary Table 1). Codon usage bias has been shown to correlate with GC composition (Wan et al., 2004). The mean GC contents were examined for the whole coding sequences of the four species and the higher GC content partly explained the bias for G/C-ending codons seen in T. aestivum and T. urartu. PR2-bias plot analysis showed a non-proportional use of A/T and C/G in the four genomes, with T3>A3 and C3>G3 (Supplementary Figure 1). Neutrality plot results showed that the correlation between GC12 and GC3 was statistically significant (p<0.001) in T. aestivum (r=0.650), Ae. tauschii (r=0.596), T. turgidum (r=0.551), and T. urartu (r=0.517). The slopes of the regression line were close to zero (0.16–0.20, Supplementary Figure 2), indicating that selective pressure should be the dominant factor shaping codon usage pattern in T. aestivum and its ancestor species.

In the present study, COA based on RSCU was expected to further identify the major factors affecting codon usage frequencies and synonymous codon preferences observed in the coding sequences of the four genomes. Genes from the four species were widely distributed along the two explanatory axes, suggesting a strong codon usage bias in these coding sequences (Figure 2). The first and second explanatory axis accounted for 40.2/4.3, 39.1/3.6, 47.9/4.1, and 40.6%/3.7% in T. urartu, T. turgidum, T. aestivum, and Ae. tauschii, respectively. Thus, the first axis reflects the primary factor that explains the differences in codon usage among all genes and genes with high GC and high AT were separated along this axis. The corresponding distribution of synonymous codons shows the separation of C/G-ending codons (except AGG that encodes arginine) and A/U-ending codons along this same axis. High G/C genes from T. urartu and T. turgidum were distributed on the right side of axis 1. The separation of genes on the second axis appears to be largely due to frequency differences in C-ending and G-ending codons among the GC rich genes; however, genes with more A-ending and T-ending codons tended to cluster at the center of axis 2. Correlation analysis showed that the third base preference of coding genes positioned on axis 1 of the COA were significantly related (Supplementary Table 2), indicating that base composition is the main factor affecting codon usage and the COA could be further used in optimal codon identification.

The optimal codons of common wheat and its ancestor species were identified based on the RSCU value (Supplementary Table 3) and the number of optimal codons in T. urartu, T. turgidum, T. aestivum, and Ae. tauschii were 21, 14, 22, and 16, respectively. All the optimal codons showed G/C at the third codon position, indicating that wheat and its ancestor species have similar optimal codon preferences for codons with G/C-endings instead of A/U-endings.

To investigate the effect of coding sequences length on codon usage, coding sequences were classified into four groups: <1,000bp, 1,000–2,000bp, 2,000–3,000bp, and >3,000bp, and the average GC3s values computed for genes from each species (Figure 3). Bias for G/C-ending synonymous codon usage was much stronger in genes coding for short vs. long proteins (Supplementary Table 4). Among coding sequences longer than 1,000bp, more GC-ending synonymous codons were observed in T. aestivum than in its ancestor species. However, among sequences <1,000bp, T. urartu sequences showed the greatest bias for GC-ending codons. Thus, in longer genes T. aestivum showed the greatest tendency to use G/C at the third codon position.

In humans, GC distribution and codon preference have been proved to be related to gene density (Versteeg et al., 2003). In the present study, substantially higher gene density occurred in the distal regions of chromosome arms and genes distributed in these regions exhibited higher GC3s values. In contrast, the centromere regions showed lower gene density and the GC3s values of genes in these regions were low (Figure 4).

Biased usage of synonymous codons can also be constrained by the amount of isoacceptor tRNA to increase translation efficiency (Duret and Mouchiroud, 1999). The frequency of each synonymous codon encoding 21 amino acids was estimated, and we further investigated the amount of corresponding tRNA in the four species (Supplementary Table 5). Significant differences were found among the species for codon frequency and corresponding tRNA abundance. For example, T. urartu and T. aestivum showed high correlations (r=0.667 and 0.724, p<0.01) for the amino acids with two synonymous codons (Tyr, Phe, Lys, His, Glu, Gln, Cys, Asp, and Asn; Figures 5A,B). In contrast, amino acid with more than four synonymous codons (Val, Thr, Ser, Pro, Leu, Gly, Arg, and Ala), T. turgidum and Ae. tauschii showed high correlations (r=0.511 and 0.488, p<0.01; Figures 5C,D).

Based on GO categories, the functional enrichment of codon usage biased genes from the four species was analyzed, and the top 10 GO terms with values of p<0.05 for genes from each species were identified (Figure 6). G/C-ending genes performed different functions compared to A/T-ending genes. The GC3s biased genes from T. aestivum and its ancestor species were commonly enriched in genes related to photosynthesis (Figure 6A). Also enriched were transmembrane transport (72/625), and response to acid chemical (71/625) in T. urartu, photosynthesis (46/802) and photosynthesis light reaction (43/802) in T. turgidum, response to acid chemical (158/1,214) and response to inorganic substance (113/1,214) in T. aestivum, and photosynthesis (46/751) and photosynthesis light reaction (40/751) in Ae. tauschii (Supplementary Table 6). AT3s biased genes likely perform similar functions in T. aestivum and T. urartu (Figure 6B) as these genes were enriched in GO terms such as cellular location, macromolecule localization, and protein location. A/T ending genes from Ae. tauschii were mainly enriched in response to extracellular stimulus (50/836) and response to nutrient levels (47/836; Supplementary Table 7). The GO enrichment (Supplementary Tables 8 and 9) demonstrated that the evolutionary features reliably reflect the relative function of different codon biased genes.

The present study shows that synonymous codons had bias at the whole genome level in T. aestivum and its ancestor species. We next sought to investigate if codon usage preference changed during evolution for gene families with important functions, including grain size (GS) genes, the CCT gene family regulating flowering, and the quality related prolamin gene family. Based on 36 grain size related genes identified in previous studies in Triticum (Li and Yang, 2017), similar GS genes from the A genome of T. aestivum, T. turgidum, and T. urartu were isolated and the GC3s variation of GS genes during polyploidization was calculated (Supplementary Table 10). A difference in the GC3s value between the homologous genes of >0.1 was classified as being codon biased. Triticum urartu is the progenitor of the A subgenome of tetraploid T. turgidum and hexaploid T. aestivum (Ling et al., 2018). From T. urartu to T. aestivum, seven GS genes were identified with different preferences in codon usage patterns (Table 1). Five of the seven genes, TRIUR3_05970-T1, TRIUR3_33526-T1, TRIUR3_34310-T1, TRIUR3_08952-T1, and TRIUR3_09477-T1, showed codon usage bias from the time T. urartu evolved into T. turgidum (Table 2).

The RSCU of each codon was measured for CCT and prolamin genes to investigate the variation in codon usage bias of these genes across the four genomes. Consistent with the results at the genome level, the G/C-ending codons in most of the CCT family genes were overrepresented (RSCU>1), except CCT24 which showed obvious A/U bias in the third position of synonymous codons. Among the species, the codon usage in T. aestivum CCT genes was most like that for genes in the D genome ancestor Ae. tauschii (Figure 7). Euclidean clustering of RSCU for codons in each gene showed that among the 26 CCT genes in T. aestivum, 15 were like Ae. tauschii, eight were like T. turgidum, and the last three were like T. urartu (Figure 8). For the prolamin genes, some optimal codons identified at the genome level (Supplementary Table 5) were underrepresented (RSCU<1) including GUG encoding Val, UAC encoding Try, CCG encoding Pro, UUC encoding Phe, CUC encoding Leu, CAG encoding Glu, GAC encoding Asp, and CGC encoding Arg (Figure 9).

Codon bias is determined by a balance between mutation, drift, and selection for optimal translational efficiency and/or accuracy (De La Torre et al., 2015). In many cases, codon usage correlates with the cognate tRNA abundance with highly used codons generally having higher intracellular tRNA concentrations than low-use codons (Zalucki et al., 2009). This fact indicates that tRNA abundance may be the main selection pressure for synonymous codon usage (Mohanta et al., 2020). In the present study of Triticum, more codon-anticodon interaction occurred in amino acids with two synonymous codons in T. aestivum and T. urartu and in amino acids with four or six synonymous codons in T. turgidum and Ae. tauschii. This result also provides evidence for the similar codon usage bias in T. aestivum and T. urartu. In contrast, for codons encoding glutamine from T. aestivum and T. urartu, corresponding tRNA copy numbers for the optimal codon (CAG) were much less than for the non-optimal codon (CAA), suggesting that the primary factor influencing codon-usage for these amino acids is not translation efficiency (Kanaya et al., 2001). Simultaneously, we considered the influence of codon swinging on tRNA abundance and codon usage preference. Crick (1966) proposed the wobble hypothesis based on stereochemistry to explain this phenomenon: when the first anticodon is A or C, only one codon can be recognized; when the first anticodon is G or U, two codons can be recognized, where G can identify C and U, where U can identify A and G. Due to the wobble phenomenon, a tRNA anticodon can be combined with more than one mRNA codon. We conjecture that UUG oscillates when recognizing CAA, so that more UUG can recognize the optimal codon CAG. We also found the same phenomenon in Ser. Wobble rules say that tRNAs with a GCU anticodon can pair with AGU, which explains the apparent discrepancy in the ability to translate AGU codons to Ser with extremely few tRNAs with ACU anticodons (Supplementary Table 5). In a recent research of 128 species of the plants (Mohanta et al., 2020), high and low frequencies of wobble base-pairing were occurred at the G:U base pairing to meet translational demand.

Duret and Mouchiroud (1999) found that codon usage bias is negatively correlated with protein length in Caenorhabditis elegans, D. melanogaster, and Arabidopsis thaliana. Gene length showed a negative correlation with GC3 content in non-grass monocot species (Mazumdar et al., 2017). The same phenomenon was found in the three progenitor species, as coding sequences shorter than 1,000bp had the highest GC3 codon bias (Figure 3). In T. aestivum, though, synonymous codon bias increased among coding sequences of 1,000–2,000bp compared to short sequences (<1,000bp). In the T. aestivum genome, we found that most coding sequences were between 1,000 and 2,000bp (N=53, 051), which is consistent with the result for T. aestivum transcript sequences reported by Wang et al. (2019). All the optimal codons identified in the Triticeae in the present study preferred G/C-endings at the genome wide level, and the GC3 content of coding sequences shorter than 2,000bp was higher than that for longer genes. This result supports the idea that in T. aestivum sequences of less than 2,000bp are the optimal length for translation, especially for important functional genes, and that the use of optimal codons increases their translational efficiency and accuracy. The optimal gene length of T. aestivum may be longer than that of its ancestral species, and perhaps the process of evolving into a complex hexaploid species enabled T. aestivum to translate longer proteins with minimal mismatches.

Genes for particular physiological processes, such as photosynthesis, response to salt stress/acid chemical/cold stress, and seed germination, were biased to GC-ending codons. In contrast, A/U-ending biased genes tended to be associated with the function of basic biological processes, like RNA processing, protein localization, protein transport, and membrane organization (Supplementary Tables 8 and 9). Some important genes always maintained conserved codon usage bias. For example, we identified 1,755 disease resistance genes in Ae. tauschii (Supplementary Table 11; Jia et al., 2013) and compared the GC3s contents with the corresponding homoeologous genes in D genome of T. aestivum. Only 19 genes had a GC3 difference of more than 0.1, an exceptionally low level (Supplementary Table 12).

CCT family genes are important for the flowering process, which enables the adaptation for reproductive success in various geographical environments (Yang et al., 2013). RSCU values of 59 synonymous codons in the CCT gene family from the four species were investigated and most CCT genes from T. aestivum (15/26) were clustered together with genes from Ae. tauschii, indicating that CCT genes from T. aestivum and Ae. tauschii were similar in codon usage frequency. We further measured the GC3 values of CCT genes and found that T. aestivum and Ae. tauschii CCTs showed similar G/C preference at the third position of codons (Supplementary Figure 3A; Supplementary Table 13), which was consistent with the RSCU results. Previous studies showed that the rapid evolution of the CCT gene family promotes the adaptability of Ae. tauschii (Zheng et al., 2017), the similar codon usage bias implies that CCT genes from T. aestivum maintained the functions from its D progenitor, Ae. tauschii, presumably to enhance environment adaptability. The prolamin sequences of the four species are highly similar, but there are obvious differences in the usage of codons, especially the optimal codons (Supplementary Figure 3B; Supplementary Table 14). As a donor of the D genome of wheat, Ae. tauschii provided wheat with better adaptability and quality traits (Jia et al., 2013).

For codon optimization, not only must the codon preference and tRNA abundance be considered, but also other factors including local mRNA folding, codon coordination, and codon correlation. By both natural evolution and genetic engineering, these properties can be altered to gain additional transcriptional regulatory ability for appropriate gene expression under specific cell conditions.