Pairwise comparisons across species are problematic when analyzing functional genomic data

KMRR (16) sought to test the ortholog conjecture. The ortholog conjecture (17) is the proposition that orthologs (genes that diverged from each other because of a speciation event) have more similar attributes than paralogs (genes that diverged from each other because of a gene duplication event). The ortholog conjecture can and has been applied to diverse attributes, including molecular sequence, biochemical function, and as in the study considered here, expression. It has important biological and technical implications. It shapes our understanding of the functional diversity of gene families. It is also used to relate findings from well-studied genes to homologous genes that have not been investigated in detail. While the ortholog conjecture describes a specific pattern of functional diversity across genes, it is also articulated as a hypothesis about the process of evolution—that there is greater evolutionary change in gene attributes after a duplication event than a speciation event.

Despite its importance, there have been relatively few tests of the ortholog conjecture. Previous work has shown that ontology annotations are not sufficient to test the ortholog conjecture (22, 23). Analyses of domain structure were consistent with the ortholog conjecture (24). There have been few tests of the ortholog conjecture with regards to gene expression (22), and the study by KMRR (16) is one of the most detailed to date.

KMRR (16) considered several publicly available datasets of gene expression across tissues and species. Their expression summary statistic is tau (25, 26), an indicator of tissue specificity of gene expression. Tau can range from a value of zero, which indicates no specificity (i.e., uniform expression across tissues), to a value of one, which indicates high specificity (i.e., expression in only one tissue). Tau is convenient in that it is a single number of defined range for each gene, although of course, since the original expression is multidimensional, this means much information is discarded. This includes information about the tissue to which expression is specific. For example, if one gene has expression specific to the brain and another expression specific to the kidney, both would have a tau of one.

The analyses by KMRR (16) are based on pairwise comparisons (Fig. 1A) between tau within each gene family. Rather than make every pairwise comparison within each gene tree, they considered only a subset of pairwise comparisons in each particular analysis. They first selected a focal species, which varied from analysis to analysis. Ortholog comparisons were limited to pairs that include this species, and the only paralogs considered were those with the highest expression in this species. Note that this subset of pairwise comparisons still samples the same changes multiple times.

They found the correlation coefficient of tau for orthologs to be significantly greater than the correlation coefficient of tau for paralogs (i.e., orthologs tend to have more similar expression than paralogs). From this, they concluded that their analyses support the ortholog conjecture. They also concluded that this pattern provides support for a particular evolutionary process: that “tissue-specificity evolves very slowly in the absence of duplication, while immediately after duplication the new gene copy differs” (16).

We reanalyzed the study by KMRR (16) using phylogenetic comparative methods. We focused on one of the datasets included in their analyses: that of Brawand et al. (27). This dataset is the best sampled in their analyses. It has gene expression data for six organs across 10 species (nine mammals and one bird), 8 of which were analyzed by KMRR (16) and are further considered here.

Many phylogenetic comparative methods are now available for addressing a wide range of questions relevant to comparative functional genomics. Some recent methods have been developed specifically for analyses of expression. For example, the Expression Variance and Evolution model compares the ratio of within-species expression variance with between-species evolutionary expression variance in an explicit phylogenetic framework (10). It can test a variety of hypotheses, including lineage-specific expression level shifts. This method considers strict orthologs that all have the same gene tree, which is the same as the phylogeny of the species under consideration. However, for the reanalysis presented here, we address different questions that consider a broad diversity of gene trees that include both speciation and different patterns of duplication. Phylogenetic independent contrasts (PICs) (2), the most widely used phylogenetic comparative method, are particularly well-suited to this challenge.

