Improved robustness to gene tree incompleteness, estimation errors, and systematic homology errors with weighted TREE-QMC

Summary methods are widely used to reconstruct species trees from gene trees while1 accounting for incomplete lineage sorting; however, it is increasingly recognized that their2 accuracy can be negatively impacted by incomplete and/or error-ridden gene trees. To3 address the latter, Zhang and Mirarab (2022) leverage gene tree branch lengths and4 support values to weight quartets within the popular summary method ASTRAL.5 Although these quartet weighting schemes improved the robustness of ASTRAL to gene6 tree estimation error, implementing the weighting schemes presented computational7 challenges, resulting in the authors abandoning ASTRAL’s original search algorithm (i.e.,8 computing an exact solution within a constrained search space) in favor of search heuristics9 (i.e., hill climbing with nearest neighbor interchange moves from a starting tree constructed10 via randomized taxon addition). Here, we show that these quartet weighting schemes can11 be leveraged within the Quartet Max Cut framework of Snir and Rao (2010), with only a12 small increase in time complexity compared to the unweighted algorithm, which behaves13 more like a constant factor in our simulation study. Moreover, our new algorithm,14 implemented within the TREE-QMC software, was highly competitive with weighted15 ASTRAL, even outperforming it in terms of species tree accuracy on some challenging16 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 2 HAN AND MOLLOY model conditions, such as large numbers of taxa. In comparing unweighted and weighted17 summary methods on two avian data sets, we found that weighting quartets by gene tree18 branch lengths improves their robustness to systematic homology errors and is as effective19 as removing the impacted taxa from individual gene trees or removing the impacted gene20 trees entirely. Lastly, our study revealed that TREE-QMC is highly robust to high rates of21 missing data and is promising as a supertree method. TREE-QMC is written in C++ and22 is publicly available on Github: https://github.com/molloy-lab/TREE-QMC23 Key words: species trees; coalescence; incomplete lineage sorting; summary methods;24 quartets; gene tree error; missing data; partitioned coalescence support25 26 Over the last decade summary methods have been widely adopted in species tree27 estimation pipelines, starting from the identification of orthologus genomic regions from28 target capture sequencing data or assembled genomes (Steenwyk et al., 2023, 2024). Next,29 an evolutionary history for each region (referred to as a gene tree) is reconstructed30 typically in two steps: first a multiple sequence alignment (MSA) is built from the31 unaligned sequences and second a phylogeny is estimated from the MSA via maximum32 likelihood under standard models of molecular sequence evolution (e.g., the GTR model;33 Tavar´ e, 1986). Finally, the estimated gene trees are given as input to summary methods,34 many of which are designed to reconstruct the species tree while appropriately accounting35 for gene tree discordance due to incomplete lineage sorting (ILS), unlike concatenation36 (Roch and Steel, 2015). Examples of such summary methods include MP-EST (Liu et al.,37 2010) and STELAR (Islam et al., 2020), which are based on triplets implied by rooted38 gene trees, NJst (Liu and Yu, 2011) and ASTRID (Vachaspati and Warnow, 2015), which39 are based on pairwise internode distances implied by unrooted gene trees, and ASTRAL40 (Mirarab et al., 2014) and TREE-QMC (Han and Molloy, 2023), which are based on41 quartets implied by unrooted gene trees.42 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 3 Despite the wide spread use of summary methods, it is increasingly recognized that43 their accuracy depends on gene tree quality. In simulations, summary methods typically44 decrease in accuracy as gene tree estimation error (GTEE) or incompleteness increases (Xi45 et al., 2015; Mirarab and Warnow, 2015; Molloy and Warnow, 2018; Nute et al., 2018).46 Moreover, a growing number of systematic studies show that homology errors can impact47 species tree estimates produced by different summary methods (Simmons et al., 2022;48 Springer and Gatesy, 2024).49 The dominant approach for combating poor quality gene trees is gene tree filtering,50 in which entire gene trees are removed from the data set typically based on some threshold51 for the proportion of missing taxa (Hosner et al., 2016; Cunha et al., 2021) and/or proxies52 for GTEE (Simmons et al., 2016). However, simulation studies have suggested that53 filtering can reduce the accuracy of summary methods and that only gene trees with very54 high GTEE should be removed entirely (Molloy and Warnow, 2018). This finding55 motivates the use of less aggressive approaches, for example contracting very low support56 branches in gene trees (Zhang et al., 2018; Simmons and Gatesy, 2021) and/or removing57 specific taxa from individual gene trees. Taxa are typically removed based on the fraction58 of gapped sites in the gene MSA (referred to as fragmentary missing data by Sayyari et al.,59 2017 and the One Thousand Plant Transcriptomes Initiative, 2019) or based on (long)60 outlier branch lengths in the gene tree (Mai and Mirarab, 2018; One Thousand Plant61 Transcriptomes Initiative, 2019). Critically, these approaches require substantial effort and62 decision-making on the part of researchers, who must choose how to quantify error as well63 as select thresholds for filtering gene trees, contracting branches, and/or removing64 individual taxa. Often a range of parameter settings are explored, leading to many65 different estimates of the species tree (e.g., Hosner et al., 2016).66 To address GTEE, Zhang and Mirarab (2022) proposed to leverage gene tree67 branch lengths and support values within the popular summary method ASTRAL.68 ASTRAL seeks a species tree that maximizes the number of quartets (unrooted four-leaf69 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 4 HAN AND MOLLOY trees) implied by the input gene trees. This optimization problem is NP-hard (Lafond and70 Scornavacca, 2019), so ASTRAL solves the problem within a constrained search space71 built from the input gene trees (Mirarab et al., 2014; Mirarab and Warnow, 2015; Zhang72 et al., 2018). The goal of weighted ASTRAL is to down-weight quartets with low support on73 the internal branch (an indicator that the quartet cannot be confidently resolved due to74 insufficient phylogenetic signal) and/or long lengths on the terminal branches (an indicator75 that the quartet may be impacted by long branch attraction). However, these quartet76 weighting schemes increase the complexity of ASTRAL’s original search algorithm. To77 achieve greater computational efficiency, weighted ASTRAL employs heuristic search,78 specifically hill climbing with nearest neighbor interchange moves from a starting tree79 constructed via randomized taxon addition. This new heuristic, often referred to as80 ASTER, has been found to be more robust to incomplete gene trees than the original81 ASTRAL algorithm (Zhang and Mirarab, 2022; Morel et al., 2023).82 Following the success of weighted ASTRAL/ASTER, Liu and Warnow (2023)83 incorporated weighting schemes based on gene tree branch lengths or support values but84 not into the distance method ASTRID. Overall, leveraging gene tree branch lengths and85 support values is a promising path forward but weighted summary methods are still in86 their infancy. Comparatively little is known about how they perform on real data, and87 there are significant challenges to developing efficient algorithms.88 In this paper, we introduce a new weighted summary method, called weighted89 TREE-QMC, which integrates the quartet weighting schemes of Zhang and Mirarab (2022)90 into the Quartet Max Cut framework (QMC) of Snir and Rao (2010) while also efficiently91 addressing the issue of quartet weight normalization for incomplete gene trees. A previous92 challenge in using QMC as a summary method was that it required researchers to93 explicitly extract quartets from gene trees prior to species tree estimation; we recently94 addressed this issue, introducing TREE-QMC Han and Molloy (2023). Weighting quartets95 based on gene tree branch lengths and support values presents computational challenges,96 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 5 and perhaps surprisingly, we achieve an algorithm with only a small increase in time97 complexity and no increase in storage complexity compared to the unweighted98 TREE-QMC method. Moreover, our simulation study showed that weighted TREE-QMC99 was fast and highly competitive with weighted ASTRAL/ASTER in terms of species tree100 accuracy, even outperforming weighted ASTRAL/ASTER on some challenging conditions,101 such as large numbers of taxa.102 We then evaluated weighted summary methods on three recently published data103 sets: one for plants with high rates of missing taxa (Morel et al., 2023) and two for birds,104 Palaeognathae (Cloutier et al., 2019) and Neognathae (Wu et al., 2024a), both of which105 have been found to contain homology errors (Simmons et al., 2022; Springer and Gatesy,106 2024). The homology errors for Palaeognathae were biased towards two species (Chicken107 and White-Throated Tinamou), which clustered together off a long branch in the 105108 affected ultraconserved element (UCE) gene trees (Simmons et al., 2022). We found that109 TREE-QMC and ASTRAL/ASETER with hybrid or length weighted recovered the same110 species tree topology as the unweighted methods after removing Chicken and111 White-Throated Tinamou from the 105 affected UCE trees or after filtering the 105112 affected UCE trees entirely. Otherwise, an alternative resolution of the focal branch was113 favored by the unweighted summary methods. To evaluate whether biased homology errors114 had a systematic impact on the reconstructed species tree, we implemented Partitioned115 Coalsecence Support (PCS) (Gatesy et al., 2017, 2019) within TREE-QMC and applied it116 to the focal branch in Palaeognathae. The PCS analysis showed that the 105 UCE trees117 with homology errors consistently supported an alternative resolution of the focal branch;118 however, the length and hybrid weighting schemes mitigated this issue. Quartet weighting119 schemes did not have a big impact on TREE-QMC for the other avian data set, where120 homology errors were not clearly biased towards particular species. For this data set and121 the plant data set, TREE-QMC produced trees that were closer to the concatenation tree122 and reference taxonomies than those produced by ASTRAL/ASTER, which produced in a123 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 6 HAN AND MOLLOY highly imbalanced (i.e., caterpillar-like) branching pattern not seen in the other estimated124 species trees. Overall, our results suggest that TREE-QMC is efficient, easy-to-use, and125 improves robustness to gene tree incompleteness, estimation error, and homology errors.126

I. Weighted TREE-QMC127 Overview of unweighted TREE-QMC128 To introduce weighted TREE-QMC, we briefly review the original (unweighted)129

