Background
==========

The plant seed is not only an organ of propagation and dispersal but also the major plant tissue harvested and used either directly as part of the human diet or as feed for animals. At the present time there is concern over long term food security and the impact of the move towards meat-based diets that will lead to a significant increase in the demand for plant protein for animal feed \[[@B1]\]. The amount of protein present in plant seeds varies from \~10% of the dry weight in most monocot (e.g. *O. sativa*, *S. bicolor*, *S. italica*, *Z. mays*and *B. distachyon*) to more than 30% in most dicots (e.g. *G. max*, *R. communis*, *C. sativus*and *A. thaliana*), and forms a major source of dietary protein \[[@B2]-[@B7]\]. To determine whether differences in evolutionary patterns may explain the phenotypic differences observed, a comparative investigation of evolutionary divergence in genes underlying protein synthesis in these two groups of plants is thus warranted.

Seed storage proteins can be classified into four groups: albumins, globulins, prolamins and glutelins \[[@B8]\]. Albumins and globulins comprise the storage proteins of dicots, whereas prolamins and glutelins are the major proteins in monocots \[[@B4],[@B9],[@B10]\]. 2S albumins, a major class of dicot seed storage proteins, have been most widely studied in the Cruciferae, notably *B. napus*and *A. thaliana*\[[@B9],[@B10]\]. Prolamins, the major endosperm storage proteins of all cereal grains, with the exceptions of oats and rice, can be classified into many subgroups, e.g. sulphur-rich (S-rich), sulphur-poor (S-poor) and high molecular weight (HMW) prolamins \[[@B4]\]. The globulins, the most widely distributed group of storage proteins, are part of the cupin superfamily \[[@B11]\] and are evolved from bacterial enzymes. The globulins are present not only in dicots but also in monocots \[[@B9]\] and can be divided into 7S vicilin-type and 11S legumin-type globulins according to their sedimentation coefficients. It should be noted that the genes encoding the 11S-type globulins in monocots are the same gene family, 11S globulin family, as those in dicots; whereas the genes encoding the 7S-type globulins in monocots are not evolutionarily related to those in dicots \[[@B4]\]. Thus, we focus on the genes encoding the 11S-type globulins in this study.

Many efforts have been made to describe the gene families encoding seed storage proteins. For albumins, they are encoded by multi-gene families in many dicots (e.g., *A. thaliana*and *B. napus*); and evolutionary research into the gene families suggests that the albumin genes were duplicated prior to the Brassiceae-Sysimbrieae split, and gene duplication has played a role in their evolution \[[@B12],[@B13]\]. For prolamins, they are the major seed storage proteins in most grass species (e.g., *Z. mays*and *S. bicolor*); and studies have suggested that the prolamin gene families have undergone many rounds of gene duplication \[[@B14],[@B15]\]. For globulins, eight, four, eleven and fourteen gene have been identified and classified in *G. max*, *A. thaliana*, *R. communis*and *O. sativa*, respectively \[[@B16]-[@B24]\].

From the research described above, it can be concluded that the seed storage protein gene families have expanded in a lineage-specific manner through gene duplication. As a major process in the evolution, gene duplications can provide raw material for evolution by producing new copies. The human globin gene family is a representative example: several globin genes have arisen from a single ancestral precursor, thus making individual genes available to take on specialized roles, with some genes becoming active during embryonic and fetal development, and others becoming active in the adult organism \[[@B25]\]. Gene duplications may also affect phenotype by altering gene dosage: the amount of protein synthesized is often proportional to the number of gene copies present, so extra genes can lead to excess proteins. This applies to many kinds of genes, such as rRNAs, tRNA and histones \[[@B26],[@B27]\]. A critical question thus can be asked: does gene duplication contribute to the higher levels of protein synthesis in dicots than in monocots? On the other hand, gene duplication usually brings variation in evolutionary rate \[[@B26]\], and such variation has been predicted to be associated with phenotypic differences. For example, Hunt et al. \[[@B28]\] investigated the evolution of genes associated with phenotypically plastic castes, sexes, and developmental stages of the fire ant *Solenopsis invicta*, and argued that an elevated rate is a precursor to the evolution of phenotypic differences. Thus, we are also interested in another question: does an accelerated evolutionary rate play a role in the evolution of storage protein content?

