(tissues in case of Makova and Li). However, most genes, particularly the TF genes, are likely to be functionally relevant only in a few tissues or experimental conditions[32]. Thus, a context-specific analysis of expression divergence is more infor- mative. Our finding in human is consistent with that of Makova and Li even though we measure expression divergence (hE) in a context-specific manner. In fact when we replace tissue-specific hE value with an overall measure of expression divergence using Pearson’s correlation coefficient between the expression profiles across the 79 (or 23) tissues, we do not observe a significant correlation with DNA binding site motif similarity as shown in Figure 1 in supplementary Data file S1. This phenomenon is not surprising because it is the absolute expression level in specific contexts that is relevant to evolution. While it is well-established that duplicated genes have undergone rapid divergence in expression immediately after duplication[29], our work refines this general observation by showing that the expression divergence may be influenced by divergence in other aspects of gene’s function, for instance, DNA binding site motif of the TF protein. Given the multitude and complexity of attributes affecting evolution, most of which are not completely understood, we expect that most individual correlations will be weak. Strong correlations are indeed scarce in biological literature. Liao et al. studied the correlations between evolutionary rates of mouse paralogs and several genic parameters and remarked that ‘‘because there are potentially many rate determinants …… it is not unexpected that the observed correlation coefficients are not very high’’[33]. There are some concerns with using the expression data in investigating functional diversification, or indeed in any functional study. First, the expression level serves only as a proxy for the amount of active transcription factor protein. However, an accurate genome-wide quantification of active transcription factors in various tissues is currently not feasible. Another concern is that when two paralogs have very similar expression, their corresponding probes on the microarray are likely to cross-hybridize, thus confounding the results. However, the possibility of a cross-hybridization serves to make our finding more conservative. Expression pattern and the DNA binding site motif represent two distinct aspects of TF function. There are others, such as interaction partners, and these distinct functional aspects are often encoded in distinct sequence domains. For instance, the DNA binding site motif is encoded largely within the DNA binding domain of the gene. In the human Forkhead (FKH) domain containing family of TFs, the sequence similarity in DNA binding FKH domain is significantly correlated with the DNA binding site motif similarity (R = 0.47, P = 9.6e-05; data not shown). Because selection pressure operates at the level of function, which is encoded in distinct domains, it is reasonable to study the relationship between diversification in distinct domains. Again, in the human FKH family, where the DNA binding domain is highly conserved across family members, the rest of the protein domains exhibit an accelerated divergence relative to the DNA binding domain immediately after duplication (data not shown). A specific comparison of divergence in the TF’s interaction domain with either the divergence in DNA binding similarity or the expression divergence or both, would be a natural extension of the current work. However, the interaction domains of TF proteins are currently not as well characterized as the DNA binding domains. As a proxy for divergence in interaction domain, one can consider the overlap in known interaction partners. However, such an investigation is currently limited by the availability of complete interaction data. Our work generalizes the previous studies of diversification of paralogs that have focused on divergence in protein coding sequence in its entirety without distinguishing Diversification of TF Paralogs among functional protein domains and without quantifying the functional consequence of mutations in the coding sequence. As a pragmatic concern, current genomic approaches to analyzing transcriptional regulation are confounded by the fact that multiple TFs, typically closely related paralogs, bind to similar binding sites. For instance, previous analysis of motifs enriched in promoters of genes that are differentially expressed in adult failing hearts identified the FKH family. However, to implicate specific members of the Forkhead family of TFs, directed experiments such as PCR and immunohistochemistry in specific cell types are needed[34]. Our results here not only explain in part the apparent redundancy in the DNA binding specificity of transcription factors – TFs with similar binding have dissimilar expression pattern – but also underscore the importance of incorporating the expression of the transcription factor genes in the analysis of transcriptional regulation. Materials and Methods TF proteins and their DNA binding specificities A total of 390 distinct human TFs were obtained from TRANSFAC 10.2[35]. Using the annotation table provided by TRANSFAC, PWM identifiers were mapped to 297 unique UniProt accession identifiers. These UniProt identifiers were then mapped to 295 unique Ensembl gene identifiers using Ensembl-47 (www.ensembl.org). For each of the TFs, the corresponding positional weight matrix (PWM) representing their DNA binding site motif was obtained from TRANSFAC. Identifying paralogous protein pairs As in several previous studies, a BLAST-based approach was used as one of the possible methods for determining paralogy[36]. BLAST was performed by querying human TF protein sequences from a database of all human protein sequences obtained from Ensembl-47. Two TFs were deemed to be paralogous if (i) the BLAST E value#1.0e-5, (ii) the BLAST-alignable region is $70% of the longer protein or the length of the longest High Scoring Pair (HSP) is .150 aa, and (iii) the identity (I) is $ 5% if the alignable region is longer than 150 aa or I$0.06n+4.8L20.32[1+exp(2L/1000)] where L is the length of alignable region. A similar set of criteria was used in [21]. One difference is that Makova and Li use 80% as the threshold for the coverage of alignable region while we use 70%. In addition, we also use a more stringent E value for identifying significant HSPs (E#1.0e-5 vs. E#1 in Makova nd Li). Based on these criteria we identified 95 paralogous TF gene pairs. Our criteria for paralogs, as for the previous works, are stringent and detect relatively recent paralogs. There are two difficulties in analyzing distant paralogs. First, inference of paralogy will be less accurate and second, the divergence in expression and binding site motifs are likely to have reached saturation and thus, we do not expect to detect a significant correlation between the two quantities. In fact our criteria is slightly less stringent that that of Makova and Li. When we use identical criteria, despite a reduction in the number of paralogous pairs, the results do not change appreciably. Similarity in DNA binding site motif for TF pairs The likelihood of two DNA binding site motifs, or PWMs, having identical DNA sites can be approximated using a number of pair- wise PWM similarity measures. Mahoney et al. have reported a detailed benchmarking study of several measures of PWM-PWM similarity[23]. We utilize the best performing Pearson correlation (PCC) based measure (provided by Shaun Mahoney). Very briefly, this particular method measures the similarity between two columns of the two PWMs using the Pearson correlation coefficient. Then using the PCC as the column ‘match’ score, the method aligns the PLoS ONE | www.plosone.org 10 June 2008 | Volume 3 | Issue 6 | e2345