For each internal node in each gene tree, we calculated the PIC (2) of tau. This is the difference in values of tau for descendant nodes scaled by the expected variance, which is largely determined by the lengths of the two branches that connect the node to its two descendants (Fig. 1B). These contrasts were then annotated by whether each is made across a speciation or duplication event. The original description of independent contrasts (2) focused on assessing covariance between changes in two traits. Our use of contrasts is a bit different—we look for differences in evolutionary changes of one trait (differential expression) between two categories of nodes (speciation and duplication). This is similar to previous applications of contrasts to examine shifts in evolutionary rates between clades (28). Rather than compare the distributions of contrasts for nodes in different clades, as this previous work did, we compare the distributions of contrasts for nodes with different annotations (speciation or duplication).

We mapped the tau values calculated by KMRR (16) for the dataset by Brawand et al. (27) onto 21,124 gene trees parsed from ENSEMBL Compara (29). These are the same precomputed trees on which the orthology/paralogy annotations that KMRR (16) used are based; 8,854 gene trees passed taxon sampling criteria (four genes) after removing tips without tau values and had at least one speciation event. Of these, 8,516 were successfully time calibrated. These calibrated trees were used to calculate PICs for 21,017 duplication nodes and 67,845 speciation nodes (Tables S1 and S2). One of these trees is presented in Fig. S1 to show the analysis.

It is essential to have a null hypothesis that makes a distinct prediction from the prediction of the hypothesis under consideration. A suitable null hypothesis in this case is that there is no difference in the evolution of expression after speciation or duplication events (30). Under this hypothesis, we would predict that contrasts across speciation nodes and duplication nodes are drawn from the same distribution. Under the alternative hypothesis specified by the ortholog conjecture (that there is a higher rate of change after duplication events than speciation events), we would expect to see the distribution of duplication contrasts shifted to higher values relative to the speciation contrasts.

We did not find increased evolutionary change in expression after duplication events compared with speciation events (Fig. 2B). The Wilcoxon rank test does not reject the null hypothesis that the rate of evolution after duplications is the same as or less than the rate after speciation ( $P=1$). Our phylogenetic comparative analysis, unlike the previously published pairwise comparative analysis (16), therefore finds no support for the ortholog conjecture in this system. This result holds when we exclude the genes with the lowest phylogenetic signal according to Blomberg’s $K$ (31) (Fig. S2) and is robust to model selection (32–34) (SI Text).

We examined the possibility that ascertainment biases were differentially impacting the inference of expression evolution after duplication and speciation events. Such a bias might obscure support for the ortholog conjecture. We focused on two possible sources of bias: node depth and branch length. We found no evidence that either affected our results (Fig. S3). We also examined the sensitivity of the results to the calibration times applied to speciation events on the gene trees. This is important, because it is expected that genes from separate species have a common ancestor older than the time at which the species diverged from each other (35, 36). There is also uncertainty associated with the timing of these speciation events. We added random noise to the calibration times in replicate analyses, and all still failed to reject the null hypothesis (SI Text).

An additional concern is that within-species expression variation is not considered when making comparisons of tau across homologous genes. Neglecting within-species variation is expected to mislead phylogenetic comparative analyses in some cases (10, 37, 38). Since branch length describes expected trait variation (2), within-species variation can be modeled by extending the terminal branches of each gene phylogeny (37). We modeled within-species variation by extending each terminal branch by 7 My. This value was selected, because it corresponds to the age of the clade Hominini (defined by the most recent common ancestor of humans and chimpanzees), the shallowest node in the phylogeny of examined species. This provides an opportunity to assess the impact that extreme within-species variation (i.e., as large as the variance between humans and chimpanzees and their most recent common ancestor) would have on our findings. These analyses still do not reject the null hypothesis and do not provide support for the ortholog conjecture (Fig. S4).

To better understand why our phylogenetic analysis supports a different conclusion (i.e., no support for the ortholog conjecture) than the published analysis by KMRR (16) (i.e., strong support for the ortholog conjecture), we first checked to make sure that we could reproduce their result based on pairwise analyses. This is important, since we are only looking at a subset of the data that they considered: the dataset by Brawand et al. (27) for gene trees that could be successfully time calibrated. In figure 1 in ref. 16, they present a higher tau correlation coefficient between ortholog pairs than between paralog pairs. We find the same here, with correlation coefficients of 0.75 ( $P<2\times {10}^{-16}$) for orthologs and 0.36 ( $P<2\times {10}^{-16}$) for paralogs.