To shed light on the two questions above, we investigated the process of molecular evolution of the 11S globulin gene family, which is widely distributed in dicot and monocot species, by comparing the differing evolutionary patterns in the two groups. Our analyses suggested that gene duplication and an accelerated evolutionary rate in 11S globulin genes may be associated with higher protein synthesis in dicots than in monocots.

Results
=======

Sequences retrieval and phylogenetic analysis
---------------------------------------------

We collect the sequences of 11S globulin genes through a COG method. This procedure is based a simple notion that, if any proteins from distant genomes are more similar to each other than to any other proteins from the same genomes, they are most likely to belong to an orthologous family. In such a family, there are two kinds of relationship between a pair of sequences, namely symmetrical and asymmetrical BeTs (the Best Hits). If the symmetrical and asymmetrical BeTs are linked respectively by solid and broken lines, the orthologous family would forms a network; and thus all other members in the network can be identified when one member was investigated \[[@B29]\]. The 11S globulin genes from the five dicots and the five monocots form a COG containing 56 sequences (Figure [1](#F1){ref-type="fig"}). Of these genes, four genes encoding 12S globulins \[[@B20]\] come from *A. thaliana*; six genes, *Gy1*-*Gy5*and *Gy7*\[[@B16]\], come from *G. max*; twelve genes, all the functional glutelin genes, come from *O. sativa*\[[@B24]\]; eleven genes described in \[[@B7]\] come from *R. communis*; and six, seven, six, one, one and two genes come from *C. sativus*, *P. trichocarpa, B. distachyon, Z. mays*, *S. bicolor*and *S. italica*, respectively.

![**The COG of 11S globulin gene family**. Solid lines show symmetrical BeTs (the Best Hits) and broken lines show asymmetrical BeTs. Genes from the same species are adjacent. Gene ID is indicated and the prefix \"Rc\" denotes IDs from *Ricinus communis*. Among these IDs, At1g03880, At1g03890, At4g28520 and At5g44120 are known to encode *CRB, CRU2, CRC* and *CRA1*, respectively; Glyma03g32030, Glyma03g32020, Glyma19g34780, Glyma10g04280, Glyma13g18450 and Glyma19g34770 to encode *Gy1-5*and *Gy7*, respectively; Rc29600.m000561, Rc29600.m000564, Rc30005.m001289, Rc30005.m001290, Rc29611.m000223, Rc29200.m000169, Rc29629.m001355, Rc29709.m001187, Rc29716.m000305, Rc29200.m000167 and Rc30005.m001288 to encode *RcLEG1-1*to *RcLEG1-5*and *RcLEG2-1*to *RcLEG2-6*, respectively; and Os01g55690, Os10g26060, Os03g31360, Os02g15169, Os02g15178, Os02g15150, Os02g16820, Os02g16830, Os02g14600, Os02g15070, Os02g25640 and Os02g15090 to encode *GluA-1*, *GluA-2*, *GluA-3, GluB-1a*, *GluB-1b*, *GluB-2*, *GluB-5*, *GluB-4*, *GluB-7*, *GluB-6*, *GluC-1*and *GluD*, respectively.](1471-2148-12-15-1){#F1}

In order to analyze the phylogenetic relationship of the 11S globulin genes from the above ten species, the NJ and Bayesian methods were used to reconstruct the phylogenetic tree. Similar results were achieved (data not shown), and the NJ tree is shown in Figure [2](#F2){ref-type="fig"}, in which there are three subfamily clades: dicot subfamily 1, dicot subfamily 2 and monocot subfamily. The dicot subfamily 1 contains all the genes from *A. thaliana*and *G. max*, two genes from *C. sativus*, two genes from *P. trichocarpa*and *RcLEG1*of *R. communis*; the dicot subfamily 2 contains the other genes from *C. sativus*, *P. trichocarpa*and *R. communis*; and the monocot subfamily contains all the genes from the five monocot species. In the three subfamilies, with the exception of rice, genes from the same species form monophyletic groups; and the phylogenetic relationship of the monophyletic groups is largely concordant with the species tree described by \[[@B30]\].

![**Phylogenetic relationships of sequences within the 11S globulin gene family by neighbor joining (NJ) method with bootstrap support above 50% shown at the nodes**. Letter *A*-*E*indicates the branches used in analysis of evolutionary rate and positive selection.](1471-2148-12-15-2){#F2}

Gene duplication in 11S globulin gene family
--------------------------------------------

In the 11S globulin gene family, there are four or more genes from *C. sativus*, *P. trichocarpa*, *R. communis*, *A. thaliana*, *G. max*, *B. distachyon*and *O. sativa*, suggesting that the 11S globulin genes in these species have undergone two or more rounds of duplication. Among these duplications, tandem duplication is a major type, and occurred in all the species above, e.g. Rc30005.m001289-Rc30005.m001290 and At1g03880-At1g03890 and Os02g16820-Os02g16830. Furthermore, most of the 11S globulin genes from the dicot species are located on the chromosome segments that share conserved synteny with each other \[[@B18]\], suggesting that segment duplication or whole genome duplication (WGD) is the major origin of duplicate genes.

Evolutionary rate of 11S globulin gene family
---------------------------------------------

To determine whether there were different evolutionary patterns across different subfamilies of the 11S globulin genes in monocots and dicots, the *ω*values for these genes were calculated by a branch-specific model.

In the branch-specific model, three cases were considered in this study. One is a two-ratio model that suggests there are distinct monocot (*ω*~0~) and dicot subfamilies (*ω*~1~), one is a three-ratio model that suggests there is one monocot (*ω*~0~) and two dicot subfamilies 1 (*ω*~1~) and 2 (*ω*~2~), and one is a six-ratio model that suggests there are subclades under the branches *A*-*E*(*ω*~1~\~ *ω*~5~) and the other (*ω*~0~). All the above models were favored over the one-ratio model by the likelihood ratio test (P \< 0.05, Table [1](#T1){ref-type="table"}) and in the two-, three- and six-ratio models, the *ω*estimate for the dicot subfamily is considerably higher than that for the monocot subfamily. These results suggest that the dicot 11S globulin genes are under reduced evolutionary constraints, and thus evolve at a higher evolutionary rate. We should point out that in the six-ratio model, the genes from *S. italica*, *Z. mays*and *S. bicolor*were not taken into account, because they may have a different evolutionary pattern due to the fact that 11S globulins are minor components \[[@B4],[@B31]\] and there are only one or two 11S globulin gene copies in these species; and that branches *A*-*E*were chosen because five branches can be grouped for the phylogenetic tree.

###### 

Evolutionary rate analysis of 11S globulin family using branch-specific model of PAML

  Model         *ω*setting                     \-*ln*L    Estimated parameters   Likelihood ratio test
  ------------- ------------------------------ ---------- ---------------------- -----------------------------------------
  One-ratio     entire tree: *ω*~0~            38214.64   *ω*~0~= 0.234          
  Two-ratio     monocot subfamily: *ω*~0~      38210.14   *ω*~0~= 0.212          two ratio *vs*. one ratio: *P*\< 0.01
                dicot subfamily 1& 2: *ω*~1~              *ω*~1~= 0.249          
  Three-ratio   monocot subfamily: *ω*~0~      38209.98   *ω*~0~= 0.212          three ratio *vs*. one ratio: *P*\< 0.01
                dicot subfamily 1: *ω*~1~                 *ω*~1~= 0.252          three ratio *vs*. two ratio: *P*\> 0.05
                dicot subfamily 2: *ω*~2~                 *ω*~2~= 0.243          
  Six-ratio     branch *A*: *ω*~1~             38182.58   *ω*~1~= 0.225          six ratio *vs*. one ratio: *P*\< 0.01
                branch *B*: *ω*~2~                        *ω*~2~= 0.373          six ratio *vs*. two ratio: *P*\< 0.01
                branch *C*: *ω*~3~                        *ω*~3~= 0.198          six ratio *vs*. three ratio: *P*\< 0.01
                branch *D*: *ω*~4~                        *ω*~4~= 0.242          
                branch *E*: *ω*~5~                        *ω*~5~= 0.187          
                other branches: *ω*~0~                    *ω*~0~= 0.283          

Positive selection in the 11S globulin genes
--------------------------------------------

To test for positive selection in the 11S globulin gene family, branches *A*, *B*, *C*, *D*and *E*were independently defined as the foreground branch in the branch-site model (Figure [2](#F2){ref-type="fig"}). When branches *A*, *B*, *C*or *D*were defined as the foreground branch, the null hypothesis is rejected (*P \<*0.01); and the estimated parameters of the alternative hypotheses indicate that about 5%-11% sites on these branches were under positive selection with a *ω*value (*ω*~2~) larger than one. When branch *E*was brought to the foreground, the null hypothesis cannot be rejected, and thus no significant positive selection was detected (Table [2](#T2){ref-type="table"}). This result suggests that positive selection mainly occurred in the 11S globulin genes in dicots.

###### 

Summary statistics for detecting selection using branch-site models of PAML

  Foreground branch   Null hypothesis   Alternative hypothesis                                             
  ------------------- ----------------- --------------------------------------------------- -------------- ---------------------------------------------------
  **Branch *A***      46859.30          *p*~0~= 0.66, *p*~1~= 0.28 (*p*~2~+*p*~3~= 0.06)    46851.6\*\*    *p*~0~= 0.69, *p*~1~= 0.27 (*p*~2~+*p*~3~= 0.04)
                                        *ω*~0~= 0.23, *ω*~1~= *ω*~2~= 1                                    *ω*~0~= 0.23, *ω*~1~= 1.00, *ω*~2~= 999
  **Branch *B***      46854.00          *p*~0~= 0.53, *p*~1~= 0.22 (*p*~2~+*p*~3~= 0.25)    46837.66\*\*   *p*~0~= 0.63, *p*~1~= 0.26 (*p*~2~+*p*~3~= 0.11)
                                        *ω*~0~= 0.23, *ω*~1~= *ω*~2~= 1                                    *ω*~0~= 0.23, *ω*~1~= 1.00, *ω*~2~= 56.08
  **Branch *C***      46852.38          *p*~0~= 0.53, *p*~1~= 0.21 (*p*~2~+*p*~3~= 0.36)    46831.76\*\*   *p*~0~= 0.64, *p*~1~= 0.26 (*p*~2~+*p*~3~= 0.10)
                                        *ω*~0~= 0.23, *ω*~1~= *ω*~2~= 1                                    *ω*~0~= 0.23, *ω*~1~= 1.00, *ω*~2~= 526.55
  **Branch *D***      46859.11          *p*~0~= 0.61, *p*~1~= 0.26 (*p*~2~+*p*~3~= 0.13)    46842.30\*\*   *p*~0~= 0.68, *p*~1~= 0.27 (*p*~2~+*p*~3~= 0.05)
                                        *ω*~0~= 0.23, *ω*~1~= *ω*~2~= 1                                    *ω*~0~= 0.23, *ω*~1~= 1.00, *ω*~2~= 999
  **Branch *E***      46859.33          *p*~0~= 0.70, *p*~1~= 0.30 (*p*~2~+*p*~3~\< 0.01)   46859.33       *p*~0~= 0.70, *p*~1~= 0.30 (*p*~2~+*p*~3~\< 0.01)
                                        *ω*~0~= 0.23, *ω*~1~= *ω*~2~= 1                                    *ω*~0~= 0.23, *ω*~1~= *ω*~2~= 1

Note: *\*\* P \<*0.01, calculated from LRT.

Discussion and Conclusions
==========================

Gene duplication may be associated with the higher 11S globulin content in dicots
---------------------------------------------------------------------------------

In our duplication analyses, we found four or more 11S globulin genes in each of the five dicot species analyzed. It appears that higher number of duplicates is a feature of the dicot 11S globulin genes, rather than being randomly produced by the hitchhiking effect of genome duplication, because:

i\) The copy number of the dicot 11S globulin gene is higher than the average copy number of the genome; for example, in the genome of *A. thaliana*, 80% of genes recovered to single copy through gene loss in a short period after the duplications, resulting an average of 2.3 copies per family \[[@B32]\], which is lower than the 4 duplicates in the 11S globulin gene family; in the genome of *G. max*, 74.1% and 56.6% genes were lost following the early and the recent WGD, respectively, resulting an average of about 3 copies per family \[[@B33],[@B34]\], which is lower than the 6 duplicates in the globulin gene family;

ii\) In each of the five dicots, there are 11S globulin genes that arose from tandem duplications; and

iii\) The genomes of *S. italica*, *Z. mays*and *S. bicolor*are thought to have undergone several rounds of duplications \[[@B35],[@B36]\], but they contain only one or two 11S globulin genes, suggesting that gene losses are common in the family, and thus implying that the duplicates of the 11S globulin genes are preferentially retained in the dicot species.

We hypothesize that the higher number of duplicates may be associated with higher 11S globulin content in dicots. The reasons are as follows.

First, the seed 11S globulin content in the species with a higher copy number of 11S globulin genes is greater than that in those with a lower copy number (Table [3](#T3){ref-type="table"}). In the dicot species *A. thaliana*, *R. communis*, *G. max*and *C. sativus*, the copy number of 11S globulin genes is four or more, and the 11S globulins are predominant among the seed storage proteins \[[@B2],[@B5],[@B20],[@B37]-[@B41]\]. In the monocot species *S. italica*, *Z. mays*and *S. bicolor*, the copy number of 11S globulin genes is one or two, and 11S globulins are minor \[[@B4],[@B31]\]. However, there are two exceptions, *O. sativa*and *B. distachyon*(monocots). In these two species, the copy number of 11S globulin genes is six or more and the 11S globulins are predominant among the seed storage proteins \[[@B6],[@B42]-[@B44]\], although the overall protein content is low.

###### 

The content of seed 11S globulins and the copy number of their genes

  Group     Species                     Seed storage protein^1^   11S globulins^2^   No. of genes
  --------- --------------------------- ------------------------- ------------------ --------------
  dicot     *Arabidopsis thaliana*      30-40%                    major component    4
            *Glycine max*               \~40%                     \~40%              6
            *Cucumis sativus*           \~35%                     major component    6
            *Populus trichocarpa*       major component           unknown            7
            *Ricinus communis*          \~40%                     \~75%              11
  monocot   *Brachypodium distachyon*   \< 10%                    \~60%              6
            *Oryza sativa*              \< 10%                    \~70%              12
            *Setaria italica*           \~10%                     minor component    2
            *Zea mays*                  \~7%                      minor component    1
            *Sorghum bicolor*           \< 10%                    minor component    1

Note: ^1^percentage to the seed dry weight; ^2^percentage to the seed storage protein. The references are listed in the main text.

Second, the presence of duplicate genes leads to the production of an extra amount of protein, because extra mRNA can be produced \[[@B26]\]. There is some evidence for this viewpoint: i) the absence of, or preferential expression of, the 11S globulin genes in *G. max*leads to glycinin deficiency or greater accumulation, respectively \[[@B45],[@B46]\]; ii) of the glycinin gene groups I (*Gy1*-*Gy3*), IIa (*Gy4*) and IIb (*Gy5*) in *G. max*, a mutation with the absence of one or two groups of genes leads to a decrease in glycinin content \[[@B47]\]; and simultaneous mutation of all the genes leads to the lack of all the glycinin polypeptides \[[@B48]\]. Third, enlarging the 2S albumin gene family of *A. thaliana*by introducing an extra member leads to an increase in transcript production \[[@B49]\].

Finally, in the investigation of ancestral whole genome duplication in seed plants and angiosperms, Jiao *et al.*\[[@B50]\] proposed that there were two WGDs in ancestral lineages shortly before the diversification of extant seed plants and extant angiosperms respectively, and argued that these ancestral WGDs i) resulted in the diversification of regulatory genes important to seed and flower development; ii) were involved in major innovations that ultimately contributed to the rise and eventual dominance of seed plants and angiosperms; and iii) enabled flowering plants to enjoy a distinct evolutionary advantage that allowed them to survive harsh climatic changes and even mass extinctions. The seed storage proteins provide essential nutrition for seed germination and development, and thus are vital for species survival and adaptation. Therefore, the genes governing the seed storage proteins are expected to have a higher number of duplications, leading to improved phenotypic robustness and an evolutionary advantage.

Higher evolutionary rate may be associated with the higher 11S globulin content in dicots
-----------------------------------------------------------------------------------------

In the evolutionary rate analyses, a consistent conclusion, that the dicot 11S globulin genes evolve more rapidly, was achieved from the three branch-specific models, i.e. two-, three- and six-ratio models (Table [1](#T1){ref-type="table"}). We hypothesize that accelerated evolutionary rate may also be associated with the higher 11S globulin content in dicots. The reasons are as follows.

Positive selection leads to an accelerated evolutionary rate of the 11S globulin genes in dicots, and may also lead to a higher ability of dicots to produce 11S globulins. Positive selection is a major factor affecting the evolutionary rate of a gene \[[@B51]\]. To investigate evidence for positive selection in the 11S globulin genes, we analyzed five branches, i.e. branches *A*-*E*(Figure [2](#F2){ref-type="fig"}), which represent the major origin events in the evolution of the 11S globulin gene families, e.g. branch *B*represents the origin of *G. max*11S globulin gene family. Of these branches, branches *A*-*D*, being the dicot 11S globulin genes, were proved to have undergone positive selection, whereas branch *E*, being the monocot 11S globulin genes, did not. This result provides evidence that positive selection may lead to a higher ability of dicots to produce 11S globulins.

To shed more light on the role of gene duplication and accelerated rates of evolution in producing the observed patterns of divergence in protein synthesis between dicots and monocots, future studies investigating other types of seed storage proteins and a broader range of plant species will be needed.

Methods
=======

Sequences Retrieval and Comparisons
-----------------------------------

Sequences retrieval and comparisons were performed using the method described by Tatusov *et al.*\[[@B29]\] with a slight modification. Briefly, our method included the following steps:

1\) Amino acid sequences of proteins were downloaded from JGI (<http://www.phytozome.net>/), the maize sequence from <http://www.maizesequence.org>, and used to construct a local BLAST database using BLAST 2.2.24. The species are listed in Table [4](#T4){ref-type="table"}.

###### 

Source of the 11S globulin genes used in this study

  Group     Common name          Species name                Version   No of genes   Reference
  --------- -------------------- --------------------------- --------- ------------- ------------
  dicot     thale cress          *Arabidopsis thaliana*      TAIR 9    33,410        \[[@B62]\]
            soybean              *Glycine max*               1.0       75,778        \[[@B34]\]
            cucumber             *Cucumis sativus*           122       32,509        \[[@B63]\]
            black cottonwood     *Populus trichocarpa*       2.0       45,778        \[[@B64]\]
            castor bean          *Ricinus communis*          0.1       31,221        \[[@B65]\]
  monocot   purple false brome   *Brachypodium distachyon*   1.0       32,255        \[[@B66]\]
            rice                 *Oryza sativa*              MSU 6.0   67,393        \[[@B67]\]
            foxtail millet       *Setaria italica*           2.1       38,038        
            maize                *Zea mays*                  5a        53,764        \[[@B68]\]
            sorghum              *Sorghum bicolor*           1.0       36,338        \[[@B69]\]

2\) An all-against-all protein sequence comparison was carried out.

3\) In all the comparisons produced in the step 2, the ones with *G. maxGy1*\[[@B16]\] as a query were identified, and the obvious paralogs were collapsed.

4\) All interspecies Best Hits (BeTs) of *Gy1*and their paralogs were detected.

5\) Steps 3) and 4) were repeated, with the resulting sequences as secondary BLASTp queries until no new sequence was found.

6\) The protocol of Tatusov *et al.*\[[@B29]\] was applied to all the sequences from the analysis above and used to form a Clusters of Orthologous Groups (COG).