(Han and Molloy, 2023), which leverages the Quartet Max Cut (QMC)130 framework (Snir and Rao, 2010, 2012; Avni et al., 2014). The idea behind QMC is to131 reconstruct the species tree in a divide-and-conquer fashion from a set of input quartets.132 At each step of the divide phase, a graph is built with vertices representing taxa. The edges133 between vertices (taxa) are weighted by the relevant input quartets. A cut of the quartet134 graph partitions the taxa into two disjoint subsets, corresponding to a bipartition135 (associated with an internal branch) in the unrooted species tree. The algorithm continues136 by recursion on subproblems defined by the taxa on each side of the bipartition.137 Termination occurs when there are three or fewer taxa in the subproblem as the solution138 (i.e., a phylogenetic tree on the taxon set) is trivial (Fig. 6a in Appendix A). Importantly,139 artificial taxa are added to subproblems at each step in the divide phase to represent the140 taxa on the other side of the bipartition (Fig. 6a in Appendix A). During the conquer141 phase, subproblem trees are connected at artificial taxa, ultimately producing an unrooted,142 binary species tree on the complete set of taxa (Fig. S1 in Han and Molloy, 2023).143 The original (unweighted) TREE-QMC method improves upon QMC in two ways.144 First, QMC operates on quartets, which must be extracted from gene trees prior to species145 tree estimation. Given k gene trees on n taxa, the preprocessing phase alone has time146 complexity Ω(n4k) and storage complexity Ω( n4) (note that we do not include the storage147 required for input gene trees). In contrast, TREE-QMC builds the quartet graph directly148 from the input gene trees, having time complexity O(n3k), with some assumptions on149 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 7 subproblem decomposition, and storage complexity O(n2 + b2 maxn), where bmax is the150 maximum number of artificial taxa for any subproblem. To summarize, the time151 complexity of TREE-QMC is lower than the quartet extraction step required for QMC,152 and the reduction in storage complexity should further improve efficiency by reducing153 cache misses. More importantly, our experimental evaluation showed TREE-QMC scaled154 to large data sets (with 1000 taxa and 1000 genes on which we could not run QMC) and155 achieved faster runtimes than ASTRAL-III, the leading method at the time (Fig. 2B in156 Molloy et al., 2022).157 Second, TREE-QMC normalizes quartet weights with respect to artificial taxa so158 that each gene tree gets one vote for every subset of four taxa present in the tree, although159 the vote could be split across the three possible quartets when gene trees are multi-labeled160 by artificial taxa introduced from the QMC framework, multiple gene copies (Legried161 et al., 2021), or multiple individuals/alleles per species (Rabiee et al., 2019). TREE-QMC162 can be run without normalization (n0 ), with uniform normalization ( n1), or with163 non-uniform normalization (n2), which leverages the structure implied by the introduction164 of artificial taxa to downweight quartets that are less relevant for resolving a particular165 subproblem (i.e., internal branch in the species tree) (Fig. 6a–d in Appendix A). In our166 simulation study, normalization was critical for accurate species tree estimation under167 challenging model conditions, such as high mean GTEE, with the n2 normalization scheme168 yielding the best accuracy (Fig. 4A in Molloy et al., 2022). Notably, TREE-QMC (with n2169 normalization) outperformed ASTRAL-III (Zhang et al., 2018) on challenging model170 conditions, ranging from large numbers of taxa (200–1000) to very high levels of ILS (60%)171 and mean GTEE (60%) (Figs. 2 and 4A in Han and Molloy, 2023, respectively).172 Based on the promising performance, we were motivated to further evaluate and173 improve upon the TREE-QMC method, especially with respect to gene tree estimation174 error and incompleteness, which are well-established obstacles to phylogenomics.175 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 8 HAN AND MOLLOY Normalizing Quartet Graph for Incomplete Gene Trees176 Our first improvement addresses the issue of incomplete gene trees, specifically their177 impact on quartet weight normalization. As previously mentioned, we say the weights are178 normalized if each gene tree votes once for each subset of four taxa present in the tree. We179 show that the normalization factors derived for the complete taxon set do not constitute a180 valid normalization scheme for incomplete gene trees (Fig. 6d–f in Appendix A), and181 sharing these normalization factors across all gene trees reduces species tree accuracy, as182 we show in our simulation study. In Appendix A, we introduce an efficient algorithm for183 computing the correct normalization values for each gene tree on the fly. This new184 approach does not increase the storage or time complexity of quartet graph construction185 and reduces the number of memory allocations and thus cache misses compared to earlier186 versions of the TREE-QMC code. Lastly, TREE-QMC now outputs a polytomy when187 there are no quartets for the subproblem; this can occur when the majority of taxa are188 missing from gene trees, for example the real and synthetic data sets assembled by Morel189 et al. (2023), which we use for benchmarking.190 Weighting Quartet Graph based on Gene Tree Branch Lengths and Support Values191 Our second improvement addresses the issue of GTEE by leveraging the quartet192 weighting schemes of Zhang and Mirarab (2022).193 Recall that the vertices in the quartet graph represent taxa from the subproblem and edges between taxa are weighted based on quartets implied by the input gene trees, with leaves mapped to taxa by the function L. Then, each quartet q = x, y|z, w implied by some gene tree T contributes two bad edges between each of the sibling pairs (i.e., (L(x), L(y)) and L(z), L(w)) as well as four good edges between each of the non-sibling pairs (i.e., ( L(x), L(z)), (L(x), L(w)), (L(y), L(z)) and ( L(y), L(w))), provided that the leaves are uniquely labeled (i.e., L(x) ̸= L(y) ̸= L(z) ̸= L(w)). In the unweighted algorithm, all good and bad edges have weight one; otherwise, good and bad edges inherit (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 9 the weight of the quartet that contributed them. In the length, support, and hybrid weighting schemes of Zhang and Mirarab (2022), the weights of the quartet q implied by gene tree T are defined as Wl(q) = exp  X e∈x,y→u, z,w→v −l(e)  , (0.1) Ws(q) = 1 − Y e∈u→v (1 − s(e)), and (0.2) Wh(q) = Wl(q) · Ws(q) (0.3) respectively, where l(e) denotes the length of branch e in gene tree T , s(e) denotes the194 support value of branch e in T , and u and v denote the anchor vertices of quartet195 q = x, y|z, w in T , with u closer to sibling pair x, y and v closer to sibling pair z, w. Note196 that u → v is the set of edges on the path between u and v, which corresponds to the197 internal branch of quartet q implied by T (Fig. 1a–b; also see Fig. 1 in Zhang and198 Mirarab, 2022). Similarly, x, y → u and z, w → v is the set of edges on the paths from each199 of the four leaf vertices to their respective anchor, which corresponds to the four terminal200 branches of quartet q implied by T .201 Figure 1 shows the impact of quartet weighting schemes on the quartet graph202 constructed by TREE-QMC, taking two gene trees from the Palaeognathae data set203 (Cloutier et al., 2019) as input, one with and one without identified homology errors204 (Simmons et al., 2022). For simplicity, we examine the quartet graph on the four Tinamou205 species, finding that the three weighting schemes successfully down-weight the quartet206 implied by the gene tree with homology errors compared to the gene tree without207 homology errors, which in turn impacts the best cut(s).208 The naive way to construct the quartet graph is to extract all weighted quartets209 from the gene trees (Fig. 1 in the Supplementary Materials); however, this preprocessing210 step is computationally intensive and requires a large amount of storage even for211 unweighted quartets, as previously mentioned. The challenge is building the quartet graph212 efficiently (i.e., directly from gene trees without extracting all possible quartets) while213 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 10 HAN AND MOLLOY a) UCE_chr7_5344_s gene tree, no homology error b) UCE_chr7_5446_s gene tree, homology error for Chicken and White-throated Tinamou c) MSA for UCE_chr7_5446_s indicating homology error 1 1 1 1 2 2 1 1 1 1 0 0 T WT EC C 0.033 0.844 0.844 0.033 0.877 0.877 0.844 0.033 0.033 0.844 0.000 0.000 Thicket Tinamou (T) White-throated Tinamou (WT) Elegant Crested Tinamou (EC) Chilean Tinamou (C) 0.509 1.0 0.915 0.334 1.0 1.0 1.0 0.998 1.0 0.966 1.0 1.0 Little Bush Moa Greater Rhea Lesser Rhea Chicken Ostrich Southern Cassowary Emu North Island Brown Kiwi Okarito Brown Kiwi Great Spotted Kiwi Little Spotted Kiwi u v Thicket Tinamou (T) Elegant Crested Tinamou (EC) Chilean Tinamou (C) 1.0 1.0 0.98 0.384 1.0 0.929 0.893 0.444 1.0 0.999 Little Bush Moa Greater Rhea Ostrich Southern Cassowary Emu North Island Brown Kiwi Okarito Brown Kiwi Great Spotted Kiwi Little Spotted Kiwi White-throated Tinamou (WT) 1.0 Chicken u v WT C EC T C T EC WT C WT EC T 6 4 6 3.442 1.754 1.8202 4 2 0.066 1.754 1.688 Cut corresponding to non-trivial bipartition Weighted Cut Values Unweighted Cut Values Good GoodBad Bad 0.0314 0.0354 0.0430 0.0599 1.000, 0.998 WT C EC T u v d) Quartets on four tinamous extracted from gene trees North Island Brown Kiwi Great Spotted Kiwi Little Spotted Kiwi Okarito Brown Kiwi Souther Cassowary Emu Greater Rhea Ostrich Little Bush Moa T EC C WT Chicken Consensus e) Quartets graphs with good and bad edges q1 not affected homology error wh(q1)=0.877 wl(q1)=0.877 ws(q1)=1.000 q2 affected by homology error wh(q2)=0.033 wl(q2)=0.074 ws(q2)=0.444 unweighted graph weighted graph f) Cuts of quartet graphs Goal: cut good but not bad edges EC C WT T u v 0.0776 0.0754 1.3887 0.444 0.0709 T WT EC C q2 q1 Fig. 1. Impact of hybrid weights on quartet graph and its best cut. a) Gene tree estimated on UCE chr7 5344 s by Cloutier et al. (2019), where u and v denote the anchor vertices for the quartet on the four Tinamou species: Thicket Tinamou (T), White-throated Tinamou (WT), Elegant Crested Tinamou (EC), and Chilean Tinamou (C). b) Gene tree estimated on UCE chr7 5446 s with long branch to Chicken and White-Throated Tinamou. c) MSA for UCE chr7 5446 s suggests homology errors for Chicken and White-Throated Tinamou, as noted by Simmons et al. (2022) (visualized with JalView; Waterhouse et al., 2009). d) Quartets q1 and q2 on the four Tinamous species extracted from gene trees and chr7 5344 s and chr7 5446 s, respectively, along with their hybrid, length, and support weights. e) Quartet graphs built from both q1 and q2 without and with hybrid weighting. f) A cut of the quartet graph partitions the species into two disjoint sets. There are only three cuts that yield non-trivial bipartitions on four species, corresponding to quartets. Cut quality is related to the ratio of good and bad edges removed (shown in table); the goal is to remove good edges but not bad edges thus preserving sibling relationships (Snir and Rao, 2010). For the unweighted graph, there are two best cuts (highlighted). For the weighted graph, there is one best cut corresponding to q1 (highlighted), which seems reasonable as q2 is affected by homology errors. weighting quartets based on gene tree branch lengths and support values. In Appendix B,214 we introduce efficient algorithms to construct the weighted quartet graph for a subproblem215 directly from gene trees in O((a2b + ab2n + b3n)k) time and O(a2 + ab + b2n) space, where216 a is the number of non-artificial taxa from the subproblem (referred to as singletons), b is217 the number of artificial taxa from the subproblem, and n is the total number of species,218 which we assume corresponds to the maximum number of leaves in any gene tree for219 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 11 simplicity (Theorem 4 in the Supplementary Materials).220 New features in weighted TREE-QMC221 Lastly, we implement several new features within weighted TREE-QMC.222 Bioconda. TREE-QMC is now installable via bioconda.223 Multi-labeled and non-binary gene trees. TREE-QMC now extends its framework for224 handling artificial taxa to accommodate multi-labeled gene trees as well as non-binary225 gene trees being given as input.226 Characters and “BP” mode. TREE-QMC now allows character data to be given as input;227 see the --chars and --char2tree options. TREE-QMC treats each character (i.e., site) in228 the data set as an unrooted tree with internal branches separating taxa assigned the same229 state from all other taxa. This functionality is useful because quartet-based summary230

are statistically consistent estimators of species trees for binary (0/1) characters231 evolving under the neutral Wright-Fisher model (Mendes and Hahn, 2017; Molloy et al.,232 2022) even when there is error and missingness provided it is unbiased (Han and Molloy,233 2024). Springer et al., 2020 refer to this approach as summary methods in “bipartition234 mode” or “BP mode”. If the --bp option is used with TREE-QMC instead of the --chars235 option, branch lengths in coalescent units will be computed under the neutral236 Wright-Fisher model with a fast maximum likelihood estimator, assuming a constant237 mutation rate (Molloy et al., 2022); otherwise, branch lengths are not computed for238 character data.239 Quadrapartition Quartet Support (QQS). TREE-QMC now enables users to output240 Quadrapartition Quartet Support (QQS) for the three quartet topologies “around”241 each internal branch (quadripartition) for unweighted or weighted quartets (Sayyari and242 Mirarab, 2016; Mirarab et al., 2024); see the --support or --supportonly options.243 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 12 HAN AND MOLLOY Additionally, maximum likelihood branch lengths in coalescent units are reported, along244 with the effective number (EN) of gene trees that provide information (i.e., quartets) for245 resolving the branch. Typically a branch is considered well supported when EN is246 sufficiently large and discordance follows signatures of ILS (i.e., the QQS value correspond247 to the branch in the estimated species tree is higher than the QQS values for the two248 alternative branch resolutions, which should be roughly equal to each other). Note that249 QQS and EN values are used by ASTRAL to compute branch support as the local250 posterior probability (Local PP) (Sayyari and Mirarab, 2016).251 Partitioned Coalescence Support (PCS). TREE-QMC now enables users to compute252 Partitioned Coalescence Support (PCS) for a specified focal branch in an input253 species tree (Gatesy et al., 2017, 2019). Specifically, the --pcsonly option returns the254 QQS values for each gene tree for each of the three possible resolutions of the focal branch255 in the species tree, enabling users to identify outlier gene trees that strongly or weakly256 favor a particular resolution of the focal branch compared to other gene trees. If applied to257 characters or gene trees sorted by their position in the genome, the output QQS plots can258 highlight interesting regions of the genome, for example a region of suppressed259 recombination in the avian tree of life (Fig. 2 in Mirarab et al., 2024) as well as enrichment260 for deep coalescence at the MHC locus in primates (Fig. PanGenomeS1 in Yoo et al.,261 2024).262

II. Simulation Study263 We now describe our performance study for evaluating the utility of weighted264 TREE-QMC on simulated data sets from prior studies focusing on gene tree estimation265 error and incompleteness, specifically the ASTRAL-II (Mirarab and Warnow, 2015),266 ASTRAL-III (Zhang et al., 2018), and Asteroid (Morel et al., 2023) data sets. Scripts are267 available at268 https://github.com/molloy-lab/tree-qmc-study/tree/main/han2024wtreeqmc.269 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 13 Species Tree Estimation Methods270 A large number of summary methods have been developed over the last decade. To271 make the simulation study manageable, we focused on weighted summary methods:272 weighted ASTRAL (i.e., ASTER) and weighted ASTRID. We also included Asteroid273 (Morel et al., 2023), as it is similar to (unweighted) ASTRID but with improvements for274 incomplete gene trees. All summary methods were run in default mode with the parameter275 settings summarized below. Detailed repository information, version/commit numbers, and276 software commands are available in the Supplementary Materials.277 For the Asteroid data sets, summary methods were run in unweighted mode, as the278 estimated gene trees did not include branch support and no alignments were available.279 Additionally, the Asteroid data sets had high rates of missing taxa, so we ran TREE-QMC280 on simulated data sets with and without updating the (n2) normalization values for281 incomplete gene trees (Fig. 6). The latter is denoted as n2 shared to indicate the same282 normalization values are shared across all gene trees regardless of their taxon sets.283 For the ASTRAL-II/III data sets, summary methods were run with their best284 weighting option. We ran ASTRAL with hybrid weights (denoted ASTER-wh), as this285 weighting scheme yielded the best accuracy in Zhang and Mirarab (2022). Motivated by286 these prior results, we ran TREE-QMC with hybrid weights (denoted TREE-QMC-wh287 or TQMC-wh) as well as in unweighted mode as an comparison. ASTRID does not288 implement the hybrid weighting scheme, and thus we ran ASTRID with support weights289 (denoted ASTRID-ws), as this weighting scheme yielded the best accuracy in Liu and290 Warnow (2023). Asteroid does not implement any weighting schemes; thus we ran291 Asteroid in default (unweighted) mode. For all weighted summary methods, we set the292 minimum/maximum values for gene tree branch support values to 0.333/1 for data sets293 with abayes support and 0/100 for data sets with bootstrap support, as indicated in the294 user manuals. On the smaller (i.e., Asteroid and ASTRAL-III) simulated data sets, we ran295 TREE-QMC with each of its three quartet weight normalization schemes: n0 (none), n1296 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 14 HAN AND MOLLOY (uniform), and n2 (non-uniform). The non-uniform (n2) yielded the best accuracy, so we297 used it for the larger (i.e., ASTRAL-II) simulated data sets and for biological analyses.298 Simulated Data Sets299 We evaluated summary methods using 3719 data sets simulated in three prior300 studies. These simulations were conducted by (1) simulating species trees under the Yule301 model (Yule, 1925), (2) simulating gene trees within the species tree under the302 Multi-Species Coalescent (MSC) (Rannala and Yang, 2003), which models ILS, and (3)303 simulating sequences down gene trees under standard models of molecular sequence304 evolution (e.g., the GTR model; Tavar´ e, 1986), which produces an MSA because there are305 no insertions or deletions. Different model conditions were created by varying a one or two306 model parameters at at time while keeping the others fixed. Multiple replicate data sets307 were simulated for each model condition. To enable comparisons of model conditions across308 studies, we empirically evaluated ILS, GTEE, and missingness (MISS) for each data set,309 reporting the average values across data sets with the same model condition. ILS was310 evaluated as the normalized RF distance between the true species tree and true gene trees,311 averaged across all gene trees. GTEE was evaluated as the normalized RF distance312 between each true and estimated gene tree, averaged across all gene trees. MISS was313 evaluated as the fraction of missing taxa from each gene tree, averaged across all gene314 trees. We summarize the model conditions and empirical properties below.315 Asteroid simulated data sets. The Asteroid data sets were simulated by Morel et al. (2023)316 to evaluate Asteroid in the context of missing data. After deleting sequences from gene317 MSAs, they estimated gene trees with maximum likelihood (ML) under the318 GTR+GAMMA4 model (Yang, 1993) using ParGenes (Morel et al., 2018). The empirical319 properties of the Asteroid data sets are given in Table S1 in the Supplementary Materials.320 Overall, the data sets were characterized by very high missingness (0.76–0.84), very low321 ILS (0.05–0.08), and medium GTEE (0.30–0.45), with some exceptions depending on the322 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 15 model parameter being varied. Model conditions (50 replicates each) were simulated by323 varying the324 1. effective population size and thus ILS level: 10 (ILS: 0), 50 million (default;325 ILS: 0.06), 100 million (ILS: 0.11), 500 million (ILS: 0.39), and 1 billion (ILS: 0.56)326 2. number of taxa: 25, 50 (default), 75, 100, 125, 150327 3. number of genes: 250, 500, 1000 (default), and 2000328 4. sequence length and thus GTEE level: 50 bp (GTEE: 46), 100 bp (GTEE: 0.35;329 default), 200 bp (GTEE: 0.25), and 500 bp (GTEE: 0.16)330 5. gene tree branch length scalar and thus GTEE level: 0.05 (GTEE: 0.55), 0.1331 (GTEE: 0.47), 1 (GTEE: 0.35; default) , 10 (GTEE: 0.42), 100 (GTEE: 0.70), and332 200 (GTEE: 0.77)333 6. missingness parameter: 0.5 (MISS: 0.74), 0.55 (MISS: 0.78), 0.6 (MISS: 0.81;334 default), 0.65 (MISS: 0.84), 0.7 (MISS: 0.86), and 0.75 (MISS: 0.88)335 where default indicates the value used in all other model conditions. There are 26 unique336 model conditions, after removing duplicates from the default parameters, and 1300337 replicate data sets in total.338 ASTRAL-III (S100) simulated data sets. The ASTRAL-III datasets with 100 taxa were339 simulated by Zhang et al. (2018) to evaluate ASTRAL-III in the context of GTEE. They340 estimated gene trees with maximum likelihood under the GTR+GAMMA model using341 FastTree-2 (Price et al., 2010). After which, branch support was estimated via342 bootstrapping with 100 replicates; later Zhang and Mirarab (2022) estimated abayes343 branch support (Anisimova et al., 2011) with IQ-TREE-2 (Minh et al., 2020). The344 empirical properties of the ASTRAL-III (S100) data sets are given in Table S2 in the345 Supplementary Materials. Overall, the data sets were characterized by medium-high ILS346 (0.46). Four model conditions (50 replicates each) were generated by varying the sequence347 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 16 HAN AND MOLLOY length and thus the GTEE level: 200 bp (GTEE: 0.56), 400 bp (GTEE: 0.42), 800 bp348 (GTEE: 0.31), and 1600 bp (GTEE: 0.23). The number of genes given to methods as input349 varied from 50, 200, 500, and 1000, yielding 16 model conditions and 800 replicate data350 sets in total.351 ASTRAL-II simulated data sets. The ASTER-II data sets (no missingness) were simulated352 by Mirarab and Warnow (2015). We estimated abayes branch support for these data sets353 (see Supplementary Materials for software command), following the recommendation by354 Zhang and Mirarab (2022). The empirical properties of the ASTRAL-III data sets are in355 Table S3 of the Supplementary Materials. Model conditions (50 replicates each) were356 simulated by varying the species tree height and speciation rate (deep/shallow) with a357 fixed number of taxa (200), which in turn varied ILS and GTEE. At the shortest species358 tree height (called 0.5 ×), the ILS and GTEE were 0.68/0.69 and 0.44/0.44, respectively.359 At the longest species tree height (called 5 ×), the ILS and GTEE were 0.09/0.21 and360 0.28/0.21, respectively. Conversely, the speciation rate and species tree height was fixed361 and the number of taxa was varied: 50, 100, 200, 500, 1000 (ILS: 0.31–0.35; GTEE:362 0.26–0.30). The number of genes given to methods as input varied from 50, 200, and 1000363 genes, yielding 33 model conditions and 1640 replicate data sets (note that 33 data sets364 were excluded based on criteria described in the Supplementary Materials).365 Evaluation Metrics366 We compared summary methods in terms of species tree error, runtime, quartet367 score, and mean branch support.368 Species tree error. We reported the normalized Robinson-Foulds (RF) error rate369 (Robinson and Foulds, 1981), which is equivalent to the false negative (FN) or false positve370 (FP) error rates when both the true and estimated species trees were binary (note that the371 FP/FN error rate are defined as the number of branches in the estimated/true species tree372 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 17 that are missing from the true/estimated species tree divided by the number of branches in373 the estimated/true species tree). All methods returned binary trees, except TREE-QMC,374 which returned non-binary trees for the Asteroid data sets when there were no quartets for375 resolving one or more subproblems. To ensure the RF distance was comparable to the376 other methods, we refined polytomies arbitrary. Additionally, we reported false positive377 (FP) and false negative (FN) error rates separately in the Supplementary Materials.378 For each model condition (with replicate data sets simulated under the same379 conditions), we tested for significant differences between some pairs of methods using380 two-sided, paired Wilcoxon signed-rank tests, as implemented in the R coin library. To be381 conservative, we corrected for ties between methods in two different ways: the Wilcoxon382

