suggests that for all families with greater than two members and for all 23 non-redundant tissues, in general, exactly one family member is expressed highly, among the members with similar DNA binding site motif. Decrease versus increase in tissue-specific expression We next investigated the relative prevalence of expression increase and expression decrease in a tissue among members of a paralogous family, irrespective of their DNA binding site motif similarities. Using the 11 families with greater than two members identified above, for each TF family and each tissue we enumerated the cases for which all but one family member have high expression and one member has low expression. High expression is defined as an expression value that is one Median Absolute Deviation or more than the median expression for that particular tissue. Low expression is the complement of high expression. We say that a family has a selective decrease in expression if all but one of the family members has high expression. Conversely, a family is deemed to have a selective increase in expression if all but one of the family members has low expression. We found that 38% of all cases were of expression increase and 20% cases of expression decrease. As a control, if we randomize the expression data, we found 48% of the cases were of expression increase and 16% were of expression decrease. The one-tailed Fisher exact test comparing decrease versus increase of expression for the true and random sets yields a p-value of 0.02, indicating that there is an excess of expression decrease in the actual data. If we repeat the analysis using families derived from KOG data, we find that 35% of the cases were of expression increase and 40% were of expression decrease. By comparison in the cases generated by randomizing the expression data, we found that 40% of the cases were of expression increase and 24% of the expression cases were of expression decrease. The one-tailed Fisher exact test comparing decrease versus increase of expression for the actual and random sets yields a p-value of 8.5e-7. This result further reinforces our findings that there is a significant tendency towards expression decrease in the families of TFs. To further quantify the relationship of expression decrease versus expression increase we performed the following analysis. For each family of TF (F1, F2, .., Fk) and each tissue T, we computed the ratio RI = (Emax-Emax2)/(Emax-Emin), where Emax and Emin represent the maximum and the minimum expression levels and Emax2 represents the second highest expression among the family members in tissue T. The ratio RI captures the increase in expression, i.e., high value of RG represents an expression increase. Similarly, define RD = (Emin2-Emin)/(Emax-Emin), where Emin2 repre- sents the second lowest expression among the family members in tissue T. The ratio RD captures decrease in expression, i.e., high value of RL represents an expression decrease. We computed RI and RD values for all TF family-tissue pairs and did the same for 10000 randomly generated families of TFs of the same size as the true set of families for families defined using the BLAST criterion and families defined using KOG data. We found that the RI values (increase) in the actual data were significantly smaller than those for the randomized families (Mann-Whitney one-sided test p- value = 0.0023 for paralogs defined using BLAST criterion and p- value = 5.2e-4 for the families defined using KOG data). In contrast, the RD values (decrease) in the actual data were greater than those for the randomized families. However, the Mann- Whitney one-sided test p-value was at the cusp of significance (p- value = 0.06) for families defined using BLAST criterion. For families defined using KOG data, we found that RD values (decrease) were significantly less than those for randomized families (Mann-Whitney, one-sided test p-value = 3.4e-8. Thus, it Diversification of TF Paralogs appears that relative to background expectation, a decrease in expression appears to be more common and of a greater magnitude than an increase in expression in gene families. Discussion Gene duplication followed by functional diversification of the duplicated genes is a major driver of evolution[1–4]. However, there are evolutionary scenarios where duplicated genes are maintained if increased dosage is advantageous, e.g. in the case of immunity related genes[12,13]. Gene paralogs may functionally diversify through mutations in their coding sequences, thereby changing the activity of the gene product, or through mutations in the regulatory sequences, thereby altering the gene expression[28], or through both processes. The pathways of functional diversifi- cation and their interrelationships are of fundamental interest from the perspective of understanding adaptive expansion of gene families[2,29,30]. While mutations may alter either the coding or the regulatory sequence of a gene, it is the resulting functional change – in expression or in protein activity – that is likely to determine the gene’s evolutionary fate. Although it is often not possible to quantify the functional consequence of a mutation, TF genes offer a unique opportunity in this respect. By explicitly quantifying one aspect of functional diversification, viz., DNA binding site motifs, we have demonstrated that paralogs with conserved DNA binding site motifs tend to diverge in their tissue- specific expression and conversely, paralogs that are highly expressed in a tissue tend to have dissimilar DNA binding site motifs and thus, different target genes. Although we have observed an overall significant trend, the inverse relationship between the two modes of diversification is likely to apply to a subset of paralogs, and several other evolutionary processes are likely to be active. We have taken a number of precautions to ensure the robustness of our conclusions. For instance, Pearson’s correlation coefficient is a standard measure for quantifying correlations between two data samples. However, PCC is applicable only when the data are normally-distributed. Consequently, we have chosen the more appropriate but conservative Kendall’s tau to measure correlation between the tissue-specific expression divergence (hE) and DNA binding site motif similarity (§B). We preferred Kendall’s tau to another non-parametric correlation measure – Spearman’s rho, because Kendall’s tau has better statistical properties and is more directly interpretable [31] (see methods). However, using either Spearman’s rho or PCC does not affect our results appreciably and does not affect our conclusions. Also, because multiple paralogous pairs involve the same TF gene, the individual hE and §B values cannot be assumed to be independent. To avoid biases caused by within-sample depen- dence, we estimated the significance of Kendall’s tau by randomly permuting the expression values among the TFs and our results remain statistically significant. Thus, our overall results remain significant across a comprehensive set of combinations of experimental design parameters - (i) multiple ways of probe aggregation for a gene or a transcript (ii) multiple PWM aggregation for a gene or a transcript, (iii) multiple measures of correlation, and (iv) multiple ways of estimating significance. Prior investigations of functional diversification among paralogs have only considered sequence divergence in the entire coding portion of the gene[2,29,30]. Makova and Li found a negative correlation between expression divergence and protein sequence divergence in human paralogs[21]. To quantify expression divergence between paralogs, the prior works have used Pearson’s correlation coefficient or Euclidian distance (Jordan, et al 2005) using expression profiles of two genes across multiple contexts PLoS ONE | www.plosone.org 9 June 2008 | Volume 3 | Issue 6 | e2345