7\) A case-by-case analysis of the COG was conducted. This analysis served to eliminate false-positives and to ensure all homologs were included.

This approach was based on the consistency between genome-specific best hits, rather than the absolute level of similarity; it therefore allows the detection of orthologs among both slowly and quickly evolving genes.

Phylogenetic Analyses
---------------------

The cDNA sequences were aligned using the codon model of program PRANK (100701 version) using the default options \[[@B52]\], and were then translated into amino acid sequences. Phylogenetic tree reconstruction was carried out using both Neighbor-Joining (NJ) and Bayesian approaches based on the aligned amino acid sequences. In the NJ method, the phylogenetic analyses were conducted using the MEGA 5 program \[[@B53]\]. The parameter setups were as follows: model, Jones-Taylor-Thornton (JTT) \[[@B54]\]; bootstrap, 1000 replicates; and gap/missing data, pairwise-deletion.

In the Bayesian method, the analyses were conducted using MrBayes v3.1 \[[@B55]\]. The parameter setups were as follows: JTT substitution model, four chains, one million generations, two runs, sampling every 100 generations and discarding a burn-in of 250,000 generations.

Estimation of d~N~/d~S~Ratios
-----------------------------

The coding sequences were aligned using the codon model of PRANK software (100701 version) using the default options \[[@B52]\], and alignment gaps were deleted manually.