(Wilcoxon, 1949) and the Pratt method (Pratt, 1959), taking the higher of the two383 p-values. Lastly, we corrected for multiple comparisons based on Bonferroni correction,384 dividing the p-value by the number of tests performed for the experiment (note that we385 treated data sets simulated for the ASTRAL-II, ASTRAL-III (S100), and Asteroid studies386 as three different experiments). If the p-value was less than 0.05 after Bonferroni correction,387 we said there was a significant difference between two methods for the model condition.388 Runtime. We recorded the wall clock time (i.e., the total time for the method to run on389 an input data set from start to finish). When collecting runtime data, we gave all methods390 exclusive access to the compute node (architecture: AMD EPYC-7313 with 32 CPUs;391 maximum RAM: 64 GB). ASTRAL/ASTER did not complete within our maximum wall392 clock time of 20 hours for three data sets when given a single-thread, so we also ran it with393 16 threads. Otherwise, methods were run with a single thread. Branch support calculations394 were turned off to enable fair runtime comparisons.395 Quartet Score and Branch Support Values. For the quartet-based methods (i.e.,396 TREE-QMC and ASTRAL/ASTER), we reported the total number of quartets in the397 input gene trees that were satisfied by each estimated species tree as well as mean branch398 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 18 HAN AND MOLLOY support. Branch support was computed using ASTRAL’s Local PP (Sayyari and Mirarab,399 2016), which is based on QQS values from either unweighted or weighted quartets. Both400 total quartet score and mean branch support were computed with the same quartet401 weighting scheme as was used for species tree estimation.402

on Simulated Data Sets403 We now report the results of benchmarking (weighted) summary methods on three404 collections of data sets simulated for prior studies.405 Asteroid Data Sets. The Asteroid data sets were used to evaluate unweighted summary406

under very high levels of missing taxa (see Figure 2a–c as well as Figures S2–S7407 and Tables S4–S11 in the Supplementary Materials). We found that TREE-QMC (n2),408 without shared normalization factors, achieved lowered species tree (RF) error in 86% of409 data sets compared to TREE-QMC (n2 shared) (Fig. 2a). Moreover, TREE-QMC (n2)410 achieved significantly lower FN and FP error than no normalization (n0) for 17 and 16 of411 the 26 model conditions, respectively (Table S11). There were no significant differences412 between non-uniform (n2) and uniform (n1) normalization (Table S10); thus, we compared413 TREE-QMC (n2) against the other summary methods.414 Asteroid was the only method that achieved significantly lower FN error than415 TREE-QMC (n2), which occurred for 3 of the 26 model conditions (Table S6). Asteroid416 was not significantly better than TREE-QMC when considering FP error. The two417

tied on 20% of data sets in terms of RF error. Asteroid was better on 46% of data418 sets and TREE-QMC (n2) was better on the remaining 34%.419 In comparison to ASTRID and ASTRAL/ASTER, TREE-QMC (n2) outperformed420 or tied with them in terms of RF error on more than 75% of data sets (it was better on421 74% and 62% of data sets and tied on 9% and 14% of data sets, respectively) (Fig. 2b–c).422 TREE-QMC (n2) achieved significantly lower FN and FP error compared to ASTRID for423 18 of the 26 model conditions (Table S6–S8). Although TREE-QMC (n2) achieved424 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 19 significantly lower FP error than ASTRAL/ASTER on 17 of the 26 model conditions,425 respectively (Table S6–S8), this dropped to 3 model conditions when considering at FN426 error. Interestingly, TREE-QMC and ASTRAL/ASTER performed similarly in terms of427 mean branch support (Fig. 2d), but ASTRAL/ASTER achieved a higher total quartet428 score on 97% of data sets.429 Data sets from Morel et al. (2023) 0 50 100 TREE-QMC (n2) 0 50 100(shared n2) a) 0 50 100 TREE-QMC (n2) 0 50 100ASTRID b) 0 50 100 TREE-QMC (n2) 0 50 100ASTRAL c) 0.50 0.75 1.00 TREE-QMC (n2) 0.50 0.75 1.00ASTRAL d) Data sets from Zhang et al. (2018) w/ abayes support 0 20 40 TREE-QMC-wh (n2) 0 20 40unweighted (n2) e) 0 20 40 TREE-QMC-wh (n2) 0 20 40ASTRID-ws f) 0 20 40 TREE-QMC-wh (n2) 0 20 40ASTRAL-wh g) 0.8 0.9 1.0 TREE-QMC-wh (n2) 0.8 0.9 1.0ASTRAL-wh h) Data sets from Mirarab et al. (2015) w/ abayes support 0 20 40 TREE-QMC-wh (n2) % RF Error 0 20 40unweighted (n2) i) 0 20 40 TREE-QMC-wh (n2) % RF Error 0 20 40ASTRID-ws j) 0 20 40 TREE-QMC-wh (n2) % RF Error 0 20 40ASTER-wh k) 0.8 0.9 1.0 TREE-QMC-wh (n2) Local PP 0.8 0.9 1.0ASTER-wh l) Fig. 2. Rows of subplots show results on data sets simulated from three prior studies. Each dot in the scatter plot is a data set, with its position based on the evaluation metric (indicated below columns of subplots) for the method on the x-axis and the method on the y-axis. Red lines indicate there is no difference between methods. Data points are colored based on density, with lighter colors indicating higher density and darker colors indicating lower density. For % RF error, higher density of data points above the red line indicates that the version of TREE-QMC on the x-axis achieves lower error (and thus is more accurate) than the method on the y-axis. ASTRAL-III (S100) data sets430 The ASTRAL-III (S100) data sets were used to evaluate weighted summary431

under varying levels of GTEE (see Figure 2d–f as well as Figure S8 and432 Tables S12–S15 in the Supplementary Materials). We compared weighted methods given433 gene trees with abayes support values, as bootstrap support values resulted in worse434 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 20 HAN AND MOLLOY species tree accuracy for all methods, as previously observed by Zhang and Mirarab435 (2022). We found that TREE-QMC-wh (n2) lowered or tied for species tree (RF) error on436 57% and 27% of data sets, compared to TREE-QMC (n2), without hybrid weighting437 (Fig. 2e). Moreover, TREE-QMC-wh (n2) typically achieved the lowest mean species tree438 (RF) error per model condition compared to TREE-QMC’s other normalization schemes439 as well as the other summary methods (Fig. S8 and Table S12). TREE-QMC-wh (n2)440 achieved significantly lower RF error than ASTRID-ws and ASTRAL/ASTER-wh in 4 and441 3 of the 16 model conditions, respectively (Tables S13–S14), whereas ASTRID-ws and442 ASTRAL/ASTER-wh were not significantly better than TREE-QMC-wh (n2) for any443 model conditions. Overall, TREE-QMC-wh (n2) outperformed or tied with ASTRID-ws444 and ASTRAL/ASTER-wh on more than 75% of data sets (it was better on 53% and 44%445 of data sets and tied on 26% and 35% of data sets, respectively) (Fig. 2f–g). However,446 TREE-QMC-wh (n2) never achieved a higher quartet score than ASTRAL/ASTER-wh447 and had better mean branch support on just 5% of data sets (Fig. 2h). Note Asteroid448 should perform the same as (unweighted) ASTRID, as there are no missing taxa and that449 ASTRID performed worse than ASTRID-ws on these data sets in the study by Liu and450 Warnow (2023).451 ASTRAL-II data sets452 The ASTRAL-II data sets were used to evaluate weighted summary methods under453 varying levels of ILS and numbers of taxa (see Figure 2g–h as well as Figures S9–S10 and454 Tables S16–S23 in the Supplementary Materials). Similar to S100 data sets, we found that455 TREE-QMC-wh (n2) lowered or tied for species tree (RF) error of TREE-QMC (n2) in456 58% and 24% of data sets, compared to TREE-QMC (n2), without hybrid weighting457 (Fig. 2i). In comparison to ASTRID-ws and ASTRAL/ASTER-wh, TREE-QMC-wh (n2)458 outperformed or tied with them on more than 75% of data sets (it was better on 54% and459 46% of data sets and it tied on 30% and 37%) (Fig. 2j–k). Moreover, TREE-QMC (n2)460 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 21 achieved significantly lower RF error than ASTRID-ws and ASTRAL/ASTER-wh on 16461 and 10 of the 33 model conditions, respectively (Tables S19–S20 and S22–S23). The462 difference between these pairs of methods was the most pronounced for model conditions463 with 500 and 1000 taxa, where TREE-QMC had a clear advantage. However,464 TREE-QMC-wh (n2) only achieved a higher quartet score than ASTRAL/ASTER-wh on465 just 1 of the 1619 data sets and had better mean branch support on just 5% of data sets466 (Fig. 2l).467 Mirarab and Warnow (2015) also published the concatenation (CA-ML) trees468 (except the 1000-taxon, 1000-gene model condition). In comparison to CA-ML469 TREE-QMC (n2) was significantly better in terms of RF error on 19 of the 33 model470 conditions (Tables S18 and S21). The only condition CA-ML had a notable advantage over471 TREE-QMC (n2) was when the species tree height was longer (in generations) and472 speciation occurred closer towards the root (Table S18), although it is worth noting that473 the ILS level for these data sets was only 9%.474 Runtime475 Lastly, we evaluated the runtime of summary methods on the ASTRAL-II data sets476 with varying numbers of taxa and 1000 gene trees. At 1000 taxa, the fastest method was477 ASTRID-ws, followed by ASTRAL/ASTER-wh (16 threads), TREE-QMC (n2),478 TREE-QMC-wh (n2), and then ASTRAL/ASTER-wh (1 thread) at 0.0, 0.35, 0.36, 1.24,479 and 8.00 hours, respectively (Fig. 3a). As the number of taxa increased from 100 to 1000480 (so by a factor of 10), the runtime ratio between the weighted and unweighted TREE-QMC481 algorithms only increased from 2.9 to 3.4 (so by a factor of 1.2) (Fig. 3b). Thus, although482 TREE-QMC-wh (n2) runtime was slower than the original (i.e., unweighted) algorithm, it483 still scaled to large numbers of taxa, where it had the greatest advantage in terms of484 species tree accuracy compared to ASTRID-ws and ASTRAL/ASTER-wh.485 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 22 HAN AND MOLLOY 10 50 100 200 500 1000 0 2 4 6 8Runtime (h) Solid lines = 1 thread Dashed line = 16 threads a) 10 50 100 200 500 1000 Number of Taxa 0 1 2 3 4Weighted / unweighted b) ASTER-wh TREE-QMC-wh (n2) TREE-QMC (n2) ASTRID-ws Fig. 3. Runtime for ASTRAL-II data with varying numbers of taxa and 1000 gene trees with abayes support. Subfigure (a) shows the average runtime in hours for each of summary method (shaded region is standard error). Note that unweighted TREE-QMC (n2) is directly beneath ASTER-wh (16 threads). Subfigure (b) shows the ratio between the runtime for hybrid weighted TREE-QMC (n2) and the original TREE-QMC (n2) method. Biological Analyses486 We now describe our re-analysis of three published biological data sets using487 summary methods Asteroid, ASTRID, TREE-QMC, and ASTRAL/ASTER. Software488 commands are similar to those used in the simulation study. Scripts used to analyze489 biological data sets are available at490 https://github.com/molloy-lab/tree-qmc-study/tree/main/han2024wtreeqmc.491 Plant data set from Morel et al., (2023)492 We first re-analyzed the 81-taxa, 6176-gene plant data set curated by Morel et al.493 (2023) to evaluate their method Asteroid. This data set was characterized by extreme494 missing data, with the majority of gene trees have 9 or fewer taxa. Specifically, the number495 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 23 of gene trees with 30–35 taxa (MISS: 57–63%), 20-29 taxa (MISS: 64–75%), 10-19 taxa496 (MISS: 77–88%), and 4–9 taxa (MISS: 89–95%) was 4, 95, 490, and 5587, respectively.497 Overall, the mean percentage of missing taxa across the 6176 gene trees was 95%. The498 high rates of incompleteness make the plant data set more similar to a supertree analysis499 (rather than a coalescent analysis) so we ran summary methods in unweighted mode only.500 After species tree estimation, we estimated branch support using ASTRAL’s Local PP501 (Sayyari and Mirarab, 2016).502 We compared the estimated species trees to the (non-binary) NCBI reference503 taxonomy as well as the (binary) concatenation tree from Morel et al. (2023), which was504 missing 7 out of 56 (12.5%) internal branches in the NCBI reference tree. Overall,505 Asteroid, TREE-QMC, ASTRAL/ASTER, and ASTRID returned binary trees missing 14506 (18%), 15 (19%), 20 (26%), and 29 (37%) of branches in the concatenation tree,507 respectively. In comparison to the reference tree, Asteroid, TREE-QMC,508 ASTRAL/ASTER, and ASTRID were missing 6 (12.5%), 6 (12.15%), 17 (30%), and 24509 (43%) branches. All methods recovered Basal Angiosperm, although only two taxa were510 sampled from this group. All methods but ASTRID recovered Superastrids and Red Alga.511 Asteroid and TREE-QMC recovered Monocots as well as Superrosids, unlike512 ASTRAL/ASTER. Only ASTRAL/ASTER recovered Green Alga but this clade was513 separated from Red Alga by all of Monocots, which formed a highly imbalanced514 (caterpillar-like) topology. The branching ordered for Monocots differed from all other515

