BACKGROUND
following speciation, there are many possible molecular events that can drive the divergence of species.
three of the most important mechanisms include changes in the coding sequence of proteins that alter protein function, changes in regulatory regions that affect gene expression, and changes in regulatory regions that affect mrna splicing.
the evolution of protein-coding sequences has been studied systematically in the adaptive evolution database , where such sequences were grouped into gene families  <cit> .
within these gene families, the ratio of nonsynonymous to synonymous nucleotide substitution rates  was used to detect an excess of nonsynonymous substitution, with positive selection as a proxy for potential functional change.
all cases of positive selection were mapped together from the gene tree to the species tree.
no systematic approach has been taken to examine relative gene expression or mrna splicing in the same way, partly because both appropriate methods and datasets are lacking.
one approach to examine the evolution of gene expression is to examine the substitution rate in promoters and look for lineages with excess substitution, analogous to ka/ks for protein coding sequences  <cit> .
this can then be correlated with relative expression levels.
an alternative approach is to reconstruct ancestral gene expression states and to examine lineages that show a significant change.
this has recently been implemented using a maximum likelihood approach for gene expression data  <cit> .
another evolutionary approach to reconstructing ancestral states is minimum evolution.
similar to the principle of parsimony, minimum evolution assumes that significant lineage-specific changes of gene expression through evolution are rarer than conservation of gene expression patterns.
naive comparison of gene expression values of genes across very closely related species does not discourage this assumption.
the use of minimum evolution methods to evaluate gene expression ancestral states can be likened to the use of parsimony to evaluate sequence ancestral states and the algorithm is modeled after fitch parsimony  <cit> .
as shown in figure  <dig>  examining evolution along branches improves the signal to noise ratio, compared with examining changes between extant sequences resulting in an analysis with more power to detect causative substitutions.
as with all methods related to parsimony, the method is expected to be most accurate with short branch lengths and well articulated  trees.
using the principle of minimum evolution, a general fast method has been developed that explicitly reconstructs the ancestral state of continuous character traits, like gene expression and mrna splicing.
the speed of this method will enable application to large datasets with many species and readily enables a subsequent mapping of data from gene expression trees to species trees.
another limitation towards extending taed-like approaches is the lack of applicable datasets.
for mrna splicing, comparisons of quantitative expressed sequence tag  data and genomic sequence data are used to evaluate relative splicing levels, but existing cross-species comparisons include very long branches  <cit> .
for gene expression, several datasets now exist including closely related species or isolates of yeast  <cit>  and primates  <cit> .
while these datasets are preliminary, they are a starting point to enable testing of methods.
here, we present our minimum evolution method, which is available as free software to download at  and test its performance on the cross-species primate dataset of enard et al.
<cit> .
methods
gene expression and sequence data
gene expression data was collected from enard et al.
<cit>  and contains samples from brain and liver of human, chimpanzee and orangutan.
sequence data was collected from ensembl  <cit>  and consists of the sequence  <dig> bp upstream of the gene transcription start site of the genes in the gene expression dataset.
the reference species tree was taken from arnason et al., an accepted phylogeny in the field  <cit> .
reconstruction of ancestral gene expression states
the reconstruction of continuous characters was done using a minimum evolution approach.
a range of values was obtained by running up and down a phylogeny and determining intervals consistent with minimum total evolution over a tree.
once the values converged on final intervals, the mid-point of the range was selected.
intervals for the first iteration are taken from those in the descendant leaves or nodes.
when they overlap the algorithm calculates the intersection between the intervals, but when they don't the algorithm constructs an interval range between the descendant intervals, as seen in figure  <dig>  the second iteration ran through the phylogeny to minimize the allowed intervals.
this was done by checking the upper and lower limit of the parent node interval and the interval of the node itself.
if the upper limit of the parent node was lower than the upper limit of the node itself or the lower limit of the parent node was higher than the lower limit of the node itself, then the limits of the node were changed to be the same as the parent node limits.
this is depicted as an example in figure  <dig>  with a corresponding sequence reconstruction depicted in figure  <dig>
ancestral sequence reconstruction
clustalw  <cit>  was used for multiple sequence alignment of promoter sequences and baseml from the paml package  <cit>  was used for ancestral sequence reconstruction.
along branch analysis
tess-transcription element search software  <cit>  was used to search the transfac database.
tess takes a candidate sequence as input and searches transfac for transcription factor binding sites that can be locally aligned with regions of the input sequence.
the output from tess consists of a list of transcription factor binding sites that match the input sequence and the position and length of the transcription factor.
these lists were manually controlled for each input sequence and correlated with the substitutional information calculated from paml.
promoters with more than 5% pairwise substitution between human and chimpanzee were discarded.
if a substitution occurred within a transcription factor binding site along any branch, it was annotated.
distributions were generated of the amount of along branch change in gene expression and ultimately the number of transcription factor binding sites from tess that were mutated along any branch.
this was normalized by the total amount of substitution to generate an enrichment value.
RESULTS
a software package that utilizes a minimum evolution algorithm to reconstruct ancestral states of continuous character data, like relative gene expression or alternative splicing levels and parse the amount of change to each branch of a phylogenetic tree is presented.
this software package is available for download at .
enard et al.
present an analysis of gene expression in a set of genes in brain and liver from human, chimpanzee with orangutan as an outgroup  <cit> .
using this dataset, we reconstructed ancestral gene expression values at the last common ancestor of human and chimpanzee.
the promoter sequences  for these genes from human, chimpanzee, and mouse as an outgroup were downloaded from ensembl  <cit> , aligned, and the last common ancestor sequence from human and chimpanzee was reconstructed using baseml from the paml package  <cit> , as described in the methods section.
while enhancers can regulate gene expression over long distances and can be critical to changes in gene expression, many important regulators of transcription are located in the  <dig> bp immediately upstream of the gene start site  <cit> .
while our knowledge of enhancer function does not permit a fully systematic analysis, analysis of promoter regions can be used to identify a non-exhaustive set of candidates.
the distribution of gene expression changes across branches is shown in figures  <dig> and  <dig> for the human and chimpanzee lineages, respectively.
the strong central peak was expected, given the conservative properties of the method.
the asymmetry of the distributions was not expected and may reflect problems with the original dataset.
if chimpanzee genes are hybridized to human sequences and then normalized to correct for substitution rates, this type of bias may be expected.
however, despite the unexpected shape of the distributions, there is still signal in the data, reflected in the significant enrichment values obtained.
the  <dig> % of genes along each lineage that were most upregulated and downregulated were collected with the middle 1% as a control.
this corresponded to  <dig> and  <dig> genes respectively.
approximately 35% of sequences were eliminated from the analysis because of >5% divergence between human and chimpanzee promoter sequences .
for the remaining sequences, the promoter sequence substitutions that occurred along each lineage were examined for known transcription factor binding sites using tess  <cit>  and the number of substitutions in such sites evaluated in each group, as shown in table  <dig>  the total substitution rate was also calculated and used to calculate an enrichment of substitutions in transcription factor binding sites along the lineages with the largest shifts in gene expression.
this enrichment, was seen in all four categories, with the human branches showing a better signal to noise ratio than the chimpanzee branches, which may have been expected given the distributions observed in figures  <dig> and  <dig> and the methodology for generating the data.
the supplementary materials  show the actual genes that have been implicated by this analysis, including the prospective transcription factor binding sites that have undergone substitution.
the dataset of genes is too small to pick out significant gene function signal from the upregulated and downregulated genes along each lineage.
along the human lineage, there were fewer substitutions predicted to destroy transcription factor binding sites on the up-regulated gene lineages compared with the control, while other lineage data were not different from the control.
because it is not clear which destroyed binding sites are normally occupied by transcriptional activators, it is difficult to interpret the biological significance of this result.
while the binding sites predicted may be candidates for playing an important role in the lineage-specific divergence of human and chimpanzee and warrant further testing for their activity in regulating expression from the respective promoters, little experimental data is currently available to further validate the study beyond the statistical validation seen in the enrichment values.
however, evolutionary approaches that consider along branch change as opposed to pairwise comparison of extant sequences  do hold promise in pinpointing substitutions that cause the divergence of gene expression during species diversification.
CONCLUSIONS
all together, a method  are made available for analysis of gene expression and alternative splicing shifts in a phylogenetic context and for detecting substitutions responsible for driving such shifts.
given some of the approximations made , the method performs surprisingly well and is a valuable starting point for this type of analysis, as well as being subject to future improvements.
ultimately, it will be valuable in comparative genomics to compare lineage-specific changes in gene content and in coding sequences, with changes in gene expression and alternative splicing to get a fuller picture of evolution.
authors' contributions
this project was conceived by dal.
all programming and analysis was performed by rr, under the supervision and with the technical support of dal and ie.
the paper was written by dal.
availability and requirements
project name: phyrex
project home page:
operating systems: linux
programming language: java
other requirements: java  <dig> .2
license: none
any restrictions to use by non-academics: none

 following speciation, there are many possible molecular events that can drive the divergence of species. three of the most important mechanisms include changes in the coding sequence of proteins that alter protein function, changes in regulatory regions that affect gene expression, and changes in regulatory regions that affect mrna splicing. the evolution of protein-coding sequences has been studied systematically in the adaptive evolution database, where such sequences were grouped into gene families. within these gene families, the ratio of nonsynonymous to synonymous nucleotide substitution rates was used to detect an excess of nonsynonymous substitution, with positive selection as a proxy for potential functional change. all cases of positive selection were mapped together from the gene tree to the species tree. no systematic approach has been taken to examine relative gene expression or mrna splicing in the same way, partly because both appropriate methods and datasets are lacking. one approach to examine the evolution of gene expression is to examine the substitution rate in promoters and look for lineages with excess substitution, analogous to ka/ks for protein coding sequences. this can then be correlated with relative expression levels. an alternative approach is to reconstruct ancestral gene expression states and to examine lineages that show a significant change. this has recently been implemented using a maximum likelihood approach for gene expression data. another evolutionary approach to reconstructing ancestral states is minimum evolution. similar to the principle of parsimony, minimum evolution assumes that significant lineage-specific changes of gene expression through evolution are rarer than conservation of gene expression patterns. naive comparison of gene expression values of genes across very closely related species does not discourage this assumption. the use of minimum evolution methods to evaluate gene expression ancestral states can be likened to the use of parsimony to evaluate sequence ancestral states and the algorithm is modeled after fitch parsimony. as shown in figure examining evolution along branches improves the signal to noise ratio, compared with examining changes between extant sequences resulting in an analysis with more power to detect causative substitutions. as with all methods related to parsimony, the method is expected to be most accurate with short branch lengths and well articulated trees. using the principle of minimum evolution, a general fast method has been developed that explicitly reconstructs the ancestral state of continuous character traits, like gene expression and mrna splicing. the speed of this method will enable application to large datasets with many species and readily enables a subsequent mapping of data from gene expression trees to species trees. another limitation towards extending taed-like approaches is the lack of applicable datasets. for mrna splicing, comparisons of quantitative expressed sequence tag data and genomic sequence data are used to evaluate relative splicing levels, but existing cross-species comparisons include very long branches. for gene expression, several datasets now exist including closely related species or isolates of yeast and primates. while these datasets are preliminary, they are a starting point to enable testing of methods. here, we present our minimum evolution method, which is available as free software to download at and test its performance on the cross-species primate dataset of enard et al.. methods gene expression and sequence data gene expression data was collected from enard et al. and contains samples from brain and liver of human, chimpanzee and orangutan. sequence data was collected from ensembl and consists of the sequence bp upstream of the gene transcription start site of the genes in the gene expression dataset. the reference species tree was taken from arnason et al., an accepted phylogeny in the field. reconstruction of ancestral gene expression states the reconstruction of continuous characters was done using a minimum evolution approach. a range of values was obtained by running up and down a phylogeny and determining intervals consistent with minimum total evolution over a tree. once the values converged on final intervals, the mid-point of the range was selected. intervals for the first iteration are taken from those in the descendant leaves or nodes. when they overlap the algorithm calculates the intersection between the intervals, but when they don't the algorithm constructs an interval range between the descendant intervals, as seen in figure the second iteration ran through the phylogeny to minimize the allowed intervals. this was done by checking the upper and lower limit of the parent node interval and the interval of the node itself. if the upper limit of the parent node was lower than the upper limit of the node itself or the lower limit of the parent node was higher than the lower limit of the node itself, then the limits of the node were changed to be the same as the parent node limits. this is depicted as an example in figure with a corresponding sequence reconstruction depicted in figure ancestral sequence reconstruction clustalw was used for multiple sequence alignment of promoter sequences and baseml from the paml package was used for ancestral sequence reconstruction. along branch analysis tess-transcription element search software was used to search the transfac database. tess takes a candidate sequence as input and searches transfac for transcription factor binding sites that can be locally aligned with regions of the input sequence. the output from tess consists of a list of transcription factor binding sites that match the input sequence and the position and length of the transcription factor. these lists were manually controlled for each input sequence and correlated with the substitutional information calculated from paml. promoters with more than 5% pairwise substitution between human and chimpanzee were discarded. if a substitution occurred within a transcription factor binding site along any branch, it was annotated. distributions were generated of the amount of along branch change in gene expression and ultimately the number of transcription factor binding sites from tess that were mutated along any branch. this was normalized by the total amount of substitution to generate an enrichment value.  a software package that utilizes a minimum evolution algorithm to reconstruct ancestral states of continuous character data, like relative gene expression or alternative splicing levels and parse the amount of change to each branch of a phylogenetic tree is presented. this software package is available for download at. enard et al. present an analysis of gene expression in a set of genes in brain and liver from human, chimpanzee with orangutan as an outgroup. using this dataset, we reconstructed ancestral gene expression values at the last common ancestor of human and chimpanzee. the promoter sequences for these genes from human, chimpanzee, and mouse as an outgroup were downloaded from ensembl, aligned, and the last common ancestor sequence from human and chimpanzee was reconstructed using baseml from the paml package, as described in the methods section. while enhancers can regulate gene expression over long distances and can be critical to changes in gene expression, many important regulators of transcription are located in the bp immediately upstream of the gene start site. while our knowledge of enhancer function does not permit a fully systematic analysis, analysis of promoter regions can be used to identify a non-exhaustive set of candidates. the distribution of gene expression changes across branches is shown in figures and for the human and chimpanzee lineages, respectively. the strong central peak was expected, given the conservative properties of the method. the asymmetry of the distributions was not expected and may reflect problems with the original dataset. if chimpanzee genes are hybridized to human sequences and then normalized to correct for substitution rates, this type of bias may be expected. however, despite the unexpected shape of the distributions, there is still signal in the data, reflected in the significant enrichment values obtained. the % of genes along each lineage that were most upregulated and downregulated were collected with the middle 1% as a control. this corresponded to and genes respectively. approximately 35% of sequences were eliminated from the analysis because of >5% divergence between human and chimpanzee promoter sequences. for the remaining sequences, the promoter sequence substitutions that occurred along each lineage were examined for known transcription factor binding sites using tess and the number of substitutions in such sites evaluated in each group, as shown in table the total substitution rate was also calculated and used to calculate an enrichment of substitutions in transcription factor binding sites along the lineages with the largest shifts in gene expression. this enrichment, was seen in all four categories, with the human branches showing a better signal to noise ratio than the chimpanzee branches, which may have been expected given the distributions observed in figures and and the methodology for generating the data. the supplementary materials show the actual genes that have been implicated by this analysis, including the prospective transcription factor binding sites that have undergone substitution. the dataset of genes is too small to pick out significant gene function signal from the upregulated and downregulated genes along each lineage. along the human lineage, there were fewer substitutions predicted to destroy transcription factor binding sites on the up-regulated gene lineages compared with the control, while other lineage data were not different from the control. because it is not clear which destroyed binding sites are normally occupied by transcriptional activators, it is difficult to interpret the biological significance of this result. while the binding sites predicted may be candidates for playing an important role in the lineage-specific divergence of human and chimpanzee and warrant further testing for their activity in regulating expression from the respective promoters, little experimental data is currently available to further validate the study beyond the statistical validation seen in the enrichment values. however, evolutionary approaches that consider along branch change as opposed to pairwise comparison of extant sequences do hold promise in pinpointing substitutions that cause the divergence of gene expression during species diversification.  all together, a method are made available for analysis of gene expression and alternative splicing shifts in a phylogenetic context and for detecting substitutions responsible for driving such shifts. given some of the approximations made, the method performs surprisingly well and is a valuable starting point for this type of analysis, as well as being subject to future improvements. ultimately, it will be valuable in comparative genomics to compare lineage-specific changes in gene content and in coding sequences, with changes in gene expression and alternative splicing to get a fuller picture of evolution. authors' contributions this project was conceived by dal. all programming and analysis was performed by rr, under the supervision and with the technical support of dal and ie. the paper was written by dal. availability and requirements project name: phyrex project home page: operating systems: linux programming language: java other requirements: java.2 license: none any restrictions to use by non-academics: none 
