Conserved Genes, Sampling Error, and Phylogenomic Inference

Disagreement or conflict among phylogenetic hypotheses obtained by analysis of large, genome-wide data sets has incited debate over potential benefits, pitfalls, and best practices associated with phylogenomic approaches ((Jeffroy et al. 2006; Philippe et al. 2009; Philippe et al. 2011). Phylogenetic analyses of genomic data typically reveal high levels of topological incongruence among gene trees. Incongruence can be real, where different genes have different genealogies due to the complex dynamics of the evolutionary processes (e.g., paralogy, horizontal transfer, or incomplete lineage sorting; Maddison 1997; Edwards et al. 2007; Rosenberg and Tao 2008; Degnan and Rosenberg 2009; Edwards 2009; Kubatko et al. 2009; Liu et al. 2009; Huang et al. 2010), or it can be the consequence of methodological artifacts. Methodological artifacts include spurious results due to sampling error (low signal-to-noise ratio), or a failure of evolutionary models to accommodate convoluted patterns in the data (e.g., López et al. 2002; Pagel and Meade 2004; Collins et al. 2005; Lockhart et al. 2006; Phillips et al. 2006; Rasmussen and Kellis 2007; Galtier and Daubin 2008; Sheffield et al. 2009; Roure and Philippe 2011; Betancur-R. et al. 2013; Romiguier et al. 2013). The ultimate effect of gene-tree discordances on inferences of organismal phylogeny will thus depend on character sampling, biological processes, and the methods and models used to attempt to correct for errors caused by one or more of these factors. In a recent article, Salichos and Rokas (2013; S&R) examine several elements that account for methodological artifacts in phylogenomics. They assert that selection of genes with “strong genealogical signals” (as measured by bootstrap support of their inferred trees) can lead to increased phylogenetic accuracy. They also imply that this increase in accuracy can be obtained without explicitly accounting for other factors, particularly coalescent variance. In that study, S&R compare 23 yeast genomes and observed that genealogies obtained for 1070 orthologous genes are all different from each other and also differ from the topology obtained either by concatenating all genes or by an extended consensus phylogeny of all gene trees. They develop a new measure of incongruence, which they term “internode certainty”, to gauge the level of conflict inherent in the data supporting specific internodes of the phylogeny. Based on this measure, S&R conclude that “the use of slowly evolving genes and conserved sites increased incongruence across many internodes of the yeast phylogeny” (S&R, p. 328). This result is shown in their Figure 2, where slowly evolving genes yield lower internode certainty scores than do fast evolving genes. They also find that “high bootstrap support is a good indicator of a gene's phylogenetic usefulness” and that “concatenating genes with high bootstrap support reduces incongruence and improves resolution” (S&R, p. 329). Based on these results they propose the adoption of three “fundamental revisions to current practices” one of which is to use conditional combination approaches selecting genes with strong genealogical signal based on bootstrap resampling scores. The notion that slowly evolving genes are a bad choice to resolve basal nodes at deep phylogenetic levels is contrary to widespread practice in recent studies (e.g., Li et al. 2007; Jian et al. 2008; Li et al. 2008; Regier et al. 2008; Zhang et al. 2012; Lang et al. 2013). We challenge S&R's interpretations and recommendations with new analyses of their yeast data. We demonstrate that the high phylogenetic incongruence among conserved genes observed by S&R is an artifact due to sampling error. We also note that S&R's recommendation to focus on genes with strong signal make the task of choosing genes upon which to base phylogenetic inference difficult, as most measures of signal can only be determined after data have been collected and complete genome sequences do not yet exist for most taxa under study. Therefore, S&R's recommendation has little practical value for defining sets of target loci for increasingly popular phylogenomic approaches based on target enrichment (Faircloth et al. 2012; Lemmon et al. 2012; Li et al. 2013). S&R use gene support frequency and their new incongruence index to measure tree certainty and to assess the relative phylogenetic performance of different data partitions in the yeast genome. They identify 131 genes with high tree certainty scores and bootstrap support scores above 80% (reported in S&R Fig. 2) and report that the lowest scores are obtained for a set of 100 slowly evolving genes. They explore two alternative species-tree topologies derived from the data sets of either all 1070 concatenated genes (S&R's Fig. 1a; hereafter concatenation tree) or the 131 genes with strongest genealogical signal (S&R's Supplementary Fig. S11c,f; hereafter strong-signal tree). While the two topologies are highly similar, they differ in the relative placement of three lineages (Candida glabrata, Saccharomyces castellii, and the Saccharomyces sensu stricto group). It is apparent that S&R favor the strong-signal topology on the basis of shared chromosomal rearrangements previously reported (Scannell et al. 2006), which in turn seems to guide their preference for the strong-signal genes. Several relevant issues not discussed by S&R, however, could have affected their results of marker performance (S&R's Fig. 2) and may prompt alternative interpretations. Relationship between tree incongruence and sequence length based on the 1070 yeast gene alignments. Incongruence is measured as the topological distance (PH85) from each gene-tree topology with respect to (a) the concatenation or (b) the “strong signal” species-tree topologies of S&R (see S&R's Fig. 1 and Supplementary Fig. S11c,f, respectively). We first assessed the relationship between gene sequence length and gene-tree incongruence, measured as the topological distance between the 1070 gene trees and each of the two alternative species trees (concatenation and strong signal). All other factors being equal, phylogenetic analyses of short sequences are more prone to sampling error than are long sequences, simply because they contain less phylogenetic information (Rasmussen and Kellis 2007). Consistent with this prediction, we found a significant negative correlation between gene sequence length and gene-tree incongruence (r2 = 0.19 and 0.17 for concatenation and strong-signal trees, respectively; P < 0.0001; Fig. 1). We then selected two groups of gene alignments for further analysis, one that includes the 131 genes with strongest genealogical signal reported by S&R (hereafter 131S) and a second group of the 260 most conserved genes (260C), based on pairwise sequence identities (see Methods section). In contrast to S&R's analysis, for which they used the 100 most conserved genes, we based our comparisons on 260 conserved genes to increase sample size so that the concatenated alignments of both 260C and 131S sets are similar in length (ca. 127,400 sites each). This condition is important for the binning experiment described below. These two sets of genes have significantly different pairwise sequence identities (U-test P < 0.00001; Supplementary Fig. S1b, ranging from 64.5 to 95.4% for the 260C set (mean = 73.5%) and from 26.8% to 85.0% for the 131S set (mean = 49.5%). Consistent with the results obtained by S&R, gene-tree incongruence is significantly higher (18.3 vs. 11.6) in the 260C set than in the 131S set (Fig. 2a,b; U-test P < 0.00001). Likewise, genealogical signal is weaker in the 260C set, with mean bootstrap support of 62.7% for 260C versus 83.1% for 131S (Fig. 2c; U-test P < 0.00001). It is important to note, however, that genes in the 260C set not only contain fewer variable sites per unit length relative to the 131S set (62.8% vs. 83.7%, respectively; Supplementary Fig. S1c; U-test P < 0.00001), but their sequence length also is on average two times shorter (260C mean = 490.0 sites vs. 131S mean = 971.8 sites; Supplementary Fig. S1a; U-test P < 0.00001). Higher incongruence among the most conserved gene trees obtained here and by S&R could be an artifact resulting from analysing shorter fragments with less information content (Fig. 1), rather than due to conflicting signals. Assessment of incongruence (a, b) and genealogical signal (c) of yeast trees estimated for the 260 most conserved (260C) genes and the 131 genes with the strongest genealogical signal (131S), under multiple binning strategies. Incongruence is measured as the topological distance (PH85) from each gene- or supergene-tree topology with respect to (a) the concatenation or (b) the strong-signal species-tree topologies of S&R (see S&R's Fig. 1 and Supplementary Fig. S11c, respectively). Incongruence among trees from binned data decreases and genealogical signal increases as a function of the length in binned sequences. Higher incongruence for 260C over 131S initially obtained with individual gene alignments disappears for bin sets and is significantly reversed for the 12.5K bins in left panel. Weaker signal for 260C over 131S initially obtained with individual gene alignments disappears for the 7.5–12.5K bin sets. Values and thick lines inside/outside boxes indicate data means and medians, respectively. *U-test P < 0.05; ***U-test P < 0.0001. To compensate for differences in sequence length and sampling error among the 131S and 260C gene sets, we explored the effect of binning genes into “supergenes” of various fixed lengths. Data binning is a statistical technique that can account for sampling error by increasing signal-to-noise ratio, and it has been applied in phylogenetics recently (Bayzid and Warnow 2013). We concatenated genes to generate bins of 2500 (2.5K), 5000 (5K), 7500 (7.5K), 10,000 (10K), or 12,500 (12.5K) sites. Analyses of binned data show that the relative incongruence between 260C and 131S initially obtained with gene trees is no longer observed for bins (Fig. 2a,b). Likewise, relative differences in bootstrap support disappear for the 7.5–12.5K bins (Fig. 2c). Remarkably, the 12.5K bins in the 260C set show a slightly better performance than in the 131S set in terms of both incongruence (3.2–5.0 vs. 5.6–6.4, respectively; Fig. 2a,b) and bootstrap support (93.1% vs. 92.2%, respectively; Fig. 2c), although the differences are only significant in one comparison (Fig. 2a; P = 0.016) presumably due to small sample size (N = 10). Taken together, these results strongly suggest that higher levels of incongruence and weaker signal reported for the conserved genes are the result of sampling error. Because of their intrinsic properties, conserved genes are individually unlikely to yield well-resolved and highly supported trees, but collectively can be informative because they are expected to meet the assumptions of evolutionary models more closely than the complex dynamics of highly variable genes (which are otherwise prone to systematic error and model misspecification issues). Previous studies offer empirical evidence for this prediction, showing that slowly evolving genes, despite having fewer informative sites, provide higher resolution and support for deep-level radiations than rapidly evolving genes (e.g., Jian et al. 2008; Regier et al. 2008; Zhang et al. 2012; Lang et al. 2013). Thus, we question the notion that inferring ancient divergences requires selecting individual genes with strong signals, as recommended by S&R. Our re-analysis of S&R's data set reinforces previous studies that indicate that conserved genes are indeed good phylogenomic markers. Identification of conserved genes can be done effectively before assembling phylogenomic data sets, without the need of a posteriori selection based on inferred trees. For instance, recently developed approaches that target regions of the genome that are both highly conserved and shared across genome-enabled taxa, can be used to identify potential markers for phylogenetic inference (Li et al. 2007; Faircloth et al. 2012; Lemmon et al. 2012; Li et al. 2013). Reduced-representation libraries can be obtained via target-capture protocols (Gnirke et al. 2009) and subsequently sequenced using massively parallel sequencing platforms. Recommendations to focus on selective data partitions assessed by phylogenetic analysis may become relevant in the future, however, once complete genomes are available for all species of interest. We have examined some of the factors that account for methodological artifacts in a phylogenomic data set comprising orthologous genes. It would be naive to assume that selection of genes on the basis of sampling error or evolutionary rate alone guarantees phylogenetic accuracy, as both systematic biases and biological processes acting alone or in concert generate discordance. While the relative impact of coalescent variance and systematic bias is rarely examined in empirical or simulated data sets (but see Huang et al. 2010; Camargo et al. 2012; Betancur-R. et al. 2013), coalescent-based species-tree methods are increasingly gaining attention (Maddison 1997; Edwards et al. 2007; Rosenberg and Tao 2008; Degnan and Rosenberg 2009; Edwards 2009; Kubatko et al. 2009; Liu et al. 2009; Huang et al. 2010). An unsettling corollary of the binning experiment described above, however, is that gene-tree discordance due to sampling error is so extensive in this data set (Fig. 2a,b) that the efficacy of approaches that account for coalescent variance without successfully decoupling this effect or other systematic biases inherent to the deep-scale nature of the yeast data set are unclear (see also discussion in S&R). We concur with S&R's recommendation that best phylogenomic practices demand careful exploration of internal conflict in the data sets, but measures of incongruence among partitions proposed in their article (e.g., internode certainty) should be modified to consider sampling error to avoid misleading results. The most conserved gene set (260C) includes genes from the pool of 1070 gene orthologs with the highest percentage of pairwise sequence identity, as estimated with Geneious Pro v6 (Biomatters Ltd.). Following S&R, the 131S set includes genes whose trees have mean support values > 80%. Initially, 17 genes were common to both gene sets but they were finally excluded from the conserved set to avoid overlap (i.e., the 260C set includes 260 of the 277 most conserved genes). Supergenes for the 260C and 131S sets were assembled using the “naive binning” approach (Bayzid and Warnow 2013), by randomly concatenating all genes for each set and exporting alignments of fixed lengths using PAUP* v4.0b10 (Swofford 2002). The concatenated alignments of the 260C and 131S sets are both of similar length (ca. 127,400 sites each). Trees from bin sets were inferred in RAxML v7.3 (Stamatakis 2006) using the PROTGAMMAIWAGF model with 30 independent replicate runs and 300 bootstrap replicates. Incongruence is measured as the topological distance between the resulting tree topologies under each treatment and two alternative species-tree topologies (S&R's Fig. 1a and Supplementary Fig. S11c,f). Topological distances were calculated using the R package Ape (Paradis et al. 2004; R Development Core Team 2011). We used the PH85 distance, which is twice the number of internal branches defining different bipartitions of the tips (Penny and Hendy 1985) and is equal to 2x Robinson–Foulds, the most commonly used metrics of tree distance. Genealogical signal for each tree is measured as mean bootstrap support across nodes as estimated using the Ape package and custom R code. Statistical significance among treatments was assessed with the non-parametric Mann–Whitney U-test in R. All supplementary files listed here, including figures, tables, and aligned data sets are available from the Dryad data repository at http://datadryad.org, doi:10.5061/dryad.157d7. Data files and/or other supplementary information related to this article have been deposited at Dryad under doi:10.5061/dryad.157d7. This research was partially funded by the GWU Selective Excellence in Diversity of Life program to R.B.R. and by the National Science Foundation Grants DEB-1004765 (to G.O.) and DEB 1132229 (to G.J.P.N.). Protein alignments and individual gene trees were kindly provided by A. Rokas. We thank A. Pyron, for insightful comments on earlier versions of the manuscript. Associate Editor M. Holder, Editor F. Anderson, B. Faircloth, and an anonymous reviewer, provided thoughtful suggestions to improve the article. R.B.R. conceived the study and analysed the data; R.B.R., G.J.P.N., and G.O. discussed the results and wrote the article.