but was highly supported (Local PP), although the effective number (EN) of gene516 trees with information around some of these branch was quite low (e.g., 1–49), in which517 case ASTRAL recommends ignoring them. The ASTRAL/ASTER tree satisfied 3790,518 4815, and 26141 more quartets than TREE-QMC, Asteroid, and ASTRID trees (note that519 the fraction of quartets satisfied was above 0.89 for all methods).520 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 24 HAN AND MOLLOY Ananas Comosus Panicum Hallii Fil2 Oryza Barthii...Oryza Rufipogon[8] Corchorus Capsularis Citrus Clementina...Pistacia Vera[2] Chlamydomonas Reinhardtii Arabidopsis Thaliana Selaginella Moellendorffii Lupinus Angustifolius Panicum Hallii Hal2 Manihot Esculenta...Populus Trichocarpa[2] Setaria Italica...Setaria Viridis[2] Gossypium Raimondii Medicago Truncatula...Trifolium Pratense[2] Marchantia Polymorpha Theobroma Cacao Criollo...Theobroma Cacao Matina[2] Eragrostis Curvula...Eragrostis Tef[2] Daucus Carota Physcomitrella Patens Basal angiosperm[2] Cannabis Sativa Female...Rosa Chinensis[6] Leersia Perrieri Glycine Max...Vigna Radiata[4] Beta Vulgaris Sorghum Bicolor...Zea Mays[2] Chara Braunii Arabis Alpina Capsicum Annuum...Solanum Tuberosum[7] Oryza Brachyantha Cyanidioschyzon Merolae...Galdieria Sulphuraria[2] Vitis Vinifera Actinidia Chinensis Dioscorea Rotundata Brachypodium Distachyon Musa Acuminata Chondrus Crispus Cynara Cardunculus...Helianthus Annuus[2] Ostreococcus Lucimarinus Citrullus Lanatus...Cucumis Sativus[3] Brassica Napus...Brassica Rapa[3] Aegilops Tauschii...Triticum Urartu[2] Arabidopsis Halleri...Arabidopsis Lyrata[2] Hordeum Vulgare...Hordeum Vulgare Goldenpromise[2] 0.99 0.47 1.0 1.0 0.41 0.99 0.19 0.84 1.0 0.93 0.66 0.19 0.30 0.32 0.99 0.44 0.68 0.74 0.71 0.56 (0.68,0.17,0.13,27) (0.39,0.28,0.32,16) (0.79,0.11,0.09,517) (0.97,0.01,0.01,638) (0.36,0.27,0.35,21) (0.86,0.05,0.08,49) (0.24,0.34,0.41,11) (0.50,0.24,0.25,20) (0.75,0.11,0.12,300) (0.39, 0.29, 0.31, 277) (1.0,0.0,0.0,1) (0.30,0.31, 0.38,45) (0.34,0.36, 0.29,214) (0.42,0.48,0.08,15) (0.62,0.16,0.20,107) (0.36,0.35,0.27,56) (0.43,0.26,0.29,25) (0.46,0.19,0.34,26) (0.38,0.33,0.28,136) (0.50,0.03,0.46,13) Superasterids Green alga Monocots Red alga Superasterids Monocots Red alga Theobroma Cacao Criollo...Theobroma Cacao Matina[2] Eragrostis Curvula...Eragrostis Tef[2] Leersia Perrieri Daucus Carota Chara Braunii Beta Vulgaris Oryza Barthii...Oryza Rufipogon[8] Marchantia Polymorpha Chondrus Crispus Brachypodium Distachyon Musa Acuminata Panicum Hallii Fil2 Selaginella Moellendorffii Physcomitrella Patens Glycine Max...Vigna Radiata[4] Sorghum Bicolor...Zea Mays[2] Arabidopsis Thaliana Aegilops Tauschii...Triticum Urartu[2] Hordeum Vulgare...Hordeum Vulgare Goldenpromise[2] Manihot Esculenta...Populus Trichocarpa[2] Chlamydomonas Reinhardtii Cynara Cardunculus...Helianthus Annuus[2] Ostreococcus Lucimarinus Arabidopsis Halleri...Arabidopsis Lyrata[2] Dioscorea Rotundata Cyanidioschyzon Merolae...Galdieria Sulphuraria[2] Brassica Napus...Brassica Rapa[3] Panicum Hallii Hal2 Actinidia Chinensis Ananas Comosus Vitis Vinifera Lupinus Angustifolius Citrullus Lanatus...Cucumis Sativus[3] Citrus Clementina...Pistacia Vera[2] Capsicum Annuum...Solanum Tuberosum[7] Oryza Brachyantha Setaria Italica...Setaria Viridis[2] Cannabis Sativa Female...Rosa Chinensis[6] Arabis Alpina Gossypium Raimondii Medicago Truncatula...Trifolium Pratense[2] Basal angiosperm[2] Corchorus Capsularis 0.76 0.47 0.61 0.46 0.54 0.50 0.66 0.33 0.52 0.50 0.33 0.74 0.99 0.74 (0.41,0.25, 0.33,70) (0.39,0.32, 0.28,16) (0.42,0.34, 0.23,26) (0.37,0.28,0.34),30 (0.38,0.27, 0.33,37) (0.47,0.44,0.07,10) (0.42,0.30,0.26),26 (0.33,0.33,0.33,24) (0.36,0.32,0.30,90) (0.36,0.33, 0.29,68) (0.33,0.33,0.33,6) (0.41,0.32, 0.25,51) (0.41,0.27, 0.30,315) (0.46,0.34, 0.19,26) Superasterids Superrosids Monocots Red alga Malvids Fabids Medicago Truncatula...Trifolium Pratense[2] Arabidopsis Thaliana Physcomitrella Patens Gossypium Raimondii Citrullus Lanatus...Cucumis Sativus[3] Arabidopsis Halleri...Arabidopsis Lyrata[2] Chondrus Crispus Oryza Barthii...Oryza Rufipogon[8] Chara Braunii Manihot Esculenta...Populus Trichocarpa[2] Marchantia Polymorpha Beta Vulgaris Aegilops Tauschii...Triticum Urartu[2] Basal angiosperm[2] Brassica Napus...Brassica Rapa[3] Setaria Italica...Setaria Viridis[2] Daucus Carota Cynara Cardunculus...Helianthus Annuus[2] Panicum Hallii Fil2 Glycine Max...Vigna Radiata[4] Ananas Comosus Oryza Brachyantha Panicum Hallii Hal2 Chlamydomonas Reinhardtii Arabis Alpina Hordeum Vulgare...Hordeum Vulgare Goldenpromise[2] Cyanidioschyzon Merolae...Galdieria Sulphuraria[2] Eragrostis Curvula...Eragrostis Tef[2] Capsicum Annuum...Solanum Tuberosum[7] Dioscorea Rotundata Musa Acuminata Citrus Clementina...Pistacia Vera[2] Lupinus Angustifolius Selaginella Moellendorffii Cannabis Sativa Female...Rosa Chinensis[6] Brachypodium Distachyon Actinidia Chinensis Leersia Perrieri Ostreococcus Lucimarinus Corchorus Capsularis Theobroma Cacao Criollo...Theobroma Cacao Matina[2] Vitis Vinifera Sorghum Bicolor...Zea Mays[2] 0.35 0.30 0.34 0.50 0.54 0.74 0.33 0.52 0.66 0.58 0.52 0.47 0.74 0.98 0.48 (0.42,0.5,0.07,7) (0.34,0.36, 0.28,221) (0.33,0.33,0.33,24) (0.47,0.44,0.07,10) (0.38,0.27,0.33,37) (0.46,0.19,0.34,26) (0.33,0.33, 0.33,6) (0.46,0.22,0.31,6) (0.42,0.30,0.26,26) (0.37,0.32, 0.30,80) (0.35,0.31, 0.32,139) (0.39,0.32, 0.28,16) (0.41,0.25, 0.32,51) (0.41,0.32, 0.26,273) (0.41,0.34,0.24,11) Brachypodium Distachyon Oryza Barthii...Oryza Rufipogon[8] Manihot Esculenta...Populus Trichocarpa[2] Cannabis Sativa Female...Rosa Chinensis[6] Hordeum Vulgare...Hordeum Vulgare Goldenpromise[2] Eragrostis Curvula...Eragrostis Tef[2] Theobroma Cacao Criollo...Theobroma Cacao Matina[2] Lupinus Angustifolius Aegilops Tauschii...Triticum Urartu[2] Panicum Hallii Hal2 Glycine Max...Vigna Radiata[4] Arabis Alpina Vitis Vinifera Daucus Carota Physcomitrella Patens Arabidopsis Halleri...Arabidopsis Lyrata[2] Ostreococcus Lucimarinus Selaginella Moellendorffii Marchantia Polymorpha Ananas Comosus Chondrus Crispus Citrullus Lanatus...Cucumis Sativus[3] Sorghum Bicolor...Zea Mays[2] Chara Braunii Actinidia Chinensis Panicum Hallii Fil2 Medicago Truncatula...Trifolium Pratense[2] Cyanidioschyzon Merolae...Galdieria Sulphuraria[2] Corchorus Capsularis Setaria Italica...Setaria Viridis[2] Basal angiosperm[2] Capsicum Annuum...Solanum Tuberosum[7] Cynara Cardunculus...Helianthus Annuus[2] Oryza Brachyantha Arabidopsis Thaliana Beta Vulgaris Citrus Clementina...Pistacia Vera[2] Brassica Napus...Brassica Rapa[3] Gossypium Raimondii Leersia Perrieri Dioscorea Rotundata Chlamydomonas Reinhardtii Musa Acuminata 0.18 1.0 0.12 0.97 0.00 1.0 0.65 0.59 1.0 0.96 0.90 0.76 0.79 0.21 0.99 0.93 0.86 0.69 0.61 0.92 0.33 0.22 0.15 0.07 0.13 0.00 0.09 0.03 0.04 a) Asteroid d) ASTRAL/ASTER b) ASTRID c) TREE-QMC Superrosids Fig. 4. Species trees estimated from data curated by Morel et al. (2023). Subfigures (a), (b), (c), and (d) show species trees we estimated using Asteroid, ASTRID, TREE-QMC, ASTRAL/ASTER, respectively. Blue branches are not in the RAxML tree from Morel et al. (2023). Branch support (i.e., Local PP) is shown for branches that disagree with the RAxML tree. QQS values and EN are given in parentheses. All species trees are shown on the same reduced the leaf set created based on subtrees in the RAxML that were shared by all estimated species trees (leaf labels have the form X...Y [N] indicating a subtree in the RAxML tree with N taxa, starting taxon X and ending with taxon Y ). The reduction in leaf labels from 81 to 43 indicates that the estimated trees are on many of the lower order relationships. Avian (Neognathae+Palaeognathae) data set from Wu et al., (2024)521 Next, we re-analyzed the avian data set curated by Wu et al. (2024a), which522 included 5756 coding sequences (CDs), 4871 introns, and 2384 intergenic segments for 124523 avian species. Later, homology errors in CDs were identified by Springer and Gatesy524 (2024), although these errors did not appear biased towards any particular set of taxa (Fig.525 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 25 1B in Springer and Gatesy, 2024 and Fig. 1A in Wu et al., 2024b). We used IQ-TREE-2526 (version 2.3.5) to compute abayes support on the best maximum likelihood (ML) gene527 trees estimated by Wu et al. (2024a) and then estimated species trees using Asteroid,528 ASTRID, TREE-QMC, and ASTRAL/ASTER (with branch support calculations were529 turned off to enable fair runtime comparisons). Species tree estimation was restricted to530 the 10,627 CDs and introns, as these markers may be more robust to orthology errors than531 intergenic regions (personal communication with Dr. Mark Springer and Dr. John Gatesy).532 All three weighted summary methods were run with and without their different weighting533 schemes. For unweighted analyses only, we also filtered the input gene trees by removing534 outlier taxa with TreeShrink version 1.3.9 (options: -m “per-species” -q 0.05) (Mai and535 Mirarab, 2018). We confirmed that TreeShrink successfully removed the homology errors536 from the CDs highlighted in Figure 1B in Springer and Gatesy (2016). After species tree537 estimation, we estimated branch support as Local PP using the same quartet weighting538 scheme that was used for species tree estimation (Sayyari and Mirarab, 2016). For539 weighted methods, we used the original gene trees when computing Local PP and QQS540 values; for unweighted methods only, we used gene trees with outlier taxa removed.541 Comparison of unweighted versus weighted summary methods. ASTRAL/ASTER542 produced four species trees: one tree for hybrid weighting, one tree for length weighting,543 one tree for support weighting OR no weighting after taxon filtering, and one tree for no544 weighting before taxon filtering. The ASTRAL/ASTER with length weighting differed545 from unweighted ASTRAL/ASTER after taxon filtering (Fig. S11d in Supplementary546 Materials) by just one internal branch. Similarly, unweighted ASTRAL/ASTER (after547 taxon filtering) differened from ASTRAL/ASTER with hybrid weighting (Fig. 5d) by just548 three branches, but this difference jumped to six branches when taxon filtering was not549 used for the unweighted method. In contrast, TREE-QMC produced two species trees: one550 tree for hybrid and length weighting (Fig. 5c) and one tree for all other ways of running551 TREE-QMC (Fig. S11b in Supplementary Materials). The two TREE-QMC trees differed552 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 26 HAN AND MOLLOY from each other on just three branches. ASTRID also produced two species trees: one tree553 for length weighting and one tree for all other ways of running ASTRID (Fig. 5b). The two554 ASTRID trees differed on 18 internal branches, with the length weighting tree being much555 farther from the NJst tree (thus we did not consider the length weighting ASTRID tree in556 further comparisons). Lastly, Asteroid recovered the same tree as ASTRID if taxon filtering557 was used; otherwise, the Asteroid tree differed from the ASTRID tree by four branches.558 Comparison to published species trees. Wu et al. (2024a) published two main species trees:559 one from applying NJst (Liu and Yu, 2011) (Fig. 5a and Fig. 2 in Cloutier et al., 2019) and560 one from applying RAxML (Stamatakis, 2014). The NJst and RAxML trees, which561 differed by 13 internal branches, were given different collection inputs than our methods,562 as they were given all three marker types (CDs, introns, and intergenetic regions) after563 removing 30% outlier genes (i.e., genes with the highest quartet distance between the best564 ML gene tree and the estimated species tree were removed). Our estimated species trees565 differed from the NJst (and RAxML) trees by 8 (16), 8 (13), and 16 (18) branches for566 ASTRID/Asteroid, TREE-QMC-wh (n2), and weighted ASTRAL/ASTER-wh,567 respectively. Thus, the TREE-QMC tree was closer to NJst and RAxML trees than the568 ASTRAL/ASTER tree by 8 and 5 branches, respectively.569 When considering major clades, all methods tested in our study recovered570 Palaeognathae/Neognathae, Galloanserae/Neoaves, Aequorlitornithes, and Hieraves.571 ASTRAL/ASTER did not recover Aquaterraves (because it clustered the two pigeons into572 Telluraves) or Columbaves, unlike ASTRID and TREE-QMC. No methods recovered573 Listusilvae, Coraciiorphae, or Australaves, but the ASTRID and TREE-QMC trees574 required fewer edit moves (either contractions/refinement moves or subtree prune and575 regraft moves) to achieve monophyly for these groups than the ASTRAL/ASTER tree. On576 the other hand, ASTRAL/ASTER always achieved the highest quartet score and its577 branches that disagreed with NJst or RAxML had high branch support (Local PP),578 although the QQS values around these branches were all close to one third, which is579 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 27 consistent with a rapid radiation for Neoaves, as demonstrated by many prior studies580 (Jarvis et al., 2014; Wu et al., 2024a; Stiller et al., 2024). Lastly, ASTRAL/ASTER581 recovered a highly imbalanced (i.e., caterpillar-like) branching pattern for Aquaterraves582 minus the two pigeons, which was striking compared to the other estimated species trees.583 Runtime. ASTRID completed in just a few seconds, weighted TREE-QMC completed in584 11 minutes, weighted ASTRAL/ASTER-wh completed in 2 hours, and ASTRAL/ASTER585 completed in 2.6 hours (both with and without taxon filtering the input using TreeShrink).586 The wall clock time of ASTRAL/ASTER dropped to 5–8 minutes when using 16 threads587 instead of one thread.588 Avian (Palaeognathae) data set from Cloutier et al., (2019)589 Lastly, we re-analyzed the avian (Palaeognathae) data set curated by Cloutier et al.590 (2019), with 3,158 ultraconserved elements (UCEs) for 15 species. Systematic homology591 errors were identified in 105 of the UCEs by Simmons et al. (2022), with the two impacted592 taxa, White-Throated Tinamou and Chicken, clustering together. We used IQ-TREE-2 to593 compute abayes support on the best maximum likelihood (ML) gene trees estimated by594 Cloutier et al. (2019) and then estimated species trees with the four summary methods,595 comparing weighting schemes as well as filtering practices for unweighted methods. All596 analyses produced one of two trees, called “A” and “B,” which differed by a single focal597 branch (quadrapartition). In tree A, the focal branch split Tinamous plus598 Kiwis+Cassowary+Emu from Rheas plus Ostrich+Chicken, whereas in tree B, Tinamous599 swapped with Rheas so they were on the same side of the focal branch as Chicken600 (Fig. S12 in Supplementary Materials). Notably, the 105 UCEs with homology errors601 strongly favored tree B, but the hybrid quartet weighting scheme reduced the magnitude of602 support for tree B, as shown by partitioned coalescent support (PCS) analysis around the603 focal branch (Fig. S13 in Supplementary Materials).604 Weighted TREE-QMC and ASTRAL (hybrid and length) returned tree A, as did605 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 28 HAN AND MOLLOY Aequorlitor -nithes LitusilvanaeColumbavesAustra -laves Coraciimorphae Palaeognathae AquaterravesTelluraves Dalmatian pelican...Red-throated loon[11] Speckled mousebird Oilbird MacQueen's bustard...Red-crested turaco[2] Cuckoo-roller Downy woodpecker...Paradise jacamar[3] Great spotted kiwi Okarito brown kiwi Woodland kingfisher Carmine bee-eater Sunbittern...White-tailed tropicbird[3] Galloanserae[19] Little spotted kiwi Blue-throated roller Common ostrich Ruff...Bar-tailed godwit[3] Hieraves[7] Killdeer...Couble-striped thick-knee[2] White-throated tinamou Common cuckoo Brown mesite...Yellow-throated sandgrouse[2] Chimney swift...Tawny frogmouth[5] Chuck-will's-widow Great-crested grebe...American flamingo[2] Bar-tailed trogon Dark-eyed junco...Kakapo[37] Black-legged kittiwake Rhinoceros hornbill...Hoopoe[2] Pigeon...Band-tailed pigeon[2] Great potoo Saker falcon...Peregrine falcon[2] Hoatzin Grey-crowned crane...Inaccessible island rail[2] Little buttonquail Greater rhea Red-legged seriema Aequorlitor -nithes Columbaves Palaeognathae AquaterravesTelluraves Killdeer...Couble-striped thick-knee[2] Red-legged seriema Dark-eyed junco...Kakapo[37] Common cuckoo White-throated tinamou Grey-crowned crane...Inaccessible island rail[2] Speckled mousebird Downy woodpecker...Paradise jacamar[3] Little spotted kiwi Chuck-will's-widow Rhinoceros hornbill...Hoopoe[2] Great-crested grebe...American flamingo[2] Little buttonquail Great spotted kiwi Hoatzin Hieraves[7] Cuckoo-roller Galloanserae[19] Great potoo Bar-tailed trogon Black-legged kittiwake Chimney swift...Tawny frogmouth[5] Oilbird Saker falcon...Peregrine falcon[2] Woodland kingfisher Okarito brown kiwi Sunbittern...White-tailed tropicbird[3] Carmine bee-eater Ruff...Bar-tailed godwit[3] Brown mesite...Yellow-throated sandgrouse[2] Blue-throated roller Common ostrich Pigeon...Band-tailed pigeon[2] Greater rhea Dalmatian pelican...Red-throated loon[11] MacQueen's bustard...Red-crested turaco[2] 1 1 0.48 0 1 1 0 1 1 1 0.97 0.99 0.99 0 0.42 1 (0.38, 0.27,0.33) (0.36,0.32,0.31) (0.33,0.33,0.32) (0.32, 0.33,0.34) (0.39,0.35,0.25) (0.37,0.31,0.31) (0.33,0.35,0.31) (0.83,0.07,0.08) (0.51, 0.30,0.18) (0.38, 0.30,0.30) (0.34,0.31,0.33) (0.35,0.32,0.32) (0.34,0.33,0.32) (0.36,0.37,0.26) (0.33,0.32,0.34) (0.38, 0.35,0.25) Aequorlitor -nithes Palaeognathae Telluraves Common cuckoo Blue-throated roller MacQueen's bustard...Red-crested turaco[2] Little buttonquail Little spotted kiwi Woodland kingfisher Greater rhea Great potoo Okarito brown kiwi Red-legged seriema Oilbird White-throated tinamou Brown mesite...Yellow-throated sandgrouse[2] Cuckoo-roller Black-legged kittiwake Pigeon...Band-tailed pigeon[2] Sunbittern...White-tailed tropicbird[3] Hieraves[7] Hoatzin Bar-tailed trogon Great spotted kiwi Saker falcon...Peregrine falcon[2] Common ostrich Dark-eyed junco...Kakapo[37] Galloanserae[19] Rhinoceros hornbill...Hoopoe[2] Dalmatian pelican...Red-throated loon[11] Speckled mousebird Chuck-will's-widow Great-crested grebe...American flamingo[2] Carmine bee-eater Killdeer...Couble-striped thick-knee[2] Grey-crowned crane...Inaccessible island rail[2] Downy woodpecker...Paradise jacamar[3] Chimney swift...Tawny frogmouth[5] Ruff...Bar-tailed godwit[3] 0.99 0.99 0.99 0.76 0.95 1 1 1 0.99 0.97 0.99 1 0.99 0.76 1 0.99 0.99 0.81 1 0.93 (0.35,0.30,0.34) (0.35,0.33,0.31) (0.35,0.30,0.34) (0.34, 0.31, 0.34) (0.34,0.33,0.31) (0.39,0.35,0.25) (0.37, 0.36, 0.25) (0.83,0.08,0.07) (0.35,0.34,0.29) (0.34,0.31,0.33) (0.35,0.34,0.30) (0.51,0.30,0.18) (0.34,0.33,0.31) (0.34,0.32,0.33) (0.38, 0.30,0.30) (0.34,0.33,0.32) (0.35,0.29,0.34) (0.34,0.33,0.31) (0.39, 0.31, 0.28) (0.34, 0.33, 0.32) Aequorlitor -nithes Columbaves Palaeognathae AquaterravesTelluraves Cuckoo-roller Red-legged seriema Common cuckoo Brown mesite...Yellow-throated sandgrouse[2] Woodland kingfisher Little buttonquail Grey-crowned crane...Inaccessible island rail[2] Okarito brown kiwi Galloanserae[19] Saker falcon...Peregrine falcon[2] Ruff...Bar-tailed godwit[3] Speckled mousebird Common ostrich Chuck-will's-widow Hieraves[7] Great-crested grebe...American flamingo[2] Killdeer...Couble-striped thick-knee[2] Bar-tailed trogon Greater rhea Great potoo White-throated tinamou Hoatzin Rhinoceros hornbill...Hoopoe[2] Carmine bee-eater Pigeon...Band-tailed pigeon[2] Dark-eyed junco...Kakapo[37] Downy woodpecker...Paradise jacamar[3] Black-legged kittiwake Sunbittern...White-tailed tropicbird[3] Dalmatian pelican...Red-throated loon[11] Oilbird Blue-throated roller Chimney swift...Tawny frogmouth[5] Great spotted kiwi MacQueen's bustard...Red-crested turaco[2] Little spotted kiwi 0 1 0.99 0.54 1 1 1 0.28 0.47 1 1 0.92 0 0 0.99 1 0 (0.32, 0.33,0.34) (0.36,0.29,0.33) (0.35,0.33,0.31) (0.33,0.33,0.32) (0.72,0.12,0.15) (0.37, 0.30,0.31) (0.51,0.27,0.20) (0.33,0.34,0.32) (0.33,0.33,0.32) (0.40, 0.30,0.29) (0.35,0.32,0.31) (0.34,0.31,0.33) (0.33,0.35,0.31) (0.32,0.39,0.27) (0.35,0.32,0.32) (0.36,0.31,0.31) (0.33, 0.28,0.38) a) NJST from Wu et al., (2024) c) weighted TREE-QMC (hybrid) d) weighted ASTRAL/ASTER (hybrid) b) ASTRID Fig. 5. Species trees estimated from data curated by Wu et al. (2024a). Subfigure (a) is the NJst tree from Figure 2 in Wu et al. (2024a). Subfigure (b), (c), and (d) show species trees we estimated on CDs and introns using ASTRID (no weighting and support weighting scheme), weighted TREE-QMC (hybrid or length weighting scheme), weighted ASTRAL (hybrid weighting scheme), respectively. Red and blue branches are not in the NJst and RAxML trees from Wu et al. (2024a), respectively. Branch support (i.e., Local PP) is shown for branches that disagree with the NJst and/or RAxML trees. QQS values (in parentheses) are based on unweighted quartets (with taxon filtering) for subfigure (b) and hybrid weighted quartets (without taxon filtering) for subfigures (c) and (d). All species trees are shown on the same reduced leaf set created based on subtrees in the NJst that were shared by all estimated species trees. Leaf labels have the form X...Y [N] indicating the subtree in the NJst tree has N taxa, starting taxon X and ending with taxon Y . The reduction in leaf labels from 124 to 36 indicates that the estimated trees agree on most of the lower order relationships. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 29 unweighted TREE-QMC and ASTRAL after filtering the UCE data set (either removing606 the two impacted taxa from the 105 gene trees with homology errors or removing the gene607 trees entirely). Otherwise, tree B was returned. Branch support (i.e., Local PP) for the two608 resolutions of the focal branch was low ( <0.5; see Fig. S12). Filtering taxa or entire gene609 trees resulted in tree A satisfying 2507 and 2269 more quartets than tree B respectively.610 Likewise, tree A satisfied 972.5 more quartets than tree B when using hybrid weights;611 otherwise it satisfied 1783 fewer quartets. The concatenation trees from UCEs estimated612 by Cloutier et al. (2019) and Springer and Gatesy (2016) (who removed homology errors)613 resolved the focal branch in the same way as tree A but differed from tree A on one614 branch, placing Emu and Cassowary as sister to Tinamous instead of Kiwis.615 Discussion616 Recently, Zhang and Mirarab (2022) introduced quartet weighting schemes to617 improve the robustness of the popular method ASTRAL to gene tree estimation error.618 Although their quartet weighting schemes are straightforward to describe, they present619 challenges for developing efficient summary methods. Here, we show that these weighting620 schemes can be effectively integrated into the Quartet Max Cut framework of Snir and Rao621 (2010), introducing weighted TREE-QMC. Our theoretical and empirical study shows that622 weighted TREE-QMC is a highly promising approach and may even be the leading623 weighted summary method developed to date.624 Scalability. Weighted TREE-QMC only increases the time complexity of quartet graph625 construction by a factor of b compared to the original (unweighted) algorithm, where626 b = O(n) is the number of artificial taxa for the subproblem (see proof of Theorem 1 in627 Han and Molloy, 2023). Interestingly, the empirical runtime of weighted TREE-QMC did628 not grow much compared in our computational experiments (Fig. 3b). We conjecture this629 difference between our theoretical and empirical study is due to the number b of artificial630 taxa per subproblem growing sub-linearly in the number n species or behaving more like a631 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 30 HAN AND MOLLOY constant factor in practice. Under the latter assumption, both weighted and unweighted632 TREE-QMC have storage complexity O(n2) and time complexity O(n2 log(n)k) with some633 assumptions on subproblem decomposition (Theorem 5 in Supplementary Materials).634 Taken together, these theoretical and empirical results on the scalability of weighted635 TREE-QMC are promising.636 It is worth noting that weighted ASTRAL/ASTER (with 16 threads) was faster637 than weighted TREE-QMC (with one thread) in our computational experiments (Fig. 3a).638 Weighted ASTRAL uses multi-threading to speed up the quartet score calculation,639 parallelizing the computation across gene trees, which contribute independently to the640 quartet score. Parallelism across gene trees could also be exploited by weighted641 TREE-QMC’s during graph construction algorithm, and future versions of the software642 should take advantage of multi-threading and vector instructions available on modern643 processors.644 Missingness and supertrees. Quartet weight normalization was an important algorithmic645 development from the original TREE-QMC method. Here, we showed that incomplete646 gene trees require their own normalization factors to be computed, which means that647 normalization cannot be correctly performed on quartets that have been extracted from648 gene trees, as required by QMC or related algorithms. We introduced a more efficient649 algorithm for computing the correct normalization factors for each gene tree and found in650 simulations that correct normalization is critical for species tree accuracy when gene trees651 are incomplete (Fig. 2a). Because the simulated data sets were characterized by extreme652 missingness and relatively low ILS, our results suggest that TREE-QMC is a promising653 supertree method. We found similar results for the plant data set, for which TREE-QMC654 successfully recovered many established clades in the NCBI taxonomy similar to Asteroid655 and concatenation. The ASTRAL/ASTER tree, in contrast, was farther from the NCBI656 taxonomy and concatenation tree, displaying a highly imbalanced (caterpillar-like)657 branching pattern towards the root. Although this topology was better supported (Local658 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 31 PP) and achieved a higher quartet score, our results on simulated data sets do not suggest659 that higher mean branch support and/or total quartet score are always indicators of660 greater species tree accuracy.661 Species tree accuracy and robustness to gene tree errors. As in Zhang and Mirarab (2022),662 we found that the hybrid quartet weighting scheme typically improved the accuracy of663 TREE-QMC compared to the unweighted version (Figs. 2e and 2f). A limitation of the664 simulated data sets is that error is due to insufficient phylogenetic signal or other issues665 arising during gene tree estimation. Such errors may follow the “MSC+Error+Support”666 model, under which Zhang and Mirarab (2022) prove weighted ASTRAL/ASTER is667 statistical consistent, unlike the unweighted method. However, these errors do not fully668 encompass those observed on real data sets due to incorrect orthology inference, poor669 quality multiple sequence alignments, and violations in the model of molecular sequence670 evolution. Although the effectiveness of methods on biological data is difficult to discern, it671 is promising that the hybrid and length weighting schemes made TREE-QMC and672 ASTRAL/ASTER more robust to systematic homology errors in the Palaeognathae data673 set (Cloutier et al., 2019). However, we caution that there was still insufficient signal in674 this data set to resolve this part of the avian evolutionary history (also see Simmons et al.,675 2016).676 Species tree topological stability. We also reanalyzed a recent avian data set curated by677 Wu et al. (2024a). A major finding of this study is the separation of Neoaves into two678 major clades: Telluraves (land birds) and Aquaterraves (waterbirds and relatives). Springer679 and Gatesy (2024) later identified homology errors in the data set and showed that680 ASTRAL-III on the cleaned data set did not recover the Telluraves and Aquaterraves split.681 Subsequently, Wu et al. (2024a) responded that NJst on cleaned data set still recovered682 the split and suggested that ASTRAL-III was less stable than NJst.683 Our reanalysis also suggests that ASTRAL/ASTER is less stable than the other684 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 32 HAN AND MOLLOY methods, as ASTRAL/ASTER returned different trees when outlier taxa were removed685 with TreeShrink, unlike ASTRID and TREE-QMC, which returned the same tree686 regardless of taxon filtering. When considering the five ways of running weighted summary687

