BACKGROUND
with the expansion of genome projects, phylogenomics - the use of numerous genes to infer phylogenetic trees - is becoming a common way to resolve controversial relationships .
since large datasets increase the amount of phylogenetic signal included in the analysis, phylogenomics is less subject to stochastic errors than single gene phylogenies.
nevertheless, some nodes remain unresolved even at the genome-scale level  <cit> .
this can either be due to intrinsic properties of the data,  or to inadequate inference methods  <cit> .
in fact, systematic errors may be more pronounced in phylogenomics: in some cases, the gain in phylogenetic signal is masked by an increased level of systematic error, which can attain the same order of magnitude  <cit> .
in the worst case, this leads to erroneous phylogenies with a high statistical support  <cit> .
molecular sequence evolution exhibits a high complexity that is not fully accounted for in current models of sequence evolution.
since the first substitution model  <cit>  several simplifying assumptions have been relaxed; the evolutionary process is considered as heterogeneous  between character states,  over time and  along the alignment.
models that relaxed these assumptions  improve fit to the data and phylogenetic accuracy, and are therefore widely used.
more recently, probabilistic models have been developed to take into account heterogeneity of the qualitative aspect of amino-acid replacements along the alignment, by assigning sites to different classes of substitutional processes  <cit> .
of particular interest is the cat model  <cit> , a mixture model that infers categories from the data without any a priori biological assumptions and takes advantage of the dirichlet process prior  <cit>  to control the number of categories through a set of hyperparameters .
in this model, the substitution process is assumed to be site-independent and is entirely defined by the equilibrium frequencies of amino acids, while their exchangeabilities are assumed equal .
the equilibrium frequencies over the twenty amino acids constitute a good proxy to represent the functional constraints acting on each position during evolution.
in the following, we will call such categories substitution profiles, or simply profiles.
the number of profiles is often several hundreds, showing that previous models lack flexibility in handling heterogeneity of the substitution process across sites and, for large datasets, the cat model has a better fit to the data than standard models based on substitution matrices   <cit>  and renders phylogenetic inference less sensitive to long branch attraction artefact  <cit> .
despite that the inferred number of categories is large, it has been shown, using a posterior predictive approach, that the number of categories estimated by the cat model is conservative  <cit> .
the models that assume different evolutionary processes across sites consider that they are homogeneous over time.
however, except for a few constant positions for which the evolutionary constraints are uniformly strong, the functional constraints of most positions are likely to have changed over their evolutionary history  <cit> .
indeed, some amino acids are replaced at sites without a significant effect on function, but these changes might modify the environment of the protein, the intra- or inter-protein interactions, and so on, leading to changes of the selective pressure at other sites  <cit> .
this results in variation of the site-specific evolutionary rate across time  <cit> , a phenomenon called heterotachy  <cit> , which is recurrent in biological sequences  <cit> .
to handle this heterogeneity, models have been proposed that allow the substitution rate to vary over time  <cit> .
in this study, we will extend to the principle of heterotachy to the heterogeneity of the substitution process over time.
we call homopecilly  the hypothesis of an identical substitution process over time at a given site.
since functional constraints acting on proteins change over time, we however expect that not only the rate but also the amino acid substitution process may vary.
in particular, the subset of acceptable amino acids or the exchangeability matrix at a given site may change throughout evolutionary time.
we will therefore test the hypothesis of homopecilly by evaluating whether the nature of the substitution process varies significantly over time at a given site.
briefly, the substitution process will be characterized by the set of stationary frequencies of amino acids, as estimated by the cat model  <cit> .
several large datasets will be divided into monophyletic taxa to test the null hypothesis of a homogeneous substitution process, that is a site should be affiliated to the same cat category  in all predefined monophyletic groups.
we demonstrate not only that this null hypothesis is significantly rejected, but also that heteropecilly might generate phylogenetic artifacts.
RESULTS
evidence for a significant qualitative heterogeneity over time
to globally estimate the presence of heteropecilly, the frequency of different profiles , the frequency of positions that are stably affiliated to two different profiles in a pair of taxonomic groups, is computed .
for each comparison , only positions showing at least two substitutions are considered, as very slowly evolving positions do not contain a signal strong enough to provide stable affiliations .
for both datasets, most of the comparisons show high values of fdp: between 40% and 80% for the mitochondrial dataset , and between 24% and 48% for the nuc <dig> dataset .
in other words, about half of the stably affiliated positions are best described by two different profiles in two different taxonomic groups.
importantly, the distribution of fdps is clearly shifted to lower values in simulations under homopecilly than in real data, demonstrating that the observed heteropecilly is not due to stochastic variations.
the significance of the fdp statistics is limited by the fact that only 15-49% of the sites are considered.
to increase the number of sites stably affiliated, it would be necessary to increase the number of substitutions, hence the number of species in predefined monophyletic groups.
however, if heteropecilly is frequent, this will have the effect of increasing the probability that a site changed evolutionary properties, and hence cannot be stably affiliated to a single profile.
to eschew this dilemma, we devised another criterion, the probability of identical profile over n clades , to estimate the heteropecilly for all sites .
a pipn value of  <dig> indicates that the site is described by different cat profiles in at least two groups, whereas a value close to  <dig> indicates that the site is always affiliated to the same profiles.
the distribution of -ln for real and simulated datasets is displayed in figure  <dig>  as expected, we observe an excess of sites with high values of -ln in real data, whereas the number of sites that show a medium value of -ln is higher for the simulated data; results are highly significant .
moreover, since the pipn criterion is highly correlated with evolutionary rate , the lower evolutionary rate in nuclear genes, relative to mitochondrial genes, decreases the power of the pipn test, making heteropecilly less marked in figure 2b than in figure 2a.
altogether, both the fdp that considers all stably affiliated sites in a pairwise comparison and the pipn that considers every site separately, while comparing all taxonomic groups simultaneously, demonstrate that qualitative time-heterogeneity of the substitution process-which we refer to as heteropecilly-is widespread in real data.
heteropecilly might be due to the comparison of paralogous genes, with different functions.
however, the genes encoded in the mitochondrial genome of animals have an extremely high probability of being truly orthologous, since no protein encoding gene duplication is known in these organelles.
similarly, the orthology was carefully checked in the nuclear supermatrix using the protocol described in  <cit> ; in fact, the proteins considered, mainly ribosomal proteins, are not subject to frequent gene duplications in animals and more importantly the duplicates do not change its function.
as a result, it is unlikely that the observed heteropecilly is due to the comparison of paralogous genes.
relationship between heteropecilly and other biological properties
we first evaluated whether heteropecilly is related to the evolutionary rate .
sites with a pipn equal to  <dig> are generally fast-, or even very fast-, evolving.
for instance, more than five-sixths of such nuclear positions have accumulated over  <dig> substitutions, whereas only  <dig> % have undergone less than  <dig> substitutions.
the relationship between heteropecilly and evolutionary rate is highly significant .
this result is expected.
a slowly evolving position is under strong functional pressure that limits not only the number of substitutions, but also the diversity of acceptable amino acids, thus reducing the probability of being affiliated to different profiles in different clades.
moreover, the statistical power is reduced.
in contrast, a fast-evolving position has more opportunities to explore different types of selective constraints, and therefore accepts substitutions in different ways over time; this leads to a small value of pipn.
it is important to note that tables  <dig> and  <dig> suggest that many fast evolving positions are indeed under rather strong selective pressure.
if they were completely free to vary, these fast evolving sites would explore all the twenty amino acids.
in such cases, substitution profiles would be the same in different clades, characterized by a small equilibrium frequency for many amino acids, and these sites would be associated with a pipn far from  <dig>  in contrast, a small pipn indicates that negative selection in a given clade is strong , but that this selection pattern changes over time; for instance, a position might accept only asp and glu in one clade, and only asp and asn in another.
a relationship between heteropecilly and rate heterogeneity over time  is possible, since both heterogeneities are due to changes in functional constraints.
accordingly, tables  <dig> and  <dig>  reveal a link of qualitative heterogeneity with heterotachy .
for instance, 80% of mitochondrial positions with a pipn of  <dig> are very heterotachous , whereas only 65% of all positions are heterotachous at that level.
nevertheless, heterotachy and heteropecilly appear anti-correlated: the proportion of heteropecillous sites is greater when sites are less heterotachous for the nuclear dataset, whereas proportion of heteropecillous sites increases with heterotachy for the mitochondrial dataset.
since the species number  is reduced and some missing data are present in the nuclear alignment, estimation of heterotachy and of heteropecilly is less accurate than in the mitochondrial alignment.
more importantly, the relationship with heterotachy is not only less marked than with the evolutionary rate, but it is in fact mainly a consequence of the correlation between heterotachy and evolutionary rate : when only the fast evolving positions  are considered, heteropecilly is almost unrelated to heterotachy, especially for the mitochondrial dataset .
further studies are therefore needed to determine whether a change of rate over time is independent, or not, of the fact that a site may be affiliated to different amino acid profiles over time.
the absence of a strong correlation between heteropecilly and heterotachy makes sense because the two heterogeneities do not apply on the same criteria: heterotachy is mainly due to loss or gain of functional constraints , whereas heteropecilly is rather due to variation in the nature of functional constraints, not in their strength.
we analyzed the correlation between heteropecilly and change in hydropathy, using the standard deviation of the site-wise profile hydrophobic score  index, which measures the diversity of hydropathy a site displays in various clades.
tables  <dig> and  <dig>  demonstrate a highly significant relationship .
if we make the assumption that a change in substitution profile reflects a variation in functional constraints, the correlation between the pipn criterion and the phs score suggests that, when functional constraints change over time, the new spectrum of acceptable amino acids becomes biochemically different from the previous one.
nevertheless, changes from one profile to another are much easier to detect with our protocol when no acceptable amino acids are common to both than when at least one amino acid remains in the set of acceptable residues.
in the first case, sites will surely be affiliated to different profiles in different clades, even if the signal is weak, whereas in the second case the sites will be likely affiliated to several, overlapping, profiles showing closer hydrophobic scores.
the use of larger datasets  is required to address this issue.
finally, we estimated whether changes between profiles showing similar physico-chemical properties  are more frequent than changes between profiles with different properties.
only sites that are stably affiliated in two different clades are considered.
about half of sites  are affiliated to profiles with different biochemical properties.
this suggests that heteropecilly is driven not only by different fine-grained functional constraints, but also by important functional changes.
analyses at the codon level  <cit>  are nevertheless required to estimate whether heteropecilly is driven by positive selection or is the result of changes in purifying selection properties.
phylogenetic structure of heteropecilly
we have shown that the substitution process, as characterized by the cat profiles, varies over time.
although the way profile affiliations change over time is not known, it is reasonable to assume that profile affiliations are often inherited from ancestors.
in other words, for orthologs, two sister clades would generally have the same substitution process, i.e.
a given site would generally be affiliated to the same profile in two sister clades.
therefore, a phylogenetic signal is expected to exist in the profile distribution.
we tested this hypothesis by recoding the multiple amino acid sequences of a clade into a single artificial sequence.
in these new sequences, the state of a given site is the profile to which this site is stably affiliated or a question mark otherwise.
only the  <dig> most frequent profiles are considered  and trees are inferred using the gtr+Γ <dig>  and the cat+Γ <dig>  models.
with the mt <dig> dataset, most of the known clades are recovered with the gtr+Γ <dig> model , many of them with high posterior probabilities , e.g.
bilateria, pancrustacea, vertebrata, amphibia, sauria, archosauria, lepidosauria, theria and eutheria.
the only problematic case is the relative position of amphibia, mammalia and sauria.
the mono/paraphyly of deuterostomia is known to be a difficult question, even using a hundred nuclear genes  <cit> .
although we don't know which statistical model has to be used for analyzing the artificial recoded sequences, recovering at least some clades is reassuring.
to test the possibility of a correct grouping of clades by chance, we applied the recoding protocol to the ten datasets simulated under a homopecillous model : except in one case, the expected monophyletic groups are not recovered or only recovered with very small posterior probabilities .
similar results were obtained for the nuc <dig> dataset , but are less clear because only five clades are available.
the congruence of the phylogeny inferred from recoded data with current taxonomy demonstrates the presence of a strong phylogenetic signal in changes of the evolutionary process, especially since only few parsimony informative positions are available .
changes in the evolutionary process are therefore relatively rare since sites remain affiliated to the same profile over hundred millions of years.
the same observation has been made for evolutionary rates, leading to local molecular clock  <cit>  or auto-correlated relaxed molecular clock models  <cit> .
nevertheless, the incongruence observed in a few cases indicates that homoplasy is present in recoded data and is not correctly handled by the cat+Γ model; not surprisingly, when a profile affiliation changes, it can either revert to the ancestral state or converge toward a state independently acquired in a distantly related clade.
the misleading effect of homoplasy is enhanced by the very heterogeneous rate of change of profile affiliations: for the mitochondrial dataset, the rate is high on the branch leading to bilateria  and for the nuc <dig> dataset, nematodes and platyhelminthes evolved several times faster than the other bilaterians .
to further characterize homoplasy in the distribution of profile affiliation, we looked at the distribution of profiles across clades.
if the process of change in the substitution process is stationary, one expects a homogeneous distribution.
the excess or lack of profile affiliation within a clade is estimated with respect to the average of the distribution among clades under consideration.
for the nuclear dataset, the profile distribution is plotted for the four clades of interest; deuterostomia -used as outgroup- are not shown .
profiles are not equally distributed among clades: some profiles are in excess for one or two clades .
this unequal distribution of profiles across clades could be studied separately within sub-groups of profiles according to their physico-chemical properties.
three sub-groups  show a large variation, whereas other sub-groups  are more homogeneous.
interestingly, fast evolving platyhelminthes are the most divergent group, and in the vast majority of cases do not have the same bias as other lophotrochozoa .
for instance, for the small and uncharged amino acid sub-group, they are associated preferentially with ags or as profiles, whereas mollusks and annelids show an affiliation preference for st or g profiles.
this bias in affiliation frequency across clades generates homoplasy that is difficult to handle, and could explain why platyhelminthes are so difficult to position .
importantly, this heterogeneity is not found in sequences simulated under the cat model, which is homogenous over time , and is significant according to a χ <dig> test .
the mt <dig> dataset yields similar results .
since the amino acid composition is influenced by the genomic nucleotide composition  <cit> , heteropecilly could be the result of the heterogeneity of the mutation process over time, i.e.
the frequency of some profiles will increase in a clade because their most frequent amino acids are becoming more frequent due to changes in the nucleotide bias.
this hypothesis predicts that profiles will be heterogeneously distributed across clades, especially for the fast evolving positions, which are most likely to reflect mutational bias.
the profiles are indeed heterogeneously distributed across clades for both mitochondrial and nuclear datasets  but, when the fastest evolving positions are considered, this remains highly significant for the mitochondrial alignment only.
for this alignment, when one compares an equal number of the most and the least heteropecillous positions, the p-value is slightly lower for the former than for the latter, even if the profiles are significantly heterogeneously distributed across clades.
this suggests that changes in the mutational process over time, albeit particularly marked in the mitochondrial genome, are likely a minor cause of heteropecilly.
in summary, the time heterogeneity of the amino acid evolutionary process  is a general phenomenon in animal evolution, present in both nuclear and mitochondrial coding genomes.
although the rate of change of the substitution process is sufficiently low to allow the recovery of a phylogenetic signal in the recoded sequences, homoplasy is present as suggested by the different rate of evolution across the tree , and by the heterogeneity of profile frequencies across clades .
we do not recommend using the recoding protocol to avoid model violation in cases where heteropecilly is suspected; even if some phylogenetic signal can be captured in the recoded sequences, the signal is probably too weak to obtain an accurate phylogeny.
heteropecilly may generate phylogenetic reconstruction artefacts
the observations presented so far demonstrate that the assumption of homopecilly, i.e.
no change of the substitution process over time, is violated, potentially leading to tree reconstruction artifacts.
as a case study, we chose the relationships between porifera, cnidaria and bilateria , which are difficult to resolve with mitochondrial genomes  <cit> .
since slow-evolving porifera and cnidaria are grouped together  <cit> , the monophyly of eumetazoa , long proposed by morphologists and recovered with nuclear genes  <cit> , is not observed for the mitochondrial data.
we make the assumption that this is due to a long-branch attraction  artifact between the fast-evolving bilateria and the distant outgroup .
this could be aggravated when using the cat model by the very long branch observed at the base of bilateria in recoded sequences , which indicates a large amount of profile affiliation changes, and hence a serious violation of one hypothesis of the cat model .
we analyzed a mitochondrial encoded dataset with  <dig> species .
when using the complete dataset , the cat+Γ <dig> model groups together cnidaria and porifera with a posterior probability of  <dig> .
we then progressively removed heteropecillous positions according to their increasing pipn value; that is, we first removed positions that most likely violate the assumption of homogeneity over time of the cat model.
for the five sub-datasets analyzed with the cat+Γ <dig> model, three kinds of topologies are observed :  with  <dig> positions, the same topology as with the complete dataset;  with  <dig> and  <dig> positions, eumetazoa  are recovered, but porifera are not monophyletic ;  with the smallest subsets, eumetazoa and porifera are both monophyletic.
with the removal of sites with heterogeneous profile affiliations across clades, support for the monophyly of eumetazoa  increases steadily, up to about one, and decreases to  <dig>  in the smallest dataset, probably due to the limited amount of data.
this result suggests that sites showing a substitutional heterogeneity of profiles over time interfere with the phylogenetic signal and may eventually result in an erroneous topology.
we computed the number of sites affiliated to the same profile in cnidaria and porifera  and in cnidaria and bilateria .
choanoflagellates are not considered because only two species are available and an accurate profile affiliation to sites is not possible.
figure  <dig> shows that, upon removal of heteropecillous sites, according to the pipn criterion, the number of sites having the same profile in porifera and cnidaria decreases much more rapidly than the number of sites having the same profile in cnidaria and bilateria.
one can reasonably argue that sites with the same profile in porifera and cnidaria generate a spurious signal that is erroneously interpreted by the cat+Γ <dig> model as synapomorphies for porifera+cnidaria, since this model assumes that profiles are identical over the whole tree.
two controls were performed.
first, since a correlation exists between evolutionary rate and heteropecilly , improvements in phylogenetic inference could be due to the removal of fast evolving sites  <cit> .
when fast evolving sites are progressively removed according to the sf method  <cit> , support for the incorrect grouping of cnidaria and porifera slightly increases .
this is in sharp contrast with the removal of heteropecillous positions .
second, since the negative effect of heteropecilly may constitute a model violation more important for the cat+Γ <dig> model than for site-homogeneous models , the gtr+Γ <dig> or mtrev+Γ <dig> models are also used.
results are completely different: the removal of heteropecillous positions does not affect phylogenetic inferences conducted with gtr+Γ <dig> or mtrev+Γ <dig> models, for which high support for cnidaria + porifera is always observed .
this result is not unexpected because, while heteropecilly constitutes a model violation for the cat model, its effect on the site-homogeneous gtr model is less clear.
to understand the different behavior of these two models, one has to distinguish two model violations, heterogeneity of the substitution process across sites and over time.
the first violation is known to seriously exacerbate lba artifacts, because the amount of homoplasy is underestimated  <cit> .
it is not expected to decrease with the removal of heteropecillous positions and probably dominates when using the site-homogeneous gtr and mtrev models.
the effect of the second model violation can only be observed with the cat model, which handles heterogeneity across sites.
removing heteropecillous positions will reduce the time-homogeneous violation and increase the accuracy of the cat model.
this hypothesis is corroborated by an evaluation of model fit by cross-validation .
the gtr+Γ <dig> model fits the complete mitochondrial mt <dig> dataset better than the cat+Γ <dig> .
this is probably due to the large number of parameters of the latter model and the limited amount of positions, since, with the larger nuclear dataset, the cat+Γ <dig> model has a better fit than the two models based on exchange matrices .
importantly, after removal of most heteropecillous sites, the cat+Γ <dig> model appears to best fit the data .
this indicates that heteropecilly constitutes a serious violation of the cat model, because the best fit is obtained despite the limited number of sites .
CONCLUSIONS
numerous heterogeneities of the evolutionary process have been discovered.
most have a clear negative impact on phylogenetic inference when not adequately handled: heterogeneity of rate across lineages  <cit> , of substitution type  <cit> , of rate across sites  <cit> , of composition across taxa  <cit>  and of substitution process across sites  <cit> .
some of them, such as heterotachy, seem to have a more limited effect on phylogenetic accuracy  <cit> .
further studies are needed to know in which category heteropecilly has to be classified.
an important prerequisite is the development of models that handle the fact that substitution properties change over time.
this could be achieved via a markov modulated cat model, similar to the covarion model  <cit> , or via the use of breakpoints  <cit> .
to choose between these two approaches, it would be important to estimate whether sites generally change their properties in a collective manner or not, since only the second approach can model this aspect.
the costs of these improved models  could be major and would be useful for phylogenetic inference only if heteropecilly turned out to seriously impair accuracy.
in any case such improved models would be helpful to advance our knowledge of protein evolution, since in most cases one can select a set of species for whose relationships are confidently known, which drastically simplifies the problem and allows the use of complex but computationally demanding models.
what are the main reasons for these shifts in profile affiliation?
the correlation between the pipn criterion and the substitution number might suggest a neutral explanation  of the variation in substitutional profile over time.
since the pipn criterion is also correlated with a change in hydropathy, this neutral explanation could be insufficient.
these changes in the site-specific substitution process could be related to functional shifts, such as adaptation of organisms to new environments  or of the protein to a new cellular environment .
to answer these questions, it would be particularly relevant to study paralogous genes in which functional shifts are well characterized.

 with the expansion of genome projects, phylogenomics - the use of numerous genes to infer phylogenetic trees - is becoming a common way to resolve controversial relationships. since large datasets increase the amount of phylogenetic signal included in the analysis, phylogenomics is less subject to stochastic errors than single gene phylogenies. nevertheless, some nodes remain unresolved even at the genome-scale level. this can either be due to intrinsic properties of the data, or to inadequate inference methods. in fact, systematic errors may be more pronounced in phylogenomics: in some cases, the gain in phylogenetic signal is masked by an increased level of systematic error, which can attain the same order of magnitude. in the worst case, this leads to erroneous phylogenies with a high statistical support. molecular sequence evolution exhibits a high complexity that is not fully accounted for in current models of sequence evolution. since the first substitution model several simplifying assumptions have been relaxed; the evolutionary process is considered as heterogeneous between character states, over time and along the alignment. models that relaxed these assumptions improve fit to the data and phylogenetic accuracy, and are therefore widely used. more recently, probabilistic models have been developed to take into account heterogeneity of the qualitative aspect of amino-acid replacements along the alignment, by assigning sites to different classes of substitutional processes. of particular interest is the cat model, a mixture model that infers categories from the data without any a priori biological assumptions and takes advantage of the dirichlet process prior to control the number of categories through a set of hyperparameters. in this model, the substitution process is assumed to be site-independent and is entirely defined by the equilibrium frequencies of amino acids, while their exchangeabilities are assumed equal. the equilibrium frequencies over the twenty amino acids constitute a good proxy to represent the functional constraints acting on each position during evolution. in the following, we will call such categories substitution profiles, or simply profiles. the number of profiles is often several hundreds, showing that previous models lack flexibility in handling heterogeneity of the substitution process across sites and, for large datasets, the cat model has a better fit to the data than standard models based on substitution matrices and renders phylogenetic inference less sensitive to long branch attraction artefact. despite that the inferred number of categories is large, it has been shown, using a posterior predictive approach, that the number of categories estimated by the cat model is conservative. the models that assume different evolutionary processes across sites consider that they are homogeneous over time. however, except for a few constant positions for which the evolutionary constraints are uniformly strong, the functional constraints of most positions are likely to have changed over their evolutionary history. indeed, some amino acids are replaced at sites without a significant effect on function, but these changes might modify the environment of the protein, the intra- or inter-protein interactions, and so on, leading to changes of the selective pressure at other sites. this results in variation of the site-specific evolutionary rate across time, a phenomenon called heterotachy, which is recurrent in biological sequences. to handle this heterogeneity, models have been proposed that allow the substitution rate to vary over time. in this study, we will extend to the principle of heterotachy to the heterogeneity of the substitution process over time. we call homopecilly the hypothesis of an identical substitution process over time at a given site. since functional constraints acting on proteins change over time, we however expect that not only the rate but also the amino acid substitution process may vary. in particular, the subset of acceptable amino acids or the exchangeability matrix at a given site may change throughout evolutionary time. we will therefore test the hypothesis of homopecilly by evaluating whether the nature of the substitution process varies significantly over time at a given site. briefly, the substitution process will be characterized by the set of stationary frequencies of amino acids, as estimated by the cat model. several large datasets will be divided into monophyletic taxa to test the null hypothesis of a homogeneous substitution process, that is a site should be affiliated to the same cat category in all predefined monophyletic groups. we demonstrate not only that this null hypothesis is significantly rejected, but also that heteropecilly might generate phylogenetic artifacts.  evidence for a significant qualitative heterogeneity over time to globally estimate the presence of heteropecilly, the frequency of different profiles, the frequency of positions that are stably affiliated to two different profiles in a pair of taxonomic groups, is computed. for each comparison, only positions showing at least two substitutions are considered, as very slowly evolving positions do not contain a signal strong enough to provide stable affiliations. for both datasets, most of the comparisons show high values of fdp: between 40% and 80% for the mitochondrial dataset, and between 24% and 48% for the nuc dataset. in other words, about half of the stably affiliated positions are best described by two different profiles in two different taxonomic groups. importantly, the distribution of fdps is clearly shifted to lower values in simulations under homopecilly than in real data, demonstrating that the observed heteropecilly is not due to stochastic variations. the significance of the fdp statistics is limited by the fact that only 15-49% of the sites are considered. to increase the number of sites stably affiliated, it would be necessary to increase the number of substitutions, hence the number of species in predefined monophyletic groups. however, if heteropecilly is frequent, this will have the effect of increasing the probability that a site changed evolutionary properties, and hence cannot be stably affiliated to a single profile. to eschew this dilemma, we devised another criterion, the probability of identical profile over n clades, to estimate the heteropecilly for all sites. a pipn value of indicates that the site is described by different cat profiles in at least two groups, whereas a value close to indicates that the site is always affiliated to the same profiles. the distribution of -ln for real and simulated datasets is displayed in figure as expected, we observe an excess of sites with high values of -ln in real data, whereas the number of sites that show a medium value of -ln is higher for the simulated data; results are highly significant. moreover, since the pipn criterion is highly correlated with evolutionary rate, the lower evolutionary rate in nuclear genes, relative to mitochondrial genes, decreases the power of the pipn test, making heteropecilly less marked in figure 2b than in figure 2a. altogether, both the fdp that considers all stably affiliated sites in a pairwise comparison and the pipn that considers every site separately, while comparing all taxonomic groups simultaneously, demonstrate that qualitative time-heterogeneity of the substitution process-which we refer to as heteropecilly-is widespread in real data. heteropecilly might be due to the comparison of paralogous genes, with different functions. however, the genes encoded in the mitochondrial genome of animals have an extremely high probability of being truly orthologous, since no protein encoding gene duplication is known in these organelles. similarly, the orthology was carefully checked in the nuclear supermatrix using the protocol described in; in fact, the proteins considered, mainly ribosomal proteins, are not subject to frequent gene duplications in animals and more importantly the duplicates do not change its function. as a result, it is unlikely that the observed heteropecilly is due to the comparison of paralogous genes. relationship between heteropecilly and other biological properties we first evaluated whether heteropecilly is related to the evolutionary rate. sites with a pipn equal to are generally fast-, or even very fast-, evolving. for instance, more than five-sixths of such nuclear positions have accumulated over substitutions, whereas only % have undergone less than substitutions. the relationship between heteropecilly and evolutionary rate is highly significant. this result is expected. a slowly evolving position is under strong functional pressure that limits not only the number of substitutions, but also the diversity of acceptable amino acids, thus reducing the probability of being affiliated to different profiles in different clades. moreover, the statistical power is reduced. in contrast, a fast-evolving position has more opportunities to explore different types of selective constraints, and therefore accepts substitutions in different ways over time; this leads to a small value of pipn. it is important to note that tables and suggest that many fast evolving positions are indeed under rather strong selective pressure. if they were completely free to vary, these fast evolving sites would explore all the twenty amino acids. in such cases, substitution profiles would be the same in different clades, characterized by a small equilibrium frequency for many amino acids, and these sites would be associated with a pipn far from in contrast, a small pipn indicates that negative selection in a given clade is strong, but that this selection pattern changes over time; for instance, a position might accept only asp and glu in one clade, and only asp and asn in another. a relationship between heteropecilly and rate heterogeneity over time is possible, since both heterogeneities are due to changes in functional constraints. accordingly, tables and reveal a link of qualitative heterogeneity with heterotachy. for instance, 80% of mitochondrial positions with a pipn of are very heterotachous, whereas only 65% of all positions are heterotachous at that level. nevertheless, heterotachy and heteropecilly appear anti-correlated: the proportion of heteropecillous sites is greater when sites are less heterotachous for the nuclear dataset, whereas proportion of heteropecillous sites increases with heterotachy for the mitochondrial dataset. since the species number is reduced and some missing data are present in the nuclear alignment, estimation of heterotachy and of heteropecilly is less accurate than in the mitochondrial alignment. more importantly, the relationship with heterotachy is not only less marked than with the evolutionary rate, but it is in fact mainly a consequence of the correlation between heterotachy and evolutionary rate: when only the fast evolving positions are considered, heteropecilly is almost unrelated to heterotachy, especially for the mitochondrial dataset. further studies are therefore needed to determine whether a change of rate over time is independent, or not, of the fact that a site may be affiliated to different amino acid profiles over time. the absence of a strong correlation between heteropecilly and heterotachy makes sense because the two heterogeneities do not apply on the same criteria: heterotachy is mainly due to loss or gain of functional constraints, whereas heteropecilly is rather due to variation in the nature of functional constraints, not in their strength. we analyzed the correlation between heteropecilly and change in hydropathy, using the standard deviation of the site-wise profile hydrophobic score index, which measures the diversity of hydropathy a site displays in various clades. tables and demonstrate a highly significant relationship. if we make the assumption that a change in substitution profile reflects a variation in functional constraints, the correlation between the pipn criterion and the phs score suggests that, when functional constraints change over time, the new spectrum of acceptable amino acids becomes biochemically different from the previous one. nevertheless, changes from one profile to another are much easier to detect with our protocol when no acceptable amino acids are common to both than when at least one amino acid remains in the set of acceptable residues. in the first case, sites will surely be affiliated to different profiles in different clades, even if the signal is weak, whereas in the second case the sites will be likely affiliated to several, overlapping, profiles showing closer hydrophobic scores. the use of larger datasets is required to address this issue. finally, we estimated whether changes between profiles showing similar physico-chemical properties are more frequent than changes between profiles with different properties. only sites that are stably affiliated in two different clades are considered. about half of sites are affiliated to profiles with different biochemical properties. this suggests that heteropecilly is driven not only by different fine-grained functional constraints, but also by important functional changes. analyses at the codon level are nevertheless required to estimate whether heteropecilly is driven by positive selection or is the result of changes in purifying selection properties. phylogenetic structure of heteropecilly we have shown that the substitution process, as characterized by the cat profiles, varies over time. although the way profile affiliations change over time is not known, it is reasonable to assume that profile affiliations are often inherited from ancestors. in other words, for orthologs, two sister clades would generally have the same substitution process, i.e. a given site would generally be affiliated to the same profile in two sister clades. therefore, a phylogenetic signal is expected to exist in the profile distribution. we tested this hypothesis by recoding the multiple amino acid sequences of a clade into a single artificial sequence. in these new sequences, the state of a given site is the profile to which this site is stably affiliated or a question mark otherwise. only the most frequent profiles are considered and trees are inferred using the gtr+Γ and the cat+Γ models. with the mt dataset, most of the known clades are recovered with the gtr+Γ model, many of them with high posterior probabilities, e.g. bilateria, pancrustacea, vertebrata, amphibia, sauria, archosauria, lepidosauria, theria and eutheria. the only problematic case is the relative position of amphibia, mammalia and sauria. the mono/paraphyly of deuterostomia is known to be a difficult question, even using a hundred nuclear genes. although we don't know which statistical model has to be used for analyzing the artificial recoded sequences, recovering at least some clades is reassuring. to test the possibility of a correct grouping of clades by chance, we applied the recoding protocol to the ten datasets simulated under a homopecillous model: except in one case, the expected monophyletic groups are not recovered or only recovered with very small posterior probabilities. similar results were obtained for the nuc dataset, but are less clear because only five clades are available. the congruence of the phylogeny inferred from recoded data with current taxonomy demonstrates the presence of a strong phylogenetic signal in changes of the evolutionary process, especially since only few parsimony informative positions are available. changes in the evolutionary process are therefore relatively rare since sites remain affiliated to the same profile over hundred millions of years. the same observation has been made for evolutionary rates, leading to local molecular clock or auto-correlated relaxed molecular clock models. nevertheless, the incongruence observed in a few cases indicates that homoplasy is present in recoded data and is not correctly handled by the cat+Γ model; not surprisingly, when a profile affiliation changes, it can either revert to the ancestral state or converge toward a state independently acquired in a distantly related clade. the misleading effect of homoplasy is enhanced by the very heterogeneous rate of change of profile affiliations: for the mitochondrial dataset, the rate is high on the branch leading to bilateria and for the nuc dataset, nematodes and platyhelminthes evolved several times faster than the other bilaterians. to further characterize homoplasy in the distribution of profile affiliation, we looked at the distribution of profiles across clades. if the process of change in the substitution process is stationary, one expects a homogeneous distribution. the excess or lack of profile affiliation within a clade is estimated with respect to the average of the distribution among clades under consideration. for the nuclear dataset, the profile distribution is plotted for the four clades of interest; deuterostomia -used as outgroup- are not shown. profiles are not equally distributed among clades: some profiles are in excess for one or two clades. this unequal distribution of profiles across clades could be studied separately within sub-groups of profiles according to their physico-chemical properties. three sub-groups show a large variation, whereas other sub-groups are more homogeneous. interestingly, fast evolving platyhelminthes are the most divergent group, and in the vast majority of cases do not have the same bias as other lophotrochozoa. for instance, for the small and uncharged amino acid sub-group, they are associated preferentially with ags or as profiles, whereas mollusks and annelids show an affiliation preference for st or g profiles. this bias in affiliation frequency across clades generates homoplasy that is difficult to handle, and could explain why platyhelminthes are so difficult to position. importantly, this heterogeneity is not found in sequences simulated under the cat model, which is homogenous over time, and is significant according to a χ test. the mt dataset yields similar results. since the amino acid composition is influenced by the genomic nucleotide composition, heteropecilly could be the result of the heterogeneity of the mutation process over time, i.e. the frequency of some profiles will increase in a clade because their most frequent amino acids are becoming more frequent due to changes in the nucleotide bias. this hypothesis predicts that profiles will be heterogeneously distributed across clades, especially for the fast evolving positions, which are most likely to reflect mutational bias. the profiles are indeed heterogeneously distributed across clades for both mitochondrial and nuclear datasets but, when the fastest evolving positions are considered, this remains highly significant for the mitochondrial alignment only. for this alignment, when one compares an equal number of the most and the least heteropecillous positions, the p-value is slightly lower for the former than for the latter, even if the profiles are significantly heterogeneously distributed across clades. this suggests that changes in the mutational process over time, albeit particularly marked in the mitochondrial genome, are likely a minor cause of heteropecilly. in summary, the time heterogeneity of the amino acid evolutionary process is a general phenomenon in animal evolution, present in both nuclear and mitochondrial coding genomes. although the rate of change of the substitution process is sufficiently low to allow the recovery of a phylogenetic signal in the recoded sequences, homoplasy is present as suggested by the different rate of evolution across the tree, and by the heterogeneity of profile frequencies across clades. we do not recommend using the recoding protocol to avoid model violation in cases where heteropecilly is suspected; even if some phylogenetic signal can be captured in the recoded sequences, the signal is probably too weak to obtain an accurate phylogeny. heteropecilly may generate phylogenetic reconstruction artefacts the observations presented so far demonstrate that the assumption of homopecilly, i.e. no change of the substitution process over time, is violated, potentially leading to tree reconstruction artifacts. as a case study, we chose the relationships between porifera, cnidaria and bilateria, which are difficult to resolve with mitochondrial genomes. since slow-evolving porifera and cnidaria are grouped together, the monophyly of eumetazoa, long proposed by morphologists and recovered with nuclear genes, is not observed for the mitochondrial data. we make the assumption that this is due to a long-branch attraction artifact between the fast-evolving bilateria and the distant outgroup. this could be aggravated when using the cat model by the very long branch observed at the base of bilateria in recoded sequences, which indicates a large amount of profile affiliation changes, and hence a serious violation of one hypothesis of the cat model. we analyzed a mitochondrial encoded dataset with species. when using the complete dataset, the cat+Γ model groups together cnidaria and porifera with a posterior probability of. we then progressively removed heteropecillous positions according to their increasing pipn value; that is, we first removed positions that most likely violate the assumption of homogeneity over time of the cat model. for the five sub-datasets analyzed with the cat+Γ model, three kinds of topologies are observed: with positions, the same topology as with the complete dataset; with and positions, eumetazoa are recovered, but porifera are not monophyletic; with the smallest subsets, eumetazoa and porifera are both monophyletic. with the removal of sites with heterogeneous profile affiliations across clades, support for the monophyly of eumetazoa increases steadily, up to about one, and decreases to in the smallest dataset, probably due to the limited amount of data. this result suggests that sites showing a substitutional heterogeneity of profiles over time interfere with the phylogenetic signal and may eventually result in an erroneous topology. we computed the number of sites affiliated to the same profile in cnidaria and porifera and in cnidaria and bilateria. choanoflagellates are not considered because only two species are available and an accurate profile affiliation to sites is not possible. figure shows that, upon removal of heteropecillous sites, according to the pipn criterion, the number of sites having the same profile in porifera and cnidaria decreases much more rapidly than the number of sites having the same profile in cnidaria and bilateria. one can reasonably argue that sites with the same profile in porifera and cnidaria generate a spurious signal that is erroneously interpreted by the cat+Γ model as synapomorphies for porifera+cnidaria, since this model assumes that profiles are identical over the whole tree. two controls were performed. first, since a correlation exists between evolutionary rate and heteropecilly, improvements in phylogenetic inference could be due to the removal of fast evolving sites. when fast evolving sites are progressively removed according to the sf method, support for the incorrect grouping of cnidaria and porifera slightly increases. this is in sharp contrast with the removal of heteropecillous positions. second, since the negative effect of heteropecilly may constitute a model violation more important for the cat+Γ model than for site-homogeneous models, the gtr+Γ or mtrev+Γ models are also used. results are completely different: the removal of heteropecillous positions does not affect phylogenetic inferences conducted with gtr+Γ or mtrev+Γ models, for which high support for cnidaria + porifera is always observed. this result is not unexpected because, while heteropecilly constitutes a model violation for the cat model, its effect on the site-homogeneous gtr model is less clear. to understand the different behavior of these two models, one has to distinguish two model violations, heterogeneity of the substitution process across sites and over time. the first violation is known to seriously exacerbate lba artifacts, because the amount of homoplasy is underestimated. it is not expected to decrease with the removal of heteropecillous positions and probably dominates when using the site-homogeneous gtr and mtrev models. the effect of the second model violation can only be observed with the cat model, which handles heterogeneity across sites. removing heteropecillous positions will reduce the time-homogeneous violation and increase the accuracy of the cat model. this hypothesis is corroborated by an evaluation of model fit by cross-validation. the gtr+Γ model fits the complete mitochondrial mt dataset better than the cat+Γ. this is probably due to the large number of parameters of the latter model and the limited amount of positions, since, with the larger nuclear dataset, the cat+Γ model has a better fit than the two models based on exchange matrices. importantly, after removal of most heteropecillous sites, the cat+Γ model appears to best fit the data. this indicates that heteropecilly constitutes a serious violation of the cat model, because the best fit is obtained despite the limited number of sites.  numerous heterogeneities of the evolutionary process have been discovered. most have a clear negative impact on phylogenetic inference when not adequately handled: heterogeneity of rate across lineages, of substitution type, of rate across sites, of composition across taxa and of substitution process across sites. some of them, such as heterotachy, seem to have a more limited effect on phylogenetic accuracy. further studies are needed to know in which category heteropecilly has to be classified. an important prerequisite is the development of models that handle the fact that substitution properties change over time. this could be achieved via a markov modulated cat model, similar to the covarion model, or via the use of breakpoints. to choose between these two approaches, it would be important to estimate whether sites generally change their properties in a collective manner or not, since only the second approach can model this aspect. the costs of these improved models could be major and would be useful for phylogenetic inference only if heteropecilly turned out to seriously impair accuracy. in any case such improved models would be helpful to advance our knowledge of protein evolution, since in most cases one can select a set of species for whose relationships are confidently known, which drastically simplifies the problem and allows the use of complex but computationally demanding models. what are the main reasons for these shifts in profile affiliation? the correlation between the pipn criterion and the substitution number might suggest a neutral explanation of the variation in substitutional profile over time. since the pipn criterion is also correlated with a change in hydropathy, this neutral explanation could be insufficient. these changes in the site-specific substitution process could be related to functional shifts, such as adaptation of organisms to new environments or of the protein to a new cellular environment. to answer these questions, it would be particularly relevant to study paralogous genes in which functional shifts are well characterized. 