Why is it that pairwise methods and phylogenetic methods lead to opposite conclusions? One reason is that multiple pairwise comparisons repeatedly sample the same evolutionary changes and in so doing, violate statistical assumptions of independence, whereas phylogenetic comparative methods make multiple independent comparisons across nonoverlapping branches of the tree. The other reason is that pairwise comparisons and phylogenetic comparative methods describe different things. Pairwise comparisons describe contemporary patterns, while phylogenetic methods infer historical processes (15). Paralogs can be more different from orthologs, even when the processes of evolution after speciation and duplication are the same. Any difference could be caused by the structure of the gene phylogenies alone. If paralogs tend to be more distantly related to each other than orthologs, then there would be more time for differences to accumulate, even if the rate of change is the same between the two. This is, in fact, the case for these data. While the mean distance (i.e., total branch length) between orthologs is 332.6 My, the mean distance between paralogs is 1,763.5 My. This is because the oldest speciation event is, by definition, the most recent common ancestor of the species included in the study, but many gene families underwent duplication before this time.

To test the hypothesis that ancient duplications that precede the oldest speciation event (Fig. S3A) impact the lower correlation of tau between paralogs than between orthologs, we removed them. When we consider only the duplication events the same age or younger than the oldest speciation event (Fig. S3B), the paralog correlation coefficient increases from 0.36 to 0.55 ( $P<2\times {10}^{-16}$). This is much closer to the ortholog correlation of 0.75.

The study by KMRR (16) did investigate the impact of node age on correlation but in a different way. In figure 2 in ref. 16, they grouped orthologs and paralogs according to the ENSEMBL node name of their most recent common ancestor and plotted the correlation of tau for each of these groups by the node age. They found that, across the investigated range of node ages, ortholog pairs have higher tau correlation than paralogs. We confirmed that we can replicate this result (Fig. 2A). There are, however, several difficulties with interpreting this plot. First, it does not just reflect the evolutionary processes that generated the data, but it is also impacted by the phylogenies along which these processes acted. The expected covariance of traits that evolve under neutral processes is, in fact, defined by the phylogeny (3). Second, the correlation for each group is based on multiple nonindependent pairwise comparisons. Third, instead of using the actual age of duplication events, the age of an adjacent named node in the species tree is used.

To better understand this plot (Fig. 2A), we performed simulations of tau on the calibrated gene trees. We did not modify the gene tree topologies or their inferred histories of duplication and loss. We simulated the evolution of tau under the null model that it evolves at the same rate after duplication and speciation events. Under the null model, the plot of correlation coefficient to node age (Fig. 2C) is very similar to that of the observed data (Fig. 2A). As in the original study, there is higher correlation coefficient across orthologs (0.73, $P<2\times {10}^{-16}$) than paralogs (0.28, $P<2\times {10}^{-16}$) when not considering node age. Phylogenetic analysis of the data simulated under the null hypothesis (Fig. 2D) does not reject the null hypothesis (Wilcoxon $P=0.999$) as expected.

We next simulated the evolution of tau under the ortholog conjecture, where the rate of evolution of tau after duplication was twofold the rate after speciation. The pairwise results of this heterogeneous model (Fig. 2E) are nearly indistinguishable from the results under the null model (Fig. 2C) and also, have a higher correlation coefficient for orthologs (0.76, $P<2\times {10}^{-16}$) than paralogs (0.24, $P<2\times {10}^{-16}$). The phylogenetic analysis of the ortholog conjecture simulation (Fig. 2F) does reject the null hypothesis (Wilcoxon $P=1.1\times {10}^{-53}$) as expected if the phylogenetic methods are working correctly.