(with different inputs and different weighting schemes), ASTRAL/ASTER688 produced four different species trees (differing by up to six branches), whereas TREE-QMC689 and ASTRID each produced two trees (differing by three branches). This may have690 interesting implications for branch support estimated by bootstrapping procedures.691 The trees produced by ASTRAL/ASTER were distinct from those produced by692 TREE-QMC and the other methods, including concatenation. Not only did693 ASTRAL/ASTER not recover the Telluraves/Aquaterraves split, it also displayed a highly694 imbalanced (caterpillar-like) branch patterning at the deeper nodes of Aquarterraves minus695 pigeons. Although the ASTRAL/ASTER trees were better supported and achieved the696 highest quartet scores, as previously mentioned, our results on simulated data sets do not697 suggest that higher mean branch support and total quartet score are always indicators of698 greater species tree accuracy. Interestingly, the TREE-QMC tree was closer to the trees699 produced by concatenation as well as the distance methods NJst and ASTRID. It is700 possible that distance methods and TREE-QMC’s divide-and-conquer framework might701 behave more similarly to each other than ASTRAL/ASTER’s heuristic search, although702 we note that TREE-QMC was more robust to missing data than ASTRID in our703 simulation study (Fig. 2b). It is also worth noting that TREE-QMC does optimize the704 same criterion score as ASTRAL/ASTER because it downweights quartets based on705 subproblem decomposition, which could result in a scoring function that is more local or706 diffuse across the species tree, akin to likelihood scores. Future work should further study707 ASTRAL/ASTER and TREE-QMC to better understand the impact of objective function708 and search heuristics on species tree accuracy.709 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 33 Appendix A: Efficient Quartet Weight Normalization given Incomplete710 Gene Trees711 In this appendix, we describe our approach for normalizing quartet weights for a712 subproblem given incomplete gene trees.713 Preliminaries714 A quartet is an unrooted, binary tree on four leaves {a, b, c, d}; thus it takes on one715 of three possible topologies, denoted a, b|c, d or a, c|b, d or a, d|b, c. Let L denote the716 mapping function from leaf vertices to taxon labels for the subproblem. To avoid717 confusion, we use lower case letters to indicate leaf vertices and upper case letters to718 indicate leaf labels, also referred to as taxa or species . This distinction is necessary because719 the introduction of artificial taxa results in multiple leaves being labeled by the same label720 (Fig. 6a–c).721 Normalization, Importance Values, and Forest Data Structure722 A key observation in the original TREE-QMC method is that gene trees can723 contribute more than one quartet per subset of four species when they are multi-labeled724 (Han and Molloy, 2023). Thus, we wish to construct the normalized quartet graph, with725 quartet weights normalized so that each gene tree votes once for each subset of four labels726 present in the tree. To implement normalization, TREE-QMC assigns an importance727 value I(x) to each leaf vertex x and then normalization factor for each quartet728 q = x, y|z, w by multiplying the importance values of its leaves:729 I(q) = I(x, y, z, w) = I(x)I(y)I(z)I(w). (0.4) Later, we will use similar notation to indicate importance values for two or three leaves730 being multiplied together (e.g., I(x, y) = I(x)I(y) and I(x, y, z) = I(x)I(y)I(z)). It is easy731 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 34 HAN AND MOLLOY to show that if732 X a∈LeafSet(T ): L(x)=X I(x) = 1 (0.5) for every taxon label X, then the importance values form a valid normalization scheme,733 specifically734 X a,b,c,d∈LeafSet(T ): {L(x),L(y),L(z),L(w)}={X,Y,Z,W } I(x, y, z, w) = 1 (0.6) for selections of four unique taxon labels {X, Y, Z, W } from the subproblem. Note that in735 the unnormalized case (n0), the importance values are set to one for all leaves, which736 does not satisfy Equation 0.5 beyond the first subproblem in which all leaves are uniquely737 labeled.738 TREE-QMC implements two normalization schemes: uniform (n1) and739 non-uniform (n2). The latter requires storing the hierarchical relationships between740 artificial taxa and real taxa for the subproblem as a collection of disjoint trees741 F = {F1, F2, . . . , Fx}, which we refer to as a forest data structure. The roots of trees in742 F are taxa in the current subproblem, whose importance values sum to one following743 Equation 0.5. Leaf vertices of the trees in F are singleton taxa for the subproblem, which744 appear at most once in any gene tree; all other vertices represent artificial taxa for the745 subproblem.746 To construct the quartet graph for the subproblem, the leaves of each gene tree are labeled by the roots of trees in F based on ancestor-descendant relationships (Fig. 6a–c). Importance values can be computed using F for both normalization schemes. In uniform normalization (n1), all leaves in the same tree Fi ∈ F are assigned the same importance value equal to the inverse of the number of leaves in Fi. In non-uniform normalization (n2), the importance values are computed via the recurrence: I(x) = ( 1 if x is a root of some tree Fi ∈ F I(x.parent) · (1/x.parent.outdegree) otherwise (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 35 where x.parent.outdegree returns the number of outgoing edges for the parent of x in Fi747 (Fig. 6c). The recurrence can be implemented via a preorder traversal of each tree in F.748 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 10 X 0 9 X’ 1 2 6 7 8 10 X Y 3 4 5 Y’ 1 2 8 10 Y Z 6 7 X Z’ 3 4 5 Y’ 6 X 7 Z’ a) Divide Phase of QMC 1,Z|2,Y 1,Y|2,Z 6 6 1/2 1/2 1,2|8,10 1,8|2,10 1,2|Z,Y 1,10|2,8 0 0 1 0 0 1 0 0 Quartet n0 n1 2/3 1/3 0 0 1 0 n2 2 4 1/6 1/3 0 0 0 0 0 0 0 0 n0 n1 2/9 2/9 0 0 0 0 n2 2 4 1/3 2/3 0 0 0 0 0 0 0 0 n0 n1 1/2 1/2 0 0 0 0 n2 T1 using F1 T2 using F1 T2 using F2 d) Quartet weights with (n1,n2) and without (n0) normalization for artificial taxa Y 3 4 5 20 1 8 Z 6 7 X 0 9 1/1 1/3 1/3 1/3 1/11/1 1/1 1/1 1/1 1/3 1/3 1/3 1/6 1/6 Data structures for normalizing quartet weights for subproblem c) Forest data structure F1 with no missing taxa f) Forest data structure F2 missing taxa 3, 6, & 8 Y 4 5 20 1 Z 7 X 0 9 1/1 1/2 1/2 1/11/1 1/1 1/1 1/2 1/2 1/4 1/4 If no artifical taxa, normalization is not needed as missing taxa removes all quartets Otherwise, normalization requires each gene tree to cast 1 vote per subset of 4 taxa - fails for highlighted columns 1 10 2 8 Y(4) Z(7) Z(X(0)) Z(X(9)) Y(5) Y(3) Z(6) 1 10 2 Y(4)Z(7) Z(X(9))Y(5) Z(X(0)) Example of complete vs. incomplete gene trees with leaves relabeled for subproblem b) Gene tree T1 with no missing taxa e) Gene tree T2 (same as T1 without 3, 6, & 8) Fig. 6. Normalization for artificial taxa. a) The divide-and-conquer algorithm starts with a set of input gene trees on 11 taxa. Taxa are split into two disjoint subsets at each subproblem generated during the divide phase. Take the subproblem on taxon set {1, 2, 8, 10, Y, Z} as an example, where X, Y , and Z are artificial taxa representing {0, 9}, {3, 4, 5}, and {6, 7, X}, respectively. b) A complete gene tree T 1 with leaves relabeled for the subproblem. c) Forest data structure F 1 for the subproblem. d) Quartets on {1, 2, Y, Z} implied by gene tree T 1 have weights summing to 1 after normalization if using forest data structure F 1 to compute importance values. e) An incomplete gene tree T 2 with leaves relabeled for the subproblem (note that T 2 is the same as T 1 but missing taxa 3, 6, and 8). Two quartets (1, 7|2, 5 and 1, 7|2, 4) in gene tree T 2 correspond to 1, Z |2, Y after relabeling for the subproblem. If using F 1 importance values for n2 normalization, I(1, 7, 2, 5) = I(1, 7, 2, 4) = 1 · 1/3 · 1 · 1/3 = 1/9. For n1 normalization, the importance values are I(1, 7, 2, 5) = I(1, 7, 2, 4) = 1 · 1/4 · 1 · 1/3 = 1/12. Thus, the total weight of quartet 1 , Z|2, Y is 1/6 and 2/9 after n1 and n2 normalization, respectively. Repeating this process for the other two topologies 1, 2|Z, Y and 1, Y |2, Z does not produce weights that sum to 1 (highlighted columns). f) Forest data structure F 2 for subproblem but updated to exclude taxa 3, 6, and 8. Quartets on {1, 2, Y, Z} implied by gene tree T 2 have weights summing to 1 after normalization if using F 2. (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 36 HAN AND MOLLOY Normalization for Incomplete Gene Trees749 If there are no missing taxa, the same importance values can be used for all gene750 trees (Fig. 6a–d). However, if a gene tree is incomplete, the importance values for subsets751 of four taxa may not sum to one, violating the requirement of normalization (Fig. 6d–e).752 To normalize correctly, importance values need to be computed for each gene tree based on753 its missing taxa (Fig. 6f). Early versions of the TREE-QMC code addressed the issue of754 incomplete gene trees by copying the forest data structure for each gene tree and then755 deleting the missing taxa while updating the importance values accordingly. This naive756 approach suffers due to frequent memory allocation and thus cache misses.757 To improve performance, the current version of the TREE-QMC computes758 importance values from the forest data structure for the subproblem for each gene tree on759 the fly in two steps (see Algorithm 1 in the Supplementary Materials). First, the outdegree760 of internal nodes in the forest is updated to account for missing data via postorder761 traversal. Second, the importance values are updated according to the recurrence via762 preorder traversal. The time complexity to get the normalization values per gene tree is763 thus O(n), which is lower than constructing the quartet graph from the gene tree (see764 Appendix B). To implement normalization, we maintain the forest data structure for each765 subproblem on the stack during the divide phase of TREE-QMC. As the storage766 complexity of each forest is O(n) and the maximum depth of recursion is O(n), the total767 storage complexity for the forest data structure is O(n2), which is the same as the quartet768 graph on the complete species set. Thus, our normalization approach does not increase the769 time or storage complexity of TREE-QMC.770 Appendix B: Efficient Normalized Quartet Graph Construction771 leveraging Gene Tree Branch Lengths and Support V alues772 In this appendix, we describe how to construct the (normalized) quartet graph773 using the weighting schemes of Zhang and Mirarab (2022) without explicitly extracting all774 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 37 weighted quartets from the gene trees (Fig. S1 in Supplementary Materials). All referenced775 Definitions and Lemmas, Theorems, and Corollaries in this section are provided in Section776 2 of the Supplementary Materials.777 Strategy for Computing Weighted Bad Edges and Good for Subproblem778 Recall that vertices in the quartet graph represent taxa from the subproblem and779 edges between taxa are weighted by quartets implied by the gene trees, whose leaves map780 to taxa. If the leaves of some quartet q = a, b|c, d implied by a gene tree T have distinct781 labels, then q contributes two bad edges (L(x), L(y)) and ( L(z), L(w)) as well as four782 good edges (L(x), L(z)), (L(x), L(w)), (L(y), L(z)), and ( L(y), L(w)) to the quartet783 graph, provided that L(x) ̸= L(y) ̸= L(z) ̸= L(w); otherwise q does not contribute any784 good or bad edges. Good and bad edges inherit the weight of the quartet that contributed785 them, with the weight set to one in the unweighted algorithm and set using786 Equations 0.1, 0.2, or 0.3 for the length, support, or hybrid weighting schemes,787 respectively. Observe that under the hybrid weighting scheme, branch support weighting788 can be turned off by setting all support values to one and similarly length weighting can be789 turned off by setting all branch lengths to zero. Moreover, non-binary gene trees can be790 handled by simply refining polytomies and setting the resulting edges to have support zero791 and length zero (all other edges have support one and length zero in the unweighted792 algorithm or inherit their respective values from the unrefined tree in the weighted793 algorithm). As non-binary trees as well as length and support weights can be handled via794 the hybrid weighting scheme, we focus on hybrid weights henceforth.795 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 38 HAN AND MOLLOY For ease of computation, we split the calculation of the hybrid weight into two parts: Wh(q) = Wl(q) · Ws(q) = Wl(q) ·  1 − Y e∈u→v (1 − s(e))  = Wl(q) − Wl(q) · Y e∈u→v (1 − s(e)) = Wl(q) − ∆W(q) (0.7) where u and v denote the anchor vertices for quartet q in gene tree T (Fig. 1). Then, the total weight of bad edges between X and Y in the quartet graph computed from gene tree T becomes B(X, Y ) = X x,y,z,w ∈LeafSet(T ) L(x)̸=L(y)̸=L(z)̸=L(w) L(x)=X,L(y)=Y q=L(x),L(y)|L(z),L(w) Wh(q)I(q) = X Wl(q)I(q)  − X ∆W(q)I(q)  = X Wl(q)I(q)  − ∆B(X, Y ) (0.8) where I(q) corresponds to the importance value for normalizing the weight of each quartet796 q (Equation 0.6) and the summations loop over all ways of selecting four unique leaves797 x, y, z, w from LeafSet(T ) such that the implied quartet has distinct labels and separates798 two leaves labeled by X and Y from the other two leaves. Observe that an algorithm for799 computing the left term ∆B(X, Y ) given gene tree T can be used to compute the right800 term after setting all branch support values in T to zero. Summing the left and right terms801 together for all pairs of taxa ( X, Y ) and repeating across all k gene trees yields the bad802 edges B in the weighted quartet graph for the subproblem. The good edges, denoted G,803 can be computed in a similar fashion.804 We therefore focus on algorithms for computing ∆ B(X, Y ) and ∆ G(X, Y ) from a805 single unrooted gene tree T henceforth. We assume that T is binary (because any806 non-binary tree can be transformed into a binary tree with support values of zero on its807 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 39 refined edges, as previously mentioned) with n leaves mapping to s = a + b taxon labels for808 the subproblem, where a denotes the number of singleton taxa and b denotes the number809 of artificial taxa (Fig. 6). For ease of computation, we root T arbitrarily, letting h denote810 its height h (note that in our complexity analysis h and n are taken to be the maximum811 values across the k gene trees). We let t.p, t.s, t.l, and t.r denote the parent, sibling, left812 child, and right child of vertex t, respectively. We let below(t) denote the set of leaves813 “below” the subtree rooted at vertex t (i.e., they are descendants of vertex t). Conversely,814 we let above(t) denote the set of leaves “above” or “outside” the subtree rooted at t (i.e.,815 they are in the set LeafSet(T) \ below(t)).816 To provide some intuition for computing ∆B(X, Y ) and ∆ G(X, Y ), consider a pair817 of leaves x and y in LeafSet(T ) such that L(x) = X and L(y) = Y . Let x → y denote the818 (unique) path between x and y in T , and let F denote the set of subtrees obtained from819 deleting edges on the path from x → y as well as their endpoints from T . To form a bad or820 good quartet for leaf pair x, y, we must select two leaves z and w from F such that821 L(z) ̸= L(w) ̸= L(x) ̸= L(y) (Fig. 6A in Han and Molloy, 2023). All pairs of z and w822 selected from the same subtree in F correspond to a “bad quartet for X, Y ” because x, y823 are siblings (one anchor vertex is on the path from x → y and the other anchor vertex is824 within the subtree from which z, w were selected). Likewise, all pairs of z and w selected825 from different subtrees in F correspond to a “good quartet for X, Y ” because x, y are NOT826 siblings (both anchor vertices of the quartet are on the path x → y).827 Building upon this idea, we propose to compute ∆B(X, Y ) and ∆ G(X, Y ) by828 summing up the weights from bad and good quartets for X, Y at their anchor vertices in829 gene tree T . Specifically, we let ∆B t(X, Y ) and ∆ Gt(X, Y ) correspond to the weights from830 bad and good quartets for X, Y such that one leaf is below the left child of t (call x), one831 leaf is below the right child of t (call y), and one of x, y is labeled by X and the other is832 labeled by Y (w.l.o.g. X is the taxa labeling x and Y is the taxon labeling y; see833 Definitions 1 and 2). Because T can be multi-labeled (e.g., due to the introduction834 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 40 HAN AND MOLLOY artificial taxa), we must consider all pairs x, y such that x ∈ {L(l) = X : l ∈ below(t.l)}835 and y ∈ {L(l) = Y : l ∈ below(t.r)}. To complete the bad or good quartet, we must836 consider all ways of selecting two additional leaves z and w from T such that837 L(z) ̸= L(w) ̸= L(x) ̸= L(y). This allows us to break the computation of ∆B t(X, Y ) and838 ∆Gt(X, Y ) into six cases:839 (a) z, w ∈ below(t.l),840 (b) z ∈ below(t.l) and w ∈ above(t.p),841 (c) z, w ∈ above(t.p),842 (d) z ∈ above(t.p) and w ∈ below(t.r) (symmetric with case (b)),843 (e) z, w ∈ below(t.r) (symmetric with case (a)),844 (f) z ∈ below(t.l) and w ∈ below(t.r).845 It is easy to see that case (c) is not relevant for the computation of ∆ G(X, Y ) because z846 and w are drawn from the same subtree off the path from any x, y under consideration,847 thus forming a bad quartet. We define equations to compute the weights from each of these848 cases and then sum the relevant cases to get ∆B t(X, Y ) and ∆ Gt(X, Y ) (Theorems 1 and849 2). Our equations rely on the precomputation of six auxiliary values, summarized below.850 • w1 X(t) and w1 X(t): singlet auxiliary values give the sum of length weights from851 leaves labeled by taxon X in the subtree below and above vertex t, respectively852 (Definition 3 and 4). Singlets are useful for computing the branch lengths of terminal853 edges to leaves labeled X in the quartets. They can be computed for all taxa X and854 for all vertices t ∈ V (T ) in O((a + b)n) time and require O((a + b)n) storage to855 access later (Lemmas 1 and 2).856 • w2 X,Y (t) and w2 X,Y (t) : doublet auxiliary values give the sum of weights857 corresponding to pairs of leaves labeled by taxa X, Y in the subtree below and above858 vertex t (Definition 5 and 6). Doublets are useful for accessing the complete branch859 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 41 lengths of terminal edges for leaf pairs labeled X, Y attached to the same in anchor860 in a quartet and computing the branch support for the internal edge. They can be861 computed for all pairs of taxa X, Y and for all vertices t ∈ V (T ) in O((a + b)2n) time862 and require O((a + b)2n) storage to access later (Lemmas 3 and 4).863 • wx|z,w X,Y (t) : bad triplet auxiliary values give the sum of weights from three leaves864 x, z, w in the subtree below t such that the triplet corresponds to a bad quartet for865 X, Y (Definition 7). Bad triplets are useful for accessing the branch lengths for the866 three terminal edges and the branch support for the internal edge in a bad quartet.867 They can be computed for all pairs of taxa X, Y and all vertices t ∈ V (T ) in868 O((a + b)4n) time and require O((a + b)2n) storage to access later (Lemma 5).869 • wx,z|w X,Y (t) : good triplet auxiliary values give the sum of weights from three leaves870 x, z, w in the subtree below t such that the triplet corresponds to a good quartet for871 X, Y (Definition 8). Good triplets can be computed for all pairs of taxa X, Y and for872 all vertices t ∈ V (T ) in O((a + b)4n) time and require O((a + b)2n) storage to access873 later (Lemma 6).874 After precomputation, ∆B t(X, Y ) and ∆ Gt(X, Y ) can be computed in constant time.875 Summing across all vertices in T enables us to build ∆ B(X, Y ) and ∆ G(X, Y ) for all X, Y876 in O((a + b)2n) time and O((a + b)2) space (Theorem 3). This does not exceed the877 requirements for precomputation, putting the final time and space complexity for878 computing ∆B and ∆G given gene tree T at O((a + b)4n) and O(a2n + abn + b2n),879 respectively.880 Efficient Algorithm for Computing Weighted Bad Edges and Good for Subproblem881 Our core contribution is reducing the final time and space complexity to882 O(a2b + ab2h + b3n) and O(a2 + ab + b2n), respectively (Theorem 4). We summarize our883 strategy below.884 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 42 HAN AND MOLLOY Reducing the time complexity of triplet auxiliary values by aggregating singletons. Our first885 goal is to reduce the time complexity of computing triplets. To achieve this goal, we886 introduce a new label 1 to represent all singleton taxa and then extend the computation of887 singlets and doublets to this new label in the natural way (Definitions 9-11 and Corollaries888 1–3), effectively aggregating the computation across all singleton taxa. We then reorganize889 the computation of triplet recurrences to leverage the aggregated singlets and doublets,890 along with some other tricks. This drops the time complexity of the triplet recurrences891 from O((a + b)2) to O(b) (Lemmas 7–10 and Corollaries 4–5 and 7–8), enabling us to892 compute triplet auxiliary values for all vertices in T and for all taxa X, Y in O((a + b)2bn)893 time (Corollaries 6 and 9) instead of O((a + b)4n) time.894 Reducing the storage and time complexity of auxiliary values for singleton taxa. Our895 second goal is to reduce the storage and time complexity of auxiliary values for singleton896 taxa. To achieve this goal, we break the computation of ∆ B(X, Y ) and ∆ G(X, Y ) into897 three cases (i)–(iii), depending on whether zero, one, or both of X and Y are singleton898 taxa. In case (i), there are multiple leaves labeled X, each with its own path to the root of899 T , and similarly for the leaves labeled Y . This means there are multiple vertices (lowest900 common ancestors of X, Y ) in gene tree T where we need to compute ∆ Bt(X, Y ) and901 ∆Gt(X, Y ), with the number of vertices is O(n). In case (ii), there are multiple leaves902 labeled Y , each with its own path to the root, but there is only one path from the unique903 leaf labeled X to the root. Again, there are multiple vertices in gene tree T where we need904 to compute ∆B t(X, Y ) and ∆ Gt(X, Y ), but the number of vertices is O(h). Lastly, in case905 (iii), there is one path from the unique leaf x labeled X to the root and one path from the906 unique leaf y labeled Y to the root. So there is just one vertex in gene tree T where we need907 to compute ∆B t(X, Y ) and ∆ Gt(X, Y ). The same analysis holds for the auxiliary values.908 We can exploit path uniqueness for singleton taxa to reduce the storage and time909 complexity of auxiliary values. Take the singlet auxiliary value w1 X as an example.910 Fortunately, the singlet w1 X is only evaluated at t.l and t.r when computing ∆ Bt(X, Y ) and911 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint WEIGHTED TREE-QMC 43 ∆Gt(X, Y ) (Theorems 1 and 2). If X is a singleton taxon, there are three possible cases at912 vertex t: (1), X labels a leaf below t.l so w1 X(t.l) exists and w1 X(t.r) is undefined, (2) X913 labels a leaf below t.r so w1 X(t.l) is undefined and w1 X(t.r) exists, and (3) X labels a leaf914 above t so both w1 X(t.l) and w1 X(t.r) are undefined. ∆B t(X, Y ) and ∆ Gt(X, Y ) are also915 undefined in case (3), so we only need to consider cases (1) and (2). In both cases, we can916 compute ∆Bt(X, Y ) and ∆ Gt(X, Y ) in a postorder traversal, so that we do not need to917 store w1 X(t.r) or w1 X(t.l) after processing vertex t. Before continuing the traversal from t, we918 need to compute and store w1 X(t) to access it at t.p. In case (1), we can compute w1 X(t)919 from w1 X(t.l), and in case (2) we can compute w1 X(t) from w1 X(t.r). Thus, if X is a singleton,920 we can simplify the recurrence for w1 X (denoted p1 X), and we can overwrite the old value at921 t.l or t.r with the new value at t. This reduces the storage and time complexity of w1 X(t)922 from O(an + bn) to O(a + bn) (Corollary 10). Similar results can be given for the other923 auxiliary quantities (Corollaries 11–14).924 Ultimately, the auxiliary values with at least one singleton taxon can be computed925 on the fly during weighted graph construction and do not affect time or storage complexity926 (proof of Theorem 4). Thus, we only need to compute and store auxiliary values for927 artificial taxa during precomputation. This reduces the time complexity of precomputation928 from O((a + b)4n) to O(b3n) and reduces the storage complexity from O((a + b)2n) to929 O(b2n) (proof of Theorem 4).930 Wrapping up. Lastly, the computation of ∆ Bt(X, Y ) and ∆ Gt(X, Y ) can be improved by931 applying the previously introduced tricks based on aggregating singleton taxa and splitting932 up the computation into three cases based on whether X and/or Y are singletons933 (Corollaries 15 and 16). After putting everything together, it takes O((a2b + ab2h + b3n)k)934 time and O((a + b)2 + b2n) = O(a2 + ab + b2n) space to construct the weighted quartet935 graph for the subproblem given k gene trees, including precomputation (Theorem 4). The936 only work remaining for the subproblem is seeking a max cut. The heuristic (Dunning937 et al., 2018) we use to seek a max cut of the weighted quartet graph does not exceed the938 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 44 HAN AND MOLLOY time complexity of graph construction (Theorem 2 in Han and Molloy, 2023).939 Conflicts of Interest940 None declared by the authors.941 Funding942 YH and EKM were supported by the State of Maryland. All computational943 experiments were performed on the Center for Bioinformatics and Computational Biology944 (CBCB) compute cluster at the University of Maryland, College Park.945 Acknowledgements946 We thank Dr. Mark Springer and Dr. John Gatesy for sharing the estimated gene947 trees used in their reply to Wu et al. (2024a) and for helpful discussions about the data.948 We thank the editor Dr. Matthew Hahn and the anonymous reviewers for detailed949 feedback that improved our study.950 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint

45 References951 Anisimova, M., M. Gil, J.-F. Dufayard, C. Dessimoz, and O. Gascuel. 2011. Survey of952 Branch Support Methods Demonstrates Accuracy, Power, and Robustness of Fast953 Likelihood-based Approximation Schemes. Systematic Biology 60:685–699.954 Avni, E., R. Cohen, and S. Snir. 2014. Weighted quartets phylogenetics. 64:233–242.955 Cloutier, A., T. B. Sackton, P. Grayson, M. Clamp, A. J. Baker, and S. V. Edwards. 2019.956 Whole-genome analyses resolve the phylogeny of flightless birds (Palaeognathae) in the957 presence of an empirical anomaly zone. Systematic Biology 68:937–955.958 Cunha, T. J., J. D. Reimer, and G. Giribet. 2021. Investigating Sources of Conflict in Deep959 Phylogenomics of Vetigastropod Snails. Systematic Biology 71:1009–1022.960 Dunning, I., S. Gupta, and J. Silberholz. 2018. What works best when? a systematic961 evaluation of heuristics for Max-Cut and QUBO. 30:421–624.962 Gatesy, J., R. W. Meredith, J. E. Janecka, M. P. Simmons, W. J. Murphy, and M. S.963 Springer. 2017. Resolution of a concatenation/coalescence kerfuffle: partitioned964 coalescence support and a robust family-level tree for mammalia. Cladistics 33:295–332.965 Gatesy, J., D. B. Sloan, J. M. Warren, R. H. Baker, M. P. Simmons, and M. S. Springer.966 2019. Partitioned coalescence support reveals biases in species-tree methods and detects967 gene trees that determine phylogenomic conflicts. Molecular Phylogenetics and968 Evolution 139:106539.969 Han, Y. and E. K. Molloy. 2023. Improving quartet graph construction for scalable and970 accurate species tree estimation from gene trees. Genome Research Pages gr–277629.971 Han, Y. and E. K. Molloy. 2024. Quartets enable statistically consistent estimation of cell972 lineage trees under an unbiased error and missingness model. Algorithms for Molecular973 Biology 18:19.974 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 46 REFERENCES Hosner, P. A., B. C. Faircloth, T. C. Glenn, E. L. Braun, and R. T. Kimball. 2016.975 Avoiding missing data biases in phylogenomic inference: An empirical study in the976 landfowl (Aves: Galliformes). Mol. Biol. Evol. 33:1110–1125.977 Islam, M., K. Sarker, T. Das, R. Reaz, and M. S. Bayzid. 2020. STELAR: a statistically978 consistent coalescent-based species tree estimation method by maximizing triplet979 consistency. BMC Genomics 21:136.980 Jarvis, E. D., S. Mirarab, et al. 2014. Whole-genome analyses resolve early branches in the981 tree of life of modern birds. Science 346:1320–1331.982 Lafond, M. and C. Scornavacca. 2019. On the weighted quartet consensus problem.983 769:1–17.984 Legried, B., E. K. Molloy, T. Warnow, and S. Roch. 2021. Polynomial-time statistical985 estimation of species trees under gene duplication and loss. 28:452–468.986 Liu, B. and T. Warnow. 2023. Weighted astrid: fast and accurate species trees from987 weighted internode distances. Algorithms for Molecular Biology 18:6.988 Liu, L. and L. Yu. 2011. Estimating species trees from unrooted gene trees. 60:661–667.989 Liu, L., L. Yu, and S. V. Edwards. 2010. A maximum pseudo-likelihood approach for990 estimating species trees under the coalescent model. BMC Evolutionary Biology 10:302.991 Mai, U. and S. Mirarab. 2018. Treeshrink: fast and accurate detection of outlier long992 branches in collections of phylogenetic trees. BMC Genomics 19:272.993 Mendes, F. K. and M. W. Hahn. 2017. Why concatenation fails near the anomaly zone.994 Systematic Biology 67:158–169.995 Minh, B. Q., H. A. Schmidt, O. Chernomor, D. Schrempf, M. D. Woodhams, A. von996 Haeseler, and R. Lanfear. 2020. IQ-TREE 2: New Models and Efficient Methods for997 Phylogenetic Inference in the Genomic Era. Molecular Biology and Evolution998 37:1530–1534.999 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint

47 Mirarab, S., R. Reaz, M. S. Bayzid, T. Zimmermann, M. S. Swenson, and T. Warnow.1000 2014. ASTRAL: genome-scale coalescent-based species tree estimation. Bioinformatics1001 30:i541–i548.1002 Mirarab, S., I. Rivas-Gonz´ alez, S. Feng, J. Stiller, Q. Fang, U. Mai, G. Hickey, G. Chen,1003 N. Brajuka, O. Fedrigo, G. Formenti, J. B. W. Wolf, K. Howe, A. Antunes, M. H.1004 Schierup, B. Paten, E. D. Jarvis, G. Zhang, and E. L. Braun. 2024. A region of1005 suppressed recombination misleads neoavian phylogenomics. Proceedings of the National1006 Academy of Sciences 121:e2319506121.1007 Mirarab, S. and T. Warnow. 2015. ASTRAL-II: coalescent-based species tree estimation1008 with many hundreds of taxa and thousands of genes. Bioinformatics 31:i44–i52.1009 Molloy, E. K., J. Gatesy, and M. S. Springer. 2022. Theoretical and practical1010 considerations when using retroelement insertions to estimate species trees in the1011 anomaly zone. Systematic Biology 71:721–740.1012 Molloy, E. K. and T. Warnow. 2018. To include or not to include: The impact of gene1013 filtering on species tree estimation methods. 67:285–303.1014 Morel, B., A. M. Kozlov, and A. Stamatakis. 2018. ParGenes: a tool for massively parallel1015 model selection and phylogenetic tree inference on thousands of genes. Bioinformatics1016 35:1771–1773.1017 Morel, B., T. A. Williams, and A. Stamatakis. 2023. Asteroid: a new algorithm to infer1018 species trees from gene trees under high proportions of missing data. Bioinformatics 29.1019 Nute, M., J. Chou, E. K. Molloy, and T. Warnow. 2018. The performance of1020 coalescent-based species tree estimation methods under models of missing data.1021 19(Suppl 5):286.1022 One Thousand Plant Transcriptomes Initiative. 2019. One thousand plant transcriptomes1023 and the phylogenomics of green plants. Nature 574:679–685.1024 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 48 REFERENCES Pratt, J. W. 1959. Remarks on zeros and ties in the wilcoxon signed rank procedures.1025 Journal of the American Statistical Association 54:655––667.1026 Price, M. N., P. S. Dehal, and A. P. Arkin. 2010. FastTree 2 – approximately1027 maximum-likelihood trees for large alignments. PLOS ONE 5:1–10.1028 Rabiee, M., E. Sayyari, and S. Mirarab. 2019. Multi-allele species reconstruction using1029 ASTRAL. 130:286–296.1030 Rannala, B. and Z. Yang. 2003. Bayes estimation of species divergence times and ancestral1031 population sizes using DNA sequences from multiple loci. Genetics 164:1645–1656.1032 Robinson, D. and L. Foulds. 1981. Comparison of phylogenetic trees. Mathematical1033 Biosciences 53:131–147.1034 Roch, S. and M. Steel. 2015. Likelihood-based tree reconstruction on a concatenation of1035 aligned sequence data sets can be statistically inconsistent. 100:56–62.1036 Sayyari, E. and S. Mirarab. 2016. Fast coalescent-based computation of local branch1037 support from quartet frequencies. 33:1654–1668.1038 Sayyari, E., J. B. Whitfield, and S. Mirarab. 2017. Fragmentary Gene Sequences1039 Negatively Impact Gene Tree and Species Tree Reconstruction. Molecular Biology and1040 Evolution 34:3279–3291.1041 Simmons, M. P. and J. Gatesy. 2021. Collapsing dubiously resolved gene-tree branches in1042 phylogenomic coalescent analyses. Molecular Phylogenetics and Evolution 158:107092.1043 Simmons, M. P., D. B. Sloan, and J. Gatesy. 2016. The effects of subsampling gene trees1044 on coalescent methods applied to ancient divergences. Molecular Phylogenetics and1045 Evolution 97:76–89.1046 Simmons, M. P., M. S. Springer, and J. Gatesy. 2022. Gene-tree misrooting drives conflicts1047 in phylogenomic coalescent analyses of palaeognath birds. Molecular Phylogenetics and1048 Evolution 167:107344.1049 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint

49 Snir, S. and S. Rao. 2010. Quartets MaxCut: A divide and conquer quartets algorithm.1050 7:704–718.1051 Snir, S. and S. Rao. 2012. Quartet MaxCut: A fast algorithm for amalgamating quartet1052 trees. 62:1–8.1053 Springer, M. S. and J. Gatesy. 2016. The gene tree delusion. 94:1–33.1054 Springer, M. S. and J. Gatesy. 2024. A new phylogeny for aves is compromised by1055 pervasive misalignment and homology problems. Proceedings of the National Academy1056 of Sciences USA 121:e2406494121.1057 Springer, M. S., E. K. Molloy, D. B. Sloan, M. P. Simmons, and J. Gatesy. 2020.1058 ILS-Aware analysis of low-homoplasy retroelement insertions: Inference of species trees1059 and introgression Using quartets. Journal of Heredity 111:147–168.1060 Stamatakis, A. 2014. RAxML version 8: a tool for phylogenetic analysis and post-analysis1061 of large phylogenies. Bioinformatics 30:1312–1313.1062 Steenwyk, J., G. Mart´ ınez-Redondo, T. Buida, E. Gluck-Thaler, X.-X. Shen, T. Gabald´ on,1063 A. Rokas, and R. Fern´ andez. 2024. Phykit: A multitool for phylogenomics. Preprints .1064 Steenwyk, J. L., Y. Li, X. Zhou, X.-X. Shen, and A. Rokas. 2023. Incongruence in the1065 phylogenomics era. Nature Reviews Genetics 24:834–850.1066 Stiller, J., S. Feng, A.-A. Chowdhury, I. Rivas-Gonz´ alez, D. A. Duchˆ ene, Q. Fang,1067 Y. Deng, A. Kozlov, A. Stamatakis, S. Claramunt, J. M. T. Nguyen, S. Y. W. Ho, B. C.1068 Faircloth, J. Haag, P. Houde, J. Cracraft, M. Balaban, U. Mai, G. Chen, R. Gao,1069 C. Zhou, Y. Xie, Z. Huang, Z. Cao, Z. Yan, H. A. Ogilvie, L. Nakhleh, B. Lindow,1070 B. Morel, J. Fjelds˚ a, P. A. Hosner, R. R. da Fonseca, B. Petersen, J. A. Tobias,1071 T. Sz´ ekely, J. D. Kennedy, A. H. Reeve, A. Liker, M. Stervander, A. Antunes, D. T.1072 Tietze, M. F. Bertelsen, F. Lei, C. Rahbek, G. R. Graves, M. H. Schierup, T. Warnow,1073 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint 50 REFERENCES E. L. Braun, M. T. P. Gilbert, E. D. Jarvis, S. Mirarab, and G. Zhang. 2024.1074 Complexity of avian evolution revealed by family-level genomes. Nature 629:851–860.1075 Tavar´ e, S. 1986. Some probabilistic and statistical problems in the analysis of DNA1076 sequences. 17:57–86.1077 Vachaspati, P. and T. Warnow. 2015. ASTRID: Accurate species trees from internode1078 distances. BMC Genomics 16:S3.1079 Waterhouse, A. M., J. B. Procter, D. M. A. Martin, M. Clamp, and G. J. Barton. 2009.1080 Jalview version 2–a multiple sequence alignment editor and analysis workbench.1081 Bioinformatics 25:1189–1191.1082 Wilcoxon, F. 1949. Some Rapid Approximate Statistical Procedures Page 6. American1083 Cyanamid Company, New York.1084 Wu, S., F. E. Rheindt, J. Zhang, J. Wang, L. Zhang, C. Quan, Z. Li, M. Wang, F. Wu,1085 Y. Qu, S. V. Edwards, Z. Zhou, and L. Liu. 2024a. Genomes, fossils, and the concurrent1086 rise of modern birds and flowering plants in the late cretaceous. Proceedings of the1087 National Academy of Sciences 121:e2319696121.1088 Wu, S., F. E. Rheindt, J. Zhang, J. Wang, L. Zhang, C. Quan, Z. Li, M. Wang, F. Wu,1089 Y. Qu, S. V. Edwards, Z. Zhou, and L. Liu. 2024b. Reply to springer and gatesy: The1090 impact of long branches and misalignments on phylogenetic analysis is minimal.1091 Proceedings of the National Academy of Sciences 121:e2409344121.1092 Xi, Z., L. Liu, and C. C. Davis. 2015. Genes with minimal phylogenetic information are1093 problematic for coalescent analyses when gene tree estimation is biased. 92:63–71.1094 Yang, Z. 1993. Maximum-likelihood estimation of phylogeny from DNA sequences when1095 substitution rates differ over sites. Molecular Biology and Evolution 10:1396–1401.1096 Yoo, D., A. Rhie, P. Hebbar, F. Antonacci, G. A. Logsdon, S. J. Solar, D. Antipov, B. D.1097 Pickett, Y. Safonova, F. Montinaro, Y. Luo, J. Malukiewicz, J. M. Storer, J. Lin, A. N.1098 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint

51 Sequeira, R. J. Mangan, G. Hickey, G. M. Anez, P. Balachandran, A. Bankevich, C. R.1099 Beck, A. Biddanda, M. Borchers, G. G. Bouffard, E. Brannan, S. Y. Brooks, L. Carbone,1100 L. Carrel, A. P. Chan, J. Crawford, M. Diekhans, E. Engelbrecht, C. Feschotte,1101 G. Formenti, G. H. Garcia, L. d. Gennaro, D. Gilbert, R. E. Green, A. Guarracino,1102 I. Gupta, D. Haddad, J. Han, R. S. Harris, G. A. Hartley, W. T. Harvey, M. Hiller,1103 K. Hoekzema, M. L. Houck, H. Jeong, K. Kamali, M. Kellis, B. Kille, C. Lee, Y. Lee,1104 W. Lees, A. P. Lewis, Q. Li, M. Loftus, Y. H. E. Loh, H. Loucks, J. Ma, Y. Mao, J. F. I.1105 Martinez, P. Masterson, R. C. McCoy, B. McGrath, S. McKinney, B. S. Meyer, K. H.1106 Miga, S. K. Mohanty, K. M. Munson, K. Pal, M. Pennell, P. A. Pevzner, D. Porubsky,1107 T. Potapova, F. R. Ringeling, J. L. Rocha, O. A. Ryder, S. Sacco, S. Saha, T. Sasaki,1108 M. C. Schatz, N. J. Schork, C. Shanks, L. Smeds, D. R. Son, C. Steiner, A. P. Sweeten,1109 M. G. Tassia, F. Thibaud-Nissen, E. Torres-Gonz´ alez, M. Trivedi, W. Wei, J. Wertz,1110 M. Yang, P. Zhang, S. Zhang, Y. Zhang, Z. Zhang, S. A. Zhao, Y. Zhu, E. D. Jarvis,1111 J. L. Gerton, I. Rivas-Gonz´ alez, B. Paten, Z. A. Szpiech, C. D. Huber, T. L. Lenz, M. K.1112 Konkel, S. V. Yi, S. Canzar, C. T. Watson, P. H. Sudmant, E. Molloy, E. Garrison,1113 C. B. Lowe, M. Ventura, R. J. O’Neill, S. Koren, K. D. Makova, A. M. Phillippy, and1114 E. E. Eichler. 2024. Complete sequencing of ape genomes. bioRxiv .1115 Yule, G. U. 1925. A mathematical theory of evolution, based on the conclusions of Dr. J.1116 C. Willis, F.R.S. Philosophical Transactions of the Royal Society of London. Series B,1117 Containing Papers of a Biological Character 213:21–87.1118 Zhang, C. and S. Mirarab. 2022. Weighting by gene tree uncertainty improves accuracy of1119 quartet-based species trees. 39:msac215.1120 Zhang, C., M. Rabiee, E. Sayyari, and S. Mirarab. 2018. ASTRAL-III: Polynomial time1121 species tree reconstruction from partially resolved gene trees. 19:153.1122 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprintthis version posted September 30, 2024. ; https://doi.org/10.1101/2024.09.27.615467doi: bioRxiv preprint