New Methods to Calculate Concordance Factors for Phylogenomic Datasets

1Measures of branch support such as the bootstrap and Bayesian posterior probability are important for making robust inferences from phylogenetic trees. However, although these measures provide useful information about the statistical support for a given branch, neither captures the topological variation present in the underlying data.

2A complementary approach involves calculating the fraction of loci consistent with a particular branch, thus capturing underlying agreement and disagreement in the data. Various approaches have been suggested, including the gene concordance factor (gCF)—sometimes also referred to as the “gene support frequency”—and internode certainty, with the former being the most widely used. The gCF describes for each branch in a reference tree the proportion of inferred single-locus trees that contain that branch. Although intuitive, the gCF suffers from three limitations: first, there are few software implementations; second, although issues with incomplete taxon sampling of gene trees have been addressed by some authors, these fixes are not available in implementations that allow for the calculation of gCF values; and third, low gCF values are hard to interpret because they may result from strongly supported discordance among individual gene trees or weak phylogenetic signal in individual loci.

3Here, we resolve these issues by implementing the calculation of gCF while accounting for unequal taxon sampling in the popular IQ-TREE package, and by introducing the site concordance factor (sCF), a measure that estimates concordance at the level of individual sites. We argue and show by detailed analyses of ten empirical datasets that concordance factors complement commonly used metrics like the bootstrap by providing additional information and insights about topological variation.

4For a given unrooted bifurcating reference tree, (e.g., an estimate of the species tree), and a set of unrooted bifurcating input trees, S = (e.g., individual gene trees), we can calculate the gCF for every internal branch in. Each gene tree,, must contain a subset of the taxa in but need not include all of the same taxa.

5Each internal branch of is associated with four clades of representing four taxon subsets,, such that the bipartition corresponds to the branch (fig. 1). The reference tree contains branch by definition, but branch may or may not be present in a gene tree of the same taxa (or a subset of them). Let be the set of taxa in that are also present in. We say is decisive for if are nonempty (we use “decisive” in the sense of, but note for clarity that use this term in a different way): that is, a gene tree is decisive as long as it contains at least one of the taxa in each of the taxon subsets. A decisive gene tree can potentially contain the branch. We say is concordant with if is decisive for and the bipartition corresponds to a branch in. In other words, a gene tree is concordant with the reference tree if it could have contained branch (i.e., it is decisive) and it does contain branch. The gCF for branch is now defined as

6In other words, the gCF is the proportion of input trees decisive for that are concordant with; that is, the proportion of all input trees that could have contained that do contain. We note that gCF is also used to calculate internode certainty and related measures.

7To further understand concordant and discordant gene trees, we categorize the discordant input trees (i.e., those that do not contain branch x) into three groups and calculate three discordance factors: gDF1, gDF2, and gDFP. The first two groups are obtained by applying a nearest neighbor interchange around branch to result in two alternative topologies (fig. 1) with two branches not appearing in: and. Accordingly, we then calculate discordance factors as and. Some input trees may not contain any of the branches x, y, or z. This will occur for input trees in which one or more of clades A, B, C, or D are not monophyletic. We define the proportion of decisive input trees that fall into this category as, where the “P” stands for “paraphyly.” The four proportions,,, and will sum to 1, as every decisive input tree must be included in one of the four categories.

8To calculate the sCF for branch, we randomly sample quartets of taxa, where are in, respectively. For each quartet, we examine the subalignment of taxa. For every site in this alignment, we call decisive for if the characters are all present and is parsimony informative when restricted to this quartet of taxa. Decisive sites can be concordant or discordant with. We say that site is concordant with if (i.e., supports the bipartition ). The concordance factor for is defined as

9Since many such quartets exist around branch (when sampling individual tips from within A, B, C, and D), we define sCF as the mean over m random quartets:

10Thus, the sCF is the average proportion of sites decisive for that are concordant with. In effect, the sCF is a measure of concordance for sites that is directly comparable to the measure of concordance for single-locus trees provided by the gCF. Unlike related quartet approaches that are based on maximum likelihood (ML) tree inference and that calculate discordance at the level of the whole alignment, the sCF uses parsimony criteria to calculate discordance at the level of individual sites. We note that the sCF is closely related to the values derived from spectral analysis. Spectral analysis is a tree-independent method of investigating phylogenetic signal. In the simple case of binary data, such that each site corresponds to a bipartition of taxa, the spectral support for a single branch counts the number of sites that correspond to the bipartition defined by that branch. In the case of binary data, this is the same as the value that we define above. The sCF in this case differs from the spectral support insofar as it is calculated by averaging over a large set of values calculated by repeatedly subsampling quartets of taxa from an alignment.

11Similarly to gene discordance factors, we also define the site discordance factors and. There is no sDFP category, because unlike the gDF values, every decisive site must contribute to one of the three proportions,, or; any quartet of taxa resolves into exactly three tree topologies shown in figure 1. In other words, the sum of,, and will always be 1.

