Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

FoxO gene family evolution in vertebrates

FoxO gene family evolution in vertebrates Background: Forkhead box, class O (FoxO) belongs to the large family of forkhead transcription factors that are characterized by a conserved forkhead box DNA-binding domain. To date, the FoxO group has four mammalian members: FoxO1, FoxO3a, FoxO4 and FoxO6, which are orthologs of DAF16, an insulin-responsive transcription factor involved in regulating longevity of worms and flies. The degree of homology between these four members is high, especially in the forkhead domain, which contains the DNA-binding interface. Yet, mouse FoxO knockouts have revealed that each FoxO gene has its unique role in the physiological process. Whether the functional divergences are primarily due to adaptive selection pressure or relaxed selective constraint remains an open question. As such, this study aims to address the evolutionary mode of FoxO, which may lead to the functional divergence. Results: Sequence similarity searches have performed in genome and scaffold data to identify homologues of FoxO in vertebrates. Phylogenetic analysis was used to characterize the family evolutionary history by identifying two duplications early in vertebrate evolution. To determine the mode of evolution in vertebrates, we performed a rigorous statistical analysis with FoxO gene sequences, including relative rate ratio tests, branch-specific d /d ratio tests, site-specific d /d N S N S ratio tests, branch-site d /d ratio tests and clade level amino acid conservation/variation patterns N S analysis. Our results suggest that FoxO is constrained by strong purifying selection except four sites in FoxO6, which have undergone positive Darwinian selection. The functional divergence in this family is best explained by either relaxed purifying selection or positive selection. Conclusion: We present a phylogeny describing the evolutionary history of the FoxO gene family and show that the genes have evolved through duplications followed by purifying selection except for four sites in FoxO6 fixed by positive selection lie mostly within the non-conserved optimal PKB motif in the C-terminal part. Relaxed selection may play important roles in the process of functional differentiation evolved through gene duplications as well. Page 1 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 parts of the adult mouse brain. Moreover, the individual Background Mammalian FoxO proteins (FoxO1, FoxO3a, FoxO4 and disruption of Foxo1, Foxo3 and Foxo4 genes in mice FoxO6) which are homologous to Caenorhabditis elegans results in different phenotypes [14,16]. While a protein DAF-16, belong to the O ('other') class of the Fox homozygous knockout of Foxo1 (FKHR) was embryonic superfamily [1,2]. FOXO1 is the first identified member of lethal due to failures in angiogenesis and vessel forma- the FoxO family of transcription factors [3] and is tion, Foxo3a-/- (FKHRL1) and Foxo4-/- (AFX) were viable involved in the transcriptional activity of alveolar rhab- and appeared to develop normally. Later in development, domyosarcomas [3]. Since then, the discovery of mamma- Foxo3a-/- females were found to be age-dependently lian FoxO genes has grown rapidly, now FoxO proteins infertile and showed abnormal ovarian follicular develop- have been identified in several different organisms, ment. As to the physiological role, each FoxO gene exhib- including zebrafish, mouse, rat and human. As transcrip- its a distinct response under a variety of conditions [20- tion factors in the nucleus, the primary function of FoxO 22]. Therefore, despite the high sequence identity shared proteins is to bind to their cognate DNA targeting by FoxO genes domain (more then 60% in humans [17]), sequences as monomers. The co-crystal structure of HNF- the physiological roles of FoxO genes are functionally 3γ with DNA shows that there are 14 protein-DNA con- diverse in mammals. tacts distributing throughout the forkhead domain, but the third α-helix (H3) plays the most important role in a Single copy genes are thought to evolve conservatively winged helix/forkhead protein's DNA-binding specificity because of strong negative selective pressure. Gene dupli- [4]. In addition, both winged loops also make important cations produce a redundant gene copy and thus release interactions with DNA [4,5]. Although the molecular one or both copies from negative selection pressure [23]. basis of the DNA-binding specificity of FoxO transcription There are a number of models for the fate of duplicate factors is poorly understood, high-affinity DNA-binding gene that predict functional differentiation of paralogs studies have identified a consensus FoxO-recognized ele- based on protein sequence or regulatory divergence ment (FRE) as (G/C) (T/A)AA(C/T)AA [6-8]. Indeed, func- [24,25]. Currently four most prominent models are neo- tional FRE sites that match this consensus sequence have functionalization [26], subfunctionalization [24], the been identified in the promoters of many genes, such as Dykhuizen-Hartl effect [27] and adaptive diversification. Fas ligand (FasL), insulin-like growth factor binding pro- Very recently, the list has been expanded by the introduc- tein 1 (IGFBP1) and the apoptotic regulator BIM[9,10]. tion of the subneofunctionalization [28] and the adaptive Additional putative FoxO-target genes and their potential radiation [29] models that predict rapid subfunctionaliza- cis-regulatory binding sites have been predicted by sys- tion after duplication followed by a prolonged period of tematic bioinformatic approaches [11]. Thus, FoxO tran- neofunctionalization and adaptive divergence of dupli- scription factors appear to be involved in various cate genes in a process analogous to species radiations, signaling pathways and control a wide range of biochem- respectively. Thus, duplications are thought to be an ical processes including cellular differentiation, tumor important precursor of functional divergence [30]. Here, suppression, metabolism, cell-cycle arrest, cell death, and we are interested in the specific role that natural selection protection from stress [1,9,10]. In the mouse, four differ- might play in the evolutionary history of this gene dupli- ent FoxO members have been identified to date: Foxo1, cation. Foxo3, Foxo4 and Foxo6 [12,13]. FoxO6 is the latest member of the FoxO family to be cloned and shares sig- The increased availability of FoxO sequences in the public nificant sequence similarity with the other members of databases allows us to explore the functional diversity the family [13]. from a phylogenetic perspective within the FoxO family in vertebrates. The study was conducted by analyzing amino FoxO1 and its close paralogus (FoxO3, FoxO4 and acid and nucleotide-based divergence data from different FoxO6) are thought to some degree of functional diversi- species covering the entire vertebrates. Our aim was to elu- fication during development [14-16] and their potential cidate the evolutionary mechanisms operating in the physiological roles might be different [14]. Indeed, a retention of these genes and evaluate the changes in selec- rapid overview of the data collected on FoxO1, 3, 4 and 6 tion pressures following duplication. We also identified highlights how these proteins may be different. First, each the sites under positive Darwinian selection. Finally, we FoxO gene showed different expression patterns in tissues tried to map the positively selected sites to the structural [7,12,17-19]. While Foxo1 was strongly expressed in the and functional regions of FoxO molecules. striatum and neuronal subsets of the hippocampus (den- tate gyrus and the ventral/posterior part of the CA Materials and methods regions), Foxo3 was more diffusely expressed throughout Sequence Data Collection the brain including all hippocampal areas, cortex and cer- The DNA sequences and amino acids sequences of ebellum, and Foxo6 expression was eminent in various FoxO genes were downloaded from NCBI's GenBank Page 2 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 http://www.ncbi.nlm.nih.gov. PSI-BLAST searches were Relative rate tests conducted against the non-redundant database of verte- The substitution rates of the FoxO genes were compared brate genomes at NCBI (e-value cutoff = 1e-24) using the between different paralogous genes that had undergone amino acid sequences of Foxo1, Foxo3, Foxo4 and Foxo6 duplication events recently, using the RRTree software of mouse (gi: 56458, gi: 56484, gi: 54601 and gi: 329934) [41]. The orthologs FoxOs (Cifoxo, BfFoxOA and as queries. Only full length coding sequences were SpFoxOl) from amphioxus (Branchiostoma floridae), included in our analysis. Jalview 2.3 [31] was used to Ciona intestinalis and Strongylocentrotus purpuratus remove the sequences with the identity higher than 95%. were used as an outgroup. The null hypothesis is that the A table with species names, abbreviations and accession rate of substitution of the tested clade is the same as that numbers are provided in supplementary materials (Addi- of the reference group. tional file 1). Estimation of substitution rates and testing natural Sequence alignment and phylogenetic analysis selection The sequences of FoxO proteins were aligned by MUSCLE We estimated the selective pressures acting on coding [32] and the resulting alignment was manually optimized regions by applying a phylogenetic-based Maximum Like- by BioEdit [33]. Incomplete sequences, and highly diver- lihood (ML) analysis. ML estimated of the relevant param- gent regions or gaps resulting in uncertain alignments eters -as branch lengths and the ratio of the were excluded from the further analysis. The final data set nonsynonymous (d ) to synonymous substitution rates included a total of 66 sequences from 19 species. The (d ), ω = d /d -that were obtained using the codeml pro- S N S amino acid alignment was subsequently transformed into gram implemented in the PAML package version 4 [42]. an aligned cds fasta file using PAL2NAL [34] which is a The ω parameter was used as a measure of the protein program to construct multiple codon alignments from selective constraints [43]. These analyses were conducted matching amino acid sequences. The nucleotide align- under different competing evolutionary hypothesis. We ment was then converted to nexus format with DnaSP first investigated whether the distribution of selective con- [35] version 4.10 for phylogenetic analysis. straints acting on the each gene fluctuated across lineages; for that, we compared the fit to the data of the "one ratio" The full alignment of 66 sequences was used to perform model (M0), which assumes a constant selective pressure the phylogenetic analysis. Tree reconstructions were done across branches, with the "free ratios" model (FR), where by the Bayesian method from the DNA alignment done in the rate parameters were estimated independently in each the MrBayes version 3.1.2 [36,37] software package, and lineage. We also examined other evolutionary scenarios; i) rooted with the BfFoxO, Cifoxo and SpFoxO from amphi- to determine which FoxO lineage had evolved at a differ- oxus (Branchiostoma floridae), Ciona intestinalis and ent rate, as compared to the rest of the phylogeny, we Strongylocentrotus purpuratus. We analyzed four inde- applied a branch-specific model to the data. Sequences pendent runs, each using the general time reversible were divided into four groups according to their phyloge- (GTR) model plus gamma distribution plus invariant sites netic analysis, and each FoxO lineage was set as the fore- model of molecular evolution (GTR+G+I), as determined ground branch. ii) to detect sites under positive selection by Modeltest version 3.7 [38]. We ran 2 million genera- in four lineages, we applied three codon-based ML substi- tion Markov Chain Monte Carlo simulations with four tution models that are site-specific (i.e., models that allow separate chains (three heated, one cold), with the first variation in the ω ratio across sites) of [44] but assume the 500,000 generations discarded as burn-in. Trees were same selection pattern for a site in all lineages; iii) to summarized for each independent run and compared to investigate the existence of sites evolving under positive check for concordant topologies. The consensus tree of all selection in only a specific lineage, we applied the modi- compatible groupings among all runs was used in all anal- fied branch-site model A of [45] in two consecutive tests yses. (test1 and test2 in [46]) to the same alignment used for the site-based models. The model allowing for positive Synonymous codon usage analyses selection is denoted model A and the lineage to be tested Codon usage bias was estimated by the effective number is the foreground lineage, whereas the remaining ones are of codons (ENC; [39]), the frequency of optimal codons the background lineages; the multiple hypothesis testing (F ; [40]) and proportion of G and C in the third codon problem [47] was taking into account using Bonferroni's OP position (G/C 3rd). For ENC, lower values indicate correction [48]. The likelihood Ratio Test (LRT) was used stronger synonymous codon usage bias, while for F to compare the fit to the data of two nested models, OP higher values indicate stronger bias. These measures were assuming that twice the log likelihood difference between calculated for all genes using the CodonW program http:/ the two models (2ΔL) follows a χ distribution with a /bioweb.pasteur.fr/seqanal/interfaces/codonw.html and number of degrees of freedom equal to the difference in used to test whether the degree of synonymous codon the number of free parameters [49]. usage biases in individual genes. Page 3 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 We used the TreeSAAP version 3.2 [50] to determine the were identified. Inclusion of these did not improve the FoxO physicochemical properties affected by natural reliability of the phylogeny, and as the aim of this study selection. This program for examining the effects of non- was to determine the evolutionary history of the FoxO synonymous substitutions on protein evolution compares gene family, only representatives from the major verte- the observed distribution of physicochemical changes brate clades were included. inferred from a phylogenetic tree with an expected distri- Phylogenetic analyses of FoxO gene lineages bution based on the assumption of completely random amino acid replacement expected under the condition of To study the molecular evolution of vertebrate FoxO selective neutrality. For all possible pairwise amino acid genes, we carried out phylogenetic inference analyses changes, the range of effect size for each of the 31 proper- based on codon alignment and inferred their evolutionary ties was determined and equally divided into 8 magnitude history using Bayesian methods. We used the Bayesian categories. Categories 1 to 3 indicate small variation in the posterior probabilities (PPs) of each node to evaluate amino acid characteristics while categories 6 to 8 repre- clades support. Figure 1 shows the consensus phylogeny sent the most radical substitutions. For all properties that obtained for FoxO gene sequences. The vertebrate FoxOs differed significantly from neutrality, Z-scores were then were assorted well to four lineages according to their FoxO calculated in each magnitude category to determine which classification, all with high PP support values (a poorly classes contributed to this deviation. The critical Z-score supported position: 0.99 PP) indicating that the forma- values for P = 0.001 are 3.09, indicating positive selection tion of the paralogous lineages occurred before the diver- on that magnitude category, and -3.09, which indicates gence of individual species, and the orthologs FoxOs negative (purifying) selection. That is, positive and nega- (Cifoxo, BfFoxOA and SpFoxOl) from amphioxus (Bran- tive Z-scores indicate positive and purifying selection, chiostoma floridae), Ciona intestinalis and Strongylocen- respectively. Radical substitutions affecting a particular trotus purpuratus were just located as an outgroup of their property that occurred more frequently than expected by assigned lineages. From Figure 1, we inferred that two chance constituted the signature of adaptive evolution major duplications had occurred early in the vertebrate [51]. lineages. The first duplication led to the emergence of two lineages which evolved into FoxO3/6 and FoxO1/4, and Testing functional divergence and structure analysis the second duplication, also early in vertebrate evolution, To study the functional divergence and structural differ- resulted in FoxO6 and FoxO3, and FoxO1 and FoxO4. ences after the gene duplication, we used the Diverge 2.0 Phylogenetic tree shows that the FoxO6 gene cluster has software to estimate the type I (θ ) and type II (θ ) func- long branches, an indication of fast-evolving lineage with I II tional divergence coefficients [52,53] among paralogous a large number of structural changes accumulating on proteins. Type I and type II refer to shifts in the evolution- them. ary rate pattern after the emergence of a new phylogenetic cluster (indicative of changes in functional constrains), Synonymous codon usage analyses and amino acid replacements completely fixed between We investigated the relationship between nucleotide con- duplicates (resulting in cluster-specific alterations of tent and codon usage by calculating different indices amino acid physiochemical properties), respectively. (Table 1) for each of the FoxO genes. We could see from Table 1 that the effective number of codons (ENC) Genes which have been predicted to subject to positive decreased with the corresponding increase of GC3. The selection were used to search for homologous sequences effective number of codons [39] is a measure of the even- in the PDB database of protein structures http:// ness of codon usage among the 61 sense codons. An www.rcsb.org/pdb/home/home.do using Blastp [54,55]. extreme case is that all codons are used equally frequently The Rasmol http://rasmol.org/ was used for all structural (given the observed frequencies of amino acids), then the manipulations and highlighting the relevant amino acid effective number of codons is 61. Reversely, only single replacements identified in the evolutionary analyses.. codon is used for each amino acid, the effective number of codons is reduced to 20. Therefore, FoxO6 gene was Results more biased than other FoxO genes as evidenced by their Sequence similarity searches and multiple alignments lower ENC values. In most cases, the observed number fell Available FoxO1, FoxO3, FoxO4 and FoxO6 sequences somewhere between the two extremes. Figure 2 shows the were retrieved from 19 species ranging from amphioxus relationship between the effective number of codons (Branchiostoma floridae) to mammals. Additional file 1 (ENC) and the GC content at the third position of each outlines the sequences (protein and DNA) used in the gene (GC3). This Figure also contains a reference line phylogenetic analyses. The highly conserved forkhead (GCref) showing the expected position of genes whose domain remained in all alignments. It should be noted codon usage is constrained solely by the nucleotide com- that additional FoxO genes for eutherians and teleosts position at the third codon position. From Figure 2, it can Page 4 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 P Figure 1 hylogenetic relationships of DNA sequences within the FoxO family Phylogenetic relationships of DNA sequences within the FoxO family. Phylogenetic tree based on the nucleotide sequence data. The numbers indicate the Bayesian probabilities for each phylogenetic clade. Shaded boxes denote the four lin- eages and one outgroup. The scale bars represent codon substitutions per site. be seen that the observed value of ENC tracks the refer- an outgroup, we evaluated the relative rates between FoxO ence line quite closely. This indicates that the nucleotide gene clusters. The analysis (Table 2) revealed that the composition at the third codon position is a major deter- FoxO6 lineage exhibited accelerated nonsynonymous minant of the effective number of codons. substitutions with respect to FoxO3 (p-value = 0.00163, Bonferroni correction) and FoxO1 (p-value = 0.0193, Relative rates of evolution of FoxO6 lineage Bonferroni correction), and that FoxO4 genes were not Using the orthologs FoxOs (Cifoxo, BfFoxOA and accelerated with respect to the other FoxO lineages. There- SpFoxOl) from amphioxus (Branchiostoma floridae), fore, evolutionary-rate changes may have occurred follow- Ciona intestinalis and Strongylocentrotus purpuratus as ing FoxO gene duplications in the evolutionary process. Table 1: Mean values of GC%, GC3%, ENC, CAI and Fop of the FoxO genes Gene GC% GC3% ENC CAI Fop FoxO1 0.5775 ± 0.0416 0.6478 ± 0.1009 51.0256 ± 5.1246 0.0661 ± 0.0119 0.3506 ± 0.0259 FoxO3 0.5797 ± 0.0538 0.6744 ± 0.1194 47.9308 ± 6.3169 0.0765 ± 0.0171 0.3851 ± 0.0283 FoxO4 0.5974 ± 0.0430 0.6082 ± 0.0459 50.9273 ± 2.5006 0.0660 ± 0.0129 0.3522 ± 0.0385 FoxO6 0.6773 ± 0.0867 0.7695 ± 0.1486 42.4620 ± 8.7015 0.0462 ± 0.0212 0.3032 ± 0.0599 Note: Mean ± Standard deviation Page 5 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 Th Figure 2 e effective number of codons (Nc) plotted for each FoxO genes The effective number of codons (Nc) plotted for each FoxO genes. The FoxO genes highlighted in blue dot. The GC(ref) line -- shown in red -- is the expected position of genes whose codon usage is only determined by the GC content at the third positions of codons (GC3s). Selective constraints and functional divergence pairs. As expected, most amino acids had very low poste- Gene duplication-specific changes in the substitution rior probability (PP) values and, therefore, they would not rates (type I functional divergence) might reflect the dif- be involved in the hypothetical functional divergence ference in evolutionary rate at amino acid sites after gene (Figure 3). Specifically, we detected 32 and 15 amino acid duplication [52,53]. We found significant evidence of positions which presumably submitted to altered func- type I functional divergence for comparisons between dif- tional constraints when the PP threshold values were set ferent gene clusters (θ = 0.23 ~0.40, P < 0.01; Table 3); to 0.87 and 0.95, respectively. Type I sites are defined as namely, there were some amino acid sites with discrepan- those with an amino acid that is conserved in one cluster cies in their evolutionary rate between these paralogous but variable in the sister cluster, implying that the site is Page 6 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 Table 2: Evolutionary Rate of the FoxO Gene Families Lineage1 Lineage2 Ka1 Ka2 dKa sd_dKa ratio_Ka P_Ka FoxO6 FoxO4 1.03324 0.997808 0.035436 0.046871 0.756029 0.449673 FoxO6 FoxO3 1.03379 0.902464 0.131326 0.041694 3.14977 0.001638 FoxO6 FoxO1 1.05943 0.952932 0.106501 0.04553 2.33915 0.019336 FoxO4 FoxO3 0.984938 0.921849 0.063089 0.044978 1.40267 0.16075 FoxO4 FoxO1 0.997224 0.942435 0.054789 0.048375 1.13257 0.257401 FoxO3 FoxO1 0.94013 0.949027 -0.0089 0.040392 -0.22028 0.825657 Note: Ka corresponds to the mean evolutionary rate measured as the number of nonsynonymous substitutions per site. dKa is the mean rate difference between the two lineages. sd_dKa is the standard deviation and ratio_Ka the ratio between dKa and sd_dKa. The P_Ka column corresponds to the P value associated to the test under structural/functional constraints in the first cluster was substantially smaller than 1 (one ratio model ω = that is absent in the variable cluster [56]. 0.084, Table 4) that indicated that purifying selection had been the predominant force acting on the evolution of Recently, a method has been developed to test for type II these vertebrate FoxOs. Omega estimates for FoxO1, functional divergence [57]. Type II sites are those that are FoxO3, FoxO4 and FoxO6 phylogenies were 0.09583, highly conserved in both clusters but are fixed for amino 0.08311, 0.14088 and 0.13464, respectively. Selective acids with different biochemical properties between sister constrains, however, are unevenly distributed across the clusters, implying these residues are responsible for the phylogeny (FR model; 2ΔL = 421.20, P < 0.001). We then functional differences between these groups. Although at ran the branch model using each FoxO lineage as the fore- least one site with evidence of type II divergence was ground branch. In this model the estimated ω was 0.0758 found for comparisons between FoxO1/FoxO3, FoxO3/ for the FoxO1, and 0.0898 for the background branches. FoxO4, and FoxO1/FoxO4 clusters, the θ values are A LR test indicated that the two-ratio model was not sig- II extremely small (θ = 0.005 ~0.074) that highlighted the nificantly different from the M0 model (2ΔL = 3.05, P > II conservation between different clusters. These results are 0.05, df = 1, Table 4). In contrast to FoxO1 analysis, the ω not unexpected given that this method calculates θ across values of the FoxO3, FoxO4 and FoxO6 lineages were dif- all sites in an alignment and thus effectively averages site- ferent from the rest of the phylogeny as the LR tests indi- wise θ values. With only ~3% of sites/cluster showing a cated that the two-ratio model fit the data better than the pattern of type II divergence in our concatenated align- M0 model for these three genes (P < 0.05). Unfortunately, ment, it is not likely that the ~9 possible type II sites have the ω estimates for FoxO3, FoxO4 and FoxO6 were not θ values high enough to compensate for the extremely indicative of positive selection, they were rather indicative II low θ values of the over 300 sites with θ effectively equal of relaxed constraint. II II to zero. Our results are similar to the analysis of Hox-gene [30]. Along with lineage heterogeneity, variations in ω across sites can also occur. Theoretically, different protein The analysis of the nonsynonymous to synonymous sub- regions with different functions may experience different stitution rate ratio can also be used to detect functional selection pressures, which can be tested by fitting the data differentiation. We estimated ω as an average over all sites to a model comprising different site classes. The results and branches from the FoxO paralogus MSA and the ratio were shown in Table 5, for each lineage, the M3 vs M0 LRT was significant, indicating that one category of ω wasn't fit Table 3: Maximum likelihood estimates of the coefficient of data well to describe the variability in selection pressure functional divergence (θ) from pairwise comparisons between across amino acid sites. The tests contrasting the models FoxO groups M1a against M2a resulted in the P value of 1 for all the a b c Comparison θ SE (θ)LRT (θ)sig. groups suggesting a lack of power and the amino acid changes within each cluster were neutral or under negative FoxO1 Vs FoxO3 0.33 0.05 50.89 P < 0.01 selection. M1a, the parameter estimates for the least FoxO1 Vs FoxO4 0.29 0.04 41.23 P < 0.01 parameter rich model describes that most sites with low ω FoxO1 Vs FoxO6 0.3 0.06 20.43 P < 0.01 estimates (indicative of strong selective constraints), that FoxO3 Vs FoxO4 0.23 0.05 19.86 P < 0.01 FoxO3 Vs FoxO6 0.4 0.05 57.92 P < 0.01 is, 82% of FoxO1 sites were under strong purifying selec- FoxO4 Vs FoxO6 0.24 0.05 21.61 P < 0.01 tion, compared to 83% for FoxO3, 74% for FoxO4 and 66% for FoxO6. The test using M7 and M8, which allows Note: θ is the coefficient of functional divergence; for beta-distributed site-specific ω ratio, detected 2 groups SE(θ) standard error; under possible positive selection at 0.05 significance level, LRT(θ) is a likelihood ratio test; Page 7 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 Typ Figure 3 e I functional divergence among the FoxO members Type I functional divergence among the FoxO members. Posterior probability (PP) profiles of the site-specific type I functional divergence. The positions with gaps involved in each paralogous comparison were not considered. Red line indicates cutoff = 0.95, while green cutoff = 0.87. one with ω = 1.36 and the other with ω = 127.02 (Table method has been developed that allows variation in ω 5). In order to test whether the estimated ω is significantly across individual codons on a specific lineage [46,58]. greater than 1, model M8 was compared with a more This model (MA) designates two categories of branches, restricted null model (M8a). For FoxO3 gene, Model M8 again foreground and background, where positive selec- did not significantly differ from model M8a (2ΔL = 0.067, tion is modeled only on the foreground branch. We then P > 0.05, df = 1). For FoxO6 gene, ω = 127.02 was signifi- applied the branch-site approach (using some pre-speci- cantly different than 1 (2ΔL = 51.92, P < 0.01, df = 1). We fied branches, i.e., foreground branches), designating also used the BEB estimation method in model M8 [46] to each FoxO gene as the foreground branch, to assess identify sites under possible positive selection. whether molecular adaptation occurred in the evolution of the FoxO genes. The results of this analysis exhibited Since positive selection will likely affect a few amino acids several positions with evidence of relaxed selection (the at specific lineages on the phylogeny, models estimating ω test 1 was significant) (Table 6). However, we could not ratios averaged by codons or by lineages are certainly reject the null hypothesis of the test 2 (ω = 1) (result not highly conservative. For this reason, a branch-sites Page 8 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 Table 4: LRTs done to detect heterogeneous selection regimes among lineages for each gene model df Parameter estimates lnL 2l p value Branch-specific models -22943.8 M0(one-ratio) ω = 0.08442 FoxO1 two-ratio vs one-ratio 1 ω = 0.0898ω = 0.0758 -22942.3 3.047108 p > 0.05 0 1 FoxO3 two-ratio vs one-ratio 1 ω = 0.0910ω = 0.0730 -22941.3 4.87108 p < 0.05 0 1 FoxO4 two-ratio vs one-ratio 1 ω = 0.0811ω = 0.1044 -22941.8 3.983174 P < 0.05 0 1 FoxO6 = 0.0786ω = 0.1358 -22935.8 16.00121 P < 0.01 two-ratio vs one-ratio 1 ω 0 1 shown); thus, these analyses do not provide any evidence expected by chance for 2 of the properties (alpha-helical for directional selection on the FoxO lineages. tendencies and compressibility). The molecular adaptation processes occurred after the Spatial distributions of possible selected FoxO6 Sites on three-dimensional structure gene duplication event were also investigated by compar- ing the magnitude of the physicochemical changes pro- Because of the evidence for possible positive selection on duced by the observed amino acid replacements with FoxO6, we predicted positively selected codon sites using those expected at random [59]. We used the program a Bayes empirical Bayes (BEB) method [45]. The sites Tree-SAAP [50] to model how 31 different physicochemi- under selection in FoxO6 are listed in Table 5. Four codon cal properties were affected by amino acid substitutions in sites were identified as positively selected at a BEB poste- each FoxO gene. Consistent overrepresentation of radical rior probability threshold of 95%. Moreover, 7 amino amino acid changes (i.e., categories 7 and 8) would sug- acid residues presumably submitted to altered functional gest repeated adaptive substitution [51]. The results indi- constraints were identified by both PAML 4 and Diverge cated that, some amino acid replacements altering these 2.0 analysis (Table 6). In order to plot positive selected physicochemical properties in the FoxO1 and FoxO3 pro- sites onto mouse (Foxo6) three-dimensional model, we teins accumulated more (or less) often than expected by first built an energy-minimized model using a homology chance (likely reflecting fitness differences) (supplemen- modeling approach [60]. The PDB entry with the highest tary materials (Additional file 1)). Moreover, for each sequence similarity -identified in the PSI-BLAST- corre- physicochemical property, the distribution of the Z-scores sponds to the human FOXO3A (PDB: 2k86). We used this across 8 magnitude classes [51] indicated that, amino acid entry as a template for the modelling. The in silico stereo- substitutions occurred less often than expected by chance chemical quality analysis [61] indicated that the generated at the most extreme magnitude-classes (supplementary model had a moderate quality (with the percentage of res- Additional file 1); these FoxO1 and FoxO3 protein prop- idues in most favored regions being no lower than the erties, therefore, were likely evolving under purifying 82.8%), with only 1.1% in disallowed regions. As selection. For FoxO4 and FoxO6 genes, less physicochem- expected, the modeled structure was roughly similar to the ical properties were affected by amino acid substitutions. template, with the three helices and two wing loops typi- The FoxO6 gene, on the contrary, seems to evolve positive cal of the Fox family in equivalent positions and with a selection, because category 8 occurs more frequently than similar predicted folding (Figure 4). Taken together, these Table 5: Site model analyses for the FoxO1, FoxO3, FoxO4 and FoxO6 phylogenies Models comparison M3 vs M0 M2a vs M1a M8 vs M7 Gene 2ΔL = (L1-L0) p-value 2ΔL = (L1-L0) p-value 2ΔL = (L1-L0) p-value ω- value Positively selected sites FoxO1 580.7848 p < 0.01 0 1 0 1 FoxO3 902.91 p < 0.01 0 1 8.23 p < 0.05 ω = 1.36 66 L (p > 0.90) FoxO4 178.24 p < 0.01 0 1 0 1 FoxO6 455.31 p < 0.01 0 1 69.66 p < 0.01 ω = 127.02 264K* 266P* 434G* 439T* Note: *P > 0.95 Page 9 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 Table 6: Parameter estimations and likelihood ratio tests for the branch-site models a b c df Parameter estimates lnL 2 P value Positive selected sites FoxO1 MA Vs M1a (test 1) 2 p = 0.70592 p = 0.18493 (p = 0.10915) -22667.231 81.89419 P < 0.01 213S** 216S* 219S* 252M* 276V** 0 1 2 w = 0.07205 (w = 1.00000) w = 285P* 296L** 306A** 340F* 360E* 0 1 2 1.00000 FoxO3 MA Vs M1a (test 1) 2 p0 = 0.77299 p1 = 0.10453 (p2 = -22646.766 122.8257 P < 0.01 5H** 25D* 26F** 33D** 34L**37N** 0.12248) w0 = 0.06958 (w1 = 1.00000) 217A* 231G** 329G* w2 = 1.00000 FoxO4 MA Vs M1a (test 1) 2 p0 = 0.67705 p1 = 0.18055 (p2 = 0.1424) -22660.108 96.14052 P < 0.01 7V** 173R** 194T** 201I** 202L** w0 = 0.07468 (w1 = 1.00000) w2 = 211F**223H* 225P** 242T* 254R** 1.00000 314S** FoxO6 MA Vs M1a (test 1) 2 p0 = 0.58516 p1 = 0.13030 (p2 = -22527.632 361.0936 P < 0.01 42Q** 46K** 155I* 164T** 165N** 0.28454) w0 = 0.06572 (w1 = 1.00000) 173R* 174E** 176E** 178L** 179F** w2 = 1.05311 180C** 188I* 189V** 203L* 207R* 223H** 230I** 231G* 232Y** 233K** 234N** 237Y** 258S** 265N* 269T** 271E** 272N** 273E** 274V** 275H** 276V** 277S** 278Q* 279G** 280L** 281H** 282P** 283S** 286N* 314S** 316V** 320H** 330Y* 366T** 367G** 368T** 369P* Note: The number of free parameters; Likelihood of the model; 2(l -l ); 1 0 * P > 0.95; ** P > 0.99; The bold amino acid residues were also found to be implicated in the functional divergence (implemented in Diverge 2.0) between FoxO data suggested that the model was stereochemically valid, performed firstly to the resolution of the evolutionary and therefore suitable for further sequence-structural relationships of these FoxOs using molecular sequence analysis. Unfortunately, we could not map any positive data. Whereas, the incorrect phylogenetic topology result- selected sites onto the surface of the 3D structure (Figure ing from mutationally saturated positions, inadequate 4A) because the crystal structures about Fox proteins are modeling of the evolutionary process and systematic bias mainly focused on the forkhead DNA-binding domain. due to variable rates of evolution among species or within Whereas the positive selected sites were mainly located in sequences [62] may make LRT generate many false posi- the region of N-terminal and C-terminal of FoxO6, which tives. Anisimova et al (2003) examined the effect of also indicated that FoxO6 underwent strong constraint on assuming a "wrong" tree [63], and he found that LRT the forkhead domain as well (Sequence logo of the fork- falsely suggested positive selection in 96% of the repli- head domain, Figure 4). cates in the M0-M3 comparison and in 86% of the repli- cates in the M7-M8 comparison at the = 5% significance Discussion level. In order to overcome this problem, we adopted a It has long been know that FoxO transcription factors play number of ways in combination. Firstly, the addition of important roles in regulating various signals, which trans- more taxa to the dataset: denser sampling of species can late various environmental stimuli into dynamic gene reduce the effect of long branch attraction (LBA) by reduc- expression programs to influence many physiological and ing the overall distances between taxa. Secondly, we used pathological processes, including cancer and aging. The the best model of DNA substitution, determined by Mod- functions of FoxO proteins are regulated at multiple lev- eltest version 3.7 [38]. And finally, our inclusion of els, which include but are not limited to phosphorylation, enough sequences in each lineage helped alleviate loss of ubiquitylation and acetylation. Interestingly, all of these LRT power from short conserved sequences. From phylo- activities affect nuclear/cytoplasmic trafficking of FoxO genetic result, we focused on the 2 main duplications proteins. The specific function of each member of this along the evolutionary history of FoxO genes, the FoxO1- family is different [14]. As the accumulation of gene FoxO4 and the FoxO3-FoxO6 duplication, which formed sequences in the database, it is feasible to explore the four gene lineages (all with the high confidence values, functional diversity from a phylogenetic perspective. We Page 10 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 The Figure 4 modeled structure of mouse Foxo6 The modeled structure of mouse Foxo6. A. The structure of the forkhead domain; B. Sequence logo of the forkhead domain and surrounding amino acids. nearly 100% posterior probability in Bayesian analysis) It is widely accepted that gene duplication can create and used for further analysis. opportunities for functional divergence in paralogues. Divergence is thought to occur where one duplicate Codon bias is largely thought to be due to weak selection retains the original protein function and the other accu- acting to optimize protein production [64-66]. Selection mulates changes, (either through redundancy or by posi- intensity for codon usage bias, therefore, is expected to tive selection) or alternatively, through the partitioning of vary among genes. Our survey of synonymous codon the functions of an unduplicated ancestor protein. What- usage in FoxO genes revealed a strong and consistent pat- ever the mechanism, if functional divergence has occurred tern of codon bias in genes with FoxO6 relative to those between duplicated genes, then it should be observable as with FoxO1, FoxO3 and FoxO4 (Table 1). At the same changes within their coding regions. time, there appears to be some conflicting results observed between F and ENC, which may be caused by The functional divergence of FoxO genes has been studied OP differences in the way that the two methods estimate by [14]. The branch length leading to the FoxO6 clade is codon bias. F is based on the frequency of a set of spe- extended relative to other FoxO genes in gene phylogeny, OP cies specific "optimal" codons, while ENC is based on the (Figure 1). This suggested that after the duplications, observed number of codons used for each amino acid. FoxO6 evolved at a faster rate than other FoxO genes. This Thus it is possible for the two methods to give different result was confirmed by significant relative rate test results estimates of codon bias. for FoxO gene lineages (Table 2). In this sense, we per- Page 11 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 formed type I functional divergence analysis, and we pared with a more restricted null model (M8a), where the detected significant type I divergence among FoxOs. The extra site class was constrained to have ω = 1. When per- comparison between the FoxO3 and FoxO6 groups forming this analysis with the sequence data from the showed the highest value for θ (0.40 ± 0.05), suggesting FoxO3 gene, model M8 did not significantly differ from that these two groups had diverged considerably more at model M8a (2ΔL = 0.067, df = 1, P > 0.05), indicating that the functional level. Next, DIVERGE was used to establish the estimated ω = 1.36 was not significantly different than the posterior probability of type I divergence at each site 1 and that there was little indication of positive selection in the alignment, employing two cut-off posterior proba- in this gene. Further, our observation of strong purifying bility values of 0.87 and 0.95. However, the cutoff value selection being the primary mode of evolution through- for residue selection is an empirical decision and is out the FoxO phylogeny is consistent with the findings of expected to depend on the intrinsic properties of the pro- a recent study about forkhead family [67,68]. tein family being analyzed. Thus, we predicted 32 candi- date functional divergence-related sites using 0.87 as a When we performed branch-site model analysis, we cutoff value (supplementary materials (Additional file found relaxed functional constraint was most consistent 1)). When we narrowed our criteria to 0.95, we got 15 can- with the molecular evolutionary analyses of the FoxO didate residues as the most likely candidate sites for type I data. Our conclusion is in contrast to a previous study functional divergence (supplementary materials (Addi- which concluded that one site was found to be under pos- tional file 1)), but we lacked a way to verify how the rate- itive selection in the FoxO3 lineage [67]. In our paper, we shift in these sites contributed to functional divergence applied the modified branch-site model A of [45] in two among the FoxO gene groups. For comparative purposes, consecutive tests (test1 and test2 in [46]) to the same the same alignment and phylogeny was submitted to a ML alignment used for the site-based models, but we could LRT, which, like the Bayesian method provided a statisti- not detected the positive selection site. To determine why cal framework where evolutionary rate shifts at particular these two articles are giving drastically different results, we protein positions could be established [46]. At last, the had a look at the sequences used for branch-site model statistically most likely positions predicted to underlie analysis in [67], and we found that only 12 sequences (4 functional divergence were agreement by both methods, for FoxO3) used for testing evolution selection. Test 2 is a particularly for the highest-ranking candidates (Table 6). more direct test for identifying positive selections in the foreground branch, and significant LRT from test 1 can be In this study we used codon substitution models imple- resulted from either positive selection or relaxed selective mented through a maximum likelihood framework to constraint in the foreground branch [46], however, the estimate the rate of evolution at silent and replacement power of this method may be limited when sample size or sites in FoxO1 and its paralogs, FoxO3, FoxO4 and divergence time is low [46]. Therefore, we concluded that FoxO6. Different models were used to investigate varia- the contradiction between our results and the previous tion in the rate of evolution between lineages of a phylog- study [67] due to the number of sequences used for anal- eny, and to estimate ω for specific lineages and sites across ysis. Moreover, our work on site-model analysis, relative phylogenies. Our objective was to determine the mode of rate test and physicochemical changes indicated that evolution on each FoxO gene lineage, and to determine FoxO6 was under positive selection. Four positive selected whether increased positive selection or decreased con- sites were identified by site-model analysis, two (264K, 337 339 and Pro ) of them straint led to the functional divergence of FoxO genes. As 266P, corresponding to mouse Gly we have demonstrated, variation among branch and sites fell into the region of non-conserved optimal PKB motif was observed in the FoxO6 phylogeny. Moreover, physic- in the C-terminal part (Thr ) [69]. The C-terminal PKB ochemical amino acid properties analysis also provided recognition sequence is not conserved in FoxO6 [13]. evidence that the entire FoxO6 gene had experienced Besides a PKB phosphorylation motif, this region contains repeated episodes of adaptive evolution. The site models a stretch of 3 additional serine residues, present in the showed that adaptation had appeared at four sites located other members of the FoxO group, but FoxO6. All these at C-terminal of FoxO6. For FoxO3 gene, Model M8 fits suggest that these serines may be functionally important data significantly better than Model M7 (2ΔL = 8.23, df = in the other sequences analyzed with the exception of 2, P < 0.05), and because 1.5% of sites are located in the FoxO6 gene, and positive selection may lead to functional positively selected site class with ω = 1.36, weak positive divergence between FoxO6 and the other members too. 545 550 selection may be indicated with this comparison. How- Another two positive selective site are Gly and Pro in ever, it has been found that a poor fit of the data to a beta the mouse Foxo6, and the functional role remains elusive. distribution may result in a high frequency of significant That is to say that the real reason for their accelerated evo- tests when comparing models M7 and M8 even in the lution is unclear. However, it should be mentioned that absence of positive selection. To take account for the ele- there is a gap in the knowledge of the relationship vated type I error rates, the original model M8 was com- between amino acid sequence and structure for full-length Page 12 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 FoxO sequence, and we are unable to speculate on the Additional material particular role of this region in these FoxO6 genes. Unfor- tunately, the shared evolutionary history and molecular Additional file 1 selection alone cannot be used as the unique criterion to Excel spreadsheet containing: A list of species, species abbreviations, and infer protein function, and the true nature of each FoxO accession numbers for sequences used in the study/A list of statistically sig- gene needs to be determined experimentally and inde- nificant physicochemical amino acid properties for each FoxO gene/A list pendently. Therefore, the positively selected site may play of the candidate residues as the most likely candidate sites for type I func- tional divergence. an important functional role and could represent an inter- Click here for file esting target site for future mutagenesis experiment thus [http://www.biomedcentral.com/content/supplementary/1471- facilitating our understanding of the structure-function 2148-9-222-S1.xls] relationships in FoxO genes. Molecular testing is required to validate this hypothesis. The result from branch-site analysis (relaxation of functional constraints) of FoxO6 also differs somewhat from previously signature of posi- Acknowledgements tive selection. We infer that the weak positive selection The authors wish to express their gratitude to the members of animal sci- and multiple branches are considered as foreground ences laboratory of Shanghai Jiao Tong University. The authors also thank branch may explain this phenomenon, because power editor for his suggestions about the manuscript and Jing Li for her help in will be reduced unless the same sites and selective con- revising the manuscript. This work is supported by the National High Tech- nology Research and Development Program of China (863 project) (grant straints are occurring along all foreground branches [46]. no. 2006AA10Z1E3, 2008AA101002), the National 973 Key Basic Research Program (grant no, 2006CB102102, 2004CB117500) and the National Nat- Conclusion ural Science Foundation of China (grant no. 30671492, 30871782). Genomic data have provided an opportunity to gain a bet- ter understanding about the evolution of FoxOs using References phylogenetic analyses. The FoxO gene family phylogeny 1. Barthel A, Schmoll D, Unterman TG: FoxO proteins in insulin showed that two duplications took place early in the evo- action and metabolism. Trends Endocrinol Metab 2005, 16(4):183-189. lution of vertebrates and triggered diversification of the 2. Kaestner KH, Knochel W, Martinez DE: Unified nomenclature for FoxO gene family into four groups. However, further the winged helix/forkhead transcription factors. Genes Dev 2000, 14(2):142-146. genome projects on a greater diversity of evolutionary lin- 3. Fredericks WJ, Galili N, Mukhopadhyay S, Rovera G, Bennicelli J, Barr eages would help better understand the gene-duplication FG, Rauscher FJ: The PAX3-FKHR fusion protein created by history. The relative rate analysis and physicochemical the t (2; 13) translocation in alveolar rhabdomyosarcomas is a more potent transcriptional activator than PAX3. Mol Cell changes indicated that FoxO6 seemed to be different from Biol 1995, 15(3):1522-1535. other members. Evolutionary rate analysis showed that 4. Clark KL, Halay ED, Lai E, Burley SK: Co-crystal structure of the molecular adaptation can also play an important role in HNF-3/fork head DNA-recognition motif resembles histone H5. Nature 1993, 364(6436):412-420. the evolution of this gene family. Indeed, positive selec- 5. Boura E, Silhan J, Herman P, Vecer J, Sulc M, Teisinger J, Obsilova V, tion was likely involved in the functional differentiation Obsil T: Both the N-terminal Loop and Wing W2 of the Fork- head Domain of Transcription Factor Foxo4 Are Important of FoxO6 gene; likewise, relaxed selection might play for DNA Binding. J Biol Chem 2007, 282(11):8265-8275. important roles over evolutionary time and shape varia- 6. Biggs Iii WH, Meisenhelder J, Hunter T, Cavenee WK, Arden KC: tion of some members of the family. Considering the evo- Protein kinase B/Akt-mediated phosphorylation promotes nuclear exclusion of the winged helix transcription factor lutionary history of the FoxO gene family, we provided FKHR1. Proc Natl Acad Sci USA 1999, 96(13):7421-7426. insight into which amino acid residues might have under- 7. Furuyama T, Nakazawa T, Nakano I, Mori N: Identification of the differential distribution patterns of mRNAs and consensus gone positive selection and could be targeted for site- binding sequences for mouse DAF-16 homologues. Biochem J directed mutagenesis. However, the identification of four 2000, 349(Pt 2):629-634. sites under positive selection requires supporting evidence 8. Gilley J, Coffer PJ, Ham J: FOXO transcription factors directly activate bim gene expression and promote apoptosis in sym- from further functional experiments to demonstrate the pathetic neurons. J Cell Biol 2003, 162(4):613-622. adaptive character of the amino acids. All these studies 9. Accili D, Arden KC: FoxOs at the Crossroads of Cellular and experiments will certainly contribute to better under- Metabolism, Differentiation, and Transformation. Cell 2004, 117(4):421-426. stand the precise role of natural selection and functional 10. Greer EL, Brunet A: FOXO transcription factors at the inter- divergence of this family. face between longevity and tumor suppression. Oncogene 2005, 24:7410-7425. 11. Xuan Z, Zhang MQ: From worm to human: bioinformatics Authors' contributions approaches to identify FOXO target genes. Mech Ageing Dev MH, Wang conceived and supervised all research. QS, 2005, 126(1):209-215. 12. Biggs Iii WH, Cavenee Karen CWK: Identification and character- Wang helped with the Type I analysis. XZ, Zhang and HB, ization of members of the FKHR (FOX O) subclass of Zhao gave the suggestion on the analysis of physicochem- winged-helix transcription factors in the mouse. Mamm Genome 2001, 12(6):416-425. ical properties. YC, Pan gave the suggestion on discussion 13. Jacobs FM, Heide LP van der, Wijchers PJ, Burbach JP, Hoekman MF, and approved the final version. Smidt MP: FoxO6, a novel member of the FoxO class of tran- Page 13 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 scription factors with distinct shuttling dynamics. J Biol Chem 39. Wright F: The'effective number of codons' used in a gene. 2003, 278(38):35959-35967. Gene 1990, 87(1):23-29. 14. Hosaka T, Biggs WH, Tieu D, Boyer AD, Varki NM, Cavenee WK, 40. Ikemura T: Correlation between the abundance of Escherichia Arden KC: Disruption of forkhead transcription factor coli transfer RNAs and the occurrence of the respective (FOXO) family members in mice reveals their functional codons in its protein genes: a proposal for a synonymous diversification. Proc Natl Acad Sci USA 2004, 101(9):2975-2980. codon choice that is optimal for the E. coli translational sys- 15. Arden KC: FoxOs in tumor suppression and stem cell mainte- tem. J Mol Biol 1981, 151(3):389-409. nance. Cell 2007, 128(2):235-237. 41. Robinson-Rechavi M, Huchon D: RRTree: Relative-Rate Tests 16. Castrillon DH, Miao L, Kollipara R, Horner JW, DePinho RA: Sup- between groups of sequences on a phylogenetic tree. Bioinfor- pression of ovarian follicle activation in mice by the tran- matics 2000, 16(3):296-297. scription factor Foxo3a. Science 2003, 301(5630):215-218. 42. Yang Z: PAML 4: phylogenetic analysis by maximum likeli- 17. Anderson MJ, Viars CS, Czekay S, Cavenee WK, Arden KC: Cloning hood. Mol Biol Evol 2007, 24(8):1586-1591. and Characterization of Three Human Forkhead Genes 43. Yang Z: Inference of selection from multiple species align- That Comprise an FKHR-like Gene Subfamily. Genomics 1998, ments. Curr Opin Genet Dev 2002, 12(6):688-694. 47(2):187-199. 44. Yang Z, Nielsen R, Goldman N, Pedersen A: Codon-substitution 18. Kitamura T, Nakae J, Kitamura Y, Kido Y, Biggs WH 3rd, Wright CV, models for heterogeneous selection pressure at amino acid White MF, Arden KC, Accili D: The forkhead transcription fac- sites. Genetics 2000, 155(1):431-449. tor Foxo1 links insulin signaling to Pdx1 regulation of pan- 45. Yang Z, Wong W, Nielsen R: Bayes empirical bayes inference of creatic beta cell growth. J Clin Invest 2002, 110(12):1839-1847. amino acid sites under positive selection. Mol Biol Evol 2005, 19. Hoekman MF, Jacobs FM, Smidt MP, Burbach JP: Spatial and tem- 22(4):1107-1118. poral expression of FoxO transcription factors in the devel- 46. Zhang J, Nielsen R, Yang Z: Evaluation of an improved branch- oping and adult murine brain. Gene Expr Patterns 2006, site likelihood method for detecting positive selection at the 6(2):134-140. molecular level. Mol Biol Evol 2005, 22(12):2472-2479. 20. Nakae J, Kitamura T, Kitamura Y, Biggs WH 3rd, Arden KC, Accili D: 47. Anisimova M, Yang Z: Multiple hypothesis testing to detect lin- The forkhead transcription factor Foxo1 regulates adipocyte eages under positive selection that affects only a few sites. differentiation. Dev Cell 2003, 4(1):119-129. Mol Biol Evol 2007, 25(4):1219-1228. 21. Richards JS, Sharma SC, Falender AE, Lo YH: Expression of FKHR, 48. Miller R: Simultaneous statistical inference. Springer-Verlag, FKHRL1, and AFX genes in the rodent ovary: evidence for New York; 1981. regulation by IGF-I, estrogen, and the gonadotropins. Mol 49. Whelan S, Goldman N: Distributions of statistics used for the Endocrinol 2002, 16(3):580-599. comparison of models of sequence evolution in phylogenet- 22. Bois PRJ, Grosveld GC: FKHR (FOXO1a) is required for myo- ics. Mol Biol Evol 1999, 16(9):1292-1299. tube fusion of primary mouse myoblasts. EMBO J 2003, 50. Woolley S, Johnson J, Smith M, Crandall K, McClellan D: TreeSAAP: 22:1147-1157. selection on amino acid properties using phylogenetic trees. 23. Hughes J, Criscuolo F: Evolutionary history of the UCP gene Bioinformatics 2003, 19(5):671-672. family: gene duplication and selection. BMC Evol Biol 2008, 51. McClellan DA, Palfreyman EJ, Smith MJ, Moss JL, Christensen RG, 8:306. Sailsbery JK: Physicochemical evolution and molecular adapta- 24. Force A, Lynch M, Pickett FB, Amores A, Yan YL, Postlethwait J: tion of the cetacean and artiodactyl cytochrome b proteins. Preservation of duplicate genes by complementary, degen- Mol Biol Evol 2005, 22(3):437-455. erative mutations. Genetics 1999, 151(4):1531-1545. 52. Gu X: Statistical methods for testing functional divergence 25. Kimura M: The neutral theory of molecular evolution. Sci Am after gene duplication. Mol Biol Evol 1999, 16(12):1664-1674. 1979, 241(5):98-100. 53. Gu X: Maximum-likelihood approach for gene family evolu- 26. Goodman M, Moore GW, Matsuda G: Darwinian evolution in the tion under functional divergence. Mol Biol Evol 2001, genealogy of haemoglobin. Nature 1975, 253(5493):603-608. 18(4):453-464. 27. Dykhuizen D, Hartl DL: Selective neutrality of 6PGD allozymes 54. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic local in E. coli and the effects of genetic background. Genetics 1980, alignment search tool. J Mol Biol 1990, 215(3):403-410. 96(4):801-817. 55. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lip- 28. He X, Zhang J: Rapid subfunctionalization accompanied by man DJ: Gapped BLAST and PSI-BLAST: a new generation of prolonged and substantial neofunctionalization in duplicate protein database search programs. Nucleic Acids Res 1997, gene evolution. Genetics 2005, 169(2):1157-1164. 25(17):3389-3402. 29. Francino MP: An adaptive radiation model for the origin of 56. Gu X: Functional divergence in protein (family) sequence new gene functions. Nat Genet 2005, 37(6):573-577. evolution. Genetica 2003, 118(2-3):2-3. 30. Lynch VJ, Roth JJ, Wagner GP: Adaptive evolution of Hox-gene 57. Gu X: A simple statistical method for estimating type-II (clus- homeodomains after cluster duplications. BMC Evol Biol 2006, ter-specific) functional divergence of protein sequences. Mol 6:86. Biol Evol 2006, 23(10):1937-1945. 31. Clamp M, Cuff J, Searle SM, Barton GJ: The Jalview Java alignment 58. Yang Z, Nielsen R: Codon-Substitution Models for Detecting editor. Bioinformatics 2004, 20(3):426-427. Molecular Adaptation at Individual Sites Along Specific Lin- 32. Edgar RC: MUSCLE: a multiple sequence alignment method eages. Mol Biol Evol 2002, 19(6):908-917. with reduced time and space complexity. BMC Bioinformatics 59. McClellan DA, McCracken KG: Estimating the influence of selec- 2004, 5:113. tion on the variable amino acid sites of the cytochrome B 33. Hall TA: BioEdit: a user-friendly biological sequence align- protein functional domains. Mol Biol Evol 2001, 18(6):917-925. ment editor and analysis program for Windows 95/98/NT. 60. Schwede T, Kopp J, Guex N, Peitsch MC: SWISS-MODEL: An 1999:95-98. automated protein homology-modeling server. Nucleic Acids 34. Suyama M, Torrents D, Bork P: PAL2NAL: robust conversion of Res 2003, 31(13):3381-3385. protein sequence alignments into the corresponding codon 61. Laskowski RA, MacArthur MW, Moss DS, Thornton JM: PRO- alignments. Nucleic Acids Res 2006:W609-612. CHECK: a program to check the stereochemical quality of 35. Rozas J, Sanchez-DelBarrio JC, Messeguer X, Rozas R: DnaSP, DNA protein structures. J Appl Cryst 1993, 26(2):283-291. polymorphism analyses by the coalescent and other meth- 62. Moreira D, Philippe H: Molecular phylogeny: pitfalls and ods. Bioinformatics 2003, 19(18):2496-2497. progress. Int Microbiol 2000, 3(1):9-16. 36. Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic 63. Anisimova M, Nielsen R, Yang Z: Effect of recombination on the inference under mixed models. Bioinformatics 2003, accuracy of the likelihood method for detecting positive 19(12):1572-1574. selection at amino acid sites. Genetics 2003, 164(3):1229-1236. 37. Huelsenbeck JP, Ronquist F: MRBAYES: Bayesian inference of 64. Akashi H: Synonymous codon usage in Drosophila mela- phylogenetic trees. Bioinformatics 2001, 17(8):754-755. nogaster: natural selection and translational accuracy. Genet- 38. Posada D, Crandall KA: MODELTEST: testing the model of ics 1994, 136(3):927-935. DNA substitution. Bioinformatics 1998, 14(9):817-818. Page 14 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 65. Akashi H: Inferring Weak Selection From Patterns of Poly- morphism and Divergence at"Silent"Sites in Drosophila DNA. Genetics 1995, 139(2):1067-1076. 66. Carlini DB, Stephan W: In Vivo Introduction of Unpreferred Synonymous Codons Into the Drosophila Adh Gene Results in Reduced Levels of ADH Protein. Genetics 2003, 163(1):239-243. 67. Christina F, Bruce R, Michael W: Identification and analysis of evolutionary selection pressures acting at the molecular level in five forkhead subfamilies. BMC Evol Biol 8:. 68. Wang M, Wang Q, Zhao H, Zhang X, Pan Y: Evolutionary selec- tion pressure of forkhead domain and functional divergence. Gene 2009, 432(1-2):19-25. 69. Heide LP van der, Jacobs FMJ, Burbach JPH, Hoekman MFM, Smidt MP: FoxO6 transcriptional activity is regulated by Thr26 and Ser184, independent of nucleo-cytoplasmic shuttling. Bio- chem J 2005, 391(Pt 3):623-629. Publish with Bio Med Central and every scientist can read your work free of charge "BioMed Central will be the most significant development for disseminating the results of biomedical researc h in our lifetime." Sir Paul Nurse, Cancer Research UK Your research papers will be: available free of charge to the entire biomedical community peer reviewed and published immediately upon acceptance cited in PubMed and archived on PubMed Central yours — you keep the copyright BioMedcentral Submit your manuscript here: http://www.biomedcentral.com/info/publishing_adv.asp Page 15 of 15 (page number not for citation purposes) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png BMC Evolutionary Biology Springer Journals

