BACKGROUND
the reliability of phylogenetic analyses on protein-coding sequences strongly depends on models designed to approximate the substitution processes of nucleotide and amino acid.
for the evolutionary analysis of protein-coding sequences, particularly phylogenetic inference, three types of substitution models can be used, provided that both dna and protein sequences are available; nucleotide
<cit> , amino acid
<cit> , and codon substitution models
<cit> .
for closely-related sequences in which most substitutions are synonymous, amino acid substitution models cannot be used, instead nucleotide substitution models may be employed.
conversely, nucleotide substitution models should not be used in the case of diverged sequences in which synonymous substitutions are almost saturated and nonsynonymous substitutions are significant.
on the other hand, codon substitution models are appropriate to both closely-related and highly-diverged sequences with the intrinsic property of detecting both synonymous and nonsynonymous substitutions.
in a reversible markov model for substitution, a substitution rate matrix must be specified to estimate the likelihood of a phylogenetic tree.
substitution tendencies between nucleotides and those between amino acids are characteristic of each species, each organelle, and each protein family.
in the case of nucleotide substitution models, full parameterization for a substitution rate matrix is possible with  <dig> parameters; the total rate is normalized to one.
however,  <dig> parameters for an amino acid substitution rate matrix are too many to be optimized for any size of a multiple sequence alignment.
thus, empirical amino acid substitution rate matrices have been estimated from a large number of substitutions inferred on phylogenetic trees of single or many protein families; the jtt
<cit> , the wag
<cit> , and the lg
<cit>  matrices from nuclear proteins, mtrev
<cit>  from vertebrate mitochondrial proteins, cprev10
<cit>  and cprev64
<cit>  from chloroplast-encoded proteins, and flu
<cit>  from influenza proteins.
even a codon substitution rate matrix, khg
<cit> , has been estimated.
these matrices significantly differ from one another, indicating that actually substitution tendencies significantly differ among these protein families.
the estimation of a substitution rate matrix requires a large number or size of alignments with intensive calculation, and therefore is not always feasible.
however, generic empirical substitution rate matrices such as jtt, wag, lg,and khg represent the average tendencies of substitutions over various protein families by sacrificing gene-level resolution
<cit> .
on the other hand, a rate matrix such as mtrev, cprev <dig>  and flu derived from a specific protein family represents substitution tendencies characteristic of the protein but often lacks generic representation of substitution tendencies enough to be applied to other protein families.
to resolve this situation, the parametrization of an amino acid substitution rate matrix has been attempted to easily generate an alignment-specific rate matrix
<cit> .
here, we propose a different approach of employing a mechanistic codon substitution model in which the biological and evolutionary mechanisms of amino acid substitutions are taken into account.
in mechanistic codon substitution models
<cit> , substitution rates are represented as the product of a codon mutation rate and the ratio of fixation, which depend on the types of codon replacement and amino acid replacement, respectively.
hence, mutational tendencies at the nucleotide level and selection at the amino acid level can be taken into account at various levels of separation.
this is a critical difference from amino acid substitution models.
as a result, the variations of codon mutation rate and selective constraint across sites can be distinguished from each other.
we formulated a codon substitution rates between codons μ and ν as
rμν≡mμνewaμbν, with a codon mutation rate mμν, a mutation equilibrium frequency
fνmut, an equilibrium frequency fν, and a selective constraint
waμbν on the substitutions between amino acids aμ and bν <cit> .
on the basis of this mechanistic codon substitution model, we estimated the wab by maximizing the likelihood of an empirical amino acid or codon substitution frequency matrix, each of jtt, wag, lg, and khg
<cit> .
it was shown
<cit>  that the mechanistic codon substitution model with a fully-parameterized codon mutation rate matrix  and a selective constraint matrix , each element of which was approximated as a linear function ,0]) of a given estimate  of selective constraint estimated from the empirical substitution frequency matrix such as jtt and lg, outperformed both nucleotide and amino acid substitution models converted into codon substitution models in the wide range of codon sequences, from closely-related to highly-diverged protein-coding sequences.
in this codon substitution model, codon mutations with multiple nucleotide changes are also taken into account, and were shown to increase significantly the likelihood of observed substitutions
<cit> .
there have been a variety of models
<cit>  for selective pressure on amino acid replacements in mechanistic codon substitution models.
for their details, please see
<cit> , in which the present mechanistic codon substitution model was discussed in comparison with other models.
in these days, dna sequences are first analyzed and protein sequences are translated from them, and consequently codon sequences are likely available for most protein sequences in databases.
however, there are protein sequences whose codon sequences are not available or not easily retrieved because of no cross link.
in such a case, analysis at the amino acid level must be forced.
here we show that even for amino acid sequences, among which synonymous substitutions cannot be identified, the present mechanistic codon substitution model with the assumption of equal codon usage outperforms amino acid substitution models using empirical substitution rate matrices.
RESULTS
on the basis of akaike information criterion
<cit>  and bayesian information criterion
<cit> , the amino acid substitution models with the empirical amino acid substitution rate matrices, cprev64
<cit> , cprev10
<cit> , mtrev
<cit> , and flu
<cit> , as well as jtt
<cit> , wag
<cit> , and lg
<cit>  that were estimated from nuclear proteins, are compared with the mechanistic codon substitution models
<cit>  with the selective constraint matrices estimated from jtt, wag, lg, and khg
<cit>  by using the  <dig> datasets: fast-evolving interspecific mitochondrial proteins  concatenating  <dig> protein-coding genes from  <dig> mammalian species
<cit> , closely-related chloroplast-encoded proteins  concatenating  <dig> protein-coding genes from  <dig> chloroplast genomes of the major angiosperm lineages
<cit> , and ha proteins of human influenza-a h1n <dig>  consisting of  <dig> sequences.
the reference tree topologies used here as the most probable one are tree- <dig> of
<cit>  for mammalian-mtprot, the one reported in
<cit>  for cpprot- <dig>  and the one inferred by the fasttree version 2
<cit>  for ha_human-flu-a-h1n <dig>  branch optimization of phylogenetic trees and their maximum log-likelihood values are calculated using phyml
<cit>  modified for the mechanistic codon substitution model.
please see the methods section for the details of the mechanistic codon substitution model and the protein sequence data used.
the naming conventions of models employed here are described in tables
<dig> and
<dig>  the aic and bic values for these  <dig> datasets are listed in tables
<dig>
<dig>  and
<dig>  respectively.
maximum log-likelihood , aic, and bic values are represented in relative to those of the reference model, which uses the empirical amino acid substitution rate matrix specific to each dataset and has the lowest aic and bic values in all the amino acid substitution models; the best amino acid substitution model is chosen here as a reference model in order to show how superior mechanistic codon substitution models are in comparison with the best amino acid substitution model.
a random effects approach , in which the discrete gamma distribution
<cit>  with m categories of unequal probabilities is assumed for a prior probability distribution, is employed for rate variation across sites in the amino acid substitution models, and also for the variation of selective constraint or mutation rate across sites in the mechanistic codon substitution models; see additional file
<dig> for details.
suffixes "-dgmr" and "-dgms" are appended to model names to denote rate and selective constraint variations across sites, respectively.
the equilibrium frequencies of amino acids are assumed to be equal to those in the aligned sequences, and equal codon usage is assumed to calculate the equilibrium codon frequencies from them; a suffix "-f" is appended to the model names.
<cit> , cprev64
<cit> , and flu
<cit> , wag
<cit> , lg
<cit> , mtrev
<cit> , cprev10
<cit> , cprev64
<cit> , and flu
the suffix "-dgmr" means that the variation of substitution rate is approximated by a discrete gamma distribution
<cit>  with m categories of unequal probabilities; see additional file
<dig> for details.
<cit>  is used to estimate wab in eq.
<dig>  the
Δε^abc and
Δε^abv represent the effects of the mean increment of contact energy between residues and of residue-volume change due to an amino acid replacement, respectively; see supporting information, text s <dig>  in
{wabjtt/wag/lg-ml91+} estimated by maximizing the likelihood of jtt/wag/lg
<cit>  in the ml-91+ model
<cit>  are used as
{wabkhg-ml200} estimated by maximizing the likelihood of the khg codon substitution matrix
<cit>  in the ml- <dig> model
<cit>  are used as
the suffix "n" means the number of parameters optimized for the substitution rate matrix.
the suffix "-f" means that equilibrium codon frequencies are assumed to be equal to codon frequencies in codon sequences; equal codon usage is assumed for amino acid sequences.
the suffix "-dgm" denotes "-dgmr", "-dgms" or "-dgmsf".
the suffixes "-dgmr" and "-dgms" mean the variation of mutation rate or selective constraint across sites, respectively, which is approximated by a discrete gamma distribution
<cit>  with m categories of unequal probabilities; see additional file
<dig> for details.
the "f" following "-dgms" means that the posterior frequencies of amino acids in each category in the first run are used in the second run as the equilibrium frequencies for each category; see the methods section.
substitution model
a
b
c
c
β
de
w
0
d
h
a"-f" means that the equilibrium frequencies are estimated to be equal to those in the alignment; equal codon usage is assumed.
"-dgmr" and "-dgms" mean discrete gamma distributions with m categories of unequal probabilities for the rate variation and the variation of selective constraint across sites, respectively.
"-dgmsf" means the equilibrium frequencies for respective categories are estimated from their posterior probabilities for sites.
the number string in the model name indicates the number of parameters optimized for the substitution rate matrix, and the remaining strings denote a rate matrix or a selective constraint matrix used.
bthe number of adjustable parameters.
cdifference from the reference state; Δℓ = ℓ +  <dig> , Δaic = aic -  <dig> , and Δbic = bic -  <dig> .
the reference tree topology is tree- <dig> in
<cit> .
dwab=min;
wabestimate is the one specified by the model name.
ethe value parenthesized means that the parameter is fixed at the value specified.
fthe average of
ewab over all amino acid pairs {a,b};
〈ewab〉≡1190∑a∑b>aewab.
gthe ratio of double to single and of triple to double nucleotide change exchangeability;
m^≡m^.
hthe ratio of mean transitional to mean transversional exchangeability;
m^tc|ag/m^.
ithe shape parameter of a discrete gamma distribution for the variation of mutation rate or selective constraint across sites.
substitution model
a
b
c
c
c
β
de
w
0
d
h
a"-f" means that the equilibrium frequencies are estimated to be equal to those in the alignment; equal codon usage is assumed.
"-dgmr" and "-dgms" mean discrete gamma distributions with m categories of unequal probabilities for the rate variation and the variation of selective constraint across sites, respectively.
"-dgmsf" means the equilibrium frequencies for respective categories are estimated from their posterior probabilities for sites.
the number string in the model name indicates the number of parameters optimized for the substitution rate matrix, and the remaining strings denote a rate matrix or a selective constraint matrix used.
bthe number of adjustable parameters.
cdifference from the reference state; Δℓ = ℓ+  <dig> , Δaic=aic-  <dig> , and Δbic=bic-  <dig> .
the reference tree topology is the one reported in
<cit> .
dwab=min;
wabestimate is the one specified by the model name.
ethe value parenthesized means that the parameter is fixed at the value specified.
fthe average of
ewab over all amino acid pairs {a,b};〈ewab〉≡1190∑a∑b>aewab.
gthe ratio of double to single and of triple to double nucleotide change exchangeability;
m^≡m^.
hthe ratio of mean transitional to mean transversional exchangeability;
m^tc|ag/m^.
ithe shape parameter of a discrete gamma distribution for the variation of mutation rate or selective constraint across sites.
substitution model
a
b
c
c
c
β
de
w
0
d
b"-f" means that the equilibrium frequencies are estimated to be equal to those in the alignment; equal codon usage is assumed.
"-dgmr" and "-dgms" mean discrete gamma distributions withmcategories of unequal probabilities for the rate variation and the variation of selective constraint across sites, respectively.
"-dgmsf" means the equilibrium frequencies for respective categories are estimated from their posterior probabilities for sites.
the number string in the model name indicates the number of parameters optimized for the substitution rate matrix, and the remaining strings denote a rate matrix or a selective constraint matrix used.
bthe number of adjustable parameters.
cdifference from the reference state; Δℓ = ℓ+ <dig> , Δaic=aic- <dig> , and Δbic = bic- <dig> .
the reference tree topology is one inferred by fasttree-2
<cit> .
dwab=min;
wabestimate is the one specified by the model name.
ethe value parenthesized means that the parameter is fixed at the value specified.
fthe average of
ewab over all amino acid pairs {a,b};
〈ewab〉≡1190∑a∑b>aewab.
gthe ratio of double to single and of triple to double nucleotide change exchangeability;
m^≡m^.
hthe ratio of mean transitional to mean transversional exchangeability;
m^tc|ag/m^.
ithe shape parameter of a discrete gamma distribution for the variation of mutation rate or selective constraint across sites.
the best models in the present amino acid substitution models for the respective datasets are cprev <dig> for cpprot- <dig>  mtrev for mammal-mtprot, and flu for ha_human-flu-a-h1n <dig>  this fact is expected because cprev <dig> was estimated
<cit>  from the full set of  <dig> protein-coding genes in the  <dig> chloroplast genomes including cpprot- <dig>  mtrev
<cit>  from the  <dig> mitochondrial proteins of  <dig> vertebrate species, and flu
<cit>  from ∼ <dig> influenza proteins including ha_human-flu-a-h1n <dig>  hence, the fact shown in the tables indicates that these matrices certainly represent substitution tendencies specific to the respective protein families.
on the other hand, cprev10
<cit>  was estimated from the smaller dataset than that for cprev <dig>  that is,  <dig> proteins in  <dig> chloroplast genomes including  <dig> land plants and the complete genome of cyanobacteria.
the size of database used for cprev <dig> may cause the cprev <dig> to perform worse than jtt for cpprot- <dig>
these tables also show that the log-likelihood values for jtt, wag, and lg are much smaller than those of the reference models and also differ largely from one another, indicating that these empirical substitution matrices represent the average tendencies of substitutions over various protein families but lack gene-level resolution.
what are the characteristics of substitutions specific to each protein family?
certainly the strength of selective constraint against amino acid substitutions depends on the type of protein, and varies across sites in a protein.
however, the dependence of selective constraint on the substituted type of amino acid may result primarily from amino acid properties.
hence, in the present mechanistic codon substitution model, we approximate the selective constraint  for a target protein family with a linear function of the particular value  that was estimated from an empirical amino acid substitution frequency matrix and represents the average strength of selective constraint against each type of amino acid substitution over various proteins.
in addition, the tendency of nucleotide mutation may differ among nuclear, mitochondrial and chloroplast dna, and selection on nucleotide substitutions at the dna/rna level may exist and depend on each gene.
in the present model, a nucleotide mutation rate matrix is fully parameterized by  <dig> parameters, and one additional parameter  is employed to represent the ratio of multiple to single nucleotide changes in a codon.
tables
<dig>
<dig>  and
<dig> clearly show that all the codon substitution models together with the respective selective constraint matrices estimated from jtt, lg, wag, and khg significantly outperform the best amino acid substitution model in all the three datasets, even though  <dig> more parameters must be estimated.
here we should notice that the best amino acid substitution rate matrices were estimated from the protein families corresponding to the respective target proteins.
in addition, it is important to notice that in agreement with common biological knowledge the mean transitional exchangeability is estimated to be far larger than the mean transversional exchangeability for all protein families in all codon substitution models.
as already claimed for codon sequences in
<cit> , the variation of selective constraint across sites  is a better model than the variation of codon mutation rate  in all the three datasets.
this fact confirms a common presumption that substitution rate variations across sites detected by amino acid substitution models are caused primarily by the variation of selective constraint rather than the mutation rate variation.
in the present analysis,  <dig> categories of unequal probabilities are employed to represent a Γ distribution.
this number of categories is chosen to be not sufficient but minimum to represent a Γ distribution.
actually, as shown in tables
<dig>
<dig>  and
<dig>  representing a Γ distribution by  <dig> categories of unequal probabilities can significantly improve the log-likelihood.
in tables
<dig>
<dig>  and
<dig>  the results for the equal constraint model are also listed; selective constraint is the same for all types of amino acid substitutions, that is, β =  <dig>  the aic and bic values of the equal constraint model are larger for the mammalian-mtprot but smaller for the cpprot and ha_human-flu-a-h1n <dig> than those of the reference amino acid substitution models.
accordingly, the estimated value of the slope β of a linear function for selective constraints, is larger for the mammalian-mtprot than for the cpprot and ha_human-flu-a-h1n <dig>  this fact indicates that the dependence  of selective constraint on amino acid type is less effective for both cpprot and ha_human-flu-a-h1n <dig>  although this may result from a property that both the datasets consist of relatively closely-related sequences and contain mostly conservative amino acid substitutions.
in the energy-increment-based  model shown in tables
<dig>
<dig>  and
<dig>  the selective constraint matrix used is the one estimated on the basis of the mean increment of contact energy and residue-volume change accompanied by an amino acid replacement
<cit> .
although the energy-increment-based selective constraints perform better than the equal constraint, it does not perform as well as the empirical selective constraints estimated from the jtt, wag, lg, and khg, indicating the good quality of their empirical values of the selective constraints.
in the present discrete gamma distribution model for the variation of selective constraint, the value of w <dig> differs among categories, but the same equilibrium frequencies of amino acids, which are estimated to be equal to those in the alignment, are employed for all the categories.
amino acid frequencies strongly depend on residue location in protein structures.
typically, non-polar residues are more and polar residues are less frequent in the interior of protein structures, where selective constraint against amino acid replacements tends to be more restrictive
<cit> , <cit> , <cit> .
there must be a correlation between the strength of selective constraint and the equilibrium frequencies of amino acids.
hence, different equilibrium frequencies should be employed in principle for each category.
substitution rate matrices that differ only in their equilibrium frequencies were employed in
<cit> , <cit> .
in
<cit> , <cit> , a different rate matrix was employed for each of gamma rate categories and biochemical and structural categories.
le et al.
<cit>  estimated from a very large alignment database and then tested four amino acid substitution rate matrices each of which corresponds to one discrete gamma category or one distribution-free category and has different exchangeabilities and equilibrium frequencies.
here, we estimate amino acid frequencies for each category from the posterior probabilities of sites being at each category and then parameters are optimized again with the estimated amino acid frequencies for each category.
the aic and bic values for this new model named with a suffix "-dg4sf" are also listed in tables
<dig>
<dig>  and
<dig>  the values of aic and bic are improved for mammal-mtprot and cpprot- <dig> but not for ha_human-flu-a-h1n <dig>  this model requiring additional  <dig> parameters needs a sufficient number of amino acid substitutions in an alignment to yield better values of aic and bic.
CONCLUSIONS
the greatest advantage of employing a mechanistic codon substitution model over amino acid substitution models resides in an intrinsic property that mutational tendencies of codons and selective constraints against amino acid changes can be separately formulated in codon substitution models.
as a result, codon mutational tendencies and the strength of selective constraint can be tailored to those in each gene with the additional  <dig> parameters, although a sufficient number of sequences and more calculations are needed for their estimation.
also, besides mutation rate variation, the variation of selective constraint across sites can be taken into account.
at the amino acid level, synonymous substitutions cannot be identified.
even so, taking account of synonymous substitutions with a codon substitution model improves its performance, as shown in the present analysis.
since mutations occur in nucleotide level, codon mutations may be well approximated by a markov process.
however, markovian properties are lost in the process of amino acid substitutions, because of redundancy of translation to amino acids.
a hidden markov model, in which codon types are hidden states changing in a markov process, is more appropriate to represent amino acid substitutions.
this may be one of reasons why the mechanistic codon substitution model outperforms any amino acid substitution model examined here.
a conclusion is drawn that the mechanistic codon substitution model is superior to amino acid substitution models even for protein sequences in evolutionary analysis.

 the reliability of phylogenetic analyses on protein-coding sequences strongly depends on models designed to approximate the substitution processes of nucleotide and amino acid. for the evolutionary analysis of protein-coding sequences, particularly phylogenetic inference, three types of substitution models can be used, provided that both dna and protein sequences are available; nucleotide, amino acid, and codon substitution models. for closely-related sequences in which most substitutions are synonymous, amino acid substitution models cannot be used, instead nucleotide substitution models may be employed. conversely, nucleotide substitution models should not be used in the case of diverged sequences in which synonymous substitutions are almost saturated and nonsynonymous substitutions are significant. on the other hand, codon substitution models are appropriate to both closely-related and highly-diverged sequences with the intrinsic property of detecting both synonymous and nonsynonymous substitutions. in a reversible markov model for substitution, a substitution rate matrix must be specified to estimate the likelihood of a phylogenetic tree. substitution tendencies between nucleotides and those between amino acids are characteristic of each species, each organelle, and each protein family. in the case of nucleotide substitution models, full parameterization for a substitution rate matrix is possible with parameters; the total rate is normalized to one. however, parameters for an amino acid substitution rate matrix are too many to be optimized for any size of a multiple sequence alignment. thus, empirical amino acid substitution rate matrices have been estimated from a large number of substitutions inferred on phylogenetic trees of single or many protein families; the jtt, the wag, and the lg matrices from nuclear proteins, mtrev from vertebrate mitochondrial proteins, cprev10 and cprev64 from chloroplast-encoded proteins, and flu from influenza proteins. even a codon substitution rate matrix, khg, has been estimated. these matrices significantly differ from one another, indicating that actually substitution tendencies significantly differ among these protein families. the estimation of a substitution rate matrix requires a large number or size of alignments with intensive calculation, and therefore is not always feasible. however, generic empirical substitution rate matrices such as jtt, wag, lg,and khg represent the average tendencies of substitutions over various protein families by sacrificing gene-level resolution. on the other hand, a rate matrix such as mtrev, cprev and flu derived from a specific protein family represents substitution tendencies characteristic of the protein but often lacks generic representation of substitution tendencies enough to be applied to other protein families. to resolve this situation, the parametrization of an amino acid substitution rate matrix has been attempted to easily generate an alignment-specific rate matrix. here, we propose a different approach of employing a mechanistic codon substitution model in which the biological and evolutionary mechanisms of amino acid substitutions are taken into account. in mechanistic codon substitution models, substitution rates are represented as the product of a codon mutation rate and the ratio of fixation, which depend on the types of codon replacement and amino acid replacement, respectively. hence, mutational tendencies at the nucleotide level and selection at the amino acid level can be taken into account at various levels of separation. this is a critical difference from amino acid substitution models. as a result, the variations of codon mutation rate and selective constraint across sites can be distinguished from each other. we formulated a codon substitution rates between codons μ and ν as rμν≡mμνewaμbν, with a codon mutation rate mμν, a mutation equilibrium frequency fνmut, an equilibrium frequency fν, and a selective constraint waμbν on the substitutions between amino acids aμ and bν. on the basis of this mechanistic codon substitution model, we estimated the wab by maximizing the likelihood of an empirical amino acid or codon substitution frequency matrix, each of jtt, wag, lg, and khg. it was shown that the mechanistic codon substitution model with a fully-parameterized codon mutation rate matrix and a selective constraint matrix, each element of which was approximated as a linear function,0]) of a given estimate of selective constraint estimated from the empirical substitution frequency matrix such as jtt and lg, outperformed both nucleotide and amino acid substitution models converted into codon substitution models in the wide range of codon sequences, from closely-related to highly-diverged protein-coding sequences. in this codon substitution model, codon mutations with multiple nucleotide changes are also taken into account, and were shown to increase significantly the likelihood of observed substitutions. there have been a variety of models for selective pressure on amino acid replacements in mechanistic codon substitution models. for their details, please see, in which the present mechanistic codon substitution model was discussed in comparison with other models. in these days, dna sequences are first analyzed and protein sequences are translated from them, and consequently codon sequences are likely available for most protein sequences in databases. however, there are protein sequences whose codon sequences are not available or not easily retrieved because of no cross link. in such a case, analysis at the amino acid level must be forced. here we show that even for amino acid sequences, among which synonymous substitutions cannot be identified, the present mechanistic codon substitution model with the assumption of equal codon usage outperforms amino acid substitution models using empirical substitution rate matrices.  on the basis of akaike information criterion and bayesian information criterion, the amino acid substitution models with the empirical amino acid substitution rate matrices, cprev64, cprev10, mtrev, and flu, as well as jtt, wag, and lg that were estimated from nuclear proteins, are compared with the mechanistic codon substitution models with the selective constraint matrices estimated from jtt, wag, lg, and khg by using the datasets: fast-evolving interspecific mitochondrial proteins concatenating protein-coding genes from mammalian species, closely-related chloroplast-encoded proteins concatenating protein-coding genes from chloroplast genomes of the major angiosperm lineages, and ha proteins of human influenza-a h1n consisting of sequences. the reference tree topologies used here as the most probable one are tree- of for mammalian-mtprot, the one reported in for cpprot- and the one inferred by the fasttree version 2 for ha_human-flu-a-h1n branch optimization of phylogenetic trees and their maximum log-likelihood values are calculated using phyml modified for the mechanistic codon substitution model. please see the methods section for the details of the mechanistic codon substitution model and the protein sequence data used. the naming conventions of models employed here are described in tables and the aic and bic values for these datasets are listed in tables and respectively. maximum log-likelihood, aic, and bic values are represented in relative to those of the reference model, which uses the empirical amino acid substitution rate matrix specific to each dataset and has the lowest aic and bic values in all the amino acid substitution models; the best amino acid substitution model is chosen here as a reference model in order to show how superior mechanistic codon substitution models are in comparison with the best amino acid substitution model. a random effects approach, in which the discrete gamma distribution with m categories of unequal probabilities is assumed for a prior probability distribution, is employed for rate variation across sites in the amino acid substitution models, and also for the variation of selective constraint or mutation rate across sites in the mechanistic codon substitution models; see additional file for details. suffixes "-dgmr" and "-dgms" are appended to model names to denote rate and selective constraint variations across sites, respectively. the equilibrium frequencies of amino acids are assumed to be equal to those in the aligned sequences, and equal codon usage is assumed to calculate the equilibrium codon frequencies from them; a suffix "-f" is appended to the model names., cprev64, and flu, wag, lg, mtrev, cprev10, cprev64, and flu the suffix "-dgmr" means that the variation of substitution rate is approximated by a discrete gamma distribution with m categories of unequal probabilities; see additional file for details. is used to estimate wab in eq. the Δε^abc and Δε^abv represent the effects of the mean increment of contact energy between residues and of residue-volume change due to an amino acid replacement, respectively; see supporting information, text s in {wabjtt/wag/lg-ml91+} estimated by maximizing the likelihood of jtt/wag/lg in the ml-91+ model are used as {wabkhg-ml200} estimated by maximizing the likelihood of the khg codon substitution matrix in the ml- model are used as the suffix "n" means the number of parameters optimized for the substitution rate matrix. the suffix "-f" means that equilibrium codon frequencies are assumed to be equal to codon frequencies in codon sequences; equal codon usage is assumed for amino acid sequences. the suffix "-dgm" denotes "-dgmr", "-dgms" or "-dgmsf". the suffixes "-dgmr" and "-dgms" mean the variation of mutation rate or selective constraint across sites, respectively, which is approximated by a discrete gamma distribution with m categories of unequal probabilities; see additional file for details. the "f" following "-dgms" means that the posterior frequencies of amino acids in each category in the first run are used in the second run as the equilibrium frequencies for each category; see the methods section. substitution model a b c c β de w 0 d h a"-f" means that the equilibrium frequencies are estimated to be equal to those in the alignment; equal codon usage is assumed. "-dgmr" and "-dgms" mean discrete gamma distributions with m categories of unequal probabilities for the rate variation and the variation of selective constraint across sites, respectively. "-dgmsf" means the equilibrium frequencies for respective categories are estimated from their posterior probabilities for sites. the number string in the model name indicates the number of parameters optimized for the substitution rate matrix, and the remaining strings denote a rate matrix or a selective constraint matrix used. bthe number of adjustable parameters. cdifference from the reference state; Δℓ = ℓ +, Δaic = aic -, and Δbic = bic -. the reference tree topology is tree- in. dwab=min; wabestimate is the one specified by the model name. ethe value parenthesized means that the parameter is fixed at the value specified. fthe average of ewab over all amino acid pairs {a,b}; 〈ewab〉≡1190∑a∑b>aewab. gthe ratio of double to single and of triple to double nucleotide change exchangeability; m^≡m^. hthe ratio of mean transitional to mean transversional exchangeability; m^tc|ag/m^. ithe shape parameter of a discrete gamma distribution for the variation of mutation rate or selective constraint across sites. substitution model a b c c c β de w 0 d h a"-f" means that the equilibrium frequencies are estimated to be equal to those in the alignment; equal codon usage is assumed. "-dgmr" and "-dgms" mean discrete gamma distributions with m categories of unequal probabilities for the rate variation and the variation of selective constraint across sites, respectively. "-dgmsf" means the equilibrium frequencies for respective categories are estimated from their posterior probabilities for sites. the number string in the model name indicates the number of parameters optimized for the substitution rate matrix, and the remaining strings denote a rate matrix or a selective constraint matrix used. bthe number of adjustable parameters. cdifference from the reference state; Δℓ = ℓ+, Δaic=aic-, and Δbic=bic-. the reference tree topology is the one reported in. dwab=min; wabestimate is the one specified by the model name. ethe value parenthesized means that the parameter is fixed at the value specified. fthe average of ewab over all amino acid pairs {a,b};〈ewab〉≡1190∑a∑b>aewab. gthe ratio of double to single and of triple to double nucleotide change exchangeability; m^≡m^. hthe ratio of mean transitional to mean transversional exchangeability; m^tc|ag/m^. ithe shape parameter of a discrete gamma distribution for the variation of mutation rate or selective constraint across sites. substitution model a b c c c β de w 0 d b"-f" means that the equilibrium frequencies are estimated to be equal to those in the alignment; equal codon usage is assumed. "-dgmr" and "-dgms" mean discrete gamma distributions withmcategories of unequal probabilities for the rate variation and the variation of selective constraint across sites, respectively. "-dgmsf" means the equilibrium frequencies for respective categories are estimated from their posterior probabilities for sites. the number string in the model name indicates the number of parameters optimized for the substitution rate matrix, and the remaining strings denote a rate matrix or a selective constraint matrix used. bthe number of adjustable parameters. cdifference from the reference state; Δℓ = ℓ+, Δaic=aic-, and Δbic = bic-. the reference tree topology is one inferred by fasttree-2. dwab=min; wabestimate is the one specified by the model name. ethe value parenthesized means that the parameter is fixed at the value specified. fthe average of ewab over all amino acid pairs {a,b}; 〈ewab〉≡1190∑a∑b>aewab. gthe ratio of double to single and of triple to double nucleotide change exchangeability; m^≡m^. hthe ratio of mean transitional to mean transversional exchangeability; m^tc|ag/m^. ithe shape parameter of a discrete gamma distribution for the variation of mutation rate or selective constraint across sites. the best models in the present amino acid substitution models for the respective datasets are cprev for cpprot- mtrev for mammal-mtprot, and flu for ha_human-flu-a-h1n this fact is expected because cprev was estimated from the full set of protein-coding genes in the chloroplast genomes including cpprot- mtrev from the mitochondrial proteins of vertebrate species, and flu from ∼ influenza proteins including ha_human-flu-a-h1n hence, the fact shown in the tables indicates that these matrices certainly represent substitution tendencies specific to the respective protein families. on the other hand, cprev10 was estimated from the smaller dataset than that for cprev that is, proteins in chloroplast genomes including land plants and the complete genome of cyanobacteria. the size of database used for cprev may cause the cprev to perform worse than jtt for cpprot- these tables also show that the log-likelihood values for jtt, wag, and lg are much smaller than those of the reference models and also differ largely from one another, indicating that these empirical substitution matrices represent the average tendencies of substitutions over various protein families but lack gene-level resolution. what are the characteristics of substitutions specific to each protein family? certainly the strength of selective constraint against amino acid substitutions depends on the type of protein, and varies across sites in a protein. however, the dependence of selective constraint on the substituted type of amino acid may result primarily from amino acid properties. hence, in the present mechanistic codon substitution model, we approximate the selective constraint for a target protein family with a linear function of the particular value that was estimated from an empirical amino acid substitution frequency matrix and represents the average strength of selective constraint against each type of amino acid substitution over various proteins. in addition, the tendency of nucleotide mutation may differ among nuclear, mitochondrial and chloroplast dna, and selection on nucleotide substitutions at the dna/rna level may exist and depend on each gene. in the present model, a nucleotide mutation rate matrix is fully parameterized by parameters, and one additional parameter is employed to represent the ratio of multiple to single nucleotide changes in a codon. tables and clearly show that all the codon substitution models together with the respective selective constraint matrices estimated from jtt, lg, wag, and khg significantly outperform the best amino acid substitution model in all the three datasets, even though more parameters must be estimated. here we should notice that the best amino acid substitution rate matrices were estimated from the protein families corresponding to the respective target proteins. in addition, it is important to notice that in agreement with common biological knowledge the mean transitional exchangeability is estimated to be far larger than the mean transversional exchangeability for all protein families in all codon substitution models. as already claimed for codon sequences in, the variation of selective constraint across sites is a better model than the variation of codon mutation rate in all the three datasets. this fact confirms a common presumption that substitution rate variations across sites detected by amino acid substitution models are caused primarily by the variation of selective constraint rather than the mutation rate variation. in the present analysis, categories of unequal probabilities are employed to represent a Γ distribution. this number of categories is chosen to be not sufficient but minimum to represent a Γ distribution. actually, as shown in tables and representing a Γ distribution by categories of unequal probabilities can significantly improve the log-likelihood. in tables and the results for the equal constraint model are also listed; selective constraint is the same for all types of amino acid substitutions, that is, β =  the aic and bic values of the equal constraint model are larger for the mammalian-mtprot but smaller for the cpprot and ha_human-flu-a-h1n than those of the reference amino acid substitution models. accordingly, the estimated value of the slope β of a linear function for selective constraints, is larger for the mammalian-mtprot than for the cpprot and ha_human-flu-a-h1n this fact indicates that the dependence of selective constraint on amino acid type is less effective for both cpprot and ha_human-flu-a-h1n although this may result from a property that both the datasets consist of relatively closely-related sequences and contain mostly conservative amino acid substitutions. in the energy-increment-based model shown in tables and the selective constraint matrix used is the one estimated on the basis of the mean increment of contact energy and residue-volume change accompanied by an amino acid replacement. although the energy-increment-based selective constraints perform better than the equal constraint, it does not perform as well as the empirical selective constraints estimated from the jtt, wag, lg, and khg, indicating the good quality of their empirical values of the selective constraints. in the present discrete gamma distribution model for the variation of selective constraint, the value of w differs among categories, but the same equilibrium frequencies of amino acids, which are estimated to be equal to those in the alignment, are employed for all the categories. amino acid frequencies strongly depend on residue location in protein structures. typically, non-polar residues are more and polar residues are less frequent in the interior of protein structures, where selective constraint against amino acid replacements tends to be more restrictive,,. there must be a correlation between the strength of selective constraint and the equilibrium frequencies of amino acids. hence, different equilibrium frequencies should be employed in principle for each category. substitution rate matrices that differ only in their equilibrium frequencies were employed in,. in,, a different rate matrix was employed for each of gamma rate categories and biochemical and structural categories. le et al. estimated from a very large alignment database and then tested four amino acid substitution rate matrices each of which corresponds to one discrete gamma category or one distribution-free category and has different exchangeabilities and equilibrium frequencies. here, we estimate amino acid frequencies for each category from the posterior probabilities of sites being at each category and then parameters are optimized again with the estimated amino acid frequencies for each category. the aic and bic values for this new model named with a suffix "-dg4sf" are also listed in tables and the values of aic and bic are improved for mammal-mtprot and cpprot- but not for ha_human-flu-a-h1n this model requiring additional parameters needs a sufficient number of amino acid substitutions in an alignment to yield better values of aic and bic.  the greatest advantage of employing a mechanistic codon substitution model over amino acid substitution models resides in an intrinsic property that mutational tendencies of codons and selective constraints against amino acid changes can be separately formulated in codon substitution models. as a result, codon mutational tendencies and the strength of selective constraint can be tailored to those in each gene with the additional parameters, although a sufficient number of sequences and more calculations are needed for their estimation. also, besides mutation rate variation, the variation of selective constraint across sites can be taken into account. at the amino acid level, synonymous substitutions cannot be identified. even so, taking account of synonymous substitutions with a codon substitution model improves its performance, as shown in the present analysis. since mutations occur in nucleotide level, codon mutations may be well approximated by a markov process. however, markovian properties are lost in the process of amino acid substitutions, because of redundancy of translation to amino acids. a hidden markov model, in which codon types are hidden states changing in a markov process, is more appropriate to represent amino acid substitutions. this may be one of reasons why the mechanistic codon substitution model outperforms any amino acid substitution model examined here. a conclusion is drawn that the mechanistic codon substitution model is superior to amino acid substitution models even for protein sequences in evolutionary analysis. 