12The above definitions of gCF and sCF apply only to unrooted trees. However, in the case that users have a rooted reference tree as well as rooted gene trees, we can extend the calculation of the gCF to allow us to calculate different gCF values on either side of the root. To do this we first add a virtual root node into the rooted reference tree, resulting in an unrooted tree. Similarly, we convert the rooted gene trees into unrooted trees by adding the same virtual root. This allows us to then follow the same procedure as above for calculating gCF values on unrooted trees, with the sole difference that the output is a rooted instead of an unrooted reference tree. In this case, it is possible to calculate gCF values on a rooted reference tree where the gCF values on either side of the root may differ. In all other cases (i.e., for sCF values, and if either the reference tree or the gene trees are unrooted), it is not possible to calculate different concordance factors on either side of the root.

13We provide two new options, –gcf and –scf, in IQ-TREE version 2 to compute gCF and sCF, respectively. A tutorial for how to use these options is provided at http://www.iqtree.org/doc/Concordance-Factor (last accessed May 13, 2020). Both options can be combined in a single run, which will calculate the gCF and sCF for every branch in the input reference tree. Although sCFs can be calculated on any alignment, the calculation of gCFs requires individual gene trees. We have therefore implemented a convenient option, -S, to specify a partition file or a directory of single-locus alignments in which IQ-TREE will infer separate trees for each partition or alignment.

14IQ-TREE provides a suite of output files to assist users in understanding and investigating gCF and sCF values. It provides tree files that can be viewed in most tree viewers, which contain information on both the proportional data (gCF, gDF1, gDF2, and gDFP; sCF, sDF1, and sDF2) and their corresponding absolute count data (gCF_N, gDF1_N, gDF2_N, and gDFP_N; sCF_N, sDF1_N, and sDF2_N), as well as the number of decisive genes and sites for each branch (gN and sN, respectively). It also provides a tree file that combines information on the gCF, sCF, and any bootstrap values that have been calculated (e.g., fig. 2). In addition to this, IQ-TREE provides a “.cf.stat” file that contains all 16 concordance and discordance values listed above for every branch in the reference tree in a machine-readable tabulated format, and through the –cf-verbose option, it provides tabulated files that detail for every branch in the reference tree whether each gene tree was concordant with that branch (the “.cf.stat_tree” file produced from a –gcf analysis), and the average number of sites in each locus that was concordant with that branch (the “.cf.stat_loci” file produced from a –scf analysis). Together these output files provide both a convenient overview of the data, and the opportunity to understand in much more detail the extent to which each locus is concordant or discordant with the reference tree.

15To demonstrate the use of gCF and sCF values, we first analyzed a dataset containing 3,220 ultra-conserved elements (UCEs) from lizards, totaling 1,301,107 bp for 43 species. To do this, we estimated a concatenated ML tree, 3,220 UCE trees, and the gCF, sCF, and bootstrap values with IQ-TREE (fig. 2).

16To investigate the relationships between dataset size, concordance factors, and bootstrap values, we also analyzed a collection of nine additional empirical phylogenomic datasets that represent a range of clades and data types (table 1). For each dataset, we performed 20 analyses across a 20-fold range of dataset sizes: the first analysis included ten randomly selected loci from the complete dataset, and each subsequent analysis added ten more randomly selected loci up to a maximum of 200 loci. For each such analysis, we calculated the gCF, sCF, and bootstrap values for every branch of the concatenated ML tree estimated from the 200-locus dataset, with the same approach as above. For the lizard dataset, we calculated both standard bootstrap (StdBoot) values and ultrafast bootstrap (UFBoot) values. As expected, we observed that StdBoot and UFBoot values were very similar (fig. 3B and supplementary fig. 8D, Supplementary Material online); therefore, because of the high computational cost of calculating StdBoot values, we calculated only UFBoot values for the remaining nine datasets. Since the results of all ten analyses are very similar, we focus here on the lizard dataset and present the results for the other nine datasets in the supplementary figures 1–10, Supplementary Material online.

17Reproducible analyses are provided in the Supplementary Material online available at https://doi.org/10.5281/zenodo.1949287.

18Our analyses highlight a number of important differences between concordance factors and bootstrap values. First, concordance factors of any type tend to be much lower than bootstrap values or other measures of statistical support for branches in the species tree (e.g., fig. 3 and supplementary figs. 1–10, Supplementary Material online). This is because bootstrap values measure the sampling variance in support of a focal branch, whereas gCF and sCF values measure the underlying variance in support of that branch at the gene- and site-level, respectively. In other words, resampled datasets may always return the same tree (i.e., 100% bootstrap support), even though incomplete lineage sorting or other processes that lead to genealogical discordance are at work (e.g., gCF and sCF values ≪100%). Of particular note is that a very high bootstrap value (e.g., of 100%) does not predict or require a similarly high concordance factor. For example, although all but one bootstrap value on the tree in figure 1 is 100%, the smallest gCF and sCF values are 4.5% and 33.2%, respectively. This pattern is repeated across all 200 analyses we performed on all ten empirical datasets (fig. 3 and supplementary figs. 1–10, Supplementary Material online).

