Evolutionary rate patterns of the Gibberellin pathway genes.

BACKGROUND: Analysis of molecular evolutionary patterns of different genes within metabolic pathways allows us to determine whether these genes are subject to equivalent evolutionary forces and how natural selection shapes the evolution of proteins in an interacting system. Although previous studies found that upstream genes in the pathway evolved more slowly than downstream genes, the correlation between evolutionary rate and position of the genes in metabolic pathways as well as its implications in molecular evolution are still less understood.
RESULTS: We sequenced and characterized 7 core structural genes of the gibberellin biosynthetic pathway from 8 representative species of the rice tribe (Oryzeae) to address alternative hypotheses regarding evolutionary rates and patterns of metabolic pathway genes. We have detected significant rate heterogeneity among 7 GA pathway genes for both synonymous and nonsynonymous sites. Such rate variation is mostly likely attributed to differences of selection intensity rather than differential mutation pressures on the genes. Unlike previous argument that downstream genes in metabolic pathways would evolve more slowly than upstream genes, the downstream genes in the GA pathway did not exhibited the elevated substitution rate and instead, the genes that encode either the enzyme at the branch point (GA20ox) or enzymes catalyzing multiple steps (KO, KAO and GA3ox) in the pathway had the lowest evolutionary rates due to strong purifying selection. Our branch and codon models failed to detect signature of positive selection for any lineage and codon of the GA pathway genes.
CONCLUSION: This study suggests that significant heterogeneity of evolutionary rate of the GA pathway genes is mainly ascribed to differential constraint relaxation rather than the positive selection and supports the pathway flux theory that predicts that natural selection primarily targets enzymes that have the greatest control on fluxes.

Background

A primary goal of molecular evolutionary studies is to elucidate the driving forces governing evolutionary change and mechanisms of molecular evolution. A general pattern arising from previous studies is that genes or proteins varied substantially in their evolutionary rates, spanning more than 3 orders of magnitude [1-4]. Although many determinants, including functional, biophysical, and fitness-related variables [3,5,6], have been proposed to explain such variation, what determines the evolutionary rate of a protein has been debated for decades and remains largely elusive despite some large-scale investigations [4,7,8]. One possible reason for this is that many studies have focused on protein interaction networks that are very heterogeneous with different types of interactions [5,8]. False correlations in protein interaction networks and noise in biological data further complicate the analyses of protein evolution [4,6]. To study evolutionary rates and patterns of proteins in well-characterized metabolic networks or pathways has thus become an alternative for understanding of molecular evolution [9-15].

Analysis of molecular evolutionary patterns of different genes within metabolic pathways allows us to determine whether these genes are subject to equivalent evolutionary forces and how natural selection shapes the evolution of proteins in an interacting system. In studies of evolutionary rates of genes in the plant anthocyanin pathway, Rausher et al. [11] and Lu and Rausher [13] demonstrated that upstream genes in the pathway evolved more slowly than downstream genes both over a broad taxonomic distance involving monocots and dicots and at the intragenic level between species within a genus. They suggested that such difference in evolutionary rates between upstream and downstream genes was due to more constraint upon the upstream genes because they participated in several different biochemical pathways. The hypothesis that earlier acting genes in genetic pathways are under strong purifying selection has been confirmed by some studies (e.g., [9,16]) but not been supported by the others (e.g., [12] and see [17]).

Another hypothesis regarding protein evolution involves a theory of pathway fluxes, indicates that natural selection would target enzymes that control metabolic fluxes and thus where selection operates in a pathway will depend on the distribution of flux control among pathway genes [17,18]. It has been suggested that pathway branch points are usually the targets of selection because they might control the expression of biochemical phenotypes disproportionately [10,14,18]. A third alternative hypothesis argued that translational selection and selection for protein folding and design might govern the rate of protein sequence evolution [7,19-21]. Despite these arguments and some other attempts [4,5,8,22], empirical studies with specific metabolic networks or pathways remain limited, particularly for plants.