it is shown that the mechanistic codon substitution model with the assumption of equal codon usage yields better values of akaike and bayesian information criteria for all three phylogenetic trees of mitochondrial, chloroplast, and influenza-a hemagglutinin proteins than the empirical amino acid substitution models with mtrev, cprev and flu, which were designed specifically for those protein families, respectively. the variation of selective constraint across sites fits the datasets significantly better than variable codon mutation rates, confirming that substitution rate variations across sites detected by amino acid substitution models are caused primarily by the variation of selective constraint against amino acid substitutions rather than the variation of codon mutation rate. the mechanistic codon substitution model is superior to amino acid substitution models even in the evolutionary analysis of protein sequences. on the other hand, in a mechanistic codon substitution model, in which each codon substitution rate is proportional to the product of a codon mutation rate and the ratio of fixation depending on the type of amino acid replacement, mutation rates and the strength of selective constraint on amino acids can be tailored to each protein family with additional parameters. selective constraints on the respective types of amino acid replacements are tailored to each gene with a linear function of a given estimate of selective constraints, which were estimated by maximizing the likelihood of an empirical amino acid or codon substitution frequency matrix, each of jtt, wag, lg, and khg. is it superior even for amino acid sequences, among which synonymous substitutions cannot be identified? amino acid substitution modelempirical amino acid substitution rate matrixmechanistic codon substitution modelstructural constraintsfunctional constraintsselective constraintsvariable selective constraint across sitesvariable mutation rate across sitesmultiple nucleotide change as a result, in the evolutionary analysis of codon sequences it outperforms codon substitution models equivalent to empirical amino acid substitution matrices. nucleotide mutations are assumed to occur independently of codon positions but multiple nucleotide changes in infinitesimal time are allowed. hence, various empirical amino acid substitution rate matrices have needed to be estimated for phylogenetic analysis: jtt, wag, and lg for nuclear proteins, mtrev for mitochondrial proteins, cprev and cprev for chloroplast-encoded proteins, and flu for influenza proteins. nucleotide and amino acid substitution tendencies are characteristic of each species, organelle, and protein family.
 nucleotide and amino acid substitution tendencies are characteristic of each species, organelle, and protein family. hence, various empirical amino acid substitution rate matrices have needed to be estimated for phylogenetic analysis: jtt, wag, and lg for nuclear proteins, mtrev for mitochondrial proteins, cprev and cprev for chloroplast-encoded proteins, and flu for influenza proteins. on the other hand, in a mechanistic codon substitution model, in which each codon substitution rate is proportional to the product of a codon mutation rate and the ratio of fixation depending on the type of amino acid replacement, mutation rates and the strength of selective constraint on amino acids can be tailored to each protein family with additional parameters. as a result, in the evolutionary analysis of codon sequences it outperforms codon substitution models equivalent to empirical amino acid substitution matrices. is it superior even for amino acid sequences, among which synonymous substitutions cannot be identified?  nucleotide mutations are assumed to occur independently of codon positions but multiple nucleotide changes in infinitesimal time are allowed. selective constraints on the respective types of amino acid replacements are tailored to each gene with a linear function of a given estimate of selective constraints, which were estimated by maximizing the likelihood of an empirical amino acid or codon substitution frequency matrix, each of jtt, wag, lg, and khg. it is shown that the mechanistic codon substitution model with the assumption of equal codon usage yields better values of akaike and bayesian information criteria for all three phylogenetic trees of mitochondrial, chloroplast, and influenza-a hemagglutinin proteins than the empirical amino acid substitution models with mtrev, cprev and flu, which were designed specifically for those protein families, respectively. the variation of selective constraint across sites fits the datasets significantly better than variable codon mutation rates, confirming that substitution rate variations across sites detected by amino acid substitution models are caused primarily by the variation of selective constraint against amino acid substitutions rather than the variation of codon mutation rate.  the mechanistic codon substitution model is superior to amino acid substitution models even in the evolutionary analysis of protein sequences. amino acid substitution modelempirical amino acid substitution rate matrixmechanistic codon substitution modelstructural constraintsfunctional constraintsselective constraintsvariable selective constraint across sitesvariable mutation rate across sitesmultiple nucleotide change 