These simulations have several implications. The pairwise comparisons used by KMRR (16) cannot distinguish between the null hypothesis and ortholog conjecture. The pairwise results are strikingly similar under both hypotheses (Fig. 2 C and E). These simulations also serve to validate the phylogenetic methods applied to this problem. As expected, our phylogenetic analysis of independent contrasts does not reject the null hypothesis when data are simulated under the null model and does reject the null hypothesis when the data are simulated under the ortholog conjecture. In contrast to pairwise methods, the phylogenetic analyses can test explicit predictions based on hypotheses about evolutionary process.

Since greater rates of expression evolution after duplication do not explain the lower correlation of tau values in the pairwise comparison plots (Fig. 2 A, C, and E), this pattern must reflect some other property that is shared across these analyses. We found that the lower correlations for duplication comparisons can largely be explained by the greater variance in the age of duplication nodes. Each point in the pairwise correlation plot (Fig. 2 A, C, and E) summarizes multiple pairwise comparisons across nodes of a given event type (speciation or duplication as indicated by color) and clade in the species tree (Theria, Mammalia, Amniota, etc. as indicated by the age of the clade along the x axis). These clade names, in turn, are derived from the node annotations in the ENSEMBL Compara trees, which apply clade names to both speciation and duplication nodes. While speciation nodes in the gene trees have a clear correspondence to clades in the species trees, duplication nodes do not, since duplication can occur at any point along any branch in the species tree. Neighboring clade ages are, therefore, a very rough approximation of the age of duplication events. This is apparent when node ages derived from calibrated trees are plotted against the age of the clade annotations for each node (Fig. S5A). Because duplication events that are all annotated with the same clade name can occur at very different times, pairwise comparisons across these nodes capture evolutionary changes in tau across very different branch lengths. They, therefore, have lower correlation than pairwise comparisons made across speciation nodes. Of the 960,481 duplication nodes considered here, 264,300 have a node age (as determined by the time-calibrated trees) within 10% of their clade annotation. If only these duplication events closest to their annotated clade nodes are retained, the pattern of lower correlation of tau evolution for duplication events disappears for the dataset simulated under the null model (Fig. S5B). This correction, however, comes at the high cost of discarding 72.5% of duplication nodes.

There has been considerable recent interest in, and controversy about, the ortholog conjecture (17, 18, 30, 39, 40). While some studies have presented support for the ortholog conjecture, our results are consistent with multiple studies that have not (17, 30, 41). It should be noted that many areas of biology restrict the term “function” to describe the biochemical function of a gene, but the literature on gene duplication often uses function in a broader sense to include any gene attribute that could impact fitness, such as biochemical function, spatial expression, or binding. Since the ortholog conjecture can be applied to such diverse attributes, it is a heterogeneous and potentially idiosyncratic topic. For example, a recent integrated analysis of multiple gene attributes after duplication found extensive divergence in protein sequence but less divergence in expression (42). In our analyses, we found that phylogenetic distance is a better predictor of the similarity of gene expression than the history of gene duplication and speciation. The ortholog conjecture does not have to be an all or nothing proposition, even as it applies to expression. It may be that the rates of tissue-specific expression evolution after duplication are greater in some organisms, gene families, and evolutionary processes (39). The lack of support that we find for the ortholog conjecture in this system as well as the mixed results of others do suggest that investigators should test for it in each situation rather than assume a priori that it is a dominant pattern in the diversity of gene expression. These tests should be done in an explicit phylogenetic framework.