Gibberellin (GA) biosynthetic pathway is involved in the production of gibberellins that control many aspects of plant growth and development, including seed germination, stem elongation, leaf expansion, and flower and seed development [23,24]. Despite a hundred GAs identified from plants, a small number of them such as GA1 and GA4 are found to function as bioactive hormones. These bioactive GAs are synthesized from trans-geranylgeranyl diphosphate (GGDP) through the 12-step conversion (Figure 1). To date 7 enzymes have been identified to be responsible for GA biosynthesis and function in multiple locations within the cell, involving the chloroplast, the ER membrane, and the cytoplasm. These 7 enzymes can be divided into 2 groups, i.e., the enzymes catalyzing early steps: ent-copalyl diphosphate synthase (CPS), ent-kaurene synthase (KS), ent-kaurene oxidase (KO), and ent-kaurenoic acid oxidase (KAO) and those for later steps: GA20-oxidase (GA20ox), GA3-oxidase (GA3ox), and GA2-oxidase (GA2ox) (Figure 1) [23]. In addition to a well-defined pathway, almost all of the genes encoding the 7 metabolic enzymes and their related mutants have been isolated and identified [23,24]. Therefore, the GA biosynthetic pathway provides an excellent system for exploring the correlation between evolutionary rate and position of the genes in metabolic pathways as well as its implications in molecular evolution. In addition, the GA biosynthetic pathway is also a good system to investigate the relative importance of selective constraint and positive selection as well as how adaptive variation is distributed among proteins in metabolic pathways. In this study, we isolated and sequenced portions of the 7 GA pathway genes from representative species of the rice tribe and analyzed their rates and patterns of molecular evolution. Our specific goals were 1) to compare evolutionary rates across the GA pathway genes and ask if evolutionary rate of the genes depends on its function or gene's position within the pathway; 2) to examine whether the GA pathway genes exhibit evidence of positive selection and particularly to understand the relative contribution of selective constraint and positive selection to variation of evolutionary rates of the GA pathway genes.

Figure 1

Simplified GA biosynthetic pathway in plants. The shaded metabolites have been demonstrated to function as bioactive hormones.

Methods

Plant materials and genes analyzed

We isolated and sequenced portions of all 7 GA pathway genes from members of the rice tribe (Oryzeae). Phylogeny of the rice tribe has been well resolved by previous studies based on sequences of plastid, mitochondrial and nuclear genes [25,26] and provided an important framework for molecular evolution study for this group (Additional file 1). To analyze the patterns of molecular evolution of the GA pathway genes, 8 diploid species were selected to represent major phylogenetic lineages of Oryzeae that are separated by a range of genetic distances [26]. They consisted of 5 Oryza species representing 5 diploid genome types, Oryza sativa (A), O. officinalis (C), O. australiensis (E), O. brachyantha (F), O. granulata (G), and one each of other 3 genera in the tribe (Chikusichloa aquatica, Luziola leiocarpa, and Rhynchoryza subulata). One species in the tribe Ehrhartoideae that is sister to Oryzeae, Ehrharta erecta, was used as an outgroup [26,27].

Of 7 enzymes in the GA biosynthetic pathway, 4 enzymes in the early steps (CPS, KS, KO and KAO) were encoded by single-copy genes (CPS1, KS1, KO2 and KAO) [23] and sampled in this study. Because the 3 enzymes in the later steps (GA20ox, GA3ox and GA2ox) were encoded by small gene families [23,24], we chose to sequence GA20ox2, GA3ox2 and GA2ox4 that were highly expressed in all organs of the wild-type rice [23]. Previous studies showed that mutants with all 7 enzymes except for GA2ox exhibited a typical phenotype of dwarfism, indicating the functional importance of these enzymes in GA biosynthesis [23]. Detailed information of the genes and the sequences and positions of the primers are listed in table 1 and Additional file 2, respectively. Species sampled and the GenBank accession numbers of the sequences obtained in this study are listed in Additional files 3 and 4, respectively.

Table 1

Information on the genes and sequences of the primers used in the study.