19Because they measure different things, concordance factors and bootstrap values are affected very differently by the addition of loci to a dataset (fig. 3 and supplementary figs. 1–10). In all of the ten empirical datasets that we analyzed, adding loci to the dataset tended to increase the bootstrap values of individual branches (e.g., fig. 3A, see also supplementary figures, Supplementary Material online) and concomitantly the proportion of branches with a bootstrap value of 100% (fig. 3B, see also Supplementary Material online). This is expected: the bootstrap is effectively measuring the standard error of the mean in a dataset (where the mean represents the tree inferred from the full dataset), and this standard error will go down with more samples. The same is not true of concordance factors, which display some estimation error when datasets are small (e.g., 30 or fewer loci in fig. 3), but subsequently remain almost completely insensitive to the addition of loci to the dataset (fig. 3 and supplementary figs. 1–10, Supplementary Material online). This is also expected: any measure of the underlying variance (or standard deviation) in a distribution will have some estimation error when the sample size is small but should not change monotonically as the sample size increases. These differences highlight the complementary information that can be gained from calculating both bootstrap values and concordance factors for individual branches. One measure is not better or worse than the other, rather, concordance factors provide useful information that bootstraps do not, and vice versa. Note that bootstrap values or posterior probabilities calculated by resampling gene trees (e.g., ) have the same behavior as the site-wise bootstrap carried out here and are also not equivalent to concordance factors.

20In principle, both gCF and sCF values can range from 0% (no genes/sites are concordant with the focal branch) to 100% (all genes/sites are concordant with the focal branch). In practice however, as exemplified in figure 3 and supplementary figures 1–10, Supplementary Material online, empirical gCF values tend to range from 0% to 100%, whereas empirical sCF values are rarely lower than 33% (represented by the dashed line in the sCF panel of fig. 3A). This is due to an important underlying difference in the way that the two values are calculated. The sCF is calculated from quartets, so a single site can only support one of three topologies (fig. 1). Because of this, if there is no consistent information in an alignment (e.g., if a long alignment were generated at random) we expect a roughly equal proportion of sites supporting each of the three trees, leading to an sCF value of ∼33% (for the same reason, sCF values for very long branches will approach 33% due to saturation). The same is not true for gCF values, because a gene tree can support not only the three possible relationships shown in figure 1 but also any other relationship in which one or more of clades A, B, C, or D are not monophyletic. The higher the number of gene trees in this latter group, the closer the gCF value will be to 0%. Because of this, we should expect gCF values to be particularly low when gene trees are estimated from alignments with limited information or where branch x is extremely short; in such cases either technical or biological processes may increase the proportion of gene trees that fail to recover the monophyly of clades A, B, C, or D found in the reference tree. Missing data may also impact gCF and sCF values, however, the relationship here is less clear. The requirement that a gene tree or a site is decisive (i.e., could in principle contain branch x, see above) should limit the impact of missing data on gCF and sCF estimates. Nevertheless, since phylogenetic estimates are known to worsen as the proportion of missing data increases (e.g.,; ), it is plausible that concordance factors may systematically decrease as the proportion of missing data increases.

21Finally, cases where the sCF value is lower than 33% may be of particular interest. These cases are, by definition, those in which maximum parsimony (MP) would favor a different resolution of a split found in the reference tree. If the reference tree was calculated from any method other than MP, there are at least two explanations for an sCF value lower than 33%. First, the branch of interest may be in an area of parameter space in which high levels of incomplete lineage sorting are known to mislead concatenated ML analyses but not MP analyses. Misleading reference trees may therefore be produced by either concatenated ML of the entire dataset or gene tree methods that use shorter sets of concatenated loci as their input (because most protein-coding genes are themselves made up of multiple topologies). Second, and more generally, there are multiple reasons why likelihood, Bayesian, or gene tree methods for producing a species tree will differ from MP resolutions. For instance, the branch of interest may be unduly affected by a small number of highly influential sites in a concatenated ML analysis. In this case, the influential sites can have an outsized influence on the ML resolution of a split because they have extreme differences in likelihood between different resolutions of that split. Because MP does not account for likelihood differences—it instead weights all sites equally—MP analyses remain unaffected by such outliers. Thus, cases in which the sCF is much lower than 33% may merit further investigation.

22We hope that the user-friendly implementation of gCF and sCF in IQ-TREE will assist researchers in gaining additional insights into their phylogenetic reconstructions. In particular, we encourage phylogeneticists to calculate both bootstrap values and concordance factors for the branches on their trees, as the two measures provide complementary information that may help to improve the accuracy of our interpretations of phylogenetic reconstructions. Indeed, the use of concordance factors may help to alleviate the commonly cited problem in phylogenomics that bootstrap values provide relatively little information when they are all 100%.

23Data are available at https://doi.org/10.5281/zenodo.1949287.