we used the cat model to define an infinite mixture of substitution processes characterized by equilibrium frequencies over the twenty amino acids, a useful proxy for qualitatively estimating the evolutionary process. its effect on phylogenetic inference has recently obtained considerable attention, which led to the development of new models of sequence evolution. however, thus far focus has been on the quantitative heterogeneity of the evolutionary process, thereby overlooking more qualitative variations. one of the properties of protein evolution is the variation of the relative rate of substitutions across sites and over time, the latter is the phenomenon called heterotachy. we studied the importance of variation of the site-specific amino-acid substitution process over time and its possible impact on phylogenetic inference. the progressive removal of the sites with the most heterogeneous cat profiles across clades led to the recovery of the monophyly of eumetazoa, suggesting that this heterogeneity can negatively influence phylogenetic inference. characterizing properties of the data that are not being correctly handled by current models is therefore of prime importance. using two large datasets, we show that qualitative changes in site-specific substitution properties over time occurred significantly.
 model violations constitute the major limitation in inferring accurate phylogenies. characterizing properties of the data that are not being correctly handled by current models is therefore of prime importance. one of the properties of protein evolution is the variation of the relative rate of substitutions across sites and over time, the latter is the phenomenon called heterotachy. its effect on phylogenetic inference has recently obtained considerable attention, which led to the development of new models of sequence evolution. however, thus far focus has been on the quantitative heterogeneity of the evolutionary process, thereby overlooking more qualitative variations.  we studied the importance of variation of the site-specific amino-acid substitution process over time and its possible impact on phylogenetic inference. we used the cat model to define an infinite mixture of substitution processes characterized by equilibrium frequencies over the twenty amino acids, a useful proxy for qualitatively estimating the evolutionary process. using two large datasets, we show that qualitative changes in site-specific substitution properties over time occurred significantly. to test whether this unaccounted qualitative variation can lead to an erroneous phylogenetic tree, we analyzed a concatenation of mitochondrial proteins in which cnidaria and porifera were erroneously grouped. the progressive removal of the sites with the most heterogeneous cat profiles across clades led to the recovery of the monophyly of eumetazoa, suggesting that this heterogeneity can negatively influence phylogenetic inference.  the time-heterogeneity of the amino-acid replacement process is therefore an important evolutionary aspect that should be incorporated in future models of sequence change. 