Gene	Locationa	Enzyme	Aligned coding length (bp)	GC (GC3)b	ENCb	Primer sequencesc
CPS1	2	ent-copalyl diphosphate synthases I	1002	0.462 (0.494)	54.6	F1: AACTTGTGGAGGTTAGCAGF2: TGTGGAGGYTAGCRGAGGR1: AAGTCGCTCAGAGGCACGR2: TAGCCCATGCAAGTCGCTC
KS1	4	ent-kaurene synthase I	1023	0.426 (0.417)	52.8	F1: TGCTGAAGCTTCCAGTTTCCF2: TCCAGTTTCCGTGAATCACR1: CTTGCACATCTTCCAGAACR2: CCTTGACGACTGCATTCAC
KO2	6	ent-kaurene oxidase II	1050	0.514 (0.639)	58.2	F: CTGTAGTTGTGCTCAATTCR1: GCCATCGTCTTGTACATGTCR2: TCAGCCTCCACYCGAACTC
KAO	6	ent-kaurenoic acid oxidase	1053	0.603 (0.820)	39.2	F: CAGGACGTTCATGTTCAGCAGR1: TCGTCGCCAAGCAGTTGTCR2: GCCAAGCAGTTGTCCAC
GA20ox2	1	GA20-oxidase II	597	0.714 (0.966)	30.1	F: ATCCCGGAGCCATTCGTSTGR: TGAAGGTGTCGCCGATGTTG
GA3ox2	1	GA3-oxidase II	786	0.710 (0.965)	30.4	F1: ACCCGCTCTRCTTCGACTTCF2: GGCGGGTGCCGGAGACGCACGR: CCATGTACTCSGGCCACGTGAC
GA2ox4	5	GA2-oxidase IV	819	0.694 (0.926)	34.0	F: GAGCAGATCTCGCTGSTGAGR1: CAGGCGGTTGTCGCYGAGR2: AGGCGGGAGAGGTAGGCAG
Total			6330

Chromosomal locations in O. sativa.

Average of 8 Oryzeae species.

Locations of the primers for each fragment are showed in Additional file 2. The internal primers were used for cycle sequencing are listed in Additional file 5.

DNA extraction, amplification and sequencing

Total DNA was isolated from silica-gel dried or fresh leaves using the method as described by Ge et al. [25]. Because the DNA fragments of each gene contained multiple exons and introns, additional internal primers were used for sequencing of four genes (CPS1, KS1, KO2 and KAO) that were about 2 kb in length. Their sequences and melting temperatures (Tm) are provided in Additional file 5. Polymerase chain reaction (PCR) was performed in a total volume of 25 μl which contained 5–50 ng of template DNA, 0.2 μM of each primer, 200 μM of each dNTP, 10 mM Tris-Cl (pH 8.3), 50 mM KCl, 1.5 mM MgCl2, and 0.75 U ExTaq DNA polymerase (TaKaRa, Shiga, Japan). Amplification was carried out in a T gradient 96 U thermocycler (Biometre, Göttingen, Germany) as follows: 3 min at 94°C, followed by 33 cycles of 30 ~ 35 s at 94°C, 30 ~ 35 s at 54°C ~ 58°C, 2 min 30 s at 72°C and a final extension at 72°C for 10 min (CPS1, KS1, KO2, and KAO). Because GC contents of GA20ox2, GA3ox2 and GA2ox4 were as high as about 70%, routine PCR amplification did not work well. We performed the amplification using LA Taq DNA polymerase, special GC Buffer I and dNTP Mixture (TaKaRa, Shiga, Japan). Except for annealing temperature at 60°C, the protocol was basically similar to the above protocol. All PCR products were separated by electrophoresis on 1.5% agarose gels stained with ethidium bromide. After purification using DNA Purification kit, most of amplification products were sequenced directly. For those PCR products that were weakly amplified and difficult to be sequenced directly, we used the pGEM-T Vector (Promega Corporation, Madison, USA) to clone and sequence at least three clones. Sequencing was carried out on a Megabase 1000 automatic DNA sequencer (Amersham Pharmacia Biotech). All sequences have been deposited in GenBank, and their accession numbers are EF577637 ~ EF577669, EU179376 ~ EU179435 (Additional file 4).