systematic analysis of lineage-specific evolution is becoming the cornerstone of comparative genomics. new methods, like phyrex, extend the capabilities of comparative genomics by tracing the evolution of additional biomolecular processes. a software package for reconstruction of continuous character traits, like relative gene expression levels or alternative splice site usage data is presented and is available for download at. as genomes evolve after speciation, gene content, coding sequence, gene expression, and splicing all diverge with time from ancestors with close relatives. this program was applied to a primate gene expression dataset to detect transcription factor binding sites that have undergone substitution, potentially having driven lineage-specific differences in gene expression.
 as genomes evolve after speciation, gene content, coding sequence, gene expression, and splicing all diverge with time from ancestors with close relatives. a minimum evolution general method for continuous character analysis in a phylogenetic perspective is presented that allows for reconstruction of ancestral character states and for measuring along branch evolution.  a software package for reconstruction of continuous character traits, like relative gene expression levels or alternative splice site usage data is presented and is available for download at. this program was applied to a primate gene expression dataset to detect transcription factor binding sites that have undergone substitution, potentially having driven lineage-specific differences in gene expression.  systematic analysis of lineage-specific evolution is becoming the cornerstone of comparative genomics. new methods, like phyrex, extend the capabilities of comparative genomics by tracing the evolution of additional biomolecular processes. 