While empirical support for the ortholog conjecture has been mixed, theoretical work has suggested several mechanisms that could lead to increased rates of evolutionary change in expression and other gene traits after gene duplication. This has been a primary motivation for the extensive discussion of the ortholog conjecture. In his seminal work, Ohno (43) outlined three outcomes of gene duplication, now often referred to as neofunctionalization (one gene copy can take on new functions, while the other retains the old functions), subfunctionalization (each gene copy has a subset of the functions present in their most recent common ancestor), and conservation (each gene copy retains the same function). Neofunctionalization was the primary focus of early work on gene duplication, although it was not found to be a widespread outcome of duplication (44). The highly influential duplication–degeneration–complementation (DDC) model (45, 46) suggested that subfunctionalization could be a better explanation for the retention of both copies of a gene. It is built on the idea that gene copies that have each lost different functions would both be needed to fulfill all of the functions present in their most recent common ancestor. Tests of the DDC model have had mixed results (47). A growing body of work suggests that retention and differentiation of genes are caused by a variety of different mechanisms in different contexts and that we should not necessarily expect one process to dominate over others (48). Our lack of support for the ortholog conjecture suggests that whatever mechanisms lead to the retention of duplicated genes in the organisms considered here do not lead to changes in the tissue specificity of gene expression as summarized by tau.

The DDC model is particularly interesting to further consider in this context, since the summary statistic tau describes tissue specificity of expression. In the context of tissue-specific expression, the ortholog conjecture predicts a greater rate of evolutionary change in tau after duplication without specifying the direction of the change. The DDC model predicts a specific type of change in expression after duplication—subfunctionalization, which would increase tissue specificity, and therefore tau, after duplication. To test for this more specific prediction of DDC, we conducted additional analyses to see if there is a greater rate of increase in tau after duplication. For each branch in each gene tree, we inferred the change in tau and standardized this change by branch lengths. We found no evidence of increased tau after duplication relative to speciation (Wilcoxon $P=0.846$). This result also holds when we discarded the 75% of genes with the lowest phylogenetic signal (Wilcoxon $P=0.715$). These results are consistent with another recent study that found little support for neofunctionalization or subfunctionalization in expression of duplicated human genes and instead, proposed that dosage compensation drives most changes in expression after duplication in mammals (49). Like the ortholog conjecture, DDC is not an all or nothing hypothesis. We fail to find support for it in these tissues in these species with the tau expression summary statistic, but it may hold in other contexts. DDC could still play a role in the retention of these genes but with respect to other attributes, such as biochemical functions, rather than tissue-specific expression.

Levin et al. (19) analyzed gene expression through the course of embryonic development for 10 animal species, each from a different clade that has been designated as having the rank of phylum. They arrived at two major conclusions. First, animal development is characterized by a well-defined middevelopmental transition that marks the transition from an early phase of gene expression to a late stage of gene expression. Second, this transition helps explain the evolution of features observed among distantly related animals. Specifically, they concluded that animals from different phyla exhibit an inverse hourglass model for the evolution of gene expression, where there is more evolutionary variance in gene expression at a midphase of development than there is at early and late phases. Closely related animals have previously been described as having an hourglass model of gene expression, where evolutionary variance in expression is greater early and late in development than at the midpoint of development (20, 50). Levin et al. (19) conclude that this contrast between distantly and closely related animals provides biological justification for the concept of phyla and may provide a definition of phyla.

Levin et al. (19) arrived at this conclusion by making multiple pairwise comparisons of ortholog expression data sampled throughout the course of embryonic development. For each species pair, they identified the orthologs shared by these species. This list of shared genes was different from species pair to species pair. They characterized each of these orthologs in each species as having expression that peaks in early, mid, or late temporal phase of development. They then calculated a similarity score for each temporal phase for each species pair based on the fraction of genes that exhibited the same patterns in each species. The distributions of similarity scores are plotted in figure 4D in ref. 19, and their Kolmogorov–Smirnov (KS) tests indicated that the early distribution and late distribution were each significantly different from middistribution ( $P<{10}^{-6}$ and $P<{10}^{-12}$, respectively). This is the support that they presented for the inverse hourglass model.