Sequence analyses

Sequences were translated into the predicted amino acid sequence and aligned using ClustalX with version 1.81 [28], followed by manual adjustment. Excluding intron regions, the stretches of amino acid residues were conserved well enough across eight species to perform unambiguous alignments. The possibility of sequence saturation was examined using DAMBE 4.5.45 [29]. A plot of the number of transitions and transversions vs. divergence offers a visual display of substitution saturation, with an asymptotic relationship indicating the presence of saturation [29]. Pairwise number of synonymous and nonsynonymous substitutions per site (dS and dN) as well as nonsynonymous/synonymous substitution rate ratio (dN/dS) between the eight species were calculated for the coding regions of all 7 genes using the improved approximate method of Yang and Nielsen [30], implemented in YN00 (CODEML) of PAML version 4 [31]. This method considers two important evolutionary features of DNA sequence: transition/transversion rate bias and codon frequency bias. Based on the pruned phylogenetic tree of the eight species (Additional file 6), we also estimated the dN/dS value of each gene using CODEML of the PAML version 4 [31].

The magnitude of codon bias is often used to reflect the degree of selective constraint in a gene and variation of synonymous substitution rates among genes may be related to codon usage [32,33]. To measure the extent of codon usage bias, we estimated effective number of codons (ENC) using DnaSP version 4.10.9 [34], which varies between 20 and 61, with the lower the value, the more biased codon usage [35].

Comparing dS and dN among genes

The pairwise synonymous and nonsynonymous genetic distances of 8 species were calculated for each gene using the method of Yang and Nielsen [30]. Of total 28 species pairs, we chose 15 species pairs involving 5 Oryza species and other three representatives of the rice tribe to represent maximum divergence time (20 myr, Additional file 1). To evaluate the reliability of estimated genetic distances of the 15 pairs, we checked the distance variance of each pair of the genes and found that the pairs involving O. sativa or O. brachyantha have high variance for GA20ox, GA3ox and GA2ox. Thus six species pairs involving either of the two species were excluded. Based on the rest 9 species pairs (O. officinalis vs. C. aquatica, O. officinalis vs. L. leiocarpa, O. officinalis vs. R. subulata, O. australiensis vs. C. aquatica, O. australiensis vs. L. leiocarpa, O. australiensis vs. R. subulata, O. granulata vs. C. aquatica, O. granulata vs. L. leiocarpa, and O. granulata vs. R. subulata), the mean distance and its standard error of each gene were calculated to show variation of dS, dN among genes. To compare dS and dN among genes, we calculated the dS and dN for each branch of the pruned tree of the six species (Additional file 1 and Additional file 6) using CODEML of the PAML version 4 [31]. Of the tree, there are total 10 evolutionary independent branches (Additional file 6). Two data matrices for seven genes were obtained for dS and dN, respectively. We used two approaches to detect whether evolutionary rates differed among the genes. First, we used multiple range test and Fisher's least significant difference (LSD) procedure to discriminate the means of the genetic distances of 9 species pairs among 7 genes. To compare dS and dN among genes, we calculated the dS and dN for each branch of the pruned tree including six species (Additional file 6) using CODEML of the PAML version 4 [31]. On the tree, there are 10 evolutionarily independent branches (Additional file 6). Then we used Wilcoxon signed rank test [36] in the pairwise comparisons of distances of the 10 branches between genes to test significant difference among the GA genes.

Comparing dN/dS among genes

Pairwise dN/dS ratios of 8 species were estimated for each gene using the method of Yang and Nielsen [30]. The mean dN/dS ratios of 9 species pairs and standard errors were calculated to show variation of dN/dS among genes. We used a method described by Lu and Rausher [13] to compare dN/dSratios for the seven genes. For this analysis, we first estimated a single value of dN/dS for each gene using the model M0 of by CODEML of PAML version 4 [31]. We then detected the significance of differences of dN/dS values between genes by comparing the likelihood of the model using the estimated value of dN/dS to the likelihoods of the same model and using dN/dS constrained to various values [13].