On the basis of the aligned coding sequences, the pairwise ratio of non-synonymous substitutions per non-synonymous site (*d*~N~) to the synonymous substitutions per synonymous site (*d*~S~) (*ω*value) of homologous genes was calculated by the maximum likelihood method in Codeml from the PAML package v4.4 \[[@B56]\]. Saturation effects were avoided by discarding the gene pairs for which *d*~S~\> 2 \[[@B57]\].

The branch-specific model, which allows the *ω*ratio to vary among the branches in the phylogeny, was used to test whether there are different *ω*values on particular lineages \[[@B58]\]. If the *ω*ratio among all the branches is a constant, the model can be changed into the one-ratio model. Thus the likelihood ratio test (LRT) was used to test whether the data fit the branch-specific model significantly better than the one-ratio model \[[@B58]\].

Detection of Positive Selection
-------------------------------

The aligned codon sequences were used to test positive selection using the branch-site model implemented in the program Codeml of PAML 4.4 \[[@B56]\]. This model allows *ω*to vary both among sites in the sequences and across branches on the tree and its purpose is to detect positive selection affecting a few sites along particular lineages (called foreground branches). The model assumes that there are four site classes in the sequence. The first class of sites is highly conserved in all lineages with a small *ω*ratio, *ω*~0~. The second class includes neutral or weakly constrained sites for which *ω*= *ω*~1~, where *ω*~1~is near or smaller than 1. In the third and fourth classes, the background lineages show *ω*~0~or *ω*~1~, but foreground branches have *ω*~2~, which may be greater than 1. In the LRT, the null hypothesis fixes *ω*~2~= 1 (neutral selection) and the alternative hypothesis constrains *ω*~2~≥ 1 (positive selection) \[[@B59],[@B60]\]. In the existence of positive selection, the posterior probabilities for the sites with positive selection were calculated by the Bayes empirical Bayes method (BEB) \[[@B61]\].

Abbreviations
=============

BEB: the Bayes empirical Bayes; BeTs: Best Hits; COG: Clusters of Orthologous Groups; HMW: high molecular weight; JTT: Jones-Taylor-Thornton; LRT: likelihood ratio test; NJ: Neighbor Joining; WGD: whole genome duplication

Authors\' contributions
=======================

YMZ designed the study, coordinated and supervised the analysis. CL and ML performed the analysis and drafted the paper. YMZ and JMD revised the manuscript. All authors read and approved the final manuscript.

Acknowledgements
================

We are grateful to Dr Hugo Zheng at McGill University for help with improvements to the English text. This work was supported by the National Basic Research Program of China (2011CB109300), the National Natural Science Foundation of China (30971848), the Fundamental Research Funds for the Central Universities (KYT201002, KJ2011003), the Specialized Research Fund for the Doctoral Program of Higher Education (20100097110035), and a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions.