FoxO gene family evolution in vertebrates

Download PDF
15 pages

Loading next page...
 
/lp/springer-journals/foxo-gene-family-evolution-in-vertebrates-d1WWvDj30b

References (71)

Publisher
Springer Journals
Copyright
Copyright © 2009 by Wang et al; licensee BioMed Central Ltd.
Subject
Life Sciences; Evolutionary Biology; Animal Systematics/Taxonomy/Biogeography; Entomology; Genetics and Population Dynamics; Life Sciences, general
eISSN
1471-2148
DOI
10.1186/1471-2148-9-222
pmid
19732467
Publisher site
See Article on Publisher Site

Abstract

Background: Forkhead box, class O (FoxO) belongs to the large family of forkhead transcription factors that are characterized by a conserved forkhead box DNA-binding domain. To date, the FoxO group has four mammalian members: FoxO1, FoxO3a, FoxO4 and FoxO6, which are orthologs of DAF16, an insulin-responsive transcription factor involved in regulating longevity of worms and flies. The degree of homology between these four members is high, especially in the forkhead domain, which contains the DNA-binding interface. Yet, mouse FoxO knockouts have revealed that each FoxO gene has its unique role in the physiological process. Whether the functional divergences are primarily due to adaptive selection pressure or relaxed selective constraint remains an open question. As such, this study aims to address the evolutionary mode of FoxO, which may lead to the functional divergence. Results: Sequence similarity searches have performed in genome and scaffold data to identify homologues of FoxO in vertebrates. Phylogenetic analysis was used to characterize the family evolutionary history by identifying two duplications early in vertebrate evolution. To determine the mode of evolution in vertebrates, we performed a rigorous statistical analysis with FoxO gene sequences, including relative rate ratio tests, branch-specific d /d ratio tests, site-specific d /d N S N S ratio tests, branch-site d /d ratio tests and clade level amino acid conservation/variation patterns N S analysis. Our results suggest that FoxO is constrained by strong purifying selection except four sites in FoxO6, which have undergone positive Darwinian selection. The functional divergence in this family is best explained by either relaxed purifying selection or positive selection. Conclusion: We present a phylogeny describing the evolutionary history of the FoxO gene family and show that the genes have evolved through duplications followed by purifying selection except for four sites in FoxO6 fixed by positive selection lie mostly within the non-conserved optimal PKB motif in the C-terminal part. Relaxed selection may play important roles in the process of functional differentiation evolved through gene duplications as well. Page 1 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 parts of the adult mouse brain. Moreover, the individual Background Mammalian FoxO proteins (FoxO1, FoxO3a, FoxO4 and disruption of Foxo1, Foxo3 and Foxo4 genes in mice FoxO6) which are homologous to Caenorhabditis elegans results in different phenotypes [14,16]. While a protein DAF-16, belong to the O ('other') class of the Fox homozygous knockout of Foxo1 (FKHR) was embryonic superfamily [1,2]. FOXO1 is the first identified member of lethal due to failures in angiogenesis and vessel forma- the FoxO family of transcription factors [3] and is tion, Foxo3a-/- (FKHRL1) and Foxo4-/- (AFX) were viable involved in the transcriptional activity of alveolar rhab- and appeared to develop normally. Later in development, domyosarcomas [3]. Since then, the discovery of mamma- Foxo3a-/- females were found to be age-dependently lian FoxO genes has grown rapidly, now FoxO proteins infertile and showed abnormal ovarian follicular develop- have been identified in several different organisms, ment. As to the physiological role, each FoxO gene exhib- including zebrafish, mouse, rat and human. As transcrip- its a distinct response under a variety of conditions [20- tion factors in the nucleus, the primary function of FoxO 22]. Therefore, despite the high sequence identity shared proteins is to bind to their cognate DNA targeting by FoxO genes domain (more then 60% in humans [17]), sequences as monomers. The co-crystal structure of HNF- the physiological roles of FoxO genes are functionally 3γ with DNA shows that there are 14 protein-DNA con- diverse in mammals. tacts distributing throughout the forkhead domain, but the third α-helix (H3) plays the most important role in a Single copy genes are thought to evolve conservatively winged helix/forkhead protein's DNA-binding specificity because of strong negative selective pressure. Gene dupli- [4]. In addition, both winged loops also make important cations produce a redundant gene copy and thus release interactions with DNA [4,5]. Although the molecular one or both copies from negative selection pressure [23]. basis of the DNA-binding specificity of FoxO transcription There are a number of models for the fate of duplicate factors is poorly understood, high-affinity DNA-binding gene that predict functional differentiation of paralogs studies have identified a consensus FoxO-recognized ele- based on protein sequence or regulatory divergence ment (FRE) as (G/C) (T/A)AA(C/T)AA [6-8]. Indeed, func- [24,25]. Currently four most prominent models are neo- tional FRE sites that match this consensus sequence have functionalization [26], subfunctionalization [24], the been identified in the promoters of many genes, such as Dykhuizen-Hartl effect [27] and adaptive diversification. Fas ligand (FasL), insulin-like growth factor binding pro- Very recently, the list has been expanded by the introduc- tein 1 (IGFBP1) and the apoptotic regulator BIM[9,10]. tion of the subneofunctionalization [28] and the adaptive Additional putative FoxO-target genes and their potential radiation [29] models that predict rapid subfunctionaliza- cis-regulatory binding sites have been predicted by sys- tion after duplication followed by a prolonged period of tematic bioinformatic approaches [11]. Thus, FoxO tran- neofunctionalization and adaptive divergence of dupli- scription factors appear to be involved in various cate genes in a process analogous to species radiations, signaling pathways and control a wide range of biochem- respectively. Thus, duplications are thought to be an ical processes including cellular differentiation, tumor important precursor of functional divergence [30]. Here, suppression, metabolism, cell-cycle arrest, cell death, and we are interested in the specific role that natural selection protection from stress [1,9,10]. In the mouse, four differ- might play in the evolutionary history of this gene dupli- ent FoxO members have been identified to date: Foxo1, cation. Foxo3, Foxo4 and Foxo6 [12,13]. FoxO6 is the latest member of the FoxO family to be cloned and shares sig- The increased availability of FoxO sequences in the public nificant sequence similarity with the other members of databases allows us to explore the functional diversity the family [13]. from a phylogenetic perspective within the FoxO family in vertebrates. The study was conducted by analyzing amino FoxO1 and its close paralogus (FoxO3, FoxO4 and acid and nucleotide-based divergence data from different FoxO6) are thought to some degree of functional diversi- species covering the entire vertebrates. Our aim was to elu- fication during development [14-16] and their potential cidate the evolutionary mechanisms operating in the physiological roles might be different [14]. Indeed, a retention of these genes and evaluate the changes in selec- rapid overview of the data collected on FoxO1, 3, 4 and 6 tion pressures following duplication. We also identified highlights how these proteins may be different. First, each the sites under positive Darwinian selection. Finally, we FoxO gene showed different expression patterns in tissues tried to map the positively selected sites to the structural [7,12,17-19]. While Foxo1 was strongly expressed in the and functional regions of FoxO molecules. striatum and neuronal subsets of the hippocampus (den- tate gyrus and the ventral/posterior part of the CA Materials and methods regions), Foxo3 was more diffusely expressed throughout Sequence Data Collection the brain including all hippocampal areas, cortex and cer- The DNA sequences and amino acids sequences of ebellum, and Foxo6 expression was eminent in various FoxO genes were downloaded from NCBI's GenBank Page 2 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 http://www.ncbi.nlm.nih.gov. PSI-BLAST searches were Relative rate tests conducted against the non-redundant database of verte- The substitution rates of the FoxO genes were compared brate genomes at NCBI (e-value cutoff = 1e-24) using the between different paralogous genes that had undergone amino acid sequences of Foxo1, Foxo3, Foxo4 and Foxo6 duplication events recently, using the RRTree software of mouse (gi: 56458, gi: 56484, gi: 54601 and gi: 329934) [41]. The orthologs FoxOs (Cifoxo, BfFoxOA and as queries. Only full length coding sequences were SpFoxOl) from amphioxus (Branchiostoma floridae), included in our analysis. Jalview 2.3 [31] was used to Ciona intestinalis and Strongylocentrotus purpuratus remove the sequences with the identity higher than 95%. were used as an outgroup. The null hypothesis is that the A table with species names, abbreviations and accession rate of substitution of the tested clade is the same as that numbers are provided in supplementary materials (Addi- of the reference group. tional file 1). Estimation of substitution rates and testing natural Sequence alignment and phylogenetic analysis selection The sequences of FoxO proteins were aligned by MUSCLE We estimated the selective pressures acting on coding [32] and the resulting alignment was manually optimized regions by applying a phylogenetic-based Maximum Like- by BioEdit [33]. Incomplete sequences, and highly diver- lihood (ML) analysis. ML estimated of the relevant param- gent regions or gaps resulting in uncertain alignments eters -as branch lengths and the ratio of the were excluded from the further analysis. The final data set nonsynonymous (d ) to synonymous substitution rates included a total of 66 sequences from 19 species. The (d ), ω = d /d -that were obtained using the codeml pro- S N S amino acid alignment was subsequently transformed into gram implemented in the PAML package version 4 [42]. an aligned cds fasta file using PAL2NAL [34] which is a The ω parameter was used as a measure of the protein program to construct multiple codon alignments from selective constraints [43]. These analyses were conducted matching amino acid sequences. The nucleotide align- under different competing evolutionary hypothesis. We ment was then converted to nexus format with DnaSP first investigated whether the distribution of selective con- [35] version 4.10 for phylogenetic analysis. straints acting on the each gene fluctuated across lineages; for that, we compared the fit to the data of the "one ratio" The full alignment of 66 sequences was used to perform model (M0), which assumes a constant selective pressure the phylogenetic analysis. Tree reconstructions were done across branches, with the "free ratios" model (FR), where by the Bayesian method from the DNA alignment done in the rate parameters were estimated independently in each the MrBayes version 3.1.2 [36,37] software package, and lineage. We also examined other evolutionary scenarios; i) rooted with the BfFoxO, Cifoxo and SpFoxO from amphi- to determine which FoxO lineage had evolved at a differ- oxus (Branchiostoma floridae), Ciona intestinalis and ent rate, as compared to the rest of the phylogeny, we Strongylocentrotus purpuratus. We analyzed four inde- applied a branch-specific model to the data. Sequences pendent runs, each using the general time reversible were divided into four groups according to their phyloge- (GTR) model plus gamma distribution plus invariant sites netic analysis, and each FoxO lineage was set as the fore- model of molecular evolution (GTR+G+I), as determined ground branch. ii) to detect sites under positive selection by Modeltest version 3.7 [38]. We ran 2 million genera- in four lineages, we applied three codon-based ML substi- tion Markov Chain Monte Carlo simulations with four tution models that are site-specific (i.e., models that allow separate chains (three heated, one cold), with the first variation in the ω ratio across sites) of [44] but assume the 500,000 generations discarded as burn-in. Trees were same selection pattern for a site in all lineages; iii) to summarized for each independent run and compared to investigate the existence of sites evolving under positive check for concordant topologies. The consensus tree of all selection in only a specific lineage, we applied the modi- compatible groupings among all runs was used in all anal- fied branch-site model A of [45] in two consecutive tests yses. (test1 and test2 in [46]) to the same alignment used for the site-based models. The model allowing for positive Synonymous codon usage analyses selection is denoted model A and the lineage to be tested Codon usage bias was estimated by the effective number is the foreground lineage, whereas the remaining ones are of codons (ENC; [39]), the frequency of optimal codons the background lineages; the multiple hypothesis testing (F ; [40]) and proportion of G and C in the third codon problem [47] was taking into account using Bonferroni's OP position (G/C 3rd). For ENC, lower values indicate correction [48]. The likelihood Ratio Test (LRT) was used stronger synonymous codon usage bias, while for F to compare the fit to the data of two nested models, OP higher values indicate stronger bias. These measures were assuming that twice the log likelihood difference between calculated for all genes using the CodonW program http:/ the two models (2ΔL) follows a χ distribution with a /bioweb.pasteur.fr/seqanal/interfaces/codonw.html and number of degrees of freedom equal to the difference in used to test whether the degree of synonymous codon the number of free parameters [49]. usage biases in individual genes. Page 3 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 We used the TreeSAAP version 3.2 [50] to determine the were identified. Inclusion of these did not improve the FoxO physicochemical properties affected by natural reliability of the phylogeny, and as the aim of this study selection. This program for examining the effects of non- was to determine the evolutionary history of the FoxO synonymous substitutions on protein evolution compares gene family, only representatives from the major verte- the observed distribution of physicochemical changes brate clades were included. inferred from a phylogenetic tree with an expected distri- Phylogenetic analyses of FoxO gene lineages bution based on the assumption of completely random amino acid replacement expected under the condition of To study the molecular evolution of vertebrate FoxO selective neutrality. For all possible pairwise amino acid genes, we carried out phylogenetic inference analyses changes, the range of effect size for each of the 31 proper- based on codon alignment and inferred their evolutionary ties was determined and equally divided into 8 magnitude history using Bayesian methods. We used the Bayesian categories. Categories 1 to 3 indicate small variation in the posterior probabilities (PPs) of each node to evaluate amino acid characteristics while categories 6 to 8 repre- clades support. Figure 1 shows the consensus phylogeny sent the most radical substitutions. For all properties that obtained for FoxO gene sequences. The vertebrate FoxOs differed significantly from neutrality, Z-scores were then were assorted well to four lineages according to their FoxO calculated in each magnitude category to determine which classification, all with high PP support values (a poorly classes contributed to this deviation. The critical Z-score supported position: 0.99 PP) indicating that the forma- values for P = 0.001 are 3.09, indicating positive selection tion of the paralogous lineages occurred before the diver- on that magnitude category, and -3.09, which indicates gence of individual species, and the orthologs FoxOs negative (purifying) selection. That is, positive and nega- (Cifoxo, BfFoxOA and SpFoxOl) from amphioxus (Bran- tive Z-scores indicate positive and purifying selection, chiostoma floridae), Ciona intestinalis and Strongylocen- respectively. Radical substitutions affecting a particular trotus purpuratus were just located as an outgroup of their property that occurred more frequently than expected by assigned lineages. From Figure 1, we inferred that two chance constituted the signature of adaptive evolution major duplications had occurred early in the vertebrate [51]. lineages. The first duplication led to the emergence of two lineages which evolved into FoxO3/6 and FoxO1/4, and Testing functional divergence and structure analysis the second duplication, also early in vertebrate evolution, To study the functional divergence and structural differ- resulted in FoxO6 and FoxO3, and FoxO1 and FoxO4. ences after the gene duplication, we used the Diverge 2.0 Phylogenetic tree shows that the FoxO6 gene cluster has software to estimate the type I (θ ) and type II (θ ) func- long branches, an indication of fast-evolving lineage with I II tional divergence coefficients [52,53] among paralogous a large number of structural changes accumulating on proteins. Type I and type II refer to shifts in the evolution- them. ary rate pattern after the emergence of a new phylogenetic cluster (indicative of changes in functional constrains), Synonymous codon usage analyses and amino acid replacements completely fixed between We investigated the relationship between nucleotide con- duplicates (resulting in cluster-specific alterations of tent and codon usage by calculating different indices amino acid physiochemical properties), respectively. (Table 1) for each of the FoxO genes. We could see from Table 1 that the effective number of codons (ENC) Genes which have been predicted to subject to positive decreased with the corresponding increase of GC3. The selection were used to search for homologous sequences effective number of codons [39] is a measure of the even- in the PDB database of protein structures http:// ness of codon usage among the 61 sense codons. An www.rcsb.org/pdb/home/home.do using Blastp [54,55]. extreme case is that all codons are used equally frequently The Rasmol http://rasmol.org/ was used for all structural (given the observed frequencies of amino acids), then the manipulations and highlighting the relevant amino acid effective number of codons is 61. Reversely, only single replacements identified in the evolutionary analyses.. codon is used for each amino acid, the effective number of codons is reduced to 20. Therefore, FoxO6 gene was Results more biased than other FoxO genes as evidenced by their Sequence similarity searches and multiple alignments lower ENC values. In most cases, the observed number fell Available FoxO1, FoxO3, FoxO4 and FoxO6 sequences somewhere between the two extremes. Figure 2 shows the were retrieved from 19 species ranging from amphioxus relationship between the effective number of codons (Branchiostoma floridae) to mammals. Additional file 1 (ENC) and the GC content at the third position of each outlines the sequences (protein and DNA) used in the gene (GC3). This Figure also contains a reference line phylogenetic analyses. The highly conserved forkhead (GCref) showing the expected position of genes whose domain remained in all alignments. It should be noted codon usage is constrained solely by the nucleotide com- that additional FoxO genes for eutherians and teleosts position at the third codon position. From Figure 2, it can Page 4 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 P Figure 1 hylogenetic relationships of DNA sequences within the FoxO family Phylogenetic relationships of DNA sequences within the FoxO family. Phylogenetic tree based on the nucleotide sequence data. The numbers indicate the Bayesian probabilities for each phylogenetic clade. Shaded boxes denote the four lin- eages and one outgroup. The scale bars represent codon substitutions per site. be seen that the observed value of ENC tracks the refer- an outgroup, we evaluated the relative rates between FoxO ence line quite closely. This indicates that the nucleotide gene clusters. The analysis (Table 2) revealed that the composition at the third codon position is a major deter- FoxO6 lineage exhibited accelerated nonsynonymous minant of the effective number of codons. substitutions with respect to FoxO3 (p-value = 0.00163, Bonferroni correction) and FoxO1 (p-value = 0.0193, Relative rates of evolution of FoxO6 lineage Bonferroni correction), and that FoxO4 genes were not Using the orthologs FoxOs (Cifoxo, BfFoxOA and accelerated with respect to the other FoxO lineages. There- SpFoxOl) from amphioxus (Branchiostoma floridae), fore, evolutionary-rate changes may have occurred follow- Ciona intestinalis and Strongylocentrotus purpuratus as ing FoxO gene duplications in the evolutionary process. Table 1: Mean values of GC%, GC3%, ENC, CAI and Fop of the FoxO genes Gene GC% GC3% ENC CAI Fop FoxO1 0.5775 ± 0.0416 0.6478 ± 0.1009 51.0256 ± 5.1246 0.0661 ± 0.0119 0.3506 ± 0.0259 FoxO3 0.5797 ± 0.0538 0.6744 ± 0.1194 47.9308 ± 6.3169 0.0765 ± 0.0171 0.3851 ± 0.0283 FoxO4 0.5974 ± 0.0430 0.6082 ± 0.0459 50.9273 ± 2.5006 0.0660 ± 0.0129 0.3522 ± 0.0385 FoxO6 0.6773 ± 0.0867 0.7695 ± 0.1486 42.4620 ± 8.7015 0.0462 ± 0.0212 0.3032 ± 0.0599 Note: Mean ± Standard deviation Page 5 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 Th Figure 2 e effective number of codons (Nc) plotted for each FoxO genes The effective number of codons (Nc) plotted for each FoxO genes. The FoxO genes highlighted in blue dot. The GC(ref) line -- shown in red -- is the expected position of genes whose codon usage is only determined by the GC content at the third positions of codons (GC3s). Selective constraints and functional divergence pairs. As expected, most amino acids had very low poste- Gene duplication-specific changes in the substitution rior probability (PP) values and, therefore, they would not rates (type I functional divergence) might reflect the dif- be involved in the hypothetical functional divergence ference in evolutionary rate at amino acid sites after gene (Figure 3). Specifically, we detected 32 and 15 amino acid duplication [52,53]. We found significant evidence of positions which presumably submitted to altered func- type I functional divergence for comparisons between dif- tional constraints when the PP threshold values were set ferent gene clusters (θ = 0.23 ~0.40, P < 0.01; Table 3); to 0.87 and 0.95, respectively. Type I sites are defined as namely, there were some amino acid sites with discrepan- those with an amino acid that is conserved in one cluster cies in their evolutionary rate between these paralogous but variable in the sister cluster, implying that the site is Page 6 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 Table 2: Evolutionary Rate of the FoxO Gene Families Lineage1 Lineage2 Ka1 Ka2 dKa sd_dKa ratio_Ka P_Ka FoxO6 FoxO4 1.03324 0.997808 0.035436 0.046871 0.756029 0.449673 FoxO6 FoxO3 1.03379 0.902464 0.131326 0.041694 3.14977 0.001638 FoxO6 FoxO1 1.05943 0.952932 0.106501 0.04553 2.33915 0.019336 FoxO4 FoxO3 0.984938 0.921849 0.063089 0.044978 1.40267 0.16075 FoxO4 FoxO1 0.997224 0.942435 0.054789 0.048375 1.13257 0.257401 FoxO3 FoxO1 0.94013 0.949027 -0.0089 0.040392 -0.22028 0.825657 Note: Ka corresponds to the mean evolutionary rate measured as the number of nonsynonymous substitutions per site. dKa is the mean rate difference between the two lineages. sd_dKa is the standard deviation and ratio_Ka the ratio between dKa and sd_dKa. The P_Ka column corresponds to the P value associated to the test under structural/functional constraints in the first cluster was substantially smaller than 1 (one ratio model ω = that is absent in the variable cluster [56]. 0.084, Table 4) that indicated that purifying selection had been the predominant force acting on the evolution of Recently, a method has been developed to test for type II these vertebrate FoxOs. Omega estimates for FoxO1, functional divergence [57]. Type II sites are those that are FoxO3, FoxO4 and FoxO6 phylogenies were 0.09583, highly conserved in both clusters but are fixed for amino 0.08311, 0.14088 and 0.13464, respectively. Selective acids with different biochemical properties between sister constrains, however, are unevenly distributed across the clusters, implying these residues are responsible for the phylogeny (FR model; 2ΔL = 421.20, P < 0.001). We then functional differences between these groups. Although at ran the branch model using each FoxO lineage as the fore- least one site with evidence of type II divergence was ground branch. In this model the estimated ω was 0.0758 found for comparisons between FoxO1/FoxO3, FoxO3/ for the FoxO1, and 0.0898 for the background branches. FoxO4, and FoxO1/FoxO4 clusters, the θ values are A LR test indicated that the two-ratio model was not sig- II extremely small (θ = 0.005 ~0.074) that highlighted the nificantly different from the M0 model (2ΔL = 3.05, P > II conservation between different clusters. These results are 0.05, df = 1, Table 4). In contrast to FoxO1 analysis, the ω not unexpected given that this method calculates θ across values of the FoxO3, FoxO4 and FoxO6 lineages were dif- all sites in an alignment and thus effectively averages site- ferent from the rest of the phylogeny as the LR tests indi- wise θ values. With only ~3% of sites/cluster showing a cated that the two-ratio model fit the data better than the pattern of type II divergence in our concatenated align- M0 model for these three genes (P < 0.05). Unfortunately, ment, it is not likely that the ~9 possible type II sites have the ω estimates for FoxO3, FoxO4 and FoxO6 were not θ values high enough to compensate for the extremely indicative of positive selection, they were rather indicative II low θ values of the over 300 sites with θ effectively equal of relaxed constraint. II II to zero. Our results are similar to the analysis of Hox-gene [30]. Along with lineage heterogeneity, variations in ω across sites can also occur. Theoretically, different protein The analysis of the nonsynonymous to synonymous sub- regions with different functions may experience different stitution rate ratio can also be used to detect functional selection pressures, which can be tested by fitting the data differentiation. We estimated ω as an average over all sites to a model comprising different site classes. The results and branches from the FoxO paralogus MSA and the ratio were shown in Table 5, for each lineage, the M3 vs M0 LRT was significant, indicating that one category of ω wasn't fit Table 3: Maximum likelihood estimates of the coefficient of data well to describe the variability in selection pressure functional divergence (θ) from pairwise comparisons between across amino acid sites. The tests contrasting the models FoxO groups M1a against M2a resulted in the P value of 1 for all the a b c Comparison θ SE (θ)LRT (θ)sig. groups suggesting a lack of power and the amino acid changes within each cluster were neutral or under negative FoxO1 Vs FoxO3 0.33 0.05 50.89 P < 0.01 selection. M1a, the parameter estimates for the least FoxO1 Vs FoxO4 0.29 0.04 41.23 P < 0.01 parameter rich model describes that most sites with low ω FoxO1 Vs FoxO6 0.3 0.06 20.43 P < 0.01 estimates (indicative of strong selective constraints), that FoxO3 Vs FoxO4 0.23 0.05 19.86 P < 0.01 FoxO3 Vs FoxO6 0.4 0.05 57.92 P < 0.01 is, 82% of FoxO1 sites were under strong purifying selec- FoxO4 Vs FoxO6 0.24 0.05 21.61 P < 0.01 tion, compared to 83% for FoxO3, 74% for FoxO4 and 66% for FoxO6. The test using M7 and M8, which allows Note: θ is the coefficient of functional divergence; for beta-distributed site-specific ω ratio, detected 2 groups SE(θ) standard error; under possible positive selection at 0.05 significance level, LRT(θ) is a likelihood ratio test; Page 7 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 Typ Figure 3 e I functional divergence among the FoxO members Type I functional divergence among the FoxO members. Posterior probability (PP) profiles of the site-specific type I functional divergence. The positions with gaps involved in each paralogous comparison were not considered. Red line indicates cutoff = 0.95, while green cutoff = 0.87. one with ω = 1.36 and the other with ω = 127.02 (Table method has been developed that allows variation in ω 5). In order to test whether the estimated ω is significantly across individual codons on a specific lineage [46,58]. greater than 1, model M8 was compared with a more This model (MA) designates two categories of branches, restricted null model (M8a). For FoxO3 gene, Model M8 again foreground and background, where positive selec- did not significantly differ from model M8a (2ΔL = 0.067, tion is modeled only on the foreground branch. We then P > 0.05, df = 1). For FoxO6 gene, ω = 127.02 was signifi- applied the branch-site approach (using some pre-speci- cantly different than 1 (2ΔL = 51.92, P < 0.01, df = 1). We fied branches, i.e., foreground branches), designating also used the BEB estimation method in model M8 [46] to each FoxO gene as the foreground branch, to assess identify sites under possible positive selection. whether molecular adaptation occurred in the evolution of the FoxO genes. The results of this analysis exhibited Since positive selection will likely affect a few amino acids several positions with evidence of relaxed selection (the at specific lineages on the phylogeny, models estimating ω test 1 was significant) (Table 6). However, we could not ratios averaged by codons or by lineages are certainly reject the null hypothesis of the test 2 (ω = 1) (result not highly conservative. For this reason, a branch-sites Page 8 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 Table 4: LRTs done to detect heterogeneous selection regimes among lineages for each gene model df Parameter estimates lnL 2l p value Branch-specific models -22943.8 M0(one-ratio) ω = 0.08442 FoxO1 two-ratio vs one-ratio 1 ω = 0.0898ω = 0.0758 -22942.3 3.047108 p > 0.05 0 1 FoxO3 two-ratio vs one-ratio 1 ω = 0.0910ω = 0.0730 -22941.3 4.87108 p < 0.05 0 1 FoxO4 two-ratio vs one-ratio 1 ω = 0.0811ω = 0.1044 -22941.8 3.983174 P < 0.05 0 1 FoxO6 = 0.0786ω = 0.1358 -22935.8 16.00121 P < 0.01 two-ratio vs one-ratio 1 ω 0 1 shown); thus, these analyses do not provide any evidence expected by chance for 2 of the properties (alpha-helical for directional selection on the FoxO lineages. tendencies and compressibility). The molecular adaptation processes occurred after the Spatial distributions of possible selected FoxO6 Sites on three-dimensional structure gene duplication event were also investigated by compar- ing the magnitude of the physicochemical changes pro- Because of the evidence for possible positive selection on duced by the observed amino acid replacements with FoxO6, we predicted positively selected codon sites using those expected at random [59]. We used the program a Bayes empirical Bayes (BEB) method [45]. The sites Tree-SAAP [50] to model how 31 different physicochemi- under selection in FoxO6 are listed in Table 5. Four codon cal properties were affected by amino acid substitutions in sites were identified as positively selected at a BEB poste- each FoxO gene. Consistent overrepresentation of radical rior probability threshold of 95%. Moreover, 7 amino amino acid changes (i.e., categories 7 and 8) would sug- acid residues presumably submitted to altered functional gest repeated adaptive substitution [51]. The results indi- constraints were identified by both PAML 4 and Diverge cated that, some amino acid replacements altering these 2.0 analysis (Table 6). In order to plot positive selected physicochemical properties in the FoxO1 and FoxO3 pro- sites onto mouse (Foxo6) three-dimensional model, we teins accumulated more (or less) often than expected by first built an energy-minimized model using a homology chance (likely reflecting fitness differences) (supplemen- modeling approach [60]. The PDB entry with the highest tary materials (Additional file 1)). Moreover, for each sequence similarity -identified in the PSI-BLAST- corre- physicochemical property, the distribution of the Z-scores sponds to the human FOXO3A (PDB: 2k86). We used this across 8 magnitude classes [51] indicated that, amino acid entry as a template for the modelling. The in silico stereo- substitutions occurred less often than expected by chance chemical quality analysis [61] indicated that the generated at the most extreme magnitude-classes (supplementary model had a moderate quality (with the percentage of res- Additional file 1); these FoxO1 and FoxO3 protein prop- idues in most favored regions being no lower than the erties, therefore, were likely evolving under purifying 82.8%), with only 1.1% in disallowed regions. As selection. For FoxO4 and FoxO6 genes, less physicochem- expected, the modeled structure was roughly similar to the ical properties were affected by amino acid substitutions. template, with the three helices and two wing loops typi- The FoxO6 gene, on the contrary, seems to evolve positive cal of the Fox family in equivalent positions and with a selection, because category 8 occurs more frequently than similar predicted folding (Figure 4). Taken together, these Table 5: Site model analyses for the FoxO1, FoxO3, FoxO4 and FoxO6 phylogenies Models comparison M3 vs M0 M2a vs M1a M8 vs M7 Gene 2ΔL = (L1-L0) p-value 2ΔL = (L1-L0) p-value 2ΔL = (L1-L0) p-value ω- value Positively selected sites FoxO1 580.7848 p < 0.01 0 1 0 1 FoxO3 902.91 p < 0.01 0 1 8.23 p < 0.05 ω = 1.36 66 L (p > 0.90) FoxO4 178.24 p < 0.01 0 1 0 1 FoxO6 455.31 p < 0.01 0 1 69.66 p < 0.01 ω = 127.02 264K* 266P* 434G* 439T* Note: *P > 0.95 Page 9 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 Table 6: Parameter estimations and likelihood ratio tests for the branch-site models a b c df Parameter estimates lnL 2 P value Positive selected sites FoxO1 MA Vs M1a (test 1) 2 p = 0.70592 p = 0.18493 (p = 0.10915) -22667.231 81.89419 P < 0.01 213S** 216S* 219S* 252M* 276V** 0 1 2 w = 0.07205 (w = 1.00000) w = 285P* 296L** 306A** 340F* 360E* 0 1 2 1.00000 FoxO3 MA Vs M1a (test 1) 2 p0 = 0.77299 p1 = 0.10453 (p2 = -22646.766 122.8257 P < 0.01 5H** 25D* 26F** 33D** 34L**37N** 0.12248) w0 = 0.06958 (w1 = 1.00000) 217A* 231G** 329G* w2 = 1.00000 FoxO4 MA Vs M1a (test 1) 2 p0 = 0.67705 p1 = 0.18055 (p2 = 0.1424) -22660.108 96.14052 P < 0.01 7V** 173R** 194T** 201I** 202L** w0 = 0.07468 (w1 = 1.00000) w2 = 211F**223H* 225P** 242T* 254R** 1.00000 314S** FoxO6 MA Vs M1a (test 1) 2 p0 = 0.58516 p1 = 0.13030 (p2 = -22527.632 361.0936 P < 0.01 42Q** 46K** 155I* 164T** 165N** 0.28454) w0 = 0.06572 (w1 = 1.00000) 173R* 174E** 176E** 178L** 179F** w2 = 1.05311 180C** 188I* 189V** 203L* 207R* 223H** 230I** 231G* 232Y** 233K** 234N** 237Y** 258S** 265N* 269T** 271E** 272N** 273E** 274V** 275H** 276V** 277S** 278Q* 279G** 280L** 281H** 282P** 283S** 286N* 314S** 316V** 320H** 330Y* 366T** 367G** 368T** 369P* Note: The number of free parameters; Likelihood of the model; 2(l -l ); 1 0 * P > 0.95; ** P > 0.99; The bold amino acid residues were also found to be implicated in the functional divergence (implemented in Diverge 2.0) between FoxO data suggested that the model was stereochemically valid, performed firstly to the resolution of the evolutionary and therefore suitable for further sequence-structural relationships of these FoxOs using molecular sequence analysis. Unfortunately, we could not map any positive data. Whereas, the incorrect phylogenetic topology result- selected sites onto the surface of the 3D structure (Figure ing from mutationally saturated positions, inadequate 4A) because the crystal structures about Fox proteins are modeling of the evolutionary process and systematic bias mainly focused on the forkhead DNA-binding domain. due to variable rates of evolution among species or within Whereas the positive selected sites were mainly located in sequences [62] may make LRT generate many false posi- the region of N-terminal and C-terminal of FoxO6, which tives. Anisimova et al (2003) examined the effect of also indicated that FoxO6 underwent strong constraint on assuming a "wrong" tree [63], and he found that LRT the forkhead domain as well (Sequence logo of the fork- falsely suggested positive selection in 96% of the repli- head domain, Figure 4). cates in the M0-M3 comparison and in 86% of the repli- cates in the M7-M8 comparison at the = 5% significance Discussion level. In order to overcome this problem, we adopted a It has long been know that FoxO transcription factors play number of ways in combination. Firstly, the addition of important roles in regulating various signals, which trans- more taxa to the dataset: denser sampling of species can late various environmental stimuli into dynamic gene reduce the effect of long branch attraction (LBA) by reduc- expression programs to influence many physiological and ing the overall distances between taxa. Secondly, we used pathological processes, including cancer and aging. The the best model of DNA substitution, determined by Mod- functions of FoxO proteins are regulated at multiple lev- eltest version 3.7 [38]. And finally, our inclusion of els, which include but are not limited to phosphorylation, enough sequences in each lineage helped alleviate loss of ubiquitylation and acetylation. Interestingly, all of these LRT power from short conserved sequences. From phylo- activities affect nuclear/cytoplasmic trafficking of FoxO genetic result, we focused on the 2 main duplications proteins. The specific function of each member of this along the evolutionary history of FoxO genes, the FoxO1- family is different [14]. As the accumulation of gene FoxO4 and the FoxO3-FoxO6 duplication, which formed sequences in the database, it is feasible to explore the four gene lineages (all with the high confidence values, functional diversity from a phylogenetic perspective. We Page 10 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 The Figure 4 modeled structure of mouse Foxo6 The modeled structure of mouse Foxo6. A. The structure of the forkhead domain; B. Sequence logo of the forkhead domain and surrounding amino acids. nearly 100% posterior probability in Bayesian analysis) It is widely accepted that gene duplication can create and used for further analysis. opportunities for functional divergence in paralogues. Divergence is thought to occur where one duplicate Codon bias is largely thought to be due to weak selection retains the original protein function and the other accu- acting to optimize protein production [64-66]. Selection mulates changes, (either through redundancy or by posi- intensity for codon usage bias, therefore, is expected to tive selection) or alternatively, through the partitioning of vary among genes. Our survey of synonymous codon the functions of an unduplicated ancestor protein. What- usage in FoxO genes revealed a strong and consistent pat- ever the mechanism, if functional divergence has occurred tern of codon bias in genes with FoxO6 relative to those between duplicated genes, then it should be observable as with FoxO1, FoxO3 and FoxO4 (Table 1). At the same changes within their coding regions. time, there appears to be some conflicting results observed between F and ENC, which may be caused by The functional divergence of FoxO genes has been studied OP differences in the way that the two methods estimate by [14]. The branch length leading to the FoxO6 clade is codon bias. F is based on the frequency of a set of spe- extended relative to other FoxO genes in gene phylogeny, OP cies specific "optimal" codons, while ENC is based on the (Figure 1). This suggested that after the duplications, observed number of codons used for each amino acid. FoxO6 evolved at a faster rate than other FoxO genes. This Thus it is possible for the two methods to give different result was confirmed by significant relative rate test results estimates of codon bias. for FoxO gene lineages (Table 2). In this sense, we per- Page 11 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 formed type I functional divergence analysis, and we pared with a more restricted null model (M8a), where the detected significant type I divergence among FoxOs. The extra site class was constrained to have ω = 1. When per- comparison between the FoxO3 and FoxO6 groups forming this analysis with the sequence data from the showed the highest value for θ (0.40 ± 0.05), suggesting FoxO3 gene, model M8 did not significantly differ from that these two groups had diverged considerably more at model M8a (2ΔL = 0.067, df = 1, P > 0.05), indicating that the functional level. Next, DIVERGE was used to establish the estimated ω = 1.36 was not significantly different than the posterior probability of type I divergence at each site 1 and that there was little indication of positive selection in the alignment, employing two cut-off posterior proba- in this gene. Further, our observation of strong purifying bility values of 0.87 and 0.95. However, the cutoff value selection being the primary mode of evolution through- for residue selection is an empirical decision and is out the FoxO phylogeny is consistent with the findings of expected to depend on the intrinsic properties of the pro- a recent study about forkhead family [67,68]. tein family being analyzed. Thus, we predicted 32 candi- date functional divergence-related sites using 0.87 as a When we performed branch-site model analysis, we cutoff value (supplementary materials (Additional file found relaxed functional constraint was most consistent 1)). When we narrowed our criteria to 0.95, we got 15 can- with the molecular evolutionary analyses of the FoxO didate residues as the most likely candidate sites for type I data. Our conclusion is in contrast to a previous study functional divergence (supplementary materials (Addi- which concluded that one site was found to be under pos- tional file 1)), but we lacked a way to verify how the rate- itive selection in the FoxO3 lineage [67]. In our paper, we shift in these sites contributed to functional divergence applied the modified branch-site model A of [45] in two among the FoxO gene groups. For comparative purposes, consecutive tests (test1 and test2 in [46]) to the same the same alignment and phylogeny was submitted to a ML alignment used for the site-based models, but we could LRT, which, like the Bayesian method provided a statisti- not detected the positive selection site. To determine why cal framework where evolutionary rate shifts at particular these two articles are giving drastically different results, we protein positions could be established [46]. At last, the had a look at the sequences used for branch-site model statistically most likely positions predicted to underlie analysis in [67], and we found that only 12 sequences (4 functional divergence were agreement by both methods, for FoxO3) used for testing evolution selection. Test 2 is a particularly for the highest-ranking candidates (Table 6). more direct test for identifying positive selections in the foreground branch, and significant LRT from test 1 can be In this study we used codon substitution models imple- resulted from either positive selection or relaxed selective mented through a maximum likelihood framework to constraint in the foreground branch [46], however, the estimate the rate of evolution at silent and replacement power of this method may be limited when sample size or sites in FoxO1 and its paralogs, FoxO3, FoxO4 and divergence time is low [46]. Therefore, we concluded that FoxO6. Different models were used to investigate varia- the contradiction between our results and the previous tion in the rate of evolution between lineages of a phylog- study [67] due to the number of sequences used for anal- eny, and to estimate ω for specific lineages and sites across ysis. Moreover, our work on site-model analysis, relative phylogenies. Our objective was to determine the mode of rate test and physicochemical changes indicated that evolution on each FoxO gene lineage, and to determine FoxO6 was under positive selection. Four positive selected whether increased positive selection or decreased con- sites were identified by site-model analysis, two (264K, 337 339 and Pro ) of them straint led to the functional divergence of FoxO genes. As 266P, corresponding to mouse Gly we have demonstrated, variation among branch and sites fell into the region of non-conserved optimal PKB motif was observed in the FoxO6 phylogeny. Moreover, physic- in the C-terminal part (Thr ) [69]. The C-terminal PKB ochemical amino acid properties analysis also provided recognition sequence is not conserved in FoxO6 [13]. evidence that the entire FoxO6 gene had experienced Besides a PKB phosphorylation motif, this region contains repeated episodes of adaptive evolution. The site models a stretch of 3 additional serine residues, present in the showed that adaptation had appeared at four sites located other members of the FoxO group, but FoxO6. All these at C-terminal of FoxO6. For FoxO3 gene, Model M8 fits suggest that these serines may be functionally important data significantly better than Model M7 (2ΔL = 8.23, df = in the other sequences analyzed with the exception of 2, P < 0.05), and because 1.5% of sites are located in the FoxO6 gene, and positive selection may lead to functional positively selected site class with ω = 1.36, weak positive divergence between FoxO6 and the other members too. 545 550 selection may be indicated with this comparison. How- Another two positive selective site are Gly and Pro in ever, it has been found that a poor fit of the data to a beta the mouse Foxo6, and the functional role remains elusive. distribution may result in a high frequency of significant That is to say that the real reason for their accelerated evo- tests when comparing models M7 and M8 even in the lution is unclear. However, it should be mentioned that absence of positive selection. To take account for the ele- there is a gap in the knowledge of the relationship vated type I error rates, the original model M8 was com- between amino acid sequence and structure for full-length Page 12 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 FoxO sequence, and we are unable to speculate on the Additional material particular role of this region in these FoxO6 genes. Unfor- tunately, the shared evolutionary history and molecular Additional file 1 selection alone cannot be used as the unique criterion to Excel spreadsheet containing: A list of species, species abbreviations, and infer protein function, and the true nature of each FoxO accession numbers for sequences used in the study/A list of statistically sig- gene needs to be determined experimentally and inde- nificant physicochemical amino acid properties for each FoxO gene/A list pendently. Therefore, the positively selected site may play of the candidate residues as the most likely candidate sites for type I func- tional divergence. an important functional role and could represent an inter- Click here for file esting target site for future mutagenesis experiment thus [http://www.biomedcentral.com/content/supplementary/1471- facilitating our understanding of the structure-function 2148-9-222-S1.xls] relationships in FoxO genes. Molecular testing is required to validate this hypothesis. The result from branch-site analysis (relaxation of functional constraints) of FoxO6 also differs somewhat from previously signature of posi- Acknowledgements tive selection. We infer that the weak positive selection The authors wish to express their gratitude to the members of animal sci- and multiple branches are considered as foreground ences laboratory of Shanghai Jiao Tong University. The authors also thank branch may explain this phenomenon, because power editor for his suggestions about the manuscript and Jing Li for her help in will be reduced unless the same sites and selective con- revising the manuscript. This work is supported by the National High Tech- nology Research and Development Program of China (863 project) (grant straints are occurring along all foreground branches [46]. no. 2006AA10Z1E3, 2008AA101002), the National 973 Key Basic Research Program (grant no, 2006CB102102, 2004CB117500) and the National Nat- Conclusion ural Science Foundation of China (grant no. 30671492, 30871782). Genomic data have provided an opportunity to gain a bet- ter understanding about the evolution of FoxOs using References phylogenetic analyses. The FoxO gene family phylogeny 1. Barthel A, Schmoll D, Unterman TG: FoxO proteins in insulin showed that two duplications took place early in the evo- action and metabolism. Trends Endocrinol Metab 2005, 16(4):183-189. lution of vertebrates and triggered diversification of the 2. Kaestner KH, Knochel W, Martinez DE: Unified nomenclature for FoxO gene family into four groups. However, further the winged helix/forkhead transcription factors. Genes Dev 2000, 14(2):142-146. genome projects on a greater diversity of evolutionary lin- 3. Fredericks WJ, Galili N, Mukhopadhyay S, Rovera G, Bennicelli J, Barr eages would help better understand the gene-duplication FG, Rauscher FJ: The PAX3-FKHR fusion protein created by history. The relative rate analysis and physicochemical the t (2; 13) translocation in alveolar rhabdomyosarcomas is a more potent transcriptional activator than PAX3. Mol Cell changes indicated that FoxO6 seemed to be different from Biol 1995, 15(3):1522-1535. other members. Evolutionary rate analysis showed that 4. Clark KL, Halay ED, Lai E, Burley SK: Co-crystal structure of the molecular adaptation can also play an important role in HNF-3/fork head DNA-recognition motif resembles histone H5. Nature 1993, 364(6436):412-420. the evolution of this gene family. Indeed, positive selec- 5. Boura E, Silhan J, Herman P, Vecer J, Sulc M, Teisinger J, Obsilova V, tion was likely involved in the functional differentiation Obsil T: Both the N-terminal Loop and Wing W2 of the Fork- head Domain of Transcription Factor Foxo4 Are Important of FoxO6 gene; likewise, relaxed selection might play for DNA Binding. J Biol Chem 2007, 282(11):8265-8275. important roles over evolutionary time and shape varia- 6. Biggs Iii WH, Meisenhelder J, Hunter T, Cavenee WK, Arden KC: tion of some members of the family. Considering the evo- Protein kinase B/Akt-mediated phosphorylation promotes nuclear exclusion of the winged helix transcription factor lutionary history of the FoxO gene family, we provided FKHR1. Proc Natl Acad Sci USA 1999, 96(13):7421-7426. insight into which amino acid residues might have under- 7. Furuyama T, Nakazawa T, Nakano I, Mori N: Identification of the differential distribution patterns of mRNAs and consensus gone positive selection and could be targeted for site- binding sequences for mouse DAF-16 homologues. Biochem J directed mutagenesis. However, the identification of four 2000, 349(Pt 2):629-634. sites under positive selection requires supporting evidence 8. Gilley J, Coffer PJ, Ham J: FOXO transcription factors directly activate bim gene expression and promote apoptosis in sym- from further functional experiments to demonstrate the pathetic neurons. J Cell Biol 2003, 162(4):613-622. adaptive character of the amino acids. All these studies 9. Accili D, Arden KC: FoxOs at the Crossroads of Cellular and experiments will certainly contribute to better under- Metabolism, Differentiation, and Transformation. Cell 2004, 117(4):421-426. stand the precise role of natural selection and functional 10. Greer EL, Brunet A: FOXO transcription factors at the inter- divergence of this family. face between longevity and tumor suppression. Oncogene 2005, 24:7410-7425. 11. Xuan Z, Zhang MQ: From worm to human: bioinformatics Authors' contributions approaches to identify FOXO target genes. Mech Ageing Dev MH, Wang conceived and supervised all research. QS, 2005, 126(1):209-215. 12. Biggs Iii WH, Cavenee Karen CWK: Identification and character- Wang helped with the Type I analysis. XZ, Zhang and HB, ization of members of the FKHR (FOX O) subclass of Zhao gave the suggestion on the analysis of physicochem- winged-helix transcription factors in the mouse. Mamm Genome 2001, 12(6):416-425. ical properties. YC, Pan gave the suggestion on discussion 13. Jacobs FM, Heide LP van der, Wijchers PJ, Burbach JP, Hoekman MF, and approved the final version. Smidt MP: FoxO6, a novel member of the FoxO class of tran- Page 13 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 scription factors with distinct shuttling dynamics. J Biol Chem 39. Wright F: The'effective number of codons' used in a gene. 2003, 278(38):35959-35967. Gene 1990, 87(1):23-29. 14. Hosaka T, Biggs WH, Tieu D, Boyer AD, Varki NM, Cavenee WK, 40. Ikemura T: Correlation between the abundance of Escherichia Arden KC: Disruption of forkhead transcription factor coli transfer RNAs and the occurrence of the respective (FOXO) family members in mice reveals their functional codons in its protein genes: a proposal for a synonymous diversification. Proc Natl Acad Sci USA 2004, 101(9):2975-2980. codon choice that is optimal for the E. coli translational sys- 15. Arden KC: FoxOs in tumor suppression and stem cell mainte- tem. J Mol Biol 1981, 151(3):389-409. nance. Cell 2007, 128(2):235-237. 41. Robinson-Rechavi M, Huchon D: RRTree: Relative-Rate Tests 16. Castrillon DH, Miao L, Kollipara R, Horner JW, DePinho RA: Sup- between groups of sequences on a phylogenetic tree. Bioinfor- pression of ovarian follicle activation in mice by the tran- matics 2000, 16(3):296-297. scription factor Foxo3a. Science 2003, 301(5630):215-218. 42. Yang Z: PAML 4: phylogenetic analysis by maximum likeli- 17. Anderson MJ, Viars CS, Czekay S, Cavenee WK, Arden KC: Cloning hood. Mol Biol Evol 2007, 24(8):1586-1591. and Characterization of Three Human Forkhead Genes 43. Yang Z: Inference of selection from multiple species align- That Comprise an FKHR-like Gene Subfamily. Genomics 1998, ments. Curr Opin Genet Dev 2002, 12(6):688-694. 47(2):187-199. 44. Yang Z, Nielsen R, Goldman N, Pedersen A: Codon-substitution 18. Kitamura T, Nakae J, Kitamura Y, Kido Y, Biggs WH 3rd, Wright CV, models for heterogeneous selection pressure at amino acid White MF, Arden KC, Accili D: The forkhead transcription fac- sites. Genetics 2000, 155(1):431-449. tor Foxo1 links insulin signaling to Pdx1 regulation of pan- 45. Yang Z, Wong W, Nielsen R: Bayes empirical bayes inference of creatic beta cell growth. J Clin Invest 2002, 110(12):1839-1847. amino acid sites under positive selection. Mol Biol Evol 2005, 19. Hoekman MF, Jacobs FM, Smidt MP, Burbach JP: Spatial and tem- 22(4):1107-1118. poral expression of FoxO transcription factors in the devel- 46. Zhang J, Nielsen R, Yang Z: Evaluation of an improved branch- oping and adult murine brain. Gene Expr Patterns 2006, site likelihood method for detecting positive selection at the 6(2):134-140. molecular level. Mol Biol Evol 2005, 22(12):2472-2479. 20. Nakae J, Kitamura T, Kitamura Y, Biggs WH 3rd, Arden KC, Accili D: 47. Anisimova M, Yang Z: Multiple hypothesis testing to detect lin- The forkhead transcription factor Foxo1 regulates adipocyte eages under positive selection that affects only a few sites. differentiation. Dev Cell 2003, 4(1):119-129. Mol Biol Evol 2007, 25(4):1219-1228. 21. Richards JS, Sharma SC, Falender AE, Lo YH: Expression of FKHR, 48. Miller R: Simultaneous statistical inference. Springer-Verlag, FKHRL1, and AFX genes in the rodent ovary: evidence for New York; 1981. regulation by IGF-I, estrogen, and the gonadotropins. Mol 49. Whelan S, Goldman N: Distributions of statistics used for the Endocrinol 2002, 16(3):580-599. comparison of models of sequence evolution in phylogenet- 22. Bois PRJ, Grosveld GC: FKHR (FOXO1a) is required for myo- ics. Mol Biol Evol 1999, 16(9):1292-1299. tube fusion of primary mouse myoblasts. EMBO J 2003, 50. Woolley S, Johnson J, Smith M, Crandall K, McClellan D: TreeSAAP: 22:1147-1157. selection on amino acid properties using phylogenetic trees. 23. Hughes J, Criscuolo F: Evolutionary history of the UCP gene Bioinformatics 2003, 19(5):671-672. family: gene duplication and selection. BMC Evol Biol 2008, 51. McClellan DA, Palfreyman EJ, Smith MJ, Moss JL, Christensen RG, 8:306. Sailsbery JK: Physicochemical evolution and molecular adapta- 24. Force A, Lynch M, Pickett FB, Amores A, Yan YL, Postlethwait J: tion of the cetacean and artiodactyl cytochrome b proteins. Preservation of duplicate genes by complementary, degen- Mol Biol Evol 2005, 22(3):437-455. erative mutations. Genetics 1999, 151(4):1531-1545. 52. Gu X: Statistical methods for testing functional divergence 25. Kimura M: The neutral theory of molecular evolution. Sci Am after gene duplication. Mol Biol Evol 1999, 16(12):1664-1674. 1979, 241(5):98-100. 53. Gu X: Maximum-likelihood approach for gene family evolu- 26. Goodman M, Moore GW, Matsuda G: Darwinian evolution in the tion under functional divergence. Mol Biol Evol 2001, genealogy of haemoglobin. Nature 1975, 253(5493):603-608. 18(4):453-464. 27. Dykhuizen D, Hartl DL: Selective neutrality of 6PGD allozymes 54. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic local in E. coli and the effects of genetic background. Genetics 1980, alignment search tool. J Mol Biol 1990, 215(3):403-410. 96(4):801-817. 55. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lip- 28. He X, Zhang J: Rapid subfunctionalization accompanied by man DJ: Gapped BLAST and PSI-BLAST: a new generation of prolonged and substantial neofunctionalization in duplicate protein database search programs. Nucleic Acids Res 1997, gene evolution. Genetics 2005, 169(2):1157-1164. 25(17):3389-3402. 29. Francino MP: An adaptive radiation model for the origin of 56. Gu X: Functional divergence in protein (family) sequence new gene functions. Nat Genet 2005, 37(6):573-577. evolution. Genetica 2003, 118(2-3):2-3. 30. Lynch VJ, Roth JJ, Wagner GP: Adaptive evolution of Hox-gene 57. Gu X: A simple statistical method for estimating type-II (clus- homeodomains after cluster duplications. BMC Evol Biol 2006, ter-specific) functional divergence of protein sequences. Mol 6:86. Biol Evol 2006, 23(10):1937-1945. 31. Clamp M, Cuff J, Searle SM, Barton GJ: The Jalview Java alignment 58. Yang Z, Nielsen R: Codon-Substitution Models for Detecting editor. Bioinformatics 2004, 20(3):426-427. Molecular Adaptation at Individual Sites Along Specific Lin- 32. Edgar RC: MUSCLE: a multiple sequence alignment method eages. Mol Biol Evol 2002, 19(6):908-917. with reduced time and space complexity. BMC Bioinformatics 59. McClellan DA, McCracken KG: Estimating the influence of selec- 2004, 5:113. tion on the variable amino acid sites of the cytochrome B 33. Hall TA: BioEdit: a user-friendly biological sequence align- protein functional domains. Mol Biol Evol 2001, 18(6):917-925. ment editor and analysis program for Windows 95/98/NT. 60. Schwede T, Kopp J, Guex N, Peitsch MC: SWISS-MODEL: An 1999:95-98. automated protein homology-modeling server. Nucleic Acids 34. Suyama M, Torrents D, Bork P: PAL2NAL: robust conversion of Res 2003, 31(13):3381-3385. protein sequence alignments into the corresponding codon 61. Laskowski RA, MacArthur MW, Moss DS, Thornton JM: PRO- alignments. Nucleic Acids Res 2006:W609-612. CHECK: a program to check the stereochemical quality of 35. Rozas J, Sanchez-DelBarrio JC, Messeguer X, Rozas R: DnaSP, DNA protein structures. J Appl Cryst 1993, 26(2):283-291. polymorphism analyses by the coalescent and other meth- 62. Moreira D, Philippe H: Molecular phylogeny: pitfalls and ods. Bioinformatics 2003, 19(18):2496-2497. progress. Int Microbiol 2000, 3(1):9-16. 36. Ronquist F, Huelsenbeck JP: MrBayes 3: Bayesian phylogenetic 63. Anisimova M, Nielsen R, Yang Z: Effect of recombination on the inference under mixed models. Bioinformatics 2003, accuracy of the likelihood method for detecting positive 19(12):1572-1574. selection at amino acid sites. Genetics 2003, 164(3):1229-1236. 37. Huelsenbeck JP, Ronquist F: MRBAYES: Bayesian inference of 64. Akashi H: Synonymous codon usage in Drosophila mela- phylogenetic trees. Bioinformatics 2001, 17(8):754-755. nogaster: natural selection and translational accuracy. Genet- 38. Posada D, Crandall KA: MODELTEST: testing the model of ics 1994, 136(3):927-935. DNA substitution. Bioinformatics 1998, 14(9):817-818. Page 14 of 15 (page number not for citation purposes) BMC Evolutionary Biology 2009, 9:222 http://www.biomedcentral.com/1471-2148/9/222 65. Akashi H: Inferring Weak Selection From Patterns of Poly- morphism and Divergence at"Silent"Sites in Drosophila DNA. Genetics 1995, 139(2):1067-1076. 66. Carlini DB, Stephan W: In Vivo Introduction of Unpreferred Synonymous Codons Into the Drosophila Adh Gene Results in Reduced Levels of ADH Protein. Genetics 2003, 163(1):239-243. 67. Christina F, Bruce R, Michael W: Identification and analysis of evolutionary selection pressures acting at the molecular level in five forkhead subfamilies. BMC Evol Biol 8:. 68. Wang M, Wang Q, Zhao H, Zhang X, Pan Y: Evolutionary selec- tion pressure of forkhead domain and functional divergence. Gene 2009, 432(1-2):19-25. 69. Heide LP van der, Jacobs FMJ, Burbach JPH, Hoekman MFM, Smidt MP: FoxO6 transcriptional activity is regulated by Thr26 and Ser184, independent of nucleo-cytoplasmic shuttling. Bio- chem J 2005, 391(Pt 3):623-629. Publish with Bio Med Central and every scientist can read your work free of charge "BioMed Central will be the most significant development for disseminating the results of biomedical researc h in our lifetime." Sir Paul Nurse, Cancer Research UK Your research papers will be: available free of charge to the entire biomedical community peer reviewed and published immediately upon acceptance cited in PubMed and archived on PubMed Central yours — you keep the copyright BioMedcentral Submit your manuscript here: http://www.biomedcentral.com/info/publishing_adv.asp Page 15 of 15 (page number not for citation purposes)

Journal

BMC Evolutionary Biology – Springer Journals

Published: Sep 7, 2009

There are no references for this article.