Detection of positive selection

The ratio of nonsynonymous to synonymous substitution rate (ω = dN/dS) provides an effective measure to detect selection or selective pressure on a gene or gene region, with ω < 1, = 1, and > 1 indicating purifying selection, neutral evolution, and positive selection, respectively [31,37]. We performed likelihood-based analyses using the CODEML program of PAML version 4 [31] to explore the selective processes acting on the GA biosynthetic genes. First, we tested whether the evolutionary rates of each gene in the pathway differed among lineages within the gene tree by using the branch models. The one ratio model (M0) assumes a single ω for all branches and all sites, whereas the free ratio model (Mf) assumes an independent ω ratio for each branch of the tree. A likelihood ratio test (LRT) was conducted to determine whether there was statistically significant difference between two models. If the LRT is significant, the null hypothesis that two models are not significantly different is rejected, and the model with higher likelihood value is assumed to be a better model [38,39].

We next used site-specific models to examine whether particular amino acid residues were subject to positive selection because the branch test often has little power to detect positive selection due to the fact that the ω ratio is seldom detected greater than 1 if all the sites are averaged [40]. The nested codon models [39,40] were performed for each gene separately. In addition to one ratio model (M0), the neutral model (M1a) classifies all the sites into 2 categories, one under strict constraint (0 < ω < 1) (purifying selection) and the other under neutral (ω = 1). Positive selection model (M2a) is based on M1a and assumes a third category under positive selection (ω > 1). Beta model (M7) assumes a beta distribution of the ω ratios, and beta&ω model (M8) extends an independent ratio estimated by the data. Models assuming positive selection M8 and M2a are compared with null models M7 and M1a, respectively. Positive selection is invoked if the LRT is significant and there is site with ω > 1 [39,40].

Results

Sequences of the 7 GA pathway genes were isolated from all sampled Oryzeae species and an outgroup, Ehrharta erecta. The sequenced regions of the 7 genes ranged from 699 bp to 2231 bp in rice and their aligned coding regions varied between 597 bp and 1053 bp for species used in this study (table 1 and Additional file 4). The GC contents for the total and 3 individual codon positions were similar for the same gene across species but varied greatly among genes, with the means ranging from 0.426 (KS1) to 0.714 (GA20ox2) (table 1 and Additional file 4). It is noted that GC contents of the 3 genes at the later steps were much higher than those of the 4 genes at the early steps in the pathway, particularly at the 3rd position (GC3) where the averages were 0.926 ~ 0.966 for the genes at the later steps in contrast to the values of 0.417 ~ 0.82 for the genes at the early steps (table 1 and Additional file 4). To estimate accurately and compare substitution rate and ω ratios of different GA pathway genes, we evaluated sequence saturation of substitution for each of the seven genes. Results of saturation plots did not reveal any mutational saturation for all the seven GA genes, indicating that the datasets were suitable for further analyses of the differences of substitution rate and ω ratios.

Comparing dN and dS among genes

