BACKGROUND
the covarion hypothesis of molecular evolution proposes that selective pressures on a given amino acid or nucleotide site are dependent on the identity of other sites in the molecule that change throughout time, resulting in changes of evolutionary rates of sites along the branches of a phylogenetic tree  <cit> .
at the sequence level, covarion-like sites are often recognizable in alignment columns as residues that are conserved among taxa in one clade but variable among taxa in one or several other clades.
changes in rates at sites in different sequences are also referred to as type-i functional divergence  <cit>  or 'heterotachy'  <cit> .
covarion-like evolution is now widely recognized as an important mode of molecular evolution in protein-coding genes, structural rna, and dna regulatory elements .
this realization has fueled the development of several kinds of phylogenetic models including:  'covarion models' that are designed to model the stochastic changes of rates at sites in sequences  <cit> ,  discrete 'rate-shift' models where sudden changes in rates at multiple sites occur at particular splits in the tree  <cit> , and  mixture of branch lengths-based heterotachy models  <cit> .
empirical studies have shown that phylogenetic estimation under the covarion models may recover different optimal topologies than when estimation is performed ignoring covarion effects .
indeed, simulation studies have shown that under some branch-length conditions, use of rates-across-sites  models that ignore covarion effects may cause long-branch repulsion biases in the resulting phylogenetic estimates  <cit> .
other studies have focused on developing computational methods to detect whether covarion-like evolution occurs in protein families  <cit> , identify covarion or heterotachous sites to analyse functional shifts in a protein family  <cit>  and detect positive selection in primate and viral genes  <cit> .
covarion models with changing rates of evolution at sites in different parts of the tree build upon the simpler ras models that assume evolutionary rates are variable among sites but constant across lineages in a gene or protein.
ras is typically modeled by a 'discretized' approximation of the gamma distribution with a series of equiprobable rate classes  <cit> .
the modeling of covarion processes is more challenging.
typically, these models allow rates at a site to vary gradually through the tree according to a stochastic process.
the gradual rate shift in a covarion context can be formulated as a markov model of rate switching between different rate classes, usually eight or less.
five specific covarion models have been proposed that differ in the complexity of the rate switching processes  <cit> .
the simplest model, proposed by tuffley and steel  <cit> , assumes that rates along a branch in a phylogenetic tree can have two states 'off' and 'on'; switching from 'off' to 'on' occurs with one rate  and from 'on' to 'off'  with another rate.
when a site is 'off', no substitutions occur and when it is 'on', substitutions occur at a constant rate.
huelsenbeck  <cit>  added additional rate classes to this model.
in the huelsenbeck model, when the site is 'on', the expected substitution rate per unit time at the site is a specific rate drawn from the discrete gamma distribution, whereas it is zero when it is 'off'.
a third covarion model was developed by galtier  <cit> , who assumed that only a subset of sites  evolve under the covarion process.
the remaining sites have a site-specific rate drawn from a discrete gamma distribution.
for sites evolving under the covarion process, rates are also drawn from a gamma distribution and the different rate classes can switch freely between each other at a single rate .
a more general model that combines features of both the tuffley-steel/huelsenbeck models and the galtier model was recently proposed  <cit> , in which a covarion site may not only switch between an 'on' and 'off' state but also can switch between different rate categories of 'on' states.
this latter model allows a variety of switching rates between the rate classes.
more recently whelan proposed a further generalized model which allows substitution rate-matrix changes as well as rate switches along the tree branches .
the first four covarion models are described in wang et al.
<cit>  which were implemented in procov for maximum likelihood  estimation of covarion parameters for a fixed phylogenetic tree and protein alignment data.
the new version of procov described herein allows ml-based tree estimation using the subtree-pruning regrafting  algorithm, under a variety of amino acid substitution models including jtt, wag and lg  <cit> .
we have also utilized several numerical libraries in procov to improve the efficiency of the likelihood calculations and thereby make computationally intensive tree searching analyses more practical.
here we demonstrate the utility of procov in performing 'deep-level' phylogenomic analyses where model misspecification can often lead to long-branch attraction.
we further explore the use of procov as a way to detect covarion sites in protein families that have structural and functional significance.
implementation
as in all common likelihood-based methods, procov implements a pruning algorithm  <cit>  for the likelihood calculation.
in conventional markov models of protein evolution, there are  <dig> amino acid states and the substitution rates of the amino acids are described by an instantaneous substitution rate matrix , such as the jtt model.
under the covarion model, character states are two dimensional, describing both the amino acid state and the substitution rate at that state at any given time.
the q matrix in a covarion model is thus a large sparse matrix.
in procov we used an algorithm introduced in  <cit>  to decompose the q matrix into a sum of two kronecker products, each consisting of two smaller matrices.
even with this efficient algorithm, the calculation of the likelihood of the data for a given tree with the general covarion model is about  <dig> to  <dig> fold slower than for an ras model with the same number of rate categories.
this is because likelihood calculations under the general covarion model have a much larger number of terms to be summed over at each ancestral node as compared to an ras model.
for instance, under the general covarion model with  <dig> rates, there is a 16-fold increase in the number of terms to be summed relative to that under the ras model.
for a given topology, ml estimates of parameters are obtained by a modified newton-raphson algorithm which requires the calculation of derivatives of the likelihood function with respect to each adjustable parameter.
as analytical derivatives are difficult to compute for the covarion parameters, numerical derivatives are computed for all three covarion switching parameters.
the derivative for the proportion of covarion sites parameter π in the general and galtier models is computed analytically as the difference of the covarion likelihood and ras likelihood across the sites.
for the tree searching function, we used the spr algorithm implemented in nhml  <cit> .
an initial tree is modified by pruning subtrees and moving them to other places.
if a rearrangement results in an increase of the likelihood, that tree is kept as a starting tree.
the algorithm iterates until no rearrangement increases the likelihood.
procov is written in ansi c, and is based on the phylogenetic inference package nhml  <cit> .
the current version of procov needs a user-supplied starting tree which should be rooted; the "retree" program of phylip  <cit>  can be used for re-rooting.
the starting tree can, for instance, be a neighbor-joining tree or a parsimony tree available from most phylogenetic packages.
compared with nhml, procov has numerous new features, including, for instance, a command-line argument for setting models, parameters, input and output data; implementing protein models and four covarion models ; new functions for matrix decomposition, matrix operations and computing derivatives.
we have also introduced the following algorithms to speed up the tree searching procedure.
since the optimization of the gamma shape parameter  and the covarion parameters takes time, during the tree searching stage, we re-optimize these parameters only when a tree with a higher likelihood than the previous best tree is found.
in this way, these parameters drift to the optimal values as the search proceeds.
furthermore, we relax the convergence condition to optimize parameters during the tree search stage; parameter optimization stops when the log-likelihood difference between two consecutive iterations is less than  <dig> .
for optimizing the final optimal tree, we impose a much stricter constraint .
although the likelihood gain from a stricter convergence threshold is usually small , according to our simulation results, this threshold yields parameter values much closer to their true values.
some of the nhml routines are particularly useful for saving tree searching time and so have been inherited by procov.
for example, if a starting tree in the newick format contains high bootstrap values that are greater than the maximum bootstrap value allowed for branch move during the spr searches , those branches will not move separately in the spr stage.
similarly it also has a function to forbid moving those branches that are longer than a user-defined value .
this branch movement restriction, resulting in partial spr searches, gives user the flexibility in choosing which internal nodes are fixed and therefore can greatly reduce tree search time if many nodes are fixed.
an extreme form of this branch movement restriction is to restrict procov to compare only a few competing topologies, as in our previously published analyses of angiosperm phylogeny .
furthermore, procov inherits from nhml a 'restart' function that can save all of the currently evaluated trees so that it will automatically bypass those topologies if the program has to be started over again.
these functions are of practical importance as ml estimation under the general covarion model will usually take several days for a moderate-sized dataset .
for compilation of the source code, we recommend the use of gcc or compatible compilers.
use of the -o <dig> and -funroll-loops for compiler optimization also significantly increases its running speed; for a small dataset we tested, this speedup can be more than two fold.
procov spends a lot of time doing matrix operations, such as matrix multiplication, matrix inversion and eigenvalue/eigenvector decomposition.
to do these kinds of calculations, phylogenetic programs including nhml commonly use c routines based on those described in numerical recipes  <cit> .
to improve speed, the current version of procov makes use of the high quality routines in basic linear algebra subprograms  implemented in automatically tuned linear algebra software  to perform basic vector and matrix operations.
this has been found to increase the speed of procov by at least three fold .
recommendations for utilizing the blas libraries other than through the free atlas  are included in the makefile.
RESULTS
comparing the speedup of procov with the new blas implementation
to compare the speedup of procov with the blas implementation versus the non-blas implementation, we tested two protein datasets  with  <dig> taxa and  <dig> sites) for fixed topologies, previously inferred with phyml  <cit>  under jtt + gamma, and optimized the parameters with jtt + the general covarion model with procov.
for acetyl-coa carboxylase, with the blas implementation, it took  <dig> minutes to finish parameter optimization and obtain the final log-likelihood score whereas the non-blas version took  <dig> hour  <dig> minutes for the same analysis.
for the hsp <dig> data set with a fixed tree, the blas and non-blas versions took  <dig> minutes and  <dig> hours  <dig> minutes, respectively.
the final likelihood scores yielded by the blas and non-blas versions of procov are the same in both cases.
to assess the performance for procov on tree searching, we simulated five datasets of  <dig> sites with seq-gen-aminocov  <cit>  based on a tree topology obtained from a 17-taxon  <dig> kda chaperonin  dataset  <cit> .
the reference tree and simulated datasets are available on the procov web site.
the simulations employed the jtt model and the ras, tuffley-steel , galtier, huelsenbeck and general covarion models, respectively.
for the models with an ras-process , four gamma rates were used in simulation.
we then used procov to estimate the topology for each dataset under the corresponding model  and   <dig> gamma rates, with a starting tree that was obtained with the neighbor-joining method by phylip for each dataset.
procov successfully recovered the same true topology in each case.
the speedups in procov with the blas versus non-blas implementations are  <dig> ,  <dig> ,  <dig>  and  <dig>  fold for the four covarion models, respectively.
there is no speedup for the ras model, as the blas libraries are not implemented for calculations under the ras model.
the above comparative results with and without the use of the blas libraries were conducted on a computer with a  <dig>  ghz intel quad core xeon processor with  <dig>  gb ram.
similar speedups were also observed on a computer with a different cpu architecture .
as the general covarion model takes much more time than the simpler covarion  and ras models in inferring a phylogeny, we asked whether the simpler models can correctly reconstruct phylogenies for data simulated under the general model.
this is of practical importance as if it is true then the simpler models would prefer to be used to reconstruct the phylogenies to save computational time.
for the cpn <dig> dataset simulated under the general model, we used each of the simpler models  to estimate the phylogenies.
each of the models were able to correctly infer the same correct topology as using the general model but the running times were very different.
the ras and ts models took  <dig> and  <dig> minutes, respectively; both the galtier and huelsenbeck models took more than  <dig> hours.
despite the same topology inferred, branch length estimates were different from those that estimated under the general model, which were the true model that were used for simulation.
figure  <dig> shows the true tree length used for simulation and estimates under the ras, ts, galtier, huelsenbeck and general models.
the tree lengths were separated as the sum of the external branch lengths and the sum of the internal branch lengths.
all simpler models, especially ras, underestimated the branch lengths.
therefore, the general model may not be replaced by simpler models for certain types of data.
previously we showed that the general model can converge to the huelsenbeck and galtier models when datasets are simulated under these models  <cit> .
here we further show that the general model can even adapt to the ras model when the data are simulated under ras.
for the cpn <dig> dataset simulated under the jtt + ras model, the general model recovered the same correct topology as the ras model.
moreover, the branch length estimates under the general model are very close to that under the ras model  for the total true tree length of  <dig> .
it turns out that the general model was able to adjust the covarion parameters  to converge to the ras model.
therefore, both the topology and the branch lengths were correctly inferred.
for the same ras-simulated dataset we also found that the huelsenbeck model was able to correctly estimate both the topology and branch lengths by adjusting parameters to mimic ras-like process .
establishing the phylogenetic position of microsporidia
with the speedups of procov made possible by the use of the blas libraries and other algorithmic improvements, it is now possible to compare topologies and, in some cases perform tree searches, for larger phylogenomic datasets  under covarion models in a reasonable time.
deep phylogenetic analyses of the eukaryote tree are often plagued with long-branch attraction  artifacts, even when large multi-gene phylogenomic data sets are used  <cit> .
one of the most famous examples of this concerns the position of microsporidia, a group of fast-evolving intracellular parasites that are now known to be relatives of fungi.
when reconstructing the phylogeny of eukaryotes rooted by archaea, if the estimation is performed with ml under an ras model, the extremely long branch leading to microsporidia is usually attracted to the long branches leading to the archaea at the base of the eukaryotes regardless of what amino acid substitution models are used.
many methods have been proposed to solve this problem, including selective taxon sampling, removal of fast-evolving proteins and saturated sites  <cit> , accounting for covarion shifts  <cit> , amino acid profile mixture modeling  <cit> , branch length mixture modeling  <cit>  and rare genomic changes of conserved amino acids-multiple substitutions  <cit> .
here we applied the general covarion model + wag to a large eukaryote phylogenomic data set made up of  <dig> proteins from  <dig> taxa and  <dig> sites  <cit>  and calculated the log likelihoods of two competing trees: the lba topology where microsporidia groups with archaea and, the correct topology where microsporidia groups with fungi.
the general covarion model clearly supports the correct microsporidia + fungi tree with a large log-likelihood gain  compared to the lba tree .
in contrast, the correct tree has a smaller log-likelihood than the lba tree under the ras + wag model .
thus, for a real example, the covarion model appears to be less susceptible to the effects of lba than the ras model.
we recently also used a site-specific class frequency mixture model implemented in qmmraxml  <cit>  to analyze the data and found the mixture model supported the correct topology, albeit with a smaller log likelihood gain .
therefore, in this particular case, the lba problem can be overcome with more realistic phylogenetic models that either account for site-specific substitution dynamics or covarion-like evolution.
detecting covarion sites of functional and structural significance
covarion models are useful not only because of improved phylogenetic estimation; they can also be used to identify patterns of sequence evolution that explain divergence in protein function or structure.
previous computational work on elongation factors  has nicely demonstrated that identifying evolutionary site-rate shifts coupled with analyses of three-dimensional structures of the protein family can pinpoint sites that are likely important in functional divergence and structural change between bacterial elongation factor tu  and eukaryotic elongation factor 1α   <cit> .
in fact a number of additional methods have been developed over the last decade to identify rate-shifted sites for the same purpose  <cit> .
most of these methods rely upon assuming that a discrete shift in rates at many sites has occurred over one branch in the protein phylogeny under examination and estimation of the phylogeny is usually performed beforehand using standard phylogenetic models.
since procov specifically models changing rates at sites during tree estimation, it may also be used to detect such rate shifts and has the added advantage that these rate shifts need not occur only on a single split in the phylogeny.
to illustrate its utility, we reanalyzed the bacterial and eukaryotic ef data  described in gaucher et al.
<cit> .
we first inferred ml phylogenetic trees for this data using procov under wag + ras and wag + the general covarion model, respectively.
we also used qmmraxml  <cit>  with the wag + ras model to obtain an ml tree.
all the three methods estimated the same tree , .
the log-likelihood difference for this tree under the general covarion model versus the ras model is  <dig>  , confirming that the ef data shows covarion-like properties  <cit> .
the parameters estimated under the general covarion model are α =  <dig> ; s <dig> =  <dig> ; s <dig> =  <dig> ; s <dig> =  <dig>  and π =  <dig> .
as described in wang et al.
<cit> , under this model site likelihoods are computed separately for both the covarion model  and the ras model  and then combined to get a total weighted likelihood of the site  × l_ras).
in order to determine which sites show a strong 'covarion' pattern, the difference in log-likelihoods between the two models at a site can be calculated as Λ = ln - ln; covarion sites are expected to have a Λ >>  <dig> as compared to sites that do not change rates across the tree.
the -l option of procov's command line arguments allows user to extract ln and ln for each site.
we used this option to get the site likelihoods for the ef data under the general covarion model and calculated Λ for each site.
two hundred and forty out of the  <dig> sites have a positive Λ  > ln) while the remaining  <dig> sites have a negative Λ. figure  <dig> shows a histogram of the distribution of Λ with a mean of  <dig>  and standard deviation  of  <dig> .
twenty one sites have a Λ > mean +  <dig> sd and  <dig> sites have a Λ > mean +  <dig> sd.
however, use of this distribution to identify covarion sites is not straightforward since the long right tail of the distribution  is likely due to the presence of many covarion sites.
thus the sd of this distribution is expected to be inflated relative to the sd of Λ distribution for non-covarion sites.
in order to get a valid cutoff value for Λ that indicates a significant likelihood difference between the two models at the site , we used seq-gen-aminocov  <cit>  to simulate  <dig>  sites based on the ef tree  under the wag + ras model with  <dig> gamma rates and α =  <dig> , which is the fitted α for the original ef data estimated with the ras model.
we then used the general covarion model to calculate the covarion and ras site likelihoods by fixing the topology and allowing all the parameters to be optimized.
we calculated Λ for each of the  <dig>  sites of the simulated ef data.
figure  <dig> shows the frequency distribution of Λ for the simulated data.
the 99th percentile of the Λ distribution is  <dig> , and can be used as a threshold for statistical significance.
note however this threshold value of  <dig>  is model  specific and therefore is only valid for the current original ef data.
using this criterion,  <dig> sites in the original ef dataset have a Λ greater than  <dig>  which can be considered covarion sites with confidence at p <  <dig> .
as the dataset has  <dig> positions, one can expect about  <dig>  sites could be covarion sites  by chance alone.
this indicates there are many more sites with significant log likelihood differences than expected in the original ef dataset.
the  <dig> covarion sites procov identified constitute  <dig> % of the total sites in the ef alignment.
this estimate is consistent with some suggestions that about 10% of the sites in this data set are covarion sites  <cit> .
however, the ef dataset used here is a relatively small one with  <dig> ef-1α and  <dig> ef-tu sequences.
inclusion of more taxa would be expected to increase the proportion of covarion sites detected  <cit> .
in any case, table  <dig> shows a list of the  <dig> identified positions in descending order of Λ. the sequence alignment columns corresponding to these sites are shown in figure  <dig>  twenty four of the sites  are the sites that were previously identified to be covarion sites of functional/structural significance by gaucher et al.
<cit> .
each of these displays a typical 'covarion-like' site pattern, i.e., a lot of different amino acid states occur in the bacterial group but no or little change occurs at that site in the eukaryotic group or vice versa.
eighteen sites  were detected as covarion sites by procov also demonstrate this typical covarion pattern but were not flagged by gaucher and colleagues' method.
as an independent test, we also used our rate-shift detection program bivar  <cit>  to estimate rate differences between the two subgroups of ef-tu and ef-1α, which recovered  <dig> sites as rate shifted with a p-value <  <dig> .
thirty one of these are the same covarion sites as picked up by procov .
eleven sites  identified by procov as covarion sites clearly show a typical covarion pattern, but these sites were not picked up by bivar .
these comparisons indicate procov may have more power to identify covarion sites than either bivar or the gaucher et al.
method.
* sequence position is based on the ef alignment  <cit> .
**Λ = ln - ln.
***c: sites were found to be covarion sites of functional or structural significance in  <cit> ; c1: sites having typical covarion site pattern but missed in  <cit> ; c2: site showing non-typical covarion site pattern.
****p-values from a bivar analysis of the ef data.
site  <dig> in table  <dig> is particularly interesting, as it ranks relatively high  among the log-likelihood differences between a covarion process and the ras process yet has a non-significant bivar p-value of  <dig> .
the method used by gaucher and colleagues also did not pick up this site as a covarion site.
inspection of the residues at this site , reveals that it does not have a typical 'covarion pattern' as the site is variable in both bacterial and eukaryotic efs.
the ef-1α subgroup is slightly more variable at this site, displaying  <dig> different amino acids that collectively can be binned into  <dig> of the six different "dayhoff" groups of amino acids  as compared to the ef-tu subgroup, which has  <dig> different amino acids from  <dig> of the dayhoff groups   <cit> .
figure  <dig> shows the amino acids at site  <dig> mapped on to the ef-tu/1α tree.
close inspection of the substitutions at this site in the ef-1α subtree reveals that a number of radical amino acid changes occur between relatively closely related sister taxa in the tree .
such radical changes are not observed in similarly closely related bacteria in the ef-tu subtree.
a subsequent analysis of the two subgroups separately with the general covarion model indicates the eukaryotic ef1α has a very big positive difference between ln and ln at site  <dig>  which suggests it could be a covarion site for the eukaryotic subset.
by contrast, Λef-tu =  <dig>  for the bacterial subset and is unlikely a covarion site, although a simulation study is needed to determine a Λ threshold for the two subtrees separately.
despite the strong support for site  <dig> being a covarion site in ef-1α, the residues at the site do not present a typical covarion pattern where variability is differentially restricted in different groups.
one possible explanation of these observations is that the covarion model is compensating for the radical substitutions between closely related taxa observed in the ef-1α subtree, which are not consistent with the wag substitution model.
a rate-switching process could accommodate such radical substitutions by in effect 'lengthening' the branches between closely related taxa.
this is in contrast to an ras model where the rates of evolution must remain constant across the tree even if radical substitutions are observed in some closely-related taxa but not in others.
to test the idea that the covarion model was compensating for this kind of substitution model misspecification at site  <dig>  we compared the likelihood of this site under a simple proportional  model  relative to the site likelihood under the wag + ras model.
as expected, for ef-1α site  <dig>  the ln = - <dig>  which is greater than ln = - <dig> , despite the fact that over all sites the wag + ras model has a greater log-likelihood  than the prop + ras model  for this subgroup.
this result suggests that it is the low exchangeability rates in the wag model corresponding to the radical amino acid changes observed at this site that lead to the poor model fit.
although unexpected, it seems that the covarion model compensates for this kind of model misspecification at sites that do not show classical covarion-type variability patterns.
in any case, although the  <dig> covarion sites listed in table  <dig> are scattered throughout the  <dig> sequence positions in the alignment, they are not randomly distributed.
for instance, there is a long sequence segment  that, with the exception of sites  <dig> and  <dig>  are all covarions.
this segment maps to a surface loop region of the ef-tu structure , that is a possible ribosome binding site in bacteria  <cit> .
the pattern at these sites indicates that the bacterial ef-tu sequences are typically variable whereas the eukaryotic ef-1α sequences are conserved, hinting a possible additional function for this loop in eukaryotes that is absent in bacteria.
discussion
we have developed procov, an ml-based phylogenetic program for modeling the covarion processes of protein evolution.
we showed that compiler optimization, especially the use of highly optimized math libraries, such as blas, can significantly speed up likelihood calculation.
although blas and related math libraries have been widely used in high performance computing software , we are not aware of other phylogenetic software that utilize these efficient libraries.
the use of the optimized math library together with some features of procov described above makes it tractable to do full tree search under the general covarion model for datasets of moderate size in a reasonable time .
for large datasets one can selectively restrict the movements of those branches and nodes that deem to be in the same group when running procov.
this partial search will considerably reduce tree search time when many nodes and branches are fixed.
for even larger phylogenomic data one can use procov to analyse several competing trees that were already established by other phylogenetic methods and see which of them is preferred by the general model.
we applied this method to the microsporidia phylogenomic dataset  <cit>  and the general model clearly supports the correct microsporidia-fungi clade tree over the lba-induced microsporidia-anchaea clan tree.
however, this may not guarantee it is the optimal tree for the general model if a tree search is conducted.
for example, a partial tree search of this microsporidia data estimated a tree of microsporidia-protist clade that had a higher likelihood than the tree of microsporidia-fungi clade.
examples in this study show that phylogenetic tree estimation under a covarion model may or may not estimate a different optimal topology than that under a non-covarion ras model.
for the simulated cpn <dig> datasets as well as the ef dataset, the ras and covarion models estimated the same optimal topologies; for the microsporidia data they differ.
our previous simulations and analytical studies explored topology and branch length conditions that the ras and covarion models will likely estimate different topologies  <cit> .
results in figure  <dig> show that even though the ras model was able to estimate the correct topology for data simulated under the general model, it would underestimate the branch lengths.
both the general and huelsenbeck models, however, will correctly infer the topology and accurately estimate branch lengths for the data simulated under the ras model.
they do so by adjusting the covarion parameters to converge to the ras model.
for real data, we do not know in advance whether the data follow covarion or ras evolutions or both.
the general model, including the ras and ts, huelsenbeck and galtier models as special cases, has the advantage of adapting to the right model in the course of parameter optimization so that it can analyse all relevant types of data appropriately, but suffers from heavy computing loads with large amounts of data.
a recent empirical test of the covarion hypothesis has shown that the frequency of covarion-sites increases with genetic distance  <cit> .
this suggests covarion-based phylogenetic inference may be useful in the estimation of the divergence time of the species spanning longer time periods.
it will therefore be interesting to revisit the estimates of dates of divergence using relaxed molecular clock methods  <cit>  in conjunction with covarion models of evolution.
in addition to the advantages of procov for phylogenetic inference under the general model, we also demonstrated that it had more power to detect covarion sites than several previous methods.
it can also be used to pinpoint those lineages where covarions are located .
like the general covarion model, covarion and ras site likelihoods are also separately calculated under the galtier model.
by contrast, the ts models is not a mixture of covarion and ras processes; the huelsenbeck model, as originally formulated, does not calculate covarion and ras site likelihoods, separately.
therefore only covarion site likelihoods are calculated for the ts and huelsenbeck models.
nevertheless, one can run two separate analyses with procov, one under either of the two models, another under the ras model, and compare their site likelihood differences to obtain Λ's for sites.
all of the four covarion models considered here are stationary time reversible models with an expectation that the proportion of variable sites  is the same in all evolutionary lineages.
however, this assumption can be overly restrictive as proportions of variable sites may vary in different lineages  <cit> .
a sequence generator for generating lineage-specific variation in the pvar is recently reported  <cit> .
a fruitful area of future development of procov may therefore be to model both changes in the proportion of variable sites and the covarion-based rate changes and switches.
furthermore, the current implemented covarion models assume that rate switching between 'on' and 'off' states and among different 'on' rates are homogenous across sites and the tree, which may not be realistic.
this is especially suspicious for large phylogenomic datasets that are from the concatenation of multiple genes of diverse functions with different functional constraints.
for instance, we previously reported that the covarion parameters, like the α parameter of the ras process, varied across different protein families .
it will be interesting to model this heterogeneity in switch rates variation across sites and lineages and implement it in procov without increasing computational load too much.
finally, the current release of procov  only handles protein sequence data.
analyses of dna substitutions under covarion models have found applications in inferring the evolutionary history of viral genes  <cit> .
future extension of procov to allow analyses of dna sequence data may be useful to investigate these kinds of data sets.
CONCLUSIONS
procov is a phylogenetic program to infer phylogeny under covarion models, which may be especially useful for problems involving estimates of deep divergences in the tree of life, where rates of evolution at sites are likely to have changed over the tree.
it can also be used to detect covarion sites, which when combined with structural bioinformatics approaches, can be a powerful method to study lineage-specific functional shifts in protein families as well as protein adaptation.
availability and requirements
* project name: procov: maximum likelihood estimation of protein phylogeny under covarion models .
* project home page: .
* operating system: any unix/linux platform.
* programming language: c.
* other requirements: gcc  or compatible compiler.
it is recommended to have the blas/atlas libraries  installed on the unix/linux system so that procov can run faster.
versions of blas and lapack, such as the generic versions from atlas, netlib, or vendor-provided libraries that work with your compiler should be installed.
the makefile should then be edited to match the type of the compiler and the path and library names of the blas and lapack libraries used.
the makefile of the procov source code gives some instances of the blas installation on a few commonly-used unix systems.
* license: gnu gpl.
* any restrictions to use by non-academics: none.
authors' contributions
all authors conceived of the study, designed the method and wrote the article.
h-cw implemented the software and ran the analyses.

 the covarion hypothesis of molecular evolution proposes that selective pressures on a given amino acid or nucleotide site are dependent on the identity of other sites in the molecule that change throughout time, resulting in changes of evolutionary rates of sites along the branches of a phylogenetic tree. at the sequence level, covarion-like sites are often recognizable in alignment columns as residues that are conserved among taxa in one clade but variable among taxa in one or several other clades. changes in rates at sites in different sequences are also referred to as type-i functional divergence or'heterotachy'. covarion-like evolution is now widely recognized as an important mode of molecular evolution in protein-coding genes, structural rna, and dna regulatory elements. this realization has fueled the development of several kinds of phylogenetic models including:'covarion models' that are designed to model the stochastic changes of rates at sites in sequences, discrete'rate-shift' models where sudden changes in rates at multiple sites occur at particular splits in the tree, and mixture of branch lengths-based heterotachy models. empirical studies have shown that phylogenetic estimation under the covarion models may recover different optimal topologies than when estimation is performed ignoring covarion effects. indeed, simulation studies have shown that under some branch-length conditions, use of rates-across-sites models that ignore covarion effects may cause long-branch repulsion biases in the resulting phylogenetic estimates. other studies have focused on developing computational methods to detect whether covarion-like evolution occurs in protein families, identify covarion or heterotachous sites to analyse functional shifts in a protein family and detect positive selection in primate and viral genes. covarion models with changing rates of evolution at sites in different parts of the tree build upon the simpler ras models that assume evolutionary rates are variable among sites but constant across lineages in a gene or protein. ras is typically modeled by a'discretized' approximation of the gamma distribution with a series of equiprobable rate classes. the modeling of covarion processes is more challenging. typically, these models allow rates at a site to vary gradually through the tree according to a stochastic process. the gradual rate shift in a covarion context can be formulated as a markov model of rate switching between different rate classes, usually eight or less. five specific covarion models have been proposed that differ in the complexity of the rate switching processes. the simplest model, proposed by tuffley and steel, assumes that rates along a branch in a phylogenetic tree can have two states'off' and'on'; switching from'off' to'on' occurs with one rate and from'on' to'off' with another rate. when a site is'off', no substitutions occur and when it is'on', substitutions occur at a constant rate. huelsenbeck added additional rate classes to this model. in the huelsenbeck model, when the site is'on', the expected substitution rate per unit time at the site is a specific rate drawn from the discrete gamma distribution, whereas it is zero when it is'off'. a third covarion model was developed by galtier, who assumed that only a subset of sites evolve under the covarion process. the remaining sites have a site-specific rate drawn from a discrete gamma distribution. for sites evolving under the covarion process, rates are also drawn from a gamma distribution and the different rate classes can switch freely between each other at a single rate. a more general model that combines features of both the tuffley-steel/huelsenbeck models and the galtier model was recently proposed, in which a covarion site may not only switch between an'on' and'off' state but also can switch between different rate categories of'on' states. this latter model allows a variety of switching rates between the rate classes. more recently whelan proposed a further generalized model which allows substitution rate-matrix changes as well as rate switches along the tree branches. the first four covarion models are described in wang et al. which were implemented in procov for maximum likelihood estimation of covarion parameters for a fixed phylogenetic tree and protein alignment data. the new version of procov described herein allows ml-based tree estimation using the subtree-pruning regrafting algorithm, under a variety of amino acid substitution models including jtt, wag and lg. we have also utilized several numerical libraries in procov to improve the efficiency of the likelihood calculations and thereby make computationally intensive tree searching analyses more practical. here we demonstrate the utility of procov in performing'deep-level' phylogenomic analyses where model misspecification can often lead to long-branch attraction. we further explore the use of procov as a way to detect covarion sites in protein families that have structural and functional significance. implementation as in all common likelihood-based methods, procov implements a pruning algorithm for the likelihood calculation. in conventional markov models of protein evolution, there are amino acid states and the substitution rates of the amino acids are described by an instantaneous substitution rate matrix, such as the jtt model. under the covarion model, character states are two dimensional, describing both the amino acid state and the substitution rate at that state at any given time. the q matrix in a covarion model is thus a large sparse matrix. in procov we used an algorithm introduced in to decompose the q matrix into a sum of two kronecker products, each consisting of two smaller matrices. even with this efficient algorithm, the calculation of the likelihood of the data for a given tree with the general covarion model is about to fold slower than for an ras model with the same number of rate categories. this is because likelihood calculations under the general covarion model have a much larger number of terms to be summed over at each ancestral node as compared to an ras model. for instance, under the general covarion model with rates, there is a 16-fold increase in the number of terms to be summed relative to that under the ras model. for a given topology, ml estimates of parameters are obtained by a modified newton-raphson algorithm which requires the calculation of derivatives of the likelihood function with respect to each adjustable parameter. as analytical derivatives are difficult to compute for the covarion parameters, numerical derivatives are computed for all three covarion switching parameters. the derivative for the proportion of covarion sites parameter π in the general and galtier models is computed analytically as the difference of the covarion likelihood and ras likelihood across the sites. for the tree searching function, we used the spr algorithm implemented in nhml. an initial tree is modified by pruning subtrees and moving them to other places. if a rearrangement results in an increase of the likelihood, that tree is kept as a starting tree. the algorithm iterates until no rearrangement increases the likelihood. procov is written in ansi c, and is based on the phylogenetic inference package nhml. the current version of procov needs a user-supplied starting tree which should be rooted; the "retree" program of phylip can be used for re-rooting. the starting tree can, for instance, be a neighbor-joining tree or a parsimony tree available from most phylogenetic packages. compared with nhml, procov has numerous new features, including, for instance, a command-line argument for setting models, parameters, input and output data; implementing protein models and four covarion models; new functions for matrix decomposition, matrix operations and computing derivatives. we have also introduced the following algorithms to speed up the tree searching procedure. since the optimization of the gamma shape parameter and the covarion parameters takes time, during the tree searching stage, we re-optimize these parameters only when a tree with a higher likelihood than the previous best tree is found. in this way, these parameters drift to the optimal values as the search proceeds. furthermore, we relax the convergence condition to optimize parameters during the tree search stage; parameter optimization stops when the log-likelihood difference between two consecutive iterations is less than. for optimizing the final optimal tree, we impose a much stricter constraint. although the likelihood gain from a stricter convergence threshold is usually small, according to our simulation results, this threshold yields parameter values much closer to their true values. some of the nhml routines are particularly useful for saving tree searching time and so have been inherited by procov. for example, if a starting tree in the newick format contains high bootstrap values that are greater than the maximum bootstrap value allowed for branch move during the spr searches, those branches will not move separately in the spr stage. similarly it also has a function to forbid moving those branches that are longer than a user-defined value. this branch movement restriction, resulting in partial spr searches, gives user the flexibility in choosing which internal nodes are fixed and therefore can greatly reduce tree search time if many nodes are fixed. an extreme form of this branch movement restriction is to restrict procov to compare only a few competing topologies, as in our previously published analyses of angiosperm phylogeny. furthermore, procov inherits from nhml a'restart' function that can save all of the currently evaluated trees so that it will automatically bypass those topologies if the program has to be started over again. these functions are of practical importance as ml estimation under the general covarion model will usually take several days for a moderate-sized dataset. for compilation of the source code, we recommend the use of gcc or compatible compilers. use of the -o and -funroll-loops for compiler optimization also significantly increases its running speed; for a small dataset we tested, this speedup can be more than two fold. procov spends a lot of time doing matrix operations, such as matrix multiplication, matrix inversion and eigenvalue/eigenvector decomposition. to do these kinds of calculations, phylogenetic programs including nhml commonly use c routines based on those described in numerical recipes. to improve speed, the current version of procov makes use of the high quality routines in basic linear algebra subprograms implemented in automatically tuned linear algebra software to perform basic vector and matrix operations. this has been found to increase the speed of procov by at least three fold. recommendations for utilizing the blas libraries other than through the free atlas are included in the makefile.  comparing the speedup of procov with the new blas implementation to compare the speedup of procov with the blas implementation versus the non-blas implementation, we tested two protein datasets with taxa and sites) for fixed topologies, previously inferred with phyml under jtt + gamma, and optimized the parameters with jtt + the general covarion model with procov. for acetyl-coa carboxylase, with the blas implementation, it took minutes to finish parameter optimization and obtain the final log-likelihood score whereas the non-blas version took hour minutes for the same analysis. for the hsp data set with a fixed tree, the blas and non-blas versions took minutes and hours minutes, respectively. the final likelihood scores yielded by the blas and non-blas versions of procov are the same in both cases. to assess the performance for procov on tree searching, we simulated five datasets of sites with seq-gen-aminocov based on a tree topology obtained from a 17-taxon kda chaperonin dataset. the reference tree and simulated datasets are available on the procov web site. the simulations employed the jtt model and the ras, tuffley-steel, galtier, huelsenbeck and general covarion models, respectively. for the models with an ras-process, four gamma rates were used in simulation. we then used procov to estimate the topology for each dataset under the corresponding model and gamma rates, with a starting tree that was obtained with the neighbor-joining method by phylip for each dataset. procov successfully recovered the same true topology in each case. the speedups in procov with the blas versus non-blas implementations are,, and fold for the four covarion models, respectively. there is no speedup for the ras model, as the blas libraries are not implemented for calculations under the ras model. the above comparative results with and without the use of the blas libraries were conducted on a computer with a ghz intel quad core xeon processor with gb ram. similar speedups were also observed on a computer with a different cpu architecture. as the general covarion model takes much more time than the simpler covarion and ras models in inferring a phylogeny, we asked whether the simpler models can correctly reconstruct phylogenies for data simulated under the general model. this is of practical importance as if it is true then the simpler models would prefer to be used to reconstruct the phylogenies to save computational time. for the cpn dataset simulated under the general model, we used each of the simpler models to estimate the phylogenies. each of the models were able to correctly infer the same correct topology as using the general model but the running times were very different. the ras and ts models took and minutes, respectively; both the galtier and huelsenbeck models took more than hours. despite the same topology inferred, branch length estimates were different from those that estimated under the general model, which were the true model that were used for simulation. figure shows the true tree length used for simulation and estimates under the ras, ts, galtier, huelsenbeck and general models. the tree lengths were separated as the sum of the external branch lengths and the sum of the internal branch lengths. all simpler models, especially ras, underestimated the branch lengths. therefore, the general model may not be replaced by simpler models for certain types of data. previously we showed that the general model can converge to the huelsenbeck and galtier models when datasets are simulated under these models. here we further show that the general model can even adapt to the ras model when the data are simulated under ras. for the cpn dataset simulated under the jtt + ras model, the general model recovered the same correct topology as the ras model. moreover, the branch length estimates under the general model are very close to that under the ras model for the total true tree length of. it turns out that the general model was able to adjust the covarion parameters to converge to the ras model. therefore, both the topology and the branch lengths were correctly inferred. for the same ras-simulated dataset we also found that the huelsenbeck model was able to correctly estimate both the topology and branch lengths by adjusting parameters to mimic ras-like process. establishing the phylogenetic position of microsporidia with the speedups of procov made possible by the use of the blas libraries and other algorithmic improvements, it is now possible to compare topologies and, in some cases perform tree searches, for larger phylogenomic datasets under covarion models in a reasonable time. deep phylogenetic analyses of the eukaryote tree are often plagued with long-branch attraction artifacts, even when large multi-gene phylogenomic data sets are used. one of the most famous examples of this concerns the position of microsporidia, a group of fast-evolving intracellular parasites that are now known to be relatives of fungi. when reconstructing the phylogeny of eukaryotes rooted by archaea, if the estimation is performed with ml under an ras model, the extremely long branch leading to microsporidia is usually attracted to the long branches leading to the archaea at the base of the eukaryotes regardless of what amino acid substitution models are used. many methods have been proposed to solve this problem, including selective taxon sampling, removal of fast-evolving proteins and saturated sites, accounting for covarion shifts, amino acid profile mixture modeling, branch length mixture modeling and rare genomic changes of conserved amino acids-multiple substitutions. here we applied the general covarion model + wag to a large eukaryote phylogenomic data set made up of proteins from taxa and sites and calculated the log likelihoods of two competing trees: the lba topology where microsporidia groups with archaea and, the correct topology where microsporidia groups with fungi. the general covarion model clearly supports the correct microsporidia + fungi tree with a large log-likelihood gain compared to the lba tree. in contrast, the correct tree has a smaller log-likelihood than the lba tree under the ras + wag model. thus, for a real example, the covarion model appears to be less susceptible to the effects of lba than the ras model. we recently also used a site-specific class frequency mixture model implemented in qmmraxml to analyze the data and found the mixture model supported the correct topology, albeit with a smaller log likelihood gain. therefore, in this particular case, the lba problem can be overcome with more realistic phylogenetic models that either account for site-specific substitution dynamics or covarion-like evolution. detecting covarion sites of functional and structural significance covarion models are useful not only because of improved phylogenetic estimation; they can also be used to identify patterns of sequence evolution that explain divergence in protein function or structure. previous computational work on elongation factors has nicely demonstrated that identifying evolutionary site-rate shifts coupled with analyses of three-dimensional structures of the protein family can pinpoint sites that are likely important in functional divergence and structural change between bacterial elongation factor tu and eukaryotic elongation factor 1α. in fact a number of additional methods have been developed over the last decade to identify rate-shifted sites for the same purpose. most of these methods rely upon assuming that a discrete shift in rates at many sites has occurred over one branch in the protein phylogeny under examination and estimation of the phylogeny is usually performed beforehand using standard phylogenetic models. since procov specifically models changing rates at sites during tree estimation, it may also be used to detect such rate shifts and has the added advantage that these rate shifts need not occur only on a single split in the phylogeny. to illustrate its utility, we reanalyzed the bacterial and eukaryotic ef data described in gaucher et al.. we first inferred ml phylogenetic trees for this data using procov under wag + ras and wag + the general covarion model, respectively. we also used qmmraxml with the wag + ras model to obtain an ml tree. all the three methods estimated the same tree,. the log-likelihood difference for this tree under the general covarion model versus the ras model is, confirming that the ef data shows covarion-like properties. the parameters estimated under the general covarion model are α =; s =; s =; s = and π =. as described in wang et al., under this model site likelihoods are computed separately for both the covarion model and the ras model and then combined to get a total weighted likelihood of the site × l_ras). in order to determine which sites show a strong'covarion' pattern, the difference in log-likelihoods between the two models at a site can be calculated as Λ = ln - ln; covarion sites are expected to have a Λ >> as compared to sites that do not change rates across the tree. the -l option of procov's command line arguments allows user to extract ln and ln for each site. we used this option to get the site likelihoods for the ef data under the general covarion model and calculated Λ for each site. two hundred and forty out of the sites have a positive Λ > ln) while the remaining sites have a negative Λ. figure shows a histogram of the distribution of Λ with a mean of and standard deviation of. twenty one sites have a Λ > mean + sd and sites have a Λ > mean + sd. however, use of this distribution to identify covarion sites is not straightforward since the long right tail of the distribution is likely due to the presence of many covarion sites. thus the sd of this distribution is expected to be inflated relative to the sd of Λ distribution for non-covarion sites. in order to get a valid cutoff value for Λ that indicates a significant likelihood difference between the two models at the site, we used seq-gen-aminocov to simulate sites based on the ef tree under the wag + ras model with gamma rates and α =, which is the fitted α for the original ef data estimated with the ras model. we then used the general covarion model to calculate the covarion and ras site likelihoods by fixing the topology and allowing all the parameters to be optimized. we calculated Λ for each of the sites of the simulated ef data. figure shows the frequency distribution of Λ for the simulated data. the 99th percentile of the Λ distribution is, and can be used as a threshold for statistical significance. note however this threshold value of is model specific and therefore is only valid for the current original ef data. using this criterion, sites in the original ef dataset have a Λ greater than which can be considered covarion sites with confidence at p <. as the dataset has positions, one can expect about sites could be covarion sites by chance alone. this indicates there are many more sites with significant log likelihood differences than expected in the original ef dataset. the covarion sites procov identified constitute % of the total sites in the ef alignment. this estimate is consistent with some suggestions that about 10% of the sites in this data set are covarion sites. however, the ef dataset used here is a relatively small one with ef-1α and ef-tu sequences. inclusion of more taxa would be expected to increase the proportion of covarion sites detected. in any case, table shows a list of the identified positions in descending order of Λ. the sequence alignment columns corresponding to these sites are shown in figure twenty four of the sites are the sites that were previously identified to be covarion sites of functional/structural significance by gaucher et al.. each of these displays a typical'covarion-like' site pattern, i.e., a lot of different amino acid states occur in the bacterial group but no or little change occurs at that site in the eukaryotic group or vice versa. eighteen sites were detected as covarion sites by procov also demonstrate this typical covarion pattern but were not flagged by gaucher and colleagues' method. as an independent test, we also used our rate-shift detection program bivar to estimate rate differences between the two subgroups of ef-tu and ef-1α, which recovered sites as rate shifted with a p-value <. thirty one of these are the same covarion sites as picked up by procov. eleven sites identified by procov as covarion sites clearly show a typical covarion pattern, but these sites were not picked up by bivar. these comparisons indicate procov may have more power to identify covarion sites than either bivar or the gaucher et al. method. * sequence position is based on the ef alignment. **Λ = ln - ln. ***c: sites were found to be covarion sites of functional or structural significance in; c1: sites having typical covarion site pattern but missed in; c2: site showing non-typical covarion site pattern. ****p-values from a bivar analysis of the ef data. site in table is particularly interesting, as it ranks relatively high among the log-likelihood differences between a covarion process and the ras process yet has a non-significant bivar p-value of. the method used by gaucher and colleagues also did not pick up this site as a covarion site. inspection of the residues at this site, reveals that it does not have a typical'covarion pattern' as the site is variable in both bacterial and eukaryotic efs. the ef-1α subgroup is slightly more variable at this site, displaying different amino acids that collectively can be binned into of the six different "dayhoff" groups of amino acids as compared to the ef-tu subgroup, which has different amino acids from of the dayhoff groups. figure shows the amino acids at site mapped on to the ef-tu/1α tree. close inspection of the substitutions at this site in the ef-1α subtree reveals that a number of radical amino acid changes occur between relatively closely related sister taxa in the tree. such radical changes are not observed in similarly closely related bacteria in the ef-tu subtree. a subsequent analysis of the two subgroups separately with the general covarion model indicates the eukaryotic ef1α has a very big positive difference between ln and ln at site which suggests it could be a covarion site for the eukaryotic subset. by contrast, Λef-tu = for the bacterial subset and is unlikely a covarion site, although a simulation study is needed to determine a Λ threshold for the two subtrees separately. despite the strong support for site being a covarion site in ef-1α, the residues at the site do not present a typical covarion pattern where variability is differentially restricted in different groups. one possible explanation of these observations is that the covarion model is compensating for the radical substitutions between closely related taxa observed in the ef-1α subtree, which are not consistent with the wag substitution model. a rate-switching process could accommodate such radical substitutions by in effect'lengthening' the branches between closely related taxa. this is in contrast to an ras model where the rates of evolution must remain constant across the tree even if radical substitutions are observed in some closely-related taxa but not in others. to test the idea that the covarion model was compensating for this kind of substitution model misspecification at site we compared the likelihood of this site under a simple proportional model relative to the site likelihood under the wag + ras model. as expected, for ef-1α site the ln = - which is greater than ln = -, despite the fact that over all sites the wag + ras model has a greater log-likelihood than the prop + ras model for this subgroup. this result suggests that it is the low exchangeability rates in the wag model corresponding to the radical amino acid changes observed at this site that lead to the poor model fit. although unexpected, it seems that the covarion model compensates for this kind of model misspecification at sites that do not show classical covarion-type variability patterns. in any case, although the covarion sites listed in table are scattered throughout the sequence positions in the alignment, they are not randomly distributed. for instance, there is a long sequence segment that, with the exception of sites and are all covarions. this segment maps to a surface loop region of the ef-tu structure, that is a possible ribosome binding site in bacteria. the pattern at these sites indicates that the bacterial ef-tu sequences are typically variable whereas the eukaryotic ef-1α sequences are conserved, hinting a possible additional function for this loop in eukaryotes that is absent in bacteria. discussion we have developed procov, an ml-based phylogenetic program for modeling the covarion processes of protein evolution. we showed that compiler optimization, especially the use of highly optimized math libraries, such as blas, can significantly speed up likelihood calculation. although blas and related math libraries have been widely used in high performance computing software, we are not aware of other phylogenetic software that utilize these efficient libraries. the use of the optimized math library together with some features of procov described above makes it tractable to do full tree search under the general covarion model for datasets of moderate size in a reasonable time. for large datasets one can selectively restrict the movements of those branches and nodes that deem to be in the same group when running procov. this partial search will considerably reduce tree search time when many nodes and branches are fixed. for even larger phylogenomic data one can use procov to analyse several competing trees that were already established by other phylogenetic methods and see which of them is preferred by the general model. we applied this method to the microsporidia phylogenomic dataset and the general model clearly supports the correct microsporidia-fungi clade tree over the lba-induced microsporidia-anchaea clan tree. however, this may not guarantee it is the optimal tree for the general model if a tree search is conducted. for example, a partial tree search of this microsporidia data estimated a tree of microsporidia-protist clade that had a higher likelihood than the tree of microsporidia-fungi clade. examples in this study show that phylogenetic tree estimation under a covarion model may or may not estimate a different optimal topology than that under a non-covarion ras model. for the simulated cpn datasets as well as the ef dataset, the ras and covarion models estimated the same optimal topologies; for the microsporidia data they differ. our previous simulations and analytical studies explored topology and branch length conditions that the ras and covarion models will likely estimate different topologies. results in figure show that even though the ras model was able to estimate the correct topology for data simulated under the general model, it would underestimate the branch lengths. both the general and huelsenbeck models, however, will correctly infer the topology and accurately estimate branch lengths for the data simulated under the ras model. they do so by adjusting the covarion parameters to converge to the ras model. for real data, we do not know in advance whether the data follow covarion or ras evolutions or both. the general model, including the ras and ts, huelsenbeck and galtier models as special cases, has the advantage of adapting to the right model in the course of parameter optimization so that it can analyse all relevant types of data appropriately, but suffers from heavy computing loads with large amounts of data. a recent empirical test of the covarion hypothesis has shown that the frequency of covarion-sites increases with genetic distance. this suggests covarion-based phylogenetic inference may be useful in the estimation of the divergence time of the species spanning longer time periods. it will therefore be interesting to revisit the estimates of dates of divergence using relaxed molecular clock methods in conjunction with covarion models of evolution. in addition to the advantages of procov for phylogenetic inference under the general model, we also demonstrated that it had more power to detect covarion sites than several previous methods. it can also be used to pinpoint those lineages where covarions are located. like the general covarion model, covarion and ras site likelihoods are also separately calculated under the galtier model. by contrast, the ts models is not a mixture of covarion and ras processes; the huelsenbeck model, as originally formulated, does not calculate covarion and ras site likelihoods, separately. therefore only covarion site likelihoods are calculated for the ts and huelsenbeck models. nevertheless, one can run two separate analyses with procov, one under either of the two models, another under the ras model, and compare their site likelihood differences to obtain Λ's for sites. all of the four covarion models considered here are stationary time reversible models with an expectation that the proportion of variable sites is the same in all evolutionary lineages. however, this assumption can be overly restrictive as proportions of variable sites may vary in different lineages. a sequence generator for generating lineage-specific variation in the pvar is recently reported. a fruitful area of future development of procov may therefore be to model both changes in the proportion of variable sites and the covarion-based rate changes and switches. furthermore, the current implemented covarion models assume that rate switching between'on' and'off' states and among different'on' rates are homogenous across sites and the tree, which may not be realistic. this is especially suspicious for large phylogenomic datasets that are from the concatenation of multiple genes of diverse functions with different functional constraints. for instance, we previously reported that the covarion parameters, like the α parameter of the ras process, varied across different protein families. it will be interesting to model this heterogeneity in switch rates variation across sites and lineages and implement it in procov without increasing computational load too much. finally, the current release of procov only handles protein sequence data. analyses of dna substitutions under covarion models have found applications in inferring the evolutionary history of viral genes. future extension of procov to allow analyses of dna sequence data may be useful to investigate these kinds of data sets.  procov is a phylogenetic program to infer phylogeny under covarion models, which may be especially useful for problems involving estimates of deep divergences in the tree of life, where rates of evolution at sites are likely to have changed over the tree. it can also be used to detect covarion sites, which when combined with structural bioinformatics approaches, can be a powerful method to study lineage-specific functional shifts in protein families as well as protein adaptation. availability and requirements * project name: procov: maximum likelihood estimation of protein phylogeny under covarion models. * project home page:. * operating system: any unix/linux platform. * programming language: c. * other requirements: gcc or compatible compiler. it is recommended to have the blas/atlas libraries installed on the unix/linux system so that procov can run faster. versions of blas and lapack, such as the generic versions from atlas, netlib, or vendor-provided libraries that work with your compiler should be installed. the makefile should then be edited to match the type of the compiler and the path and library names of the blas and lapack libraries used. the makefile of the procov source code gives some instances of the blas installation on a few commonly-used unix systems. * license: gnu gpl. * any restrictions to use by non-academics: none. authors' contributions all authors conceived of the study, designed the method and wrote the article. h-cw implemented the software and ran the analyses. 