We examined the matrix of pairwise comparisons used as the base for the KS tests and figure 4D in the work by Levin et al. (19) and thus, as support for the inverse hourglass model. We found several problems resulting from the use of multiple pairwise comparisons. The first problems are specific to this particular implementation of pairwise comparisons. We found that every data point was included twice, because both reciprocal pairwise comparisons (which have the same values) were retained. As a consequence, there are 90 entries for the 45 pairwise comparisons, and by doubling the data, the significance of the result seems stronger than it actually is. After removing the duplicate values, the P values are far less significant: 0.002 for the early to middle comparison and on the order of ${10}^{-6}$ for early to late. In addition, the test that they used (the KS test) is not appropriate for the hypothesis that they seek to evaluate. The KS test does not just evaluate whether one distribution is greater than the other; it also tests whether the shapes of the distributions are the same. The samples in this dataset are matched (i.e., for each pairwise comparison, there are early, mid, and late expression values), which the KS test does not take into account. The paired version of the Wilcoxon test is instead appropriate in this case. When applied to the deduplicated data, the $P$ value of this test is 0.02 for the early to middle comparison and on the order of ${10}^{-7}$ for early to late comparison.

After we addressed the issues above, we were able to explore more general issues that can be a problem when making multiple pairwise comparisons between species. We found that all five of the lowest values in the midphase distribution (Fig. 3A) are for pairwise comparisons that include the ctenophore (comb jelly). When the nine pairwise comparisons that include the ctenophore are removed, there is no significant difference between the early-phase and midphase distributions ( $P=0.14$ for the early to middle comparison and $P<{10}^{-5}$ for the late to middle comparison) and no support for the inverse hourglass (Fig. 3B). This highlights a well-understood property of pairwise comparisons across species (2, 51): evolutionary changes along a given branch, like those along the ctenophore branch, impact each of the multiple pairwise comparisons that include that branch. The pairwise comparisons are, therefore, not independent—different pairwise comparisons are impacted by changes along some of the same branches (Fig. 1A). This can give the impression of a general pattern across the tree that is instead specific to changes along one part of the tree. The number of comparisons impacted by each change depends on the structure of the phylogenetic tree (i.e., how the species are related to each other). Phylogenetic comparative methods were developed specifically to address this problem (2).

While we show problems with pairwise comparisons that impact the analysis by Levin et al. (19), we did not perform a phylogenetic reanalysis of this study as we did for the study by KMRR (16). This is because the similarity metric computed in the pairwise comparisons of Levin et al. (19) is not suitable for phylogenetic analysis, as it is based on different genes for different species pairs, and therefore, it cannot be modeled across the phylogeny. A full phylogenetic reanalysis would be possible using upstream analysis products to rederive new expression summary statistics.

We are not the first to apply phylogenetic comparative methods to functional genomic data. While the vast majority of comparative functional genomic studies have used standard pairwise similarity methods, a small number of comparative functional genomic studies have used phylogenetic comparative approaches (52–55). For instance, a phylogenetic ANOVA (10) of the evolution of gene expression in strict orthologs improves statistical power and drastically reduces the rate of false positives relative to pairwise approaches.

Some of the most widely used phylogenetic comparative methods (2, 3) are already directly applicable to comparative functional genomic studies. There are also interesting new challenges, such as parameterizing differential expression, so that species-specific technical biases are not mistaken for evolutionary changes in expression (56).

Branch lengths are fundamental to phylogenetic comparative methods, since the null expectation is that there are more changes along longer branches. In this study, we applied four separate modifications of branch lengths. First, we time calibrated the gene trees, so that nodes corresponding to the same speciation events had the same age across all trees. This is critical to making the comparisons of functional genomic traits equivalent across genes. Second, PICs (including the implementation that we use here) adjust internal branch lengths to accommodate the variation associated with estimating character states at internal nodes (2). Third, we randomized the calibration times to determine the sensitivity of the analyses to calibration times and the branch length estimations that they impact (SI Text). Fourth, we extended terminal branch lengths to understand the potential impact of within-species variation (Fig. S4). The first two of these modifications are technical requirements for the application of phylogenetic comparative methods. The second two explore the sensitivity of the analyses to branch length and the information that they convey and find that the results are robust to variation in branch length. This is encouraging for the future of phylogenetic comparative functional genomics given that one concern of applying these methods is the challenge of estimating branch lengths.