The averages of dS and dN of the 9 species pairs for each gene were shown in Figure 2 and table 2. Both dS and dN values varied greatly among genes, with the highest dS (0.747 for GA2ox4) being 3.5 times the lowest (0.214 for KS1) and the highest dN (0.086 for KS1) being 5.9 times the lowest (0.015 for GA20ox2). For dS, the 3 genes in the early steps of the pathway (CPS1, KS1 and KO2) have significantly lower dS values than KAO and 3 genes at the later steps of the pathway (GA20ox2, GA3ox2 and GA2ox4). For dN, GA20ox2 has the lowest value (0.015), followed by 4 genes (KO2, GA3ox2, CPS1, and KAO), and the second and last genes in the pathway, KS1 and GA2ox4, have the highest and comparable dN value (0.086 and 0.087). Unlike dS values, dN varied significantly both between and within genes in the early and late steps of the pathway. It is noted that 2 types of substitutions did not correlated (r2 = 0.009, P = 0.840), which is clearly reflected by the fact that KS1 had a highest dN value but the lowest dS (Figure 2 and table 2). Multiple range tests generated 2 and 3 homogenous groups of genes for dS and dN, respectively, with statistically significant difference between genes among groups and non-significant difference between genes within groups (table 2). To overcome the rate heterogeneity of dS and dN among lineages, we performed the Wilcoxon signed rank test for pairwise comparisons of dS and dN between genes for 10 independent branch of pruned tree consisting of 6 species. For dS, the gene KAO and three late-step genes (GA20ox2, GA3ox2, GA2ox4) showed significant rate differences from the first 3 genes in the pathway (CPS1, KS1 and KO2) (P < 0.01). No significant difference was detected between KAO and the late-step genes and among the early-step genes. For dN, none of the tests showed significant difference, involving the comparisons between and within both early- and late-step genes (table 2 and Additional file 7).

Figure 2

Evolutionary rates at sysnonymous and nonsynonymous sites of the GA pathway genes. Horizontal bars indicate the averages of dS, dN and dS/dN values. Boxes represent 95% confidence intervals.

Table 2

Tests on difference of dS, dN and dN/dS among the 7 GA genesa

Gene	dS		dN		dN/dS
	Mean ± SEb	WSRTc	Mean ± SEb	WSRTc	Mean ± SE	ωd
CPS1	0.292 ± 0.016	A	0.040 ± 0.002	B	0.139 ± 0.007	C (0.144)
KS1	0.214 ± 0.023	A	0.086 ± 0.006	C	0.420 ± 0.030	D (0.386)
KO2	0.250 ± 0.008	A	0.050 ± 0.001	B	0.201 ± 0.009	C (0.177)
KAO	0.600 ± 0.023	B	0.037 ± 0.002	B	0.062 ± 0.003	B (0.074)
GA20ox2	0.556 ± 0.029	B	0.015 ± 0.001	A	0.027 ± 0.002	A (0.022)
GA3ox2	0.699 ± 0.063	B	0.044 ± 0.004	B	0.066 ± 0.007	B (0.073)
GA2ox4	0.747 ± 0.092	B	0.087 ± 0.010	C	0.119 ± 0.006	C (0.129)

a Capital letters (A to D) indicate the homogenous groups.

b Calculated based on 9 specie pairs.

c Wilcoxon signed rank test using 10 independent branches of the pruned tree consisting of 6 species

d ω, calculated as dN/dS ratios using PAML model M0.

The dN/dS ratios of the 9 species pairs varied remarkably among the 7 genes, with the highest value (0.420 for KS1) being 15.6 times the lowest (0.027 for GA20ox2) (Figure 2 and table 2). We compared dN/dS ratios for the 7 genes and detected the significance of differences between genes using the method of Lu and Rausher [13]. The dN/dS value for each gene ranged from 0.386 (KS1) to 0.022 (GA20ox2) with a trend similar to that estimated by pairwise dN/dS ratios (table 2). Based on log-likelihood ratio comparisons, we found 4 homogenous groups, i.e., A (GA20ox2), B (KAO and GA3ox2), C (KO2, CPS1, and GA2ox4), and D (KS1). As shown in table 3, significant difference of dN/dS ratios is found between genes from different groups but not found between genes from the same groups. The results are similar to those obtained by multiple range test and Wilcoxon signed rank test, with the exception that KO2 and CPS1 were in the same group in log-likelihood ratio comparisons (table 3).

Table 3

Comparison and significant test of log-likelihood ratio of dN/dS between genes