the covarion hypothesis of molecular evolution holds that selective pressures on a given amino acid or nucleotide site are dependent on the identity of other sites in the molecule that change throughout time, resulting in changes of evolutionary rates of sites along the branches of a phylogenetic tree. at the sequence level, covarion-like evolution at a site manifests as conservation of nucleotide or amino acid states among some homologs where the states are not conserved in other homologs. covarion-like evolution has been shown to relate to changes in functions at sites in different clades, and, if ignored, can adversely affect the accuracy of phylogenetic inference. covarion models implemented in procov may be especially useful for phylogenetic estimation when ancient divergences between sequences have occurred and rates of evolution at sites are likely to have changed over the tree. procov is a software tool that implements a number of previously proposed covarion models of protein evolution for phylogenetic inference in a maximum likelihood framework. those sites with the greatest log-likelihood difference between a'covarion' and an ras process were found to be of functional or structural significance in a dataset of bacterial and eukaryotic elongation factors. procov can be used to identify covarion sites by comparing the site likelihoods under the covarion process to the corresponding site likelihoods under a rates-across-sites process. several algorithmic and implementation improvements in this tool over previous versions make computationally expensive tree searches with covarion models more efficient and analyses of large phylogenomic data sets tractable. it can also be used to study lineage-specific functional shifts in protein families that result in changes in the patterns of site variability among subtrees.
 the covarion hypothesis of molecular evolution holds that selective pressures on a given amino acid or nucleotide site are dependent on the identity of other sites in the molecule that change throughout time, resulting in changes of evolutionary rates of sites along the branches of a phylogenetic tree. at the sequence level, covarion-like evolution at a site manifests as conservation of nucleotide or amino acid states among some homologs where the states are not conserved in other homologs. covarion-like evolution has been shown to relate to changes in functions at sites in different clades, and, if ignored, can adversely affect the accuracy of phylogenetic inference.  procov is a software tool that implements a number of previously proposed covarion models of protein evolution for phylogenetic inference in a maximum likelihood framework. several algorithmic and implementation improvements in this tool over previous versions make computationally expensive tree searches with covarion models more efficient and analyses of large phylogenomic data sets tractable. procov can be used to identify covarion sites by comparing the site likelihoods under the covarion process to the corresponding site likelihoods under a rates-across-sites process. those sites with the greatest log-likelihood difference between a'covarion' and an ras process were found to be of functional or structural significance in a dataset of bacterial and eukaryotic elongation factors.  covarion models implemented in procov may be especially useful for phylogenetic estimation when ancient divergences between sequences have occurred and rates of evolution at sites are likely to have changed over the tree. it can also be used to study lineage-specific functional shifts in protein families that result in changes in the patterns of site variability among subtrees. 