Category/Gene	Estimated dN/dS	lnL	Fixed dN/dS	lnL	2ΔLa
Between groups D and C
KS1	0.386	-3341.16	0.2815	-3345.46	8.59**
KO2	0.177	-3118.44	0.2815	-3126.22	15.56***
CPS1	0.144	-3034.22	0.2815	-3049.87	31.3***
GA2ox4	0.129	-2733.05	0.2815	-2755.45	44.8***
Within group C
KO2	0.177	-3118.44	0.153	-3119.19	1.50NS
CPS1	0.144	-3034.22	0.153	-3034.35	0.26NS
GA2ox4	0.129	-2733.05	0.153	-2734.14	2.18NS
Between groups C and B
KO2	0.177	-3118.44	0.1015	-3128.93	20.98***
CPS1	0.144	-3034.22	0.1015	-3038.05	7.66**
GA2ox4	0.129	-2733.05	0.1015	-2735.16	4.21*
KAO	0.074	-3071.99	0.1015	-3075.8	7.62**
GA3ox2	0.073	-2210.46	0.1015	-2213.21	5.50*
Within group B
KAO	0.074	-3071.99	0.0735	-3071.99	0NS
GA3ox2	0.073	-2210.46	0.0735	-2210.46	0NS
Between groups B and A
KAO	0.074	-3071.99	0.0475	-3078.8	13.62***
GA3ox2	0.073	-2210.46	0.0475	-2214.72	8.52**
GA20ox2	0.022	-1373.16	0.0475	-1379.06	11.8***

a d.f. = 1, ***P < 0.001, ** P < 0.01, *P < 0.05, NS not significant.

Codon usage bias and its correlation with GC3 and substitution rates

Estimates of the codon usage indicated that the ENC for the same gene was similar across all species (χ2 test, P = 1.00) (Additional file 4), suggesting that the magnitude of codon bias underwent little change during the diversification of the rice tribe. However, substantial variation was found among genes in the degree of codon bias, with the means of ENC ranging from 30.1 (GA20ox2) to 58.2 (KO2) (table 1 and Additional file 4). Pairwise comparisons of one-way ANOVA found that differences of the ENC were significant for all pairwise comparisons (P < 0.001) except for CPS1 vs. KS1 (P = 0.090) and GA20ox2 vs. GA3ox2 (P = 0.386). Note that the 3 genes in the late steps (GA20ox2, GA3ox2, and GA2ox4) exhibited strong codon usage bias (ENC = 30.1 ~ 34.0), whereas the 3 genes in the early steps (CPS1, KS1 and KO2) possessed little codon usage bias (ENC = 52.8 ~ 58.2). KAO exhibits an intermediate degree of codon usage bias with ENC of 39.2. ANOVA showed that ENC of the late-step genes differed significantly from those of the early step genes (P < 0.001). To further investigate the correlation between evolutionary rate and expression level, we analyzed the correlation between each of three parameters (dN, dS, and dN/dS) and ENC that negatively correlates with expression level. Results showed that significantly negative correlation existed between dS and ENC (P < 0.001), but not between ENC and either dN (P = 0.540) or dN/dS (P = 0.103).

Testing for positive selection

Significant rate heterogeneity among the GA pathway genes, particularly for dN/dS ratios that span more than one order of magnitude (Figure 2), might result either from intense purifying selection on slowly evolving genes (e.g., GA20ox2) or from frequent episodes of positive selection on fast evolving genes (e.g., KS1). Therefore, we first used branch models to test whether the evolutionary rate of each gene differed among lineages within the rice tribe. For all 7 genes, free ratio model (Mf) did no have significantly higher likelihood scores than one ratio model (M0) (P > 0.05) (Additional file 8). The ω values were estimated to be 0.022 ~ 0.386 under M0 model, suggesting that purifying selection or selection constraint best explains the evolution of these genes.

Because the branch model test averages the ratio across all sites and is a very conservative test of positive selection [39,40], we applied site-specific codon models to the 7 genes separately, to determine whether there was positive selection on codon sites of these genes. Results showed that, for all 7 genes, models M2a and M8 assuming positive selection were not significantly better than the null models M1a and M7 (for M1a vs. M2a, 2ΔL = 0.0 ~